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

练习一实验一二. 熟悉简单的矩阵输入 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 1C = 4 5 6 7 8 9 1 2 34 AA = 1 2 3 4 5 6 7 8 9 BB = 9 8 7 6 5 4 3 2 1 CC = 4 5 6 7 8 9 1 2 3三. 基本序列运算1.A=1,2,3,B=4,5,6A = 1 2 3B = 4 5 6 C=A+BC = 5 7 9 D=A-BD = -3 -3 -3 E=A.*BE = 4 10 18 F=A./BF = 0.2500 0.4000 0.5000 G=A.BG = 1 32 729 stem(A) stem(B) stem(C) stem(D) stem(E) stem(F) stem(G)再举例: a=-1,-2,-3a = -1 -2 -3 b=-4,-5,-6b = -4 -5 -6 c=a+bc = -5 -7 -9 d=a-bd = 3 3 3 e=a.*be = 4 10 18 f=a./bf = 0.2500 0.4000 0.5000 g=a.bg =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,3b=2;-1;-1x=inv(A)*bA = 1 1 1 1 -2 1 1 2 3b = 2 -1 -1x = 3.0000 1.0000 -2.0000 4. A=2,3,-1;3,-2,1;1,2,1b=18;8;24x=inv(A)*bA = 2 3 -1 3 -2 1 1 2 1b = 18 8 24x = 4 6 8实验二二. 1. k=0:50x=sin(k);stem(x)xlabel(k);ylabel(sinX);title(sin(k)(k); 2. k=-25:1:25x=sin(k)+sin(pi*k);stem(k,x)xlabel(k);ylabel(f(k);title(sink+sink);3. k=3:50x=k.*sin(k);stem(k,x)xlabel(k);ylabel(f(k);title(ksink(k-3);4.%函数function y=f1(k)if k f1=1 1 1 1;f2=3 2 1;conv(f1,f2)ans = 3 5 6 6 3 13.函数定义： function r= pulse( k )if k0 r=0;else r=1;endend 运行代码for k=1:10f1(k)=pulse(k);f2(k)=(0.5k)*pulse(k);endconv(f1,f2)结果ans = Columns 1 through 100.5000 0.7500 0.8750 0.9375 0.9688 0.9844 0.9922 0.9961 0.9980 0.9990 Columns 11 through 200.9995 0.9998 0.9999 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 21 through 300.5000 0.2500 0.1250 0.0625 0.0312 0.0156 0.0078 0.0039 0.0020 0.0010 Columns 31 through 390.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.00004for i=1:10f1(i)=pulse(i);f2(i)=(-0.5)i)*pulse(i);endconv(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 390.0005 0.0002 0.0001 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000实验三2clear;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); endendF=F/10;subplot(2,1,1);plot(w,abs(F),b-);xlabel(w);ylabel(F);title(幅度频谱);gridsubplot(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.mfunction y=fun1(x)if(-pix) & (x0) y=pi+x;elseif (0x) & (xpi) y=pi-x;else y=0end%new.mclear allclcfor i=1:1000 g(i)=fun1(2/1000*i-1); w(i)=(i-1)*0.2*pi;endfor i=1001:10000 g(i)=0; w(i)=(i-1)*0.2*pi;endG=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相位频谱);2.%fun2.mfunction y=fun2(x)if x-1 y=cos(pi*x/2);else y=0;end%new2.mfor i=1:1000 g(i)=fun2(2/1000*i-1); w(i)=(i-1)*0.2*pi;endfor i=1001:10000 g(i)=0; w(i)=(i-1)*0.2*pi;endG=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.mfunction y=fun3(x)if x-1 y=1;elseif x0 & x 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-1p = -1k = 12. Ns=10Ds=1,-5,0sys2=tf(Ns,Ds)实验结果：Ns = 10Ds = 1 -5 0sys2 = 10 - s2 - 5 sz,p,k=tf2zp(10,1,-5,0)z = Empty matrix: 0-by-1p = 0 5k =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,1Ds=1,0,-1sys2=tf(Ns,Ds)实验结果：Ns = 1 1Ds = 1 0 -1sys2 = s + 1 - s2 - 1 Continuous-time transfer function.3 Ns=1,6,6,0;Ds=1,6,8;sys3=tf(Ns,Ds)实验结果：Ns = 1 6 6 0Ds = 1 6 8sys3 = s3 + 6 s2 + 6 s - s2 + 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);endmag = 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);endmag = 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);endmag = 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);endmag = 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-1p = 0k = 13. num=1,0,1;den=1,2,5;z,p,k=tf2zp(num,den)zplane(z,p)实验结果：z = 0 + 1.0000i 0 - 1.0000ip = -1.0000 + 2.0000i -1.0000 - 2.0000ik = 14. 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.1304ip = -2.3247 -0.3376 + 0.5623i -0.3376 - 0.5623ik =1.80005 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=-3P=-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.0000ip = 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点运动，同时其开环转移函数也是稳定的，由此，该线性负反馈系统也是稳定的。练习三实验三五1help windowWINDOW 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 window2.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=ABQ = 4-1 clearB=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 0b = 1 0c = 1 3d = 02.%冲激响应 clear;b=1,3;a=1,2,2;sys=tf(b,a);impulse(sys)%求零输入响应 A=1,3;1,-2;B=1;2;Q=ABQ = 1.6000 -0.2000 clearB=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 0b = 1 0c = 1 3d = 03.%冲激响应 clear;b=1,3;a=1,2,1;sys=tf(b,a);impulse(sys)%求零输入响应 A=1,3;1,-1;B=1;2;Q=ABQ = 1.7500 -0.2500 clearB=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 0b = 1 0c = 1 3d = 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); 练习三实验三五1help windowWINDOW 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 window2.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