连续时间LTI系统的时域分析

一、课程设计题目：基于 MATLAB 的连续时间LTI 系统的时域分析二、基本要求：① 掌握连续时不变信号处理的基本概念、基本理论和基本方法； ② 学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法； ③ 学会用 MATLAB 对信号进行分析和处理；④ 编程实现卷积积分或卷积和，零输入响应，零状态响应； ⑤ 撰写课程设计论文，用信号处理基本理论分析结果。

三、设计方法与步骤：一般的连续时间系统分析有以下几个步骤: ①求解系统的零输入响应; ②求解系统的零状态响应; ③求解系统的全响应; ④分析系统的卷积；⑤画出它们的图形. 下面以具体的微分方程为例说明利用MATLAB 软件分析系统的具体方法.1．连续时间系统的零输入响应描述n 阶线性时不变（LTI ）连续系统的微分方程为：已知y 及各阶导数的初始值为y(0),y (1)(0)，… y (n-1)(0)， 求系统的零输入响应。

建模当LIT 系统的输入为零时，其零输入响应为微分方程的其次解（即令微分方程的等号右端为零），其形式为（设特征根均为单根）其中p 1,p 2,…,p n 是特征方程a 1λn +a 2λn-1+…+a n λ+a n =0的根，它们可以用root(a)语句求得。

各系数 由y 及其各阶导数的初始值来确定。

对此有………………………………………………………………………………………写成矩阵形式为: P 1n-1C 1+ P 2n-1C 2+…+ P n n-1C n =D n-1y 01121111n n m n n m m n n m d y d y dy d u du a a a a y b b b u dtdt dt dt dt -++-++⋅⋅⋅⋅⋅++=+⋅⋅⋅⋅++1212()n p t p t p t n y t C e C e C e =++⋅⋅⋅⋅+120n C C C y ++⋅⋅⋅⋅+=11220n n p C p C p C Dy ++⋅⋅⋅⋅+=111111220n n n n n n p C p C p C D y----++⋅⋅⋅⋅+=1012201111120111n n n n n n n C y p p p C Dy p p p C D y ----⋅⋅⋅⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⋅⋅⋅⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⋅⋅⋅⎣⎦⎣⎦⎣⎦即 V •C=Y 0 其解为：C=V\Y 0 式中V 为范德蒙矩阵，在matlab 的特殊矩阵库中有vander 。

以下面式子为例：y(0_)=1,y(0_)=5；MATLAB 程序：a=input('输入分母系数a=[a1,a2,...]=');n=length(a)-1;Y0=input('输入初始条件向量 Y0=[y0,Dy0,D2y0,...]='); p=roots(a);V=rot90(vander(p));c=V\Y0'; dt=input('dt=');te=input('te='); t=0:dt:te;y=zeros(1,length(t)); for k=1:n y=y+c(k)*exp(p(k)*t);end plot(t,y);gridxlabel('t') ;ylabel('y'); title('零输入响应');程序运行结果：用这个通用程序来解一个三阶系统，运行此程序并输入 a=[1,5,4] Y0=[1,5] dt=0.01 te=6结果如下图：根据图可以分析零输入响应，它的起始值与输入函数无关，只与它的初始状态值有关，其起始值等于y(0_)的值。

随着时间的推移，最后零输入响应的值无限的趋近于0。

1211112111n n n n n p p p V p p p ---⋅⋅⋅⎡⎤⎢⎥⋅⋅⋅⎢⎥=⎢⎥⎢⎥⋅⋅⋅⎣⎦[]12n C C C C =⋅⋅⋅1000n C y Dy D y -⎡⎤=⋅⋅⋅⎣⎦''''()5()4()2()4()y t y t y t f t f t ++=-ty零输入响应2．卷积的计算连续时间信号1()f t 和2()f t 的卷积运算可用信号的分段求和来实现，即：1212120()()*()()()lim()()k f t f t f t f t f t d f k f t k ττ∞∞-∞∆→=-∞==-=∆-∆⋅∆∑⎰如果只求当t = n ∆（n 为整数）时f (t )的值f (n ∆) ，则上式可得：1212()()()()[()]k k f n f k f t k f k f n k ∞∞=-∞=-∞∆=∆-∆⋅∆=∆∆-∆∑∑式中的12()[()]k f k f n k ∞=-∞∆∆-∆∑实际上就是连续时间信号1()f t 和2()f t 经等时间间隔 均匀抽样的离散序列1()f k ∆和2()f k ∆的卷积和。

当∆足够小时，()f n ∆就是卷积积分的结果——连续时间信号f (t )的较好数值近似。

建模下面是利用MATLAB 实现连续信号卷积的通用程序conv()，该程序在计算出卷积积分的数值近似的同时，还绘制出f (t )的时域波形图。

应注意，程序中是如何设定f (t )的时间长度。

MATLAB 程序：f1=input('输入函数f1='); f2=input('输入函数f2='); dt=input('dt='); y=conv(f1,f2);plot(dt*([1:length(y)]-1),y);grid on; title('卷积');xlabel('t'); ylabel('f1*f2')程序运行结果： 输入以下数据：f1=sin(3*t) f2=cos(3*t+2) dt=0.01 得出图形如下：卷积tf 1*f 23．连续时间系统零状态响应的数值计算我们知道，LTI 连续系统可用如下所示的线性常系数微分方程来描述，()()0()()NMi j i j i j a yt b f t ===∑∑例如，对于以下方程:''''''''''''32103210()()()()()()()()a y t a y t a y t a y tb f t b f t b f t b f t +++=+++可用32103210[,,,],[,,,],a a a a a b b b b b ==输入函数()u f t =，得出它的冲击响应h ，再根据LTI 系统的零状态响应y （t ）是激励u （t ）与冲击响应h （t ）的卷积积分。

注意，如果微分方程的左端或右端表达式中有缺项，则其向量a 或b 中的对应元素应为零，不能省略不写，否则出错。

求函数的零状态响应及初始状态'(0)(0)0zs zs y y --==。

输入函数()sin(3*)cos(2*)f t t t =+。

建模先求出系统的冲击响应，写出其特征方程2540λλ++=求出其特征根为p 和p ，及相应的留数r ，r ；则冲击响应为''''()5()4()2()4()y t y t y t f t f t ++=-1212()p t p th t r e r e =+ 输入y （t ）可用输入u （t ）与冲击响应h （t ）的卷积求得。

MATLAB 程序：a=input('输入分母系数a=[a1,a2,...]='); b=input('输入输入信号系数b=[b1,b2,...]='); dt=input('dt=');te=input('te='); t=0:dt:te;u=input('输入函数u='); te=t(end);dt=te/(length(t)-1); [r,p,k]=residue(b,a);h=r(1)*exp(p(1)*t)+r(2)*exp(p(1)*t); subplot(2,1,1),plot(t,h);grid title('冲击函数'); y=conv(u,h)*dt; subplot(2,1,2),plot(t,y(1:length(t)));grid title('零状态响应');程序运行结果执行这个程序，取a=[1,5,4] b=[2,4] dt=0.01 te=6sin(3*)cos(2*)u t t =+得出图形如下：由于初始状态为零，所以零状态的起始值也为零，即h (t )包含了连续系统的固有特性，与系统的输入无关。

只要知道了系统的冲激响应，即可求得系统在不同输入时产生的输出。

因此，求解系统的冲激响应h 对进行连续时间系统的分析具有非常重要的意义冲击函数零状态响应4．连续时间系统的全响应计算上面通过对LTI 系统函数的描述，我们可以得知：如果在系统的初始状态不为零，在激励f （t ）的作用下，LTI 系统的响应称为全响应，它是零输入响应和零状态响应之和，即()()()zi zs y t y t y t =+故可先求出零输入响应和零状态响应，再把两者相加，得到全响应。

但简单的相加可能由于零输入与零状态的矩阵不同而不能的出正确的结果，这就需要对矩阵进行截取，使它们的阶数相同。

例如，对于以下方程：初始值为：y(0_)=1,y(0_)=5； 输入函数为：()sin(3*)cos(2*)f t t t =+求它的全响应。

建模先根据零输入响应的求法，得出零输入响应y1（t ）。

再根据零状态响应的求法，得出零状态响应y2（t ）。

最后，全响应y 等于零输入响应y1（t ）加上零状态响应y2（t ），得出全响应。

MATLAB 程序：a=input('输入分母系数a=[a1,a2,...]='); n=length(a)-1;Y0=input('输入初始条件向量 Y0=[y0,Dy0,D2y0,...]='); b=input('输入输入信号系数b=[b1,b2,...]='); u=input('输入函数u=');dt=input('dt=');te=input('te='); t=0:dt:te;p=roots(a);V=rot90(vander(p));c=V\Y0'; y1=zeros(1,length(t));for k=1:n y1=y1+c(k)*exp(p(k)*t);end te=t(end);dt=te/(length(t)-1); [r,p,k]=residue(b,a);h=r(1)*exp(p(1)*t)+r(2)*exp(p(1)*t); y2=conv(u,h)*dt;y=y1(1:length(t))+y2(1:length(t)); figure(1);subplot(3,1,1),plot(t,y1),gridxlabel('t'); ylabel('y1');title('零输入响应'); subplot(3,1,2),plot(t,y2(1:length(t)));grid xlabel('t');ylabel('y2');title('零状态响应');''''()5()4()2()4()y t y t y t f t f t ++=-subplot(3,1,3),plot(t,y),gridxlabel('t'); ylabel('y');title('全响应响应');程序运行结果执行程序，取 a=[1,5,4] Y0=[1,5] b=[1,2,4] u=sin(3*t)+cos(2*t) dt=0.01 te=6结果如下图：在零输入响应中任一时刻取值y1，在零状态响应的对应时刻取值y2，再在全响应的对应时刻取值y 。