信号与系统的SSD分析

一.实验目的学会运用MATLAB实现状态方程与系统函数之间的互换；学会运用MATLAB在变换域求解状态方程与输出方程；学会运用MATLAB在时域求解状态方程与输出方程；学会运用MATLAB数值求解系统方程。二.实验原理及实例分析1.状态方程与系统函数之间的互换连续时间系统的状态方程与输出方程用矩阵可表示为{𝑥̇(𝑡)=𝐴𝑥(𝑡)+𝐵𝑓(𝑡)𝑦(𝑡)=𝐶𝑥(𝑡)+𝐷𝑓(𝑡)离散时间LTI系统的状态方程与输出方程用矩阵可表示为{𝑥(n+1)=Ax(n)+B𝑓(n)𝑦(𝑛)=𝐶𝑥(𝑛)+𝐷𝑓(𝑛)MATLAB控制系统工具箱提供了ss2tf和tf2ss两个函数，来实现系统的状态空间(ss)表示法和系统函数(tf)表示法之间的互换。tf2ss函数是将一个系统的系统函数转化为状态空间表示法，其语句格式为[A,B,C,D]=tf2ss(num,den)其中，num和den分别表示系统函数H(s)或H(z)的分子和分母多项式的系数，tf2ss的返回值A、B、C、D分别为状态方程与出处方程的矩阵。【实例1】已知某连续系统的系统函数为H(s)=4𝑠+10𝑠3+8𝑠2+19𝑠+12试用MATLAB命令求该系统的状态方程与输出方程。解：利用tf2ss函数求出矩阵，MATLAB源程序为[A,B,C,D]=tf2ss([410],[181912])A=-8-19-12100010B=100C=0410D=0所以，系统状态方程与输出方程分别为[𝑥̇1(𝑡)𝑥̇2(𝑡)𝑥̇3(𝑡)]=[−8−19−12100010][𝑥1(𝑡)𝑥2(𝑡)𝑥3(𝑡)]+[100]𝑓(𝑡)y(t)=[0410][𝑥1(𝑡)𝑥2(𝑡)𝑥3(𝑡)]ss2tf函数是将一个系统的状态空间表示法转化为系统函数，其语句格式为[num,den]=ss2tf(A,B,C,D,iu)其中，iu表示第i个输入，当只有一个输入时可忽略。ss2tf的调用返回值为H(s)或H(z)的分子多项式的系数矩阵num和分母多项式den，表示第iu个输入信号对输出y(t)的传递函数。【实例2】已知某离散时间系统的状态方程和输出方程分别为[𝑥1(𝑛+1)𝑥2(𝑛+1)]=[01−3−4][𝑥1(𝑛)𝑥2(𝑛)]+[02]𝑓(𝑛)y(t)=[−1−2][𝑥1(𝑛)𝑥2(𝑛)]+𝑓(𝑛)试用MATLAB命令求该离散时间系统的系统函数H(z)。解：已知系统状态空间矩阵，利用ss2tf函数可求出系统函数H(z)。由于是单输入，因此参数iu可省略。MATLAB源程序为：A=[01;-3-4];B=[0;2];C=[-12];D=1;[num,den]=ss2tf(A,B,C,D)num=1.00008.00001.0000den=143Hz=tf(num,den,-1)Transferfunction:z^2+8z+1------------------------z^2+4z+3即系统函数为H(z)=1+8𝑧−1+𝑧−21+4𝑧−1+3𝑧−2【实例3】一个多输入多输出系统，其状态方程和输出方程分别为[𝑥̇1(𝑡)𝑥̇2(𝑡)]=[01−2−3][𝑥1(𝑡)𝑥2(𝑡)]+[1011][𝑓1(𝑡)𝑓2(𝑡)][𝑦1(𝑡)𝑦2(𝑡)𝑦3(𝑡)]=[101102][𝑥1(𝑡)𝑥2(𝑡)]+[001001][𝑓1(𝑡)𝑓2(𝑡)]试用MATLAB命令求该系统的系统函数。解：由于是多输入多输出系统，所求系统函数是一个矩阵。因此，先分别用ss2tf函数求出对每个输入的传递函数，MATLAB源程序为：A=[0,1;-2,-3];B=[1,0;1,1];C=[1,0;1,1;0,2];D=[0,0;1,0;0,1];[num1,den1]=ss2tf(A,B,C,D,1)num1=01.00004.00001.00005.00004.000002.0000-4.0000den1=132[num2,den2]=ss2tf(A,B,C,D,2)num2=001.000001.00001.00001.00005.00002.0000den2=132所以，系统函数为H(s)=1𝑠2+3𝑠+2[𝑠+41𝑠2+5𝑠+4𝑠+12𝑠−4𝑠2+5𝑠+2]。2.状态方程的变换域符号求解法【实例4】已知连续系统的状态方程和输出方程分别为[𝑥̇1(𝑡)𝑥̇2(𝑡)]=[−1−41−1][𝑥1(𝑡)𝑥2(𝑡)]+[0110][𝑓1(𝑡)𝑓2(𝑡)][𝑦1(𝑡)𝑦2(𝑡)]=[110−1][𝑥1(𝑡)𝑥2(𝑡)]+[1010][𝑓1(𝑡)𝑓2(𝑡)]其初始状态和信号激励分别为[𝑥1(0−)𝑥2(0−)]=[21][𝑓1(𝑡)𝑓2(𝑡)]=[𝑢(𝑡)𝑒−𝑡𝑢(𝑡)]试用MATLAB命令求该系统的状态变量和输出变量。解：分别求出系统的状态变量与输出响应。MATLAB源程序为：symssA=[-1,-4;1,-1];B=[0,1;1,0];C=[1,1;0,-1];D=[1,0;1,0];r0=[2;1];X=[1/s;1/(s+1)];phis=inv(s*eye(2)-A);rs=phis*(r0+B*X);rs=simplify(rs)rs=((14*s)/5+3/5)/(s^2+2*s+5)-4/(5*s)(s^3+5*s^2+6*s+1)/(s*(s+1)*(s^2+2*s+5))rt=ilaplace(rs);rt=simplify(rt)rt=(14*exp(-t)*(cos(2*t)-(11*sin(2*t))/28))/5-4/5exp(-t)/4+(11*exp(-t)*(cos(2*t)+(28*sin(2*t))/11))/20+1/5ys=C*phis*r0+[C*phis*B+D]*X;ys=simplify(ys)ys=(4*s^3+9*s^2+8*s+2)/(s*(s+1)*(s^2+2*s+5))(-2*s^2+s+4)/(s*(s+1)*(s^2+2*s+5))yt=ilaplace(ys);yt=simplify(yt)yt=exp(-t)/4+(67*exp(-t)*(cos(2*t)+(6*sin(2*t))/67))/20+2/54/5-(11*exp(-t)*(cos(2*t)+(28*sin(2*t))/11))/20-exp(-t)/4用ilaplace求拉普拉斯反变换时假设了信号是因果的，因此，状态变量的结果暗含着与单位阶跃信号相乘，即[𝑥1(𝑡)𝑥2(𝑡)]=[(−45+145𝑒−𝑡cos(2𝑡)−1110𝑒−𝑡sin(2𝑡))𝑢(𝑡)(15+14𝑒−𝑡+1120𝑒−𝑡cos(2𝑡)+75𝑒−𝑡sin(2𝑡))𝑢(𝑡)]输出响应为[𝑦1(𝑡)𝑦2(𝑡)]=[(25+14𝑒−𝑡+6720𝑒−𝑡cos(2𝑡)+310𝑒−𝑡sin(2𝑡))𝑢(𝑡)(45−14𝑒−𝑡−1120𝑒−𝑡cos(2𝑡)−75𝑒−𝑡sin(2𝑡))𝑢(𝑡)]根据所得结果可画出输出响应波形，MATLAB源程序为t=0:0.01:4;y1=2/5+1/4*exp(-t)+67/20*exp(-t).*cos(2*t)+3/10*exp(-t).*sin(2*t);y2=4/5-1/4*exp(-t)-11/20*exp(-t).*cos(2*t)-5/7*exp(-t).*sin(2*t);subplot(211)plot(t,y1),gridonylabel('y1(t)'),xlabel('t')subplot(212)plot(t,y2),gridonylabel('y2(t)'),xlabel('t')根据程序运行结果，绘出0≤𝑡≤4范围内输出响应的波形，如图1-1所示图1-1连续系统的输出响应【实例5】给定系统状态方程、输出方程、激励信号和系统的初始条件分别为[𝑥1(𝑛+1)𝑥2(𝑛+1)]=[01−1656][𝑥1(𝑛)𝑥2(𝑛)]+[01]𝑓(𝑛)y(n)=[−15][𝑥1(𝑛)𝑥2(𝑛)]𝑓(𝑛)=𝑢(𝑛),[𝑥1(0)𝑥2(0)]=[23]试用MATLAB命令求输出响应y(n)。解：MATLAB源程序为symszA=[0,1;-1/6,5/6];B=[0;1];C=[-1,5];D=0;r0=[2;3];X=z/(z-1);phiz=inv(z*eye(2)-A);yz=C*phiz*z*r0+(C*phiz*B+D)*X;yz=simplify(yz)yz=(6*z*(13*z^2-11*z+2))/((2*z-1)*(3*z-1)*(z-1))yn=iztrans(yz)yn=3*(1/2)^n-2*(1/3)^n+12同样，用iztrans函数求z反变换时假设了信号时因果的，因此，求得的结果暗含着与单位阶跃序列相乘，即输出响应为y(n)=[3∗(12)𝑛−2∗(13)𝑛+12]𝑢(𝑛)根据所得结果可画出输出响应序列的波形，MATLAB源程序为n=0:15;yn=3*(1/2).^n-2*(1/3).^n+12;stem(n,yn),gridonxlabel('n'),ylabel('y(n)')axis([0,15,11,14])根据程序运行结果，绘出0≤n≤15范围内该离散系统输出响应序列的波形，如图1-2所示。图1-2离散系统的输出响应3.状态方程的时域符号求解法【实例6】试用MATLAB时域求解法求【实例4】中连续时间系统的状态变量、冲激响应和输出响应。解：先利用expm函数求出矩阵指数函数eAt，然后分别求状态变量和输出响应的零输入解和零状态解。值得注意的是，在MATLAB环境中，由于符号工具箱没有提供符号卷积函数，因此，根据x(t)=eAtx(0−)+eAtB∗f(t)y(t)=CeAtx(0−)+[CeAtB+Dδ(t)]∗f(t)的结果，求零状态解仍要采用变换域法来求解。例如，可以用两信号s域乘积的拉氏变换来求时域的卷积。MATLAB源程序为symstA=[-1-4;1-1];B=[01;10];C=[11;0-1];D=[10;10];x=[heaviside(t);exp(-t)*heaviside(t)];r0=[2;1];E=expm(t*A)E=[exp(t*(-1-2*i))/2+exp(t*(-1+2*i))/2,-exp(t*(-1-2*i))*i+exp(t*(-1+2*i))*i][(exp(t*(-1-2*i))*i)/4-(exp(t*(-1+2*i))*i)/4,exp(t*(-1-2*i))/2+exp(t*(-1+2*i))/2]rzi=E*r0;rzs=ilaplace(laplace(E*B)*laplace(x));rt=simplify(rzi+rzs)rt=14/5*exp(-t)*cos(2*t)-11/10*exp(-t)*sin(2*t)-4/57/5*exp(-t)*sin(2*t)+11/20*exp(-t)*cos(2*t)+1/5+1/4*exp(-t)ht=C*E*B+D*Dirac(t)ht=[-2*exp(-t)*sin(2*t)+exp(-t)*cos(2*t)+dirac(t),exp(-t)*cos(2*t)+1/2*exp(-t)*sin(2*t)][-exp