東南大學信號與系統MATLAB實踐第一次作業(yè).doc

《東南大學信號與系統MATLAB實踐第一次作業(yè).doc》由會員分享，可在線閱讀，更多相關《東南大學信號與系統MATLAB實踐第一次作業(yè).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) 實驗結果： 五．求出下列系統的極零點，判斷系統的穩(wěn)定性。 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 由求得的極點，該系統不穩(wěn)定。 4. >>z=[-3] P=[-1,-5,-15] 所以該系統為穩(wěn)定的。 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) 所以該系統不穩(wěn)定。 七．已知反饋系統開環(huán)轉移函數如下，試作其奈奎斯特圖，并判斷系統是否穩(wěn)定。 1. >> b=[1]; a=[1,3,2]; sys=tf(b,a); nyquist(sys); 實驗結果： 由于奈奎斯特圖并未圍繞上-1點運動，同時其開環(huán)轉移函數也是穩(wěn)定的，由此，該線性負反饋系統也是穩(wěn)定的。 2． >> b=[1]; a=[1,4,4,0]; sys=tf(b,a); nyquist(sys); 實驗結果： 由于奈奎斯特圖并未圍繞上-1點運動，同時其開環(huán)轉移函數也是穩(wěn)定的，由此，該線性負反饋系統也是穩(wěn)定的。 3. >> b=[1]; a=[1,2,2]; sys=tf(b,a); nyquist(sys); 實驗結果： 由于奈奎斯特圖并未圍繞上-1點運動，同時其開環(huán)轉移函數也是穩(wěn)定的，由此，該線性負反饋系統也是穩(wěn)定的。 練習三 實驗三 五． 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