东南大学信号与系统MATLAB实践第一次作业.doc

《东南大学信号与系统MATLAB实践第一次作业.doc》由会员分享，可在线阅读，更多相关《东南大学信号与系统MATLAB实践第一次作业.doc（81页珍藏版）》请在装配图网上搜索。

<信号与系统MATLAB实践> 练习一 实验一 二. 熟悉简单的矩阵输入 1.实验代码 >>A=[1,2,3;4,5,6;7,8,9] 实验结果 A = 1 2 3 4 5 6 7 8 9 3．实验代码 >>B=[9,8,7;6,5,4;3,2,1] C=[4,5,6;7,8,9;1,2,3] 实验结果： B = 9 8 7 6 5 4 3 2 1 C = 4 5 6 7 8 9 1 2 3 4．>> A A = 1 2 3 4 5 6 7 8 9 >> B B = 9 8 7 6 5 4 3 2 1 >> C C = 4 5 6 7 8 9 1 2 3 三. 基本序列运算 1.>>A=[1,2,3],B=[4,5,6] A = 1 2 3 B = 4 5 6 >> C=A+B C = 5 7 9 >> D=A-B D = -3 -3 -3 >> E=A.*B E = 4 10 18 >> F=A./B F = 0.2500 0.4000 0.5000 >> G=A.^B G = 1 32 729 >> stem(A) >> stem(B) >> stem(C) >> stem(D) >> stem(E) >> stem(F) >> stem(G) 再举例: >> a=[-1,-2,-3] a = -1 -2 -3 >> b=[-4,-5,-6] b = -4 -5 -6 >> c=a+b c = -5 -7 -9 >> d=a-b d = 3 3 3 >> e=a.*b e = 4 10 18 >> f=a./b f = 0.2500 0.4000 0.5000 >> g=a.^b g = 1.0000 -0.0313 0.0014 >> stem(a) >> stem(b) >> stem(c) >> stem(d) >> stem(e) >> stem(f) >> stem(g) 2. >>t=0:0.001:10 f=5*exp(-t)+3*exp(-2*t); plot(t,f) ylabel('f(t)'); xlabel('t'); title('(1)'); >> t=0:0.001:3; f=(sin(3*t))./(3*t); plot(t,f) ylabel('f(t)'); xlabel('t'); title('(2)'); >> k=0:1:4; f=exp(k); stem(f) 四. 利用MATLAB求解线性方程组 2. >>A=[1,1,1;1,-2,1;1,2,3] b=[2;-1;-1] x=inv(A)*b A = 1 1 1 1 -2 1 1 2 3 b = 2 -1 -1 x = 3.0000 1.0000 -2.0000 4. >> A=[2,3,-1;3,-2,1;1,2,1] b=[18;8;24] x=inv(A)*b A = 2 3 -1 3 -2 1 1 2 1 b = 18 8 24 x = 4 6 8 实验二 二. 1. >> k=0:50 x=sin(k); stem(x) xlabel('k'); ylabel('sinX'); title('sin(k)ε(k)'); 2. >> k=-25:1:25 x=sin(k)+sin(pi*k); stem(k,x) xlabel('k'); ylabel('f(k)'); title('sink+sinπk'); 3. >> k=3:50 x=k.*sin(k); stem(k,x) xlabel('k'); ylabel('f(k)'); title('ksinkε(k-3)'); 4. %函数 function y=f1(k) if k<0 y=(-1)^k; else y=(-1)^k+(0.5)^k; end %运行代码 for k=-10:1:10; y4(k+11)=f1(k); end k=-10:1:10; stem(k,y4); xlabel('k'); ylabel('f(k)'); title('4'); 七． 2．>> f1=[1 1 1 1]; f2=[3 2 1]; conv(f1,f2) ans = 3 5 6 6 3 1 3. 函数定义： function [r]= pulse( k ) if k<0 r=0; else r=1; end end 运行代码 for k=1:10 f1(k)=pulse(k); f2(k)=(0.5^k)*pulse(k); end conv(f1,f2) 结果 ans = Columns 1 through 10 0.5000 0.7500 0.8750 0.9375 0.9688 0.9844 0.9922 0.9961 0.9980 0.9990 Columns 11 through 20 0.9995 0.9998 0.9999 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 21 through 30 0.5000 0.2500 0.1250 0.0625 0.0312 0.0156 0.0078 0.0039 0.0020 0.0010 Columns 31 through 39 0.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 4 for i=1:10 f1(i)=pulse(i); f2(i)=((-0.5)^i)*pulse(i); end conv(f1,f2) 结果 ans = Columns 1 through 10 -0.5000 -0.2500 -0.3750 -0.3125 -0.3438 -0.3281 -0.3359 -0.3320 -0.3340 -0.3330 Columns 11 through 20 -0.3325 -0.3323 -0.3322 -0.3321 -0.3321 -0.3320 -0.3320 -0.3320 -0.3320 -0.3320 Columns 21 through 30 0.1680 -0.0820 0.0430 -0.0195 0.0117 -0.0039 0.0039 -0.0000 0.0020 0.0010 Columns 31 through 39 0.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 实验三 2． clear; x=[1,2,3,4,5,6,6,5,4,3,2,1]; N=0:11; w=-pi:0.01:pi; m=length(x); n=length(w); for i=1:n F(i)=0; for k=1:m F(i)=F(i)+x(k)*exp(-1j*w(i)*k); end end F=F/10; subplot(2,1,1); plot(w,abs(F),'b-');xlabel('w');ylabel('F');title('幅度频谱');grid subplot(2,1,2); plot(w,angle(F),'b-');xlabel('w'); X=fftshift(fft(x))/10; subplot(2,1,1); hold on; plot(N*2*pi/12-pi,abs(X),'r.'); legend('DIFT算法','DFT算法'); subplot(2,1,2);hold on; plot(N*2*pi/12-pi,angle(X),'r.'); xlabel('w');ylabel('相位');title('相位频谱');grid 三． 1. %fun1.m function y=fun1(x) if((-pi-1 y=cos(pi*x/2); else y=0; end %new2.m for i=1:1000 g(i)=fun2(2/1000*i-1); w(i)=(i-1)*0.2*pi; end for i=1001:10000 g(i)=0; w(i)=(i-1)*0.2*pi; end G=fft(g)/1000; subplot(1,2,1); plot(w(1:50),abs(G(1:50))); xlabel('w');ylabel('G');title('幅度频谱'); subplot(1,2,2); plot(w(1:50),angle(G(1:50))) xlabel('w');ylabel('Fi');title('相位频谱'); 3. %fun3.m function y=fun3(x) if x<0 && x>-1 y=1; elseif x>0 && x<1 y=-1; else y=0 end %new.m for i=1:1000 g(i)=fun3(2/1000*i-1); w(i)=(i-1)*0.2*pi; end for i=1001:10000 g(i)=0; w(i)=(i-1)*0.2*pi; end G=fft(g)/1000; subplot(1,2,1); plot(w(1:50),abs(G(1:50))); xlabel('w');ylabel('G');title('DFT幅度频谱'); subplot(1,2,2); plot(w(1:50),angle(G(1:50))) xlabel('w');ylabel('Fi');title('DFT相位频谱'); 练习二 实验六 一．用MATLAB语言描述下列系统，并求出极零点、 1. >> Ns=[1]; Ds=[1,1]; sys1=tf(Ns,Ds) 实验结果： sys1 = 1 ----- s + 1 >> [z,p,k]=tf2zp([1],[1,1]) z = Empty matrix: 0-by-1 p = -1 k = 1 2. >>Ns=[10] Ds=[1,-5,0] sys2=tf(Ns,Ds) 实验结果： Ns = 10 Ds = 1 -5 0 sys2 = 10 --------- s^2 - 5 s >>[z,p,k]=tf2zp([10],[1,-5,0]) z = Empty matrix: 0-by-1 p = 0 5 k = 10 二．已知系统的系统函数如下，用MATLAB描述下列系统。 1． >> z=[0]; p=[-1,-4]; k=1; sys1=zpk(z,p,k) 实验结果： sys1 = s ----------- (s+1) (s+4) Continuous-time zero/pole/gain model. 2. >> Ns=[1,1] Ds=[1,0,-1] sys2=tf(Ns,Ds) 实验结果： Ns = 1 1 Ds = 1 0 -1 sys2 = s + 1 ------- s^2 - 1 Continuous-time transfer function. 3． >> Ns=[1,6,6,0]; Ds=[1,6,8]; sys3=tf(Ns,Ds) 实验结果： Ns = 1 6 6 0 Ds = 1 6 8 sys3 = s^3 + 6 s^2 + 6 s ----------------- s^2 + 6 s + 8 Continuous-time transfer function. 六．已知下列H(s)或H(z)，请分别画出其直角坐标系下的频率特性曲线。 1. >> clear; for n = 1:400 w(n) = (n-1)*0.05; H(n) = (1j*w(n))/(1j*w(n)+1); end mag = abs(H); phase = angle(H); subplot(2,1,1) plot(w,mag);title('幅频特性') subplot(2,1,2) plot(w,phase);title('相频特性') 实验结果： 2. >> clear; for n = 1:400 w(n) = (n-1)*0.05; H(n) = (2*j*w(n))/((1j*w(n))^2+sqrt(2)*j*w(n)+1); end mag = abs(H); phase = angle(H); subplot(2,1,1) plot(w,mag);title('幅频特性') subplot(2,1,2) plot(w,phase);title('相频特性') 实验结果： 3. >>clear; for n = 1:400 w(n) = (n-1)*0.05; H(n) = (1j*w(n)+1)^2/((1j*w(n))^2+0.61); end mag = abs(H); phase = angle(H); subplot(2,1,1) plot(w,mag);title('幅频特性') subplot(2,1,2) plot(w,phase);title('相频特性') 实验结果： 4. >>clear; for n = 1:400 w(n) = (n-1)*0.05; H(n) =3*(1j*w(n)-1)*(1j*w(n)-2)/(1j*w(n)+1)*(1j*w(n)+2); end mag = abs(H); phase = angle(H); subplot(2,1,1) plot(w,mag);title('幅频特性') subplot(2,1,2) plot(w,phase);title('相频特性') 实验结果： 实验七 三．已知下列传递函数H(s)或H(z)，求其极零点，并画出极零图。 1. >> z=[1,2]'; p=[-1,-2]'; zplane(z,p) 实验结果： 2. >> z=[1,2]; p=[-1,-2]; zplane(z,p) >> num=[1]; den=[1,0]; [z,p,k]=tf2zp(num,den); zplane(z,p) >> num=[1]; den=[1,0]; [z,p,k]=tf2zp(num,den) zplane(z,p) 实验结果： z = Empty matrix: 0-by-1 p = 0 k = 1 3. >> num=[1,0,1]; den=[1,2,5]; [z,p,k]=tf2zp(num,den) zplane(z,p) 实验结果： z = 0 + 1.0000i 0 - 1.0000i p = -1.0000 + 2.0000i -1.0000 - 2.0000i k = 1 4. >> num=[1.8,1.2,1.2,3]; den=[1,3,2,1]; [z,p,k]=tf2zp(num,den) zplane(z,p) 实验结果： z = -1.2284 0.2809 + 1.1304i 0.2809 - 1.1304i p = -2.3247 -0.3376 + 0.5623i -0.3376 - 0.5623i k = 1.8000 5． >> clear; A=[0,1,0; 0,0,1; -6,-11,-6]; B=[0;0;1]; C=[4,5,1]; D=0; sys5=ss(A,B,C,D); pzmap(sys5) 实验结果： 五．求出下列系统的极零点，判断系统的稳定性。 1. >> clear; A=[5,2,1,0; 0,4,6,0; 0,-3,-6,-1;1,-2,-1,3]; B=[1;2;3;4]; C=[1,2,5,2]; D=0; sys=ss(A,B,C,D); [z,p,k]=ss2zp(A,B,C,D,1) pzmap(sys) 实验结果： z = 4.0280 + 1.2231i 4.0280 - 1.2231i 0.2298 p = -3.4949 4.4438 + 0.1975i 4.4438 - 0.1975i 0.6074 k = 28 由求得的极点，该系统不稳定。 4. >>z=[-3] P=[-1,-5,-15] 所以该系统为稳定的。 5. >>num=100*conv([1,0],conv([1,2],conv([1,2],conv([1,3,2],[1,3,2])))); den=conv([1,1],conv([1,-1],conv([1,3,5,2],conv([1,0,2,0,4],[1,0,2,0,4])))); [z,p,k]=tf2zp(num,den) 实验结果： z = 0 -2.0005 + 0.0005i -2.0005 - 0.0005i -1.9995 + 0.0005i -1.9995 - 0.0005i -1.0000 + 0.0000i -1.0000 - 0.0000i p = 1.0000 0.7071 + 1.2247i 0.7071 - 1.2247i 0.7071 + 1.2247i 0.7071 - 1.2247i -1.2267 + 1.4677i -1.2267 - 1.4677i -0.7071 + 1.2247i -0.7071 - 1.2247i -0.7071 + 1.2247i -0.7071 - 1.2247i -1.0000 -0.5466 >>zplane(z,p) 所以该系统不稳定。 七．已知反馈系统开环转移函数如下，试作其奈奎斯特图，并判断系统是否稳定。 1. >> b=[1]; a=[1,3,2]; sys=tf(b,a); nyquist(sys); 实验结果： 由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。 2． >> b=[1]; a=[1,4,4,0]; sys=tf(b,a); nyquist(sys); 实验结果： 由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。 3. >> b=[1]; a=[1,2,2]; sys=tf(b,a); nyquist(sys); 实验结果： 由于奈奎斯特图并未围绕上-1点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。 练习三 实验三 五． 1． >>help window WINDOW Window function gateway. WINDOW(@WNAME,N) returns an N-point window of type specified by the function handle @WNAME in a column vector. @WNAME can be any valid window function name, for example: @bartlett - Bartlett window. @barthannwin - Modified Bartlett-Hanning window. @blackman - Blackman window. @blackmanharris - Minimum 4-term Blackman-Harris window. @bohmanwin - Bohman window. @chebwin - Chebyshev window. @flattopwin - Flat Top window. @gausswin - Gaussian window. @hamming - Hamming window. @hann - Hann window. @kaiser - Kaiser window. @nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window. @parzenwin - Parzen (de la Valle-Poussin) window. @rectwin - Rectangular window. @tukeywin - Tukey window. @triang - Triangular window. WINDOW(@WNAME,N,OPT) designs the window with the optional input argument specified in OPT. To see what the optional input arguments are, see the help for the individual windows, for example, KAISER or CHEBWIN. WINDOW launches the Window Design & Analysis Tool (WinTool). EXAMPLE: N = 65; w = window(@blackmanharris,N); w1 = window(@hamming,N); w2 = window(@gausswin,N,2.5); plot(1:N,[w,w1,w2]); axis([1 N 0 1]); legend('Blackman-Harris','Hamming','Gaussian'); See also bartlett, barthannwin, blackman, blackmanharris, bohmanwin, chebwin, gausswin, hamming, hann, kaiser, nuttallwin, parzenwin, rectwin, triang, tukeywin, wintool. Overloaded functions or methods (ones with the same name in other directories) help fdesign/window.m Reference page in Help browser doc window 2. >>N = 128; w = window(@rectwin,N); w1 = window(@bartlett,N); w2 = window(@hamming,N); plot(1:N,[w,w1,w2]); axis([1 N 0 1]); legend('矩形窗','Bartlett','Hamming'); 3. >>wvtool(w,w1,w2) 六． ts=0.01; N=20; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(1):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=30; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(2):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=50; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(3):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=100; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(4):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=150; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(5):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; 实验八 1． %冲激响应 >> clear; b=[1,3]; a=[1,3,2]; sys=tf(b,a); impulse(sys); 结果： %求零输入响应 >> A=[1,3;0,-2]; B=[1;2]; Q=A\B Q = 4 -1 >> clear B=[1,3]; A=[1,3,2]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[4;-1]; initial(sys,x0); grid; a = -3 -2 1 0 b = 1 0 c = 1 3 d = 0 2. %冲激响应 >> clear; b=[1,3]; a=[1,2,2]; sys=tf(b,a); impulse(sys) %求零输入响应 >> A=[1,3;1,-2]; B=[1;2]; Q=A\B Q = 1.6000 -0.2000 >> clear B=[1,3]; A=[1,2,2]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.6;-0.2]; initial(sys,x0); grid; a = -2 -2 1 0 b = 1 0 c = 1 3 d = 0 3. %冲激响应 >> clear; b=[1,3]; a=[1,2,1]; sys=tf(b,a); impulse(sys) %求零输入响应 >> A=[1,3;1,-1]; B=[1;2]; Q=A\B Q = 1.7500 -0.2500 >> clear B=[1,3]; A=[1,2,1]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.75;-0.25]; initial(sys,x0); grid; a = -2 -1 1 0 b = 1 0 c = 1 3 d = 0 二． >> clear; b=1; a=[1,1,1,0]; sys=tf(b,a); subplot(2,1,1); impulse(sys);title('冲击响应'); subplot(2,1,2); step(sys);title('阶跃响应'); t=0:0.01:20; e=sin(t); r=lsim(sys,e,t); figure; subplot(2,1,1); plot(t,e);xlabel('Time');ylabel('A');title('激励信号'); subplot(2,1,2); plot(t,r);xlabel('Time');ylabel('A');title('响应信号'); 三． 1. >> clear; b=[1,3]; a=[1,3,2]; t=0:0.08:8; e=[exp(-3*t)]; sys=tf(b,a); lsim(sys,e,t); 2. >> clear; b=[1,3]; a=[1,2,2]; t=0:0.08:8; sys=tf(b,a); step(sys) 3． >> clear; b=[1,3]; a=[1,2,1]; t=0:0.08:8; e=[exp(-2*t)]; sys=tf(b,a); lsim(sys,e,t); Doc: 1. >> clear; B=[1]; A=[1,1,1]; sys=tf(B,A,-1); n=0:200; e=5+cos(0.2*pi*n)+2*sin(0.7*pi*n); r=lsim(sys,e); stem(n,r); 2. >> clear; B=[1,1,1]; A=[1,-0.5,-0.5]; sys=tf(B,A,-1); e=[1,zeros(1,100)]; n=0:100; r=lsim(sys,e); stem(n,r); 练习三 实验三 五． 1． >>help window WINDOW Window function gateway. WINDOW(@WNAME,N) returns an N-point window of type specified by the function handle @WNAME in a column vector. @WNAME can be any valid window function name, for example: @bartlett - Bartlett window. @barthannwin - Modified Bartlett-Hanning window. @blackman - Blackman window. @blackmanharris - Minimum 4-term Blackman-Harris window. @bohmanwin - Bohman window. @chebwin - Chebyshev window. @flattopwin - Flat Top window. @gausswin - Gaussian window. @hamming - Hamming window. @hann - Hann window. @kaiser - Kaiser window. @nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window. @parzenwin - Parzen (de la Valle-Poussin) window. @rectwin - Rectangular window. @tukeywin - Tukey window. @triang - Triangular window. WINDOW(@WNAME,N,OPT) designs the window with the optional input argument specified in OPT. To see what the optional input arguments are, see the help for the individual windows, for example, KAISER or CHEBWIN. WINDOW launches the Window Design & Analysis Tool (WinTool). EXAMPLE: N = 65; w = window(@blackmanharris,N); w1 = window(@hamming,N); w2 = window(@gausswin,N,2.5); plot(1:N,[w,w1,w2]); axis([1 N 0 1]); legend('Blackman-Harris','Hamming','Gaussian'); See also bartlett, barthannwin, blackman, blackmanharris, bohmanwin, chebwin, gausswin, hamming, hann, kaiser, nuttallwin, parzenwin, rectwin, triang, tukeywin, wintool. Overloaded functions or methods (ones with the same name in other directories) help fdesign/window.m Reference page in Help browser doc window 2. >>N = 128; w = window(@rectwin,N); w1 = window(@bartlett,N); w2 = window(@hamming,N); plot(1:N,[w,w1,w2]); axis([1 N 0 1]); legend('矩形窗','Bartlett','Hamming'); 3. >>wvtool(w,w1,w2) 六． ts=0.01; N=20; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(1):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=30; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(2):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=50; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(3):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=100; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(4):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; ts=0.01; N=150; t=0:ts:(N-1)*ts; x=2*sin(4*pi*t)+5*cos(6*pi*t); g=fft(x,N); y=abs(g)/100; figure(5):plot(0:2*pi/N:2*pi*(N-1)/N,y); grid; 实验八 1． %冲激响应 >> clear; b=[1,3]; a=[1,3,2]; sys=tf(b,a); impulse(sys); 结果： %求零输入响应 >> A=[1,3;0,-2]; B=[1;2]; Q=A\B Q = 4 -1 >> clear B=[1,3]; A=[1,3,2]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[4;-1]; initial(sys,x0); grid; a = -3 -2 1 0 b = 1 0 c = 1 3 d = 0 2. %冲激响应 >> clear; b=[1,3]; a=[1,2,2]; sys=tf(b,a); impulse(sys) %求零输入响应 >> A=[1,3;1,-2]; B=[1;2]; Q=A\B Q = 1.6000 -0.2000 >> clear B=[1,3]; A=[1,2,2]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.6;-0.2]; initial(sys,x0); grid; a = -2 -2 1 0 b = 1 0 c = 1 3 d = 0 3. %冲激响应 >> clear; b=[1,3]; a=[1,2,1]; sys=tf(b,a); impulse(sys) %求零输入响应 >> A=[1,3;1,-1]; B=[1;2]; Q=A\B Q = 1.7500 -0.2500 >> clear B=[1,3]; A=[1,2,1]; [a,b,c,d]=tf2ss(B,A) sys=ss(a,b,c,d); x0=[1.75;-0.25]; initial(sys,x0); grid; a = -2 -1 1 0 b = 1 0 c = 1 3 d = 0 二． >> clear; b=1; a=[1,1,1,0]; sys=tf(b,a); subplot(2,1,1); impulse(sys);title('冲击响应'); subplot(2,1,2); step(sys);title('阶跃响应'); t=0:0.01:20; e=sin(t); r=lsim(sys,e,t); figure; subplot(2,1,1); plot(t,e);xlabel('Time');ylabel('A');title('激励信号'); subplot(2,1,2); plot(t,r);xlabel('Time');ylabel('A');title('响应信号'); 三． 1. >> clear; b=[1,3]; a=[1,3,2]; t=0:0.08:8; e=[exp(-3*t)]; sys=tf(b,a); lsim(sys,e,t); 2. >> clear; b=[1,3]; a=[1,2,2]; t=0:0.08:8; sys=tf(b,a); step(sys) 3． >> clear; b=[1,3]; a=[1,2,1]; t=0:0.08:8; e=[exp(-2*t)]; sys=tf(b,a); lsim(sys,e,t); Doc: 1