武汉科技大学Matlab演示程序.docx

1、武汉科技大学Matlab演示程序1 单位冲激信号单位冲激信号clear,close allt0=-1;tf=3;dt=0.05;t1=0;t=t0:dt:tf;st=length(t);k1=floor(t1-t0)/dt);f1=zeros(1,st);f1(k1)=1/dt;subplot(211),stairs(t,f1,m),grid onaxis(-1,2,0,22)2 单位阶跃信号单位阶跃信号clear,close allt0=-1;te=3;dt=0.05;t1=0;t=t0:dt:te;st=length(t);k1=floor(t1-t0)/dt);f2=zeros(1,k1

2、-1),ones(1,st-k1+1);subplot(211),stairs(t,f2,m),grid onaxis(-1,3,0,1.2)3 系统的零输入响应已知连续系统微分方程及初始状态，求系统的输入响应clear,close alla=1,7,14,8te=7;dt=0.1;t1=0;n=length(a)-1;y0=0,1,0;p=roots(a);V=rot90(vander(p);c=Vy0;t=0:dt:te;y=zeros(1,length(t);for k=1:n y=y+c(k)*exp(p(k)*t);endsubplot(211),plot(t,y,m),grid o

3、n4求微分方程的全解y=dsolve(D2y+3*Dy+2*y=2*t+2*t2,y(0)=1,Dy(0)=1)y =2-2*t+t2-2*exp(-2*t)+exp(-t)5 求自由响应和强迫响应 yht=dsolve(D2y+3*Dy+2*y=0) % 求齐次通解yht =C1*exp(-2*t)+C2*exp(-t) yt=dsolve(D2y+3*Dy+2*y=2*t+2*t2) % 求非齐次通解yt =2-2*t+t2+C1*exp(-2*t)+C2*exp(-t) yp=yt-yht % 求特解，即强迫响应yp =2-2*t+t2 yh=y-yp % 求齐次解，即自由响应yh =-

4、2*exp(-2*t)+exp(-t)6 求零输入响应和零状态响应 yzi=dsolve(D2y+3*Dy+2*y=0,y(0)=1,Dy(0)=1)yzi =-2*exp(-2*t)+3*exp(-t) yzs=dsolve(D2y+3*Dy+2*y=2*t+2*t2,y(0)=0,Dy(0)=0)yzs =2-2*t+t2-2*exp(-t)7用符号画图函数ezplot( )画各种响应的波形t=0:0.01:3;figure(1)ezplot(yzi,0,3);hold on;ezplot(yzs,0,3);ezplot(y,0,3)axis(0,3,-1 5), hold off;tit

5、le(全响应,零输入响应,零状态响应); figure(2)ezplot(yh,0,3);hold on;ezplot(yp,0,3);ezplot(y,0,3)axis(0,3,-1 5) , hold off;title(全响应,自由响应,强迫响应)8计算连续信号的卷积% LT2_10.mdt=0.01;t=-1:dt:5;L=length(t);tp=2*t(1):dt:2*t(L);f1=rectpuls(t-0.5);f2=0.5*rectpuls(t-2,2);y=dt*conv(f1,f2);subplot(3,1,1),plot(t,f1,linewidth,2),ylabel

6、(f1(t);axis(t(1) t(L) -0.2 1.2);grid,subplot(3,1,2),plot(t,f2,linewidth,2),ylabel(f2(t);axis(t(1) t(L) -0.2 1.2);gridsubplot(3,1,3),plot(tp,y,linewidth,2),ylabel(y(t);axis(t(1) t(L) -0.2 1);grid运行结果如下：1 方波分解方波分解为多次谐波之和clear,close allt=0:0.01:2*pi;y=zeros(10,max(size(t);x=zeros(size(t);for k=1:2:19 x

7、=x+sin(k*t)/k;y(k+1)/2,:)=x;endsubplot(211),plot(t,y(1:9,:),grid on;line(0,pi+0.5,pi/4,pi);text(pi+0.5,pi/4,pi/4);axis(0,2*pi,-1,1)halft=ceil(length(t)/2);subplot(212),mesh(t(1:halft),1:10,y(:,1:halft)例3.1方波的傅里叶级数,最高谐波次数为6, 12和34的波形比较tau_T=3/4; % 占空比3/4n_max=6 12 34; % 最高谐波次数:6,12,34N=length(n_max);

8、 % 计算N次t=-1.1:.002:1.1;omega_0=2*pi; % 基波频率for k=1:N n=; n=-n_max(k):n_max(k); L_n=length(n); F_n=zeros(1,L_n); for i=1:L_n % 计算傅里叶复系数Fn F_n(i)=tau_T*Sa(tau_T*n(i)*pi)*exp(-j*tau_T*n(i)*pi); end F=F_n*exp(j*omega_0*n*t); % 计算前几项的部分和subplot(N,1,k),plot(t,real(F),linewidth,2); % 在N幅图中的第k子图画实部波形 axis(-

9、1.1 1.1 -0.5 1.5); line(-1.1 1.1,0 0,color,r); % 画直线,表示横轴,线为红色 line(0 0,-0.5 1.5,color,r); % 画直线,表示纵轴,线为红色 bt=strcat(最高谐波次数=,num2str(n_max(k); % 字符串连接 title(bt); % 在N幅图中的第k子图上写标题end程序运行结果显示在下图中。2 周期信号的频谱半波周期信号及其频谱clear,close allT=0.01;w=200*pi;N=10;t=linspace(-T/2,T/2);dt=T/99;f=abs(cos(w*t);subplot

10、(211),plot(t,f,m),grid on;for k=0:N a(k+1)=trapz(f.*cos(k*w*t)*dt/T*2; b(k+1)=trapz(f.*sin(k*w*t)*dt/T*2; A(k+1)=sqrt(a(k+1).2+b(k+1).2);endsubplot(212),stem(0:N,A(1)/2,A(2:end),m),grid on3 非周期信号的傅里叶变换已知信号，试用Matlab计算其傅里叶变换，并画出时间函数、幅度频谱和相位频谱。% LT3_5.mt0=-2;t1=4;t=t0:0.02:t1; % 定义时间范围w0=-15;w1=15;w=w0

11、:0.02:w1; % 定义频率范围f=sym(exp(-2*t)*Heaviside(t) % 定义符号函数f(t)F=fourier(f) % 求f(t)的傅里叶变换F=simple(F) % 化简F(jw)的表达式f1=subs(f,t,t); % 将t数组代入f(t)后用f1表示fmin=min(f1)-0.2;fmax=max(f1)+0.2; % 求f1的最大和最小值Fv=subs(F,w,w); % 将w数组代入F(jw)后用Fv表示F1=abs(Fv); % 求F(jw)的模P1=angle(Fv); % 求F(jw)的相角subplot(3,1,1),plot(t,f1,li

12、newidth,2); % 在第一幅图上画f(t)grid;ylabel(f(t);axis(t0,t1,fmin,fmax);Fmin=min(F1)-0.05;Fmax=max(F1)+0.05;subplot(3,1,2),plot(w,F1,linewidth,2,color,k);% 在第二幅图上画|F(jw)|grid;ylabel(|F(jw)|);axis(w0,w1,Fmin,Fmax);subplot(3,1,3),plot(w,P1*180/pi,linewidth,2,color,k); % 在第三幅图上画相位频谱grid;ylabel(相位(度);程序运行结果为f =

13、 exp(-2*t)*Heaviside(t)F = 1/(2+i*w)程序画出的时间信号、幅度频谱及相位频谱的曲线图如下图所示。1拉普拉斯反变换拉普拉斯反变换clear,close alla=1,9,26,24;b=1,2,1;dt=0.1;te=4;c,p=residue(b,a);t=0:dt:te;f=zeros(1,length(t);for i=1:length(a)-1 f=f+c(i)*exp(p(i)*t);endsubplot(211),plot(t,f,m),grid on2 零极点图零极点图clear,close alla=1,2,3,2;b=1,0,-1;zs=roo

14、ts(b);ps=roots(a);subplot(211),plot(real(zs),imag(zs),go,real(ps),imag(ps),mx,markersize,12);grid on;legend(zero,pole);1 单位脉冲序列单位脉冲序列clear,close allk0=-1;ke=6;ks=0;k1=k0:ke;f1=(k1-ks)=0;subplot(211),stem(k1,f1,m),grid on2 单位阶跃序列单位阶跃序列clear,close allk0=-1;ke=6;ks=0;k1=k0:ke;f1=(k1-ks)=0;subplot(211),

15、stem(k1,f1,m),grid on1 离散时间傅里叶变换离散时间傅里叶变换clear,close alldt=2*pi/8;w=linspace(-2*pi,2*pi,2000)/dt;k=-2:2;f=ones(1,5);subplot(211),stem(f);F=f*exp(-j*k*w);subplot(212),plot(w,abs(F),m);grid on;2 离散傅里叶变换有限长序列f(k)，设N8，16时，分别求该序列的离散傅里叶变换clear,close allN=8;N1=16;K=4;n=0:N-1;k=0:N1-1;f1=ones(1,K),zeros(1,N

16、1-K);Fk=fft(f1,N);Fk1=fft(f1,N1);subplot(311),stem(f1,m);subplot(312),stem(n,abs(Fk),m);text(3.2,3.3,N=8),gridsubplot(313),stem(k,abs(Fk1),m);text(7.2,3.5,N=16),grid1 滤波器的频率响应已知FIR滤波器的系统函数，求频率响应clear,close alla=1;b=0.2,0.6,0.4,1,0.4,0.6,0.2;H,w=freqz(b,a);subplot(211),plot(w/pi,abs(H),m);gridtitle(FIR滤波器的幅频曲线);xlabel(w/pi),ylabel(幅度);subplot(212),plot(w/pi,angle(H),g);gridxlabel(w/pi),ylabel(相位(rad);2 系统函数的零极点图绘制给定H（z）的零极点图clear,close allA=1,2,3,3,2,2;B=1,0,3;zplane(B,A)