1、实验六：连续时间 LTI 系统的复频域分析一、实验目的1、掌握拉普拉斯变换的物理意义、基本性质及应用。2、掌握用拉普拉斯变换求解连续时间 LTI 系统的时域响应。3、掌握系统函数的概念，掌握系统函数的零、极点分布（零、极点图）与系统的稳定性、时域特性等之间的相互关系。4、掌握用 MATLAB 对系统进行变换域分析的常用函数及编程方法。二、实验原理1、连续时间 LTI 系统的复频域描述拉普拉斯变换（The Laplace transform）主要用于系统分析。描述系统的另一种数学模型就是建立在拉普拉斯变换基础上的“系统函数（System Function） ”H(s)：6.1)()( txLsX
2、yYH换系 统 激 励 信 号 的 拉 氏 变 换系 统 冲 击 响 应 的 拉 氏 变系统函数 的实质就是系统单位冲激响应（Impulse Response） 的拉普拉斯)(th变换。因此，系统函数也可以定义为：6.2dtehsHs)(所以，系统函数 的一些特点是和系统的时域响应 的特点相对应的。在教材)( )(th中，我们求系统函数的方法，除了按照拉氏变换的定义式的方法之外，更常用的是根据描述系统的线性常系数微分方程（Linear Constant-Coefficient Defrential Equation） ，经过拉氏变换之后得到系统函数 。)(sH假设描述一个连续时间 LTI 系统
3、的线性常系数微分方程为：6.3MkkNkkdtxbdtya00 )()(对式 6.3 两边做拉普拉斯变换，则有MkNk sXbsYa00 )()(即 6.4NkMksabXsYH0)(式 6.4 告诉我们，对于一个能够用线性常系数微分方程描述的连续时间 LTI 系统，它的系统函数是一个关于复变量 s 的有理多项式的分式，其分子和分母的多项式系数与系统微分方程左右两端的系数是对应的。根据这一特点，可以很容易的根据微分方程写出系统函数表达式，或者根据系统函数表达式写出系统的微分方程。系统函数 大多数情况下是复变函数，因此， 可以有多种表示形式：)(sH)(sH1、直角坐标形式： Im)Re(js2
4、、零极点形式： NiiMjjpszk1)()(3、部分分式和形式： （假设系统的 NM，且无重极点） kksAH0)(根据我们所要分析的问题的不同，可以采用不同形式的系统函数 表达式。)(sH在 MATLAB 中，表达系统函数 的方法是给出系统函数的分子多项式和分母多项式)(s的系数向量。由于系统函数的分子和分母的多项式系数与系统微分方程左右两端的系数是对应的，因此，用 MATLAB 表示系统函数，就是用系统函数的两个系数向量来表示。应用拉普拉斯变换分析系统的主要内容有：1、分析系统的稳定性；2、分析系统的频率响应。分析方法主要是通过绘制出系统函数的零极点分布图，根据零极点分布情况，判断系统的
5、稳定性。MATLAB 中有相应的复频域分析函数，下面简要介绍如下：z,p,k = tf2zp(num,den)：求系统函数的零极点，返回值 z 为零点行向量，p为极点行向量，k 为系统传递函数的零极点形式的增益。num 为系统函数分子多项式的系数向量，den 为系统函数分母多项式系数向量。H = freqs(num,den,w)：计算由 num,den 描述的系统的频率响应特性曲线。返回值 H 为频率向量规定的范围内的频率响应向量值。如果不带返回值 H，则执行此函数后，将直接在屏幕上给出系统的对数频率响应曲线（包括幅频特性取向和相频特性曲线） 。x,y = meshgrid(x1,y1)：用来
6、产生绘制平面图的区域，由 x1，y1 来确定具体的区域范围，由此产生 s 平面区域。meshgrid(x,y,fs)：绘制系统函数的零极点曲面图。H = impulse(num,den)：求系统的单位冲激响应，不带返回值，则直接绘制响应曲线，带返回值则将冲激响应值存于向量 h 之中。2、系统函数的零极点分布图系统函数的零极点图（Zero-pole diagram）能够直观地表示系统的零点和极点在 s 平面上的位置，从而比较容易分析系统函数的收敛域（Regin of convergence）和稳定性（stablity） 。下面给出一个用于绘制连续时间 LTI 系统的零极点图的扩展函数 splan
7、e(num,den)：% splane% This function is used to draw the zero-pole plot in the s-planefunction splane(num,den)p = roots(den); % Determine the polesq = roots(num); % Determine the zerosp = p; q = q; x = max(abs(p q); % Determine the range of real-axis x = x+1;y = x; % Determine the range of imaginary-a
8、xis plot(-x x,0 0,:);hold on; % Draw the real-axisplot(0 0,-y y,:);hold on; % Draw the imaginary-axisplot(real(p),imag(p),x);hold on; % Draw the polesplot(real(q),imag(q),o);hold on; % Draw the zerostitle(zero-pole plot);xlabel(Real Part);ylabel(Imaginal Part)axis(-x x -y y); % Determine the display
9、-range对于一个连续时间 LTI 系统，它的全部特性包括稳定性、因果性（Causality）和它具有何种滤波特性（Frequency-domain aspect）等完全由它的零极点在 s 平面上的位置所决定。3、拉普拉斯变换与傅里叶变换之间的关系根据课堂上所学的知识可知，拉普拉斯变换与傅里叶变换之间的关系可表述为：傅里叶变换是信号在虚轴上的拉普拉斯变换，也可用下面的数学表达式表示6.5jsHj)(上式表明，给定一个信号 h(t)，如果它的拉普拉斯变换存在的话，它的傅里叶变换不一定存在，只有当它的拉普拉斯变换的收敛域包括了整个虚轴，则表明其傅里叶变换是存在的。下面的程序可以以图形的方式，表现
10、拉普拉斯变换与傅里叶变换的这种关系。% Relation_ft_lt% This program is used to observe the relationship between the Fourier transform% and the Laplace transform of a rectangular pulse.clear, close all,a = -0:0.1:5;b = -20:0.1:20;a, b = meshgrid (a, b);c = a+i*b; %确定绘图区域c = (1-exp (-2* (c+eps)./ (c+eps);c = abs (c); %计
11、算拉普拉斯变换subplot (211)mesh (a,b,c); %绘制曲面图surf (a,b,c);view (-60,20) %调整观察视角axis (-0,5,-20,20,0,2);title (The Laplace transform of the rectangular pulse);w = -20:0.1:20;Fw = (2*sin(w+eps).*exp(i*(w+eps)./(w+eps);subplot (212); plot(w,abs(Fw)title (The Fourier transform of the rectangular pulse)xlabel
12、(frequence w)上面的程序不要求完全读懂，重点是能够从所得到的图形中，观察拉和理解普拉斯变换与傅里叶变换之间的相互关系就行。4、系统函数的极点分布与系统的稳定性和因果性之间的关系一个稳定的 LTI 系统，它的单位冲激响应 h(t)满足绝对可积条件，即6.6dth)(同时，我们还应该记得，一个信号的傅里叶变换的存在条件就是这个信号满足绝对可积条件，所以，如果系统是稳定的话，那么，该系统的频率响应也必然是存在的。又根据傅里叶变换与拉普拉斯变换之间的关系，可进一步推理出，稳定的系统，其系统函数的收敛域必然包括虚轴。稳定的因果系统，其系统函数的全部极点一定位于 s 平面的左半平面。所以，对于
13、一个给定的 LTI 系统，它的稳定性、因果性完全能够从它的零极点分布图上直观地看出。例题 6-1：已知一个因果的 LTI 系统的微分方程为 )(26)()(401)(306)(148)()(10)( 2456 txtydttydtydtytydty 编写程序，绘制出系统的零极点分布图，并说明它的稳定性如何。解：这是一个高阶系统，显然手工计算它的极点是很困难的。可以利用前面给出的扩展函数 splane()，来绘制系统的零极点分布图。范例程序如下：% Program6_1% This program plots the zero-pole diagram of an LTI system desc
14、ribed% by the linear constant-coefficient differential equationclear, close all,b = 262;a = 1 10 48 148 306 401 262;subplot (221)splane (b,a)title (The zero-pole diagram)执行该程序后，得到系统的零极点分布图如图6.1 所示。由于已知该系统是因果系统，从零极点分布图上看，它的全部极点都位于 s 平面的左半平面上，所以系统是稳定的。然后，直接在命令窗口键入 roots(a)回车后，就得到系统的极点为：ans =-0.5707 +
15、2.4716i-0.5707 - 2.4716i-2.7378 + 0.0956i-2.7378 - 0.0956i-1.6915 + 1.6014i-1.6915 - 1.6014i若题目中没有说明该系统是否是因果的，则需要做详细的分析。从零极点分布图上可以看出，该系统可能的收敛域共有四种可能，另外三种可能如下：（a） 收敛域为 Res b = 1; a = 1 3 2; r, p, k = residue (b, a)命令窗口立即给出计算结果为：r =-11p =-2-1k = 根据 r、p、k 之值，可以写出 X(s)的部分分式和的表达式为：12)(ssX然后根据不同的 ROC，可写出
16、X(s)的时域表达式 x(t)。第一种情况，ROC 为 Res -1，则 x(t)为因果信号，其数学表达式为)()()(2tuettx在这个例题中，函数 residue()仅仅完成了部分分式分解的任务，至于逆变换的数学表达式的结果是什么，还得结合收敛域的不同才能写出。如果 X(s)的分子的阶不小于分母的阶，则 k 将不等于一个空矩阵，例如，当例题图 6-3时，我们在命令窗口中键入：23)(2ssX b = 1 0 0 0; a = 1 3 2; r,p,k=residue(b,a)则：r =8-1p =-2-1k =1 -3这里的 k = 1 3，实际上是将 X(s)做了一个长除法后，得到的商
17、的多项式。所以，根据上面的 r、p、k 之值，可写出 X(s)的部分分式和的表达式为：1283)(ssX有关函数 residue()的详细用法，可通过在线帮助加以了解。三、实验内容及步骤1. 将绘制零极点图的扩展 splane 为文函数文件 splane 以件名存盘。% splane% This function is used to draw the zero-pole plot in the s-planefunction splane(num,den)p = roots(den); % Determine the polesq = roots(num); % Determine the
18、zerosp = p; q = q; x = max(abs(p q); % Determine the range of real-axis x = x+1;y = x; % Determine the range of imaginary-axis plot(-x x,0 0,:);hold on; % Draw the real-axisplot(0 0,-y y,:);hold on; % Draw the imaginary-axisplot(real(p),imag(p),x);hold on; % Draw the polesplot(real(q),imag(q),o);hol
19、d on; % Draw the zerostitle(zero-pole plot);xlabel(Real Part);ylabel(Imaginal Part)axis(-x x -y y); % Determine the display-range2. 运行程序 Relation_ft_lt，观察拉普拉斯变换与傅里叶变换之间的关系。在点击工具条上的旋转按钮，再将鼠标放在曲面图上拖动图形旋转，从各个角度观察拉普拉斯曲面图形，并同傅立叶变换的曲线图比较，加深对拉普拉斯变换与傅里叶变换之间关系的理解与记忆。3. 因果系统函数 ，绘制出系统的零极点图、系统的单位冲激响应、)2(1)(ssH系
20、统的幅度频率响应和相位频率相应的图形。上面的图形中可以看出，该系统的零点和极点分别位于：从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？从频率响应特性上看，该系统具有何种滤波特性？% Program3_1% This Program is used to compute and draw the plots of the frequency response % of a continuous-time systemb = 1; % The coefficient vector of the right side of the differential equationa = 1 3
21、 2; % The coefficient vector of the left side of the differential equationH,w = freqs(b,a); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the real part of the frequency responseHi = imag(H);
22、 % Compute the imaginary part of the frequency responsesubplot(221)plot(w,Hm), grid on, title(Magnitude response), xlabel(Frequency in rad/sec)subplot(223)plot(w,phai), grid on, title(Phase response), xlabel(Frequency in rad/sec)subplot(222)plot(w,Hr), grid on, title(Real part of frequency response)
23、, xlabel(Frequency in rad/sec)subplot(224)plot(w,Hi), grid on, title(Imaginary part of frequency response), xlabel(Frequency in rad/sec)-5 0 5-505 Magnitude responseFrequency in rad/secImaginal Part0 5 10-0.200.20.40.6 Real part of frequency responseFrequency in rad/sec0 5 10-3-2-10 Phase responseFr
24、equency in rad/sec0 5 10-0.4-0.3-0.2-0.10Imaginary part of frequency responseFrequency in rad/sec4.因果系统的系统函数12)(3ssaH此处 a 取 1，绘制出系统的零极点图、系统的单位冲激响应、系统的幅度频率响应和相位频率相应的图形。上面的图形中可以看出，该系统的零点和极点分别位于：从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？从频率响应特性上看，该系统具有何种滤波特性？改变系统函数中的 a 值，分别取 0.6、0.8、4、16 等不同的值，反复执行程序，观察系统的幅度频率响应特性曲
25、线（带宽、过渡带宽和阻带衰减等） ，贴一张 a = 4 时的图形。观察 a 取不同的值时系统的幅度频率响应特性曲线的变化（带宽、过渡带宽和阻带衰减等） ，请用一段文字说明零点位置对系统滤波特性的这些影响。% Program3_1% This Program is used to compute and draw the plots of the frequency response % of a continuous-time systemb = 1 0 1; % The coefficient vector of the right side of the differential equa
26、tiona = 1 2 2 1; % The coefficient vector of the left side of the differential equationH,w = freqs(b,a); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the real part of the frequency response
27、Hi = imag(H); % Compute the imaginary part of the frequency responsesubplot(221)plot(w,Hm), grid on, title(Magnitude response), xlabel(Frequency in rad/sec)subplot(223)plot(w,phai), grid on, title(Phase response), xlabel(Frequency in rad/sec)subplot(222)plot(w,Hr), grid on, title(Real part of freque
28、ncy response), xlabel(Frequency in rad/sec)subplot(224)plot(w,Hi), grid on, title(Imaginary part of frequency response), xlabel(Frequency in rad/sec)-5 0 5-505 Magnitude responseFrequency in rad/secImaginal Part0 5 10-0.500.51 Real part of frequency responseFrequency in rad/sec0 5 10-3-2-101 Phase r
29、esponseFrequency in rad/sec0 5 10-1-0.500.5Imaginary part of frequency responseFrequency in rad/sec5. 对于因果系统 ，已知输入信号为 ，要12)(3ssaH )8sin()i(tttx求输出信号 ，K 为一个不为零的系数，根据 4 所得到的不同 a 值时的幅度intty频率响应图形，选择一个合适的 a 值从而使本系统能够实现本题的滤波要求。你选择的 a 值为：选择 a 值的根据是：试编写一个 MATLAB 程序，仿真这个滤波过程，要求绘制出系统输入信号、系统的单位冲激响应和系统的输出信号波形。
30、你选择的 a 值为： a=1/16选择 a 值的根据是：因为 X（t）=sin（t）+sin(8t),只要把 sin（8t）滤掉就满足要求。从上面的图形看随着 a 的值增大滤波越来越好。故选择 16试编写一个 MATLAB 程序 Q5_6，仿真这个滤波过程，要求绘制出系统输入信号、系统的单位冲激响应和系统的输出信号波形。抄写程序 Q5_6 如下：b=1/16 0 1 % The coefficient vector of the right side of the differential equationa=1 2 2 1 % The coefficient vector of the le
31、ft side of the differential equationH,w = freqs(b,a); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the real part of the frequency responseHi = imag(H); % Compute the imaginary part of the f
32、requency responsesubplot(221)splane (b,a), grid on, title(The zero-pole diagram), xlabel(Real Part)%plot(w,Hm)subplot(222)plot(w,Hm), grid on, title(Magnitude response)%, xlabel(Frequency in rad/sec)%plot(w,phai)subplot(223)plot(w,Hi), grid on, title(The impulse response)subplot(224)plot(w,phai), gr
33、id on, title(Phase response)b =0.0625 0 1.0000a =1 2 2 1-5 0 5-505 The zero-pole diagramReal PartImaginal Part0 5 1000.51 Magnitude response0 5 10-1-0.500.5 The impulse response0 5 10-4-2024 Phase response6. 已知一个因果系统的系统函数为 ，作用于系统的输入信号为61)(23ssH，试用 MATLAB 帮助你求系统的响应信号 y(t)的数学表达式。)()(4tuetxf=exp(-4*t)*u(t); L=laplace(f)L =1/(s + 4)H=(s+5)/(s3+6*s2+11*s+6);Y=H*fY =(s + 5)/(s + 1)*(s3 + 6*s2 + 11*s + 6)b=1,5;a=1 7 17 17 6;r, p, k = residue (b, a)r =-0.50003.0000-2.50002.0000p =-3.0000-2.0000-1.0000-1.0000k =四. 实验报告要求 1实验目的。2实验内容与步骤。用MATLAB语言完成编程，并附上仿真结果。3简要回答问题。4记录调试运行情况及所遇问题的解决方法。
