武汉科技大学Matlab演示程序.docx
《武汉科技大学Matlab演示程序.docx》由会员分享，可在线阅读，更多相关《武汉科技大学Matlab演示程序.docx（18页珍藏版）》请在冰点文库上搜索。
武汉科技大学Matlab演示程序
1单位冲激信号
％单位冲激信号
clear,closeall
t0=-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'）,gridon
axis（[-1,2,0,22]）
2单位阶跃信号
％单位阶跃信号
clear,closeall
t0=-1;te=3;dt=0.05;t1=0;
t=[t0:
dt:
te];st=length（t）;
k1=floor（（t1-t0）/dt）;
f2=[zeros（1,k1-1）,ones（1,st-k1+1）];
subplot（211）,stairs（t,f2,'m'）,gridon
axis（[-1,3,0,1.2]）
3系统的零输入响应
％已知连续系统微分方程及初始状态，求系统的输入响应
clear,closeall
a=[1,7,14,8]
te=7;dt=0.1;t1=0;
n=length（a）-1;
y0=[0,1,0];
p=roots（a）;V=rot90（vander（p））;c=V\y0';
t=[0:
dt:
te];
y=zeros（1,length（t））;
fork=1:
n
y=y+c（k）*exp（p（k）*t）;end
subplot（211）,plot（t,y,'m'）,gridon
4求微分方程的全解
y=dsolve（'D2y+3*Dy+2*y=2*t+2*t^2','y（0）=1,Dy（0）=1'）
y=
2-2*t+t^2-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*t^2'）%求非齐次通解
yt=
2-2*t+t^2+C1*exp（-2*t）+C2*exp（-t）
>>yp=yt-yht%求特解，即强迫响应
yp=
2-2*t+t^2
>>yh=y-yp%求齐次解，即自由响应
yh=
-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*t^2','y（0）=0,Dy（0）=0'）
yzs=
2-2*t+t^2-2*exp（-t）
7用符号画图函数ezplot（）画各种响应的波形
>>t=0:
0.01:
3;figure
（1）
>>ezplot（yzi,[0,3]）;holdon;ezplot（yzs,[0,3]）;ezplot（y,[0,3]）
>>axis（[0,3,-15]）,holdoff;
>>title（'全响应,零输入响应,零状态响应'）;figure
（2）
>>ezplot（yh,[0,3]）;holdon;ezplot（yp,[0,3]）;ezplot（y,[0,3]）
>>axis（[0,3,-15]）,holdoff;
>>title（'全响应,自由响应,强迫响应'）
8计算连续信号的卷积
%LT2_10.m
dt=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（'f1（t）'）;
axis（[t
（1）t（L）-0.21.2]）;grid,
subplot（3,1,2）,plot（t,f2,'linewidth',2）,ylabel（'f2（t）'）;
axis（[t
（1）t（L）-0.21.2]）;grid
subplot（3,1,3）,plot（tp,y,'linewidth',2）,ylabel（'y（t）'）;
axis（[t
（1）t（L）-0.21]）;grid
运行结果如下：
1方波分解
％方波分解为多次谐波之和
clear,closeall
t=0:
0.01:
2*pi;
y=zeros（10,max（size（t）））;x=zeros（size（t））;
fork=1:
2:
19x=x+sin（k*t）/k;y（（k+1）/2,:
）=x;end
subplot（211）,plot（t,y（1:
9,:
））,gridon;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/4
n_max=[61234];%最高谐波次数:
6,12,34
N=length（n_max）;%计算N次
t=-1.1:
.002:
1.1;
omega_0=2*pi;%基波频率
fork=1:
N
n=[];
n=[-n_max（k）:
n_max（k）];
L_n=length（n）;
F_n=zeros（1,L_n）;
fori=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（[-1.11.1-0.51.5]）;
line（[-1.11.1],[00],'color','r'）;%画直线,表示横轴,线为红色
line（[00],[-0.51.5],'color','r'）;%画直线,表示纵轴,线为红色
bt=strcat（'最高谐波次数=',num2str（n_max（k）））;%字符串连接
title（bt）;%在N幅图中的第k子图上写标题
end
程序运行结果显示在下图中。
2周期信号的频谱
％半波周期信号及其频谱
clear,closeall
T=0.01;w=200*pi;N=10;
t=linspace（-T/2,T/2）;dt=T/99;
f=abs（cos（w*t））;subplot（211）,plot（t,f,'m'）,gridon;
fork=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）;
end
subplot（212）,stem（0:
N,[A
（1）/2,A（2:
end）],'m'）,
gridon
3非周期信号的傅里叶变换
已知信号
，试用Matlab计算其傅里叶变换，并画出时间函数
、幅度频谱和相位频谱。
%LT3_5.m
t0=-2;t1=4;t=t0:
0.02:
t1;%定义时间范围
w0=-15;w1=15;w=w0:
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,'linewidth',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=
exp（-2*t）*Heaviside（t）
F=
1/（2+i*w）
程序画出的时间信号
、幅度频谱
及相位频谱
的曲线图如下图所示。
1拉普拉斯反变换
％拉普拉斯反变换
clear,closeall
a=[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））;
fori=1:
length（a）-1
f=f+c（i）*exp（p（i）*t）;
end
subplot（211）,plot（t,f,'m'）,
gridon
2零极点图
％零极点图
clear,closeall
a=[1,2,3,2];
b=[1,0,-1];
zs=roots（b）;ps=roots（a）;
subplot（211）,plot（real（zs）,imag（zs）,'go',real（ps）,imag（ps）,'mx','markersize',12）;
gridon;legend（'zero','pole'）;
1单位脉冲序列
％单位脉冲序列
clear,closeall
k0=-1;ke=6;ks=0;
k1=k0:
ke;
f1=[（k1-ks）==0];
subplot（211）,stem（k1,f1,'m'）,gridon
2单位阶跃序列
％单位阶跃序列
clear,closeall
k0=-1;ke=6;ks=0;
k1=k0:
ke;
f1=[（k1-ks）>=0];
subplot（211）,stem（k1,f1,'m'）,gridon
1离散时间傅里叶变换
％离散时间傅里叶变换
clear,closeall
dt=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'）;
gridon;
2离散傅里叶变换
％有限长序列f（k），设N＝8，16时，分别求该序列的离散傅里叶变换
clear,closeall
N=8;N1=16;K=4;
n=0:
N-1;k=0:
N1-1;
f1=[ones（1,K）,zeros（1,N1-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'）,grid
subplot（313）,stem（k,abs（Fk1）,'m'）;
text（7.2,3.5,'N=16'）,grid
1滤波器的频率响应
％已知FIR滤波器的系统函数，求频率响应
clear,closeall
a=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'）;grid
title（'FIR滤波器的幅频曲线'）;
xlabel（'w/pi'）,ylabel（'幅度'）;
subplot（212）,plot（w/pi,angle（H）,'g'）;grid
xlabel（'w/pi'）,ylabel（'相位（rad）'）;
2系统函数的零极点图
％绘制给定H（z）的零极点图
clear,closeall
A=[1,2,3,3,2,2];B=[1,0,3];
zplane（B,A）
升级会员 