实验三连续时间LTI系统的复频域分析.docx
《实验三连续时间LTI系统的复频域分析.docx》由会员分享，可在线阅读，更多相关《实验三连续时间LTI系统的复频域分析.docx（21页珍藏版）》请在冰点文库上搜索。
实验三连续时间LTI系统的复频域分析
实验三连续时间LTI系统的复频域分析
一、实验目的
1、掌握拉普拉斯变换的物理意义、基本性质及应用；
2、掌握用拉普拉斯变换求解连续时间LTI系统的时域响应；
3、掌握系统函数的概念，掌握系统函数的零、极点分布（零、极点图）与系统的稳定性、时域特性等之间的相互关系；
4、掌握用MATLAB对系统进行变换域分析的常用函数及编程方法。
基本要求：
掌握拉普拉斯变换及其基本性质，掌握应用拉普拉斯变换求解系统的微分方程，能够自己编写程序完成对系统时域响应的求解。
掌握并理解系统函数的概念，掌握系统函数零极点与系统时域和频域特性之间的关系，能够编写程序完成对系统的一些主要特性如稳定性、因果性等的分析。
二、实验原理及方法
1、连续时间LTI系统的复频域描述
拉普拉斯变换（TheLaplacetransform）主要用于系统分析。
描述系统的另一种数学模型就是建立在拉普拉斯变换基础上的“系统函数（SystemFunction）”——H（s）：
3.1
系统函数
的实质就是系统单位冲激响应（ImpulseResponse）
的拉普拉斯变换。
因此，系统函数也可以定义为：
3.2
所以，系统函数
的一些特点是和系统的时域响应
的特点相对应的。
在教材中，我们求系统函数的方法，除了按照拉氏变换的定义式的方法之外，更常用的是根据描述系统的线性常系数微分方程（LinearConstant-CoefficientDefrentialEquation），经过拉氏变换之后得到系统函数
。
假设描述一个连续时间LTI系统的线性常系数微分方程为：
3.3
对式3.3两边做拉普拉斯变换，则有
即
3.4
式3.4告诉我们，对于一个能够用线性常系数微分方程描述的连续时间LTI系统，它的系统函数是一个关于复变量s的有理多项式的分式，其分子和分母的多项式系数与系统微分方程左右两端的系数是对应的。
根据这一特点，可以很容易的根据微分方程写出系统函数表达式，或者根据系统函数表达式写出系统的微分方程。
系统函数
大多数情况下是复变函数，因此，
可以有多种表示形式：
1、直角坐标形式：
2、零极点形式：
3、部分分式和形式：
（假设系统的N>M，且无重极点）
根据我们所要分析的问题的不同，可以采用不同形式的系统函数
表达式。
在MATLAB中，表达系统函数
的方法是给出系统函数的分子多项式和分母多项式的系数向量。
由于系统函数的分子和分母的多项式系数与系统微分方程左右两端的系数是对应的，因此，用MATLAB表示系统函数，就是用系统函数的两个系数向量来表示。
应用拉普拉斯变换分析系统的主要内容有：
1、分析系统的稳定性；
2、分析系统的频率响应。
分析方法主要是通过绘制出系统函数的零极点分布图，根据零极点分布情况，判断系统的稳定性。
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）：
用来产生绘制平面图的区域，由x1，y1来确定具体的区域范围，由此产生s平面区域。
meshgrid（x,y,fs）：
绘制系统函数的零极点曲面图。
H=impulse（num,den）：
求系统的单位冲激响应，不带返回值，则直接绘制响应曲线，带返回值则将冲激响应值存于向量h之中。
2、系统函数的零极点分布图
系统函数的零极点图（Zero-polediagram）能够直观地表示系统的零点和极点在s平面上的位置，从而比较容易分析系统函数的收敛域（Reginofconvergence）和稳定性（stablity）。
下面给出一个用于绘制连续时间LTI系统的零极点图的扩展函数splane（num,den）：
%splane
%Thisfunctionisusedtodrawthezero-poleplotinthes-plane
functionsplane（num,den）
p=roots（den）;%Determinethepoles
q=roots（num）;%Determinethezeros
p=p';q=q';
x=max（abs（[pq]））;%Determinetherangeofreal-axis
x=x+1;
y=x;%Determinetherangeofimaginary-axis
plot（[-xx],[00],':
'）;holdon;%Drawthereal-axis
plot（[00],[-yy],':
'）;holdon;%Drawtheimaginary-axis
plot（real（p）,imag（p）,'x'）;holdon;%Drawthepoles
plot（real（q）,imag（q）,'o'）;holdon;%Drawthezeros
title（'zero-poleplot'）;
xlabel（'RealPart'）;ylabel（'ImaginalPart'）
axis（[-xx-yy]）;%Determinethedisplay-range
对于一个连续时间LTI系统，它的全部特性包括稳定性、因果性（Causality）和它具有何种滤波特性（Frequency-domainaspect）等完全由它的零极点在s平面上的位置所决定。
3、拉普拉斯变换与傅里叶变换之间的关系
根据课堂上所学的知识可知，拉普拉斯变换与傅里叶变换之间的关系可表述为：
傅里叶变换是信号在虚轴上的拉普拉斯变换，也可用下面的数学表达式表示
5.5
上式表明，给定一个信号h（t），如果它的拉普拉斯变换存在的话，它的傅里叶变换不一定存在，只有当它的拉普拉斯变换的收敛域包括了整个虚轴，则表明其傅里叶变换是存在的。
下面的程序可以以图形的方式，表现拉普拉斯变换与傅里叶变换的这种关系。
%Relation_ft_lt
%ThisprogramisusedtoobservetherelationshipbetweentheFouriertransform
%andtheLaplacetransformofarectangularpulse.
clear,closeall,
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）;%计算拉普拉斯变换
subplot（211）
mesh（a,b,c）;%绘制曲面图
surf（a,b,c）;
view（-60,20）%调整观察视角
axis（[-0,5,-20,20,0,2]）;
title（'TheLaplacetransformoftherectangularpulse'）;
w=-20:
0.1:
20;
Fw=（2*sin（w+eps）.*exp（i*（w+eps）））./（w+eps）;
subplot（212）;plot（w,abs（Fw））
title（'TheFouriertransformoftherectangularpulse'）
xlabel（'frequencew'）
上面的程序不要求完全读懂，重点是能够从所得到的图形中，观察和理解普拉斯变换与傅里叶变换之间的相互关系就行。
4、系统函数的极点分布与系统的稳定性和因果性之间的关系
一个稳定的LTI系统，它的单位冲激响应h（t）满足绝对可积条件，即
3.6
同时，我们还应该记得，一个信号的傅里叶变换的存在条件就是这个信号满足绝对可积条件，所以，如果系统是稳定的话，那么，该系统的频率响应也必然是存在的。
又根据傅里叶变换与拉普拉斯变换之间的关系，可进一步推理出，稳定的系统，其系统函数的收敛域必然包括虚轴。
稳定的因果系统，其系统函数的全部极点一定位于s平面的左半平面。
所以，对于一个给定的LTI系统，它的稳定性、因果性完全能够从它的零极点分布图上直观地看出。
例题3-1：
已知一个因果的LTI系统的微分方程为
编写程序，绘制出系统的零极点分布图，并说明它的稳定性如何。
解：
这是一个高阶系统，显然手工计算它的极点是很困难的。
可以利用前面给出的扩展函数splane（），来绘制系统的零极点分布图。
范例程序如下：
%Program5_1
%Thisprogramplotsthezero-polediagramofanLTIsystemdescribed
%bythelinearconstant-coefficientdifferentialequation
clear,closeall,
b=262;
a=[11048148306401262];
subplot（221）
splane（b,a）
title（'Thezero-polediagram'）
执行该程序后，得到系统的零极点分布图如图3.1所示。
由于已知该系统是因果系统，从零极点分布图上看，它的全部极点都位于s平面的左半平面上，所以系统是稳定的。
然后，直接在命令窗口键入
>>roots（a）
回车后，就得到系统的极点为：
ans=
-0.5707+2.4716i
-0.5707-2.4716i
-2.7378+0.0956i
-2.7378-0.0956i
-1.6915+1.6014i
-1.6915-1.6014i
若题目中没有说明该系统是否是因果的，则需要做详细的分析。
从零极点分布图上可以看出，该系统可能的收敛域共有四种可能，另外三种可能如下：
（a）收敛域为Re{s}<-2.7378，此种情况说明，该系统是一个反因果系统（Anticausalsystem），由于收敛域不包含虚轴，故此系统是不稳定的。
（b）、（c）收敛域为-2.7378总之，系统的稳定性主要取决于系统函数的收敛域是否包含整个虚轴，而系统的因果性则取决于系统极点位置的分布。
需要特别强调的是，MATLAB总是把由分子和分母多项式表示任何系统都当作是因果系统。
所以，利用impulse（）函数求得的单位冲激响应总是因果信号。
5、系统函数的零极点分布与系统的滤波特性
系统具有何种滤波特性，主要取决于系统的零极点所处的位置。
没有零点的系统，通常是一个低通滤波器。
例题3-2已知一个系统的系统函数为
显然，这是一个一阶系统，无零点。
为了确定该系统具有何种滤波特性，需要把系统的频率响应特性曲线绘制出来加以判断。
绘制系统的频率特性曲线如图3.2所示。
通过编程，可以将系统的零极点分布图和系统的频率响应特性以及系统的单位冲激响应特性绘制在一个图形窗口的各个子图中，这样便于观察系统的零极点分布情况与系统的时域和频域之间的关系。
如图3.3所示
图3.3系统的零极点分布与系统的单位冲激响应、频率相应
6、拉普拉斯及其逆变换的计算
我们已经知道，直接用拉普拉斯逆变换（Inversetransform）的定义公式计算逆变换是很困难的，通常的计算拉普拉斯逆变换的方法是长除法（Longdivision）和部分分式分解法（Partialfractionexpension）。
MATLAB的内部函数residue（）可以帮助我们完成拉普拉斯逆变换的计算，此函数实际上是将F（s）函数进行部分分式展开的函数。
例题3-3已知某信号的拉普拉斯变换表达式为
求该信号的时域表达式。
解：
由于题目没有指定收敛域，所以必须考虑所有可能的情况。
为此，可以先计算出该信号的拉普拉斯变换表达式的极点。
很显然，X（s）有两个极点，分别为s=-1，s=-2。
零极点分布图如例题图5-3所示。
在MATLAB命令窗口键入：
>>b=1;
>>a=[132];
>>[r,p,k]=residue（b,a）
命令窗口立即给出计算结果为：
r=
-1
1
p=
-2
-1
k=
[]
根据r、p、k之值，可以写出X（s）的部分分式和的表达式为：
然后根据不同的ROC，可写出X（s）的时域表达式x（t）。
第一种情况，ROC为Re{s}<-2，则x（t）为反因果信号，其数学表达式为
第二种情况，ROC为-2第三种情况，ROC为Re{s}>-1，则x（t）为因果信号，其数学表达式为
在这个例题中，函数residue（）仅仅完成了部分分式分解的任务，至于逆变换的数学表达式的结果是什么，还得结合收敛域的不同才能写出。
如果X（s）的分子的阶不小于分母的阶，则k将不等于一个空矩阵，例如，当
时，我们在命令窗口中键入：
>>b=[1000];
>>a=[132];
>>[r,p,k]=residue（b,a）
则：
r=
8
-1
p=
-2
-1
k=
1-3
这里的k=[13]，实际上是将X（s）做了一个长除法后，得到的商的多项式。
所以，根据上面的r、p、k之值，可写出X（s）的部分分式和的表达式为：
有关函数residue（）的详细用法，可通过在线帮助加以了解。
MATLAB工具箱中有函数laplace（f）和ilaplace（F）函数来进行拉氏变换和拉氏反变换，这属于符号运算范畴。
详情请看教材P195。
三、实验内容及步骤
实验前，必须首先阅读本实验原理，了解所给的MATLAB相关函数，读懂所给出的全部范例程序。
实验开始时，先在计算机上运行这些范例程序，观察所得到的信号的波形图。
并结合范例程序所完成的工作，进一步分析程序中各个语句的作用，从而真正理解这些程序。
实验前，一定要针对下面的实验项目做好相应的实验准备工作，包括事先编写好相应的实验程序等事项。
Q3-1将绘制零极点图的扩展函数文件splane以splane为文件名存盘。
Q3-2运行程序Relation_ft_lt，观察拉普拉斯变换与傅里叶变换之间的关系。
在点击工具条上的旋转按钮，再将鼠标放在曲面图上拖动图形旋转，从各个角度观察拉普拉斯曲面图形，并同傅立叶变换的曲线图比较，加深对拉普拉斯变换与傅里叶变换之间关系的理解与记忆。
Q3-3系统函数为
，编写程序Q3_3，绘制出系统的零极点图、系统的单位冲激响应、系统的幅度频率响应和相位频率相应的图形。
各个子图要求按照图3.3布置
formatrat
num=[1,0];
den=[1,3,2];
[r,p]=residue（num,den）;
subplot（221）
splane（num,den）
title（'Thezero-polediagram'）
>>w=-10:
0.01:
10;b=[1,0];
a=[1,3,2];
[H,w]=freqs（b,a）;
Hm=abs（H）;
phai=angle（H）;
subplot（222）
plot（w,Hm）,gridon,title（'Magnituderesponse'）,xlabel（'Frequencyinrad/sec'）
subplot（224）
plot（w,phai）,gridon,title（'Phaseresponse'）,xlabel（'Frequencyinrad/sec'）
>>s=-10:
0.01:
10;
h=s./（s.^2+3.*s+2）;
subplot（223）
plot（s,h）
title（'Theunitimpulseresponse'）
从上面的图形中可以看出，该系统的零点和极点分别位于：
答：
零点位于原点，极点位于x轴左边
从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？
答：
H（s）的所有极点位于s的左半平面，且H（t）将随时间t衰减，当t无穷时，H（t）消失，所以为稳定系统。
从频率响应特性上看，该系统具有何种滤波特性？
答：
低通
Q3-4已知因果系统函数为：
编写程序Q3_3，绘制出系统的零极点图、系统的单位冲激响应、系统的幅度频率响应和相位频率相应的图形。
各个子图要求按照图3.3布置。
此处a取1。
图形如下：
>>formatrat
a=1;
num=[1/a,0,1];
den=[1,2,2,1];
[r,p]=residue（num,den）;
subplot（221）
splane（num,den）
title（'Thezero-polediagram'）
s=-10:
0.01:
10;
[H,s]=freqs（num,den）;
Hm=abs（H）;
phai=angle（H）;
subplot（222）
plot（s,Hm）,gridon,title（'Magnituderesponse'）,xlabel（'Frequencyinrad/sec'）
subplot（224）
plot（s,phai）,gridon,title（'Phaseresponse'）,xlabel（'Frequencyinrad/sec'）
h=s./（s.^2+3.*s+2）;
subplot（223）
plot（s,h）
title（'Theunitimpulseresponse'）
从上面的图形中可以看出，该系统的零点和极点分别位于：
答，零点位于虚轴，极点左平面。
从时域和零极点分布特征两个方面说明该系统是否是稳定的系统？
答：
H（s）的所有极点位于s的左半平面，且H（t）将随时间t衰减，当t无穷时，H（t）消失，所以为稳定系统。
从频率响应特性上看，该系统具有何种滤波特性？
答：
低通
改变系统函数中的a值，分别取0.6、0.8、4、16等不同的值，反复执行程序Q3-4，观察系统的幅度频率响应特性曲线（带宽、过渡带宽和阻带衰减等），贴一张a=4时的图形如下：
formatrat
a=4;
num=[1/a,0,1];
den=[1,2,2,1];
[r,p]=residue（num,den）;
subplot（221）
splane（num,den）
title（'Thezero-polediagram'）
s=-10:
0.01:
10;
[H,s]=freqs（num,den）;
Hm=abs（H）;
phai=angle（H）;
subplot（222）
plot（s,Hm）,gridon,title（'Magnituderesponse'）,xlabel（'Frequencyinrad/sec'）
subplot（224）
plot（s,phai）,gridon,title（'Phaseresponse'）,xlabel（'Frequencyinrad/sec'）
h=s./（s.^2+3.*s+2）;
subplot（223）
plot（s,h）
title（'Theunitimpulseresponse'）
>>
观察a取不同的值时系统的幅度频率响应特性曲线的变化（带宽、过渡带宽和阻带衰减等），请用一段文字说明零点位置对系统滤波特性的这些影响。
答：
a变大，带宽变大，过渡带宽变大，阻带衰减变慢。
Q3-5已知
，编程求其反变换。
f=sym（'（2*s+1）/（s^2+7*s+10）'）;
ft=ilaplace（f）
ft=
3/exp（5*t）-1/exp（2*t）
Q3-6已知一个因果系统的系统函数为
，作用于系统的输入信号为
，试用MATLAB帮助你求系统的响应信号y（t）的数学表达式。
请在这里抄写你用MATLAB求解的命令（结合必要的文字说明）：
x=sym（'exp（-4*t）*heaviside（t）'）;
X=laplace（x）
X=
1/（s+4）%求得x（t）的频域响应
ys=sym（'（（s+5）/（s^3+6*s^2+11*s+6））*（1/（s+4））'）;
Ys=ilaplace（ys）
这里抄写程序执行结果。
Ys=
2/（3*exp（t））-3/（2*exp（2*t））+1/exp（3*t）-1/（6*exp（4*t））
四、实验报告要求
1、按要求完整书写你所编写的全部MATLAB程序
2、详细记录实验过程中的有关信号波形图（存于自带的U盘中），图形要有明确的标题。
全部的MATLAB图形应该用打印机打印，然后贴在本实验报告中的相应位置，禁止复印件。
3、实事求是地回答相关问题，严禁抄袭。
本实验完成时间：
年月日
升级会员 