信号与系统的matlab表示.docx
《信号与系统的matlab表示.docx》由会员分享，可在线阅读，更多相关《信号与系统的matlab表示.docx（31页珍藏版）》请在冰豆网上搜索。
信号与系统的matlab表示
连续信号与系统分析
一、典型信号的matlab表示
表示连续信号，需定义自变量的范围和取样间隔，如t=0:
0.01:
3
1.实指数信号y=k*exp（a*t）
2.正弦信号k*sin（w*t+phi）k*cos（w*t+phi）
3.复指数信号y=k*exp（（a+i*b）*t）
实部real（y）虚部imag（y）模abs（y）相角angle（y）共轭conj（y）
4.抽样信号Sat=sinc（t/pi）
5.矩形脉冲信号y=rectpuls（t,width）
周期方波信号y=square（2*pi*f*t,duty）%产生频率为fHZ，占空比为duty%的方波
6.三角脉冲信号
非周期三角波y=tripuls（t,width,skew）%斜度skew，最大幅度出现在t=（width/2）*skew
周期三角波y=sawtooth（t,width）
7.单位阶跃信号functiony=uCT（t）y=（t>=0）
阶跃信号符号函数Heaviside（）y=sym（‘Heaviside（t）’）%调用时必须用sym定义
冲激信号符号函数Dirac（）
二、Matlab的符号运算
1.定义符号变量
syms变量名symsx
sym（‘变量名’）x=sym（‘x’）
sym（‘表达式’）sym（‘x+1’）
2.化简符号运算结果simple或simplify
3.绘制符号表达式图形ezplot（y,[a,b]）
三、连续信号的运算
微分和积分运算（用符号表达式来表示）
1.微分运算
Diff（function,’variable’,n）%variable为求导变量，n为求导阶数
例：
symsaxy
y=sin（a*x^2）;
dy=diff（y,’x’）
2.积分运算
int（function,’variable’,a,b）%a为积分下限，b为积分上限
3.信号的反折fliplr（x）
4.卷积计算
1）符号运算计算卷积（求解积分的方法）
例：
symsTttao
xt1=exp（-t）;
xt2=exp（-t/T）;
xt_tao=subs（xt1,t,tao）*subs（xt2,t,t-tao）;
yt=int（xt_tao,tao,0,t）;
yt=simplify（yt）;
2）数值计算法求卷积
conv（）
y=dt*conv（e,h）
例：
求e（t）=u（t）-u（t-1）和h（t）=u（t）-u（t-1）的卷积
t0=-2;t1=4;dt=0.01;
t=t0:
dt:
t1;
e=u（t）-u（t-1）;
h=u（t）-u（t-1）;
y=dt*conv（e,h）;%Computetheconvolutionofx（t）andh（t）
subplot（221）
plot（t,e）,gridon,title（'Signale（t）'）,axis（[t0,t1,-0.2,1.2]）
subplot（222）
plot（t,h）,gridon,title（'Signalh（t）'）,axis（[t0,t1,-0.2,1.2]）
subplot（212）
t=2*t0:
dt:
2*t1;%thetimerangetotheconvolutionofeandh.
plot（t,y）,gridon,title（'Theconvolutionofx（t）andh（t）'）,axis（[2*t0,2*t1,-0.1,1.2]）,
xlabel（'Timetsec'）
四、连续LTI系统的时域分析
1.系统响应的符号求解dsolve（‘eq1,eq2,…’,’cond1,cond2,…’）;
%eqi表示微分方程，condi表示初始条件
例：
eq=’D3y+2*D2y+Dy=0’;
cond=’y（0）=1,Dy（0）=1,D2y（0）=2’;
yzi=dsolve（eq,cond）;%零输入响应
simplify（yzi）;
eq1=’D3y+4*D2y+8*Dy=3*Dx+8*x’;
eq2=’x=Heaviside（t）’;
cond=’y（-0.01）=0,Dy（-0.01）=0,D2y（-0.01）=0’;
yzs=dsolve（eq1,eq2,cond）;
simplify（yzs.y）;%零状态响应
2.零状态响应的数值求解
1）y=lsim（sys,f,t）
%sys表示系统模型，由sys=tf（b,a）生成的系统函数对象
%f输入信号向量，t时间抽样点向量
例：
ts=0;te=5;dt=0.01;
sys=tf（[6],[1,5,6]）;
t=ts:
dt:
te;
f=10*sin（2*pi*t）.*UT（t）;
y=lsim（sys,f,t）;
plot（t,y）,gridon;
xlabel（‘time’）,ylabel（‘y（t）’）;
title（‘零状态响应’）;
2）y=conv（f,impul）
3.连续系统冲激响应y=impulse（sys,t）%sys表示系统模型
4.连续系统阶跃响应y=step（sys,t）
五、信号的频域分析
1.傅立叶变换
1）符号运算求法
fourier（）和ifourier（）
例：
的傅立叶变换
ft=sym（‘exp（-2*t）*Heaviside（t）’）;
fw=fourier（ft）
ezplot（abs（fw））;%或者fw_conj=conj（fw）;Gw=sqrt（fw*fw_conj）;
phase=atan（image（fw）/real（fw））;%或者angle（fw）
ezplot（phase）
的傅立叶反变换
symst
fw=sym（‘1/（1+w^2’）;
ft=ifourier（fw,t）
2）数值计算求法
例：
求
的傅立叶变换
1）数值计算
dt=0.01;
t=-4:
dt:
4;
ft=（t+4）/2.*uCT（t+4）-t.*uCT（t）+（t-4）/2.*uCT（t-4）;
N=2000;
k=-N:
N;
W=pi*k/（N*dt）;
F=dt*ft*exp（-j*t'*W）;
F=abs（F）;
plot（W,F）,gridon;
axis（[-pipi-19]）;
title（'amplitudespectrum'）;
2）符号计算
ft=sym（'（t+4）/2*Heaviside（t+4）-t*Heaviside（t）+（t-4）/2*Heaviside（t-4）'）;
Fw=simplify（fourier（ft））;
ezplot（abs（Fw）,[-pipi]）;gridon;
2.系统的频率特性
1）[H,w]=freqs（b,a）：
连续系统频率响应的函数
2）波特图：
采用对数坐标的幅频特性和相频特性曲线,可显示频响间的微小差异
bode（sys）
例：
求
的频率特性
w=0:
0.01:
8*pi;
b=[1];
a=[11];
H=freqs（b,a,w）;
subplot（211）;
plot（w,abs（H））;
subplot（212）;
plot（w,angle（H））;
figure
（2）;
sys=tf（b,a）;
bode（sys）;
3.连续时间LTI系统的频域分析
例:
，求系统的响应。
W=-6*pi:
0.01:
6*pi;
B=[5];
A=[1,5];
H1=freqs（b,a,w）;
plot（w,abs（H1））;%系统幅频特性
Hw=sym（'5/（5+i*w）'）;
xt=sym（'Heaviside（t）-Heaviside（t-1）'）;
Xw=simplify（fourier（xt））;
figure;
ezplot（abs（Xw））;
Yw=Hw*Xw;%输出信号的傅立叶变换
yt=ifourier（Yw）;
figure;
ezplot（yt,[-0.2,2]）;
例：
求稳态响应
t=0:
0.1:
20;
w1=1;w2=10;
H1=1/（-w1^2+j*3*w1+2）;
H2=1/（-w2^2+j*3*w2+2）;
f=5*cos（t）+2*cos（10*t）;
y=abs（H1）*cos（w1*t+angle（H1））+abs（H2）*cos（w2*t+angle（H2））;
subplot（211）;
plot（t,f）,gridon;
title（'输入信号波形'）;
subplot（212）;
plot（t,y）,gridon;
title（'稳态响应波形'）;
[H,W]=freqs（[1],[132]）;
figure;
plot（W,abs（H））;
5.连续系统的零极点分析
1）求多项式的根
roots（）%求多项式的根b=[1-2];zs=roots（b）;
b=[1-2];
a=[145];
zs=roots（b）;
ps=roots（a）;
plot（real（zs）,imag（zs）,'blacko',real（ps）,imag（ps）,'blackx','markersize',12）;
axis（[-33-22]）;grid
legend（'zero','pole'）;
2）画零极点分布图
pzmap（sys）%sys=tf（b,a）;表示系统的模型
b=[1-2];
a=[145];
sys=tf（b,a）;
pzmap（sys）;
axis（[-33-22]）;
3）求解零极点的值
pole（sys）;
zero（sys）;
6.零极点分布与时域特性的关系pzmap（）和impulse（）
b=[1]
a=[101];
sys=tf（b,a）;
figure
（1）;
pzmap（sys）;
axis（[-22-22]）;
figure
（2）;
impulse（sys）;
7.拉普拉斯变换
1）正变换
L=laplace（f）%f为符号表达式
例：
symsat
L=laplace（exp（-t）*sin（a*t））;
2）反变换
符号运算
f=ilaplace（L）;%L符号表达式
例：
F=sym（‘s^2/（s^2+1）’）;
ft=ilaplace（F）;
部分分式展开
[r,p,k]=residue（b,a）;%p为极点，r为部分分式系数，k为整式部分的系数
例：
的反变换
b=[1-2];
a=conv（conv（[10],[11]）,conv（[11],[11]））;
[r,p,k]=residue（b,a）
实验一连续LTI系统的时频域分析
1.图1所示为一RLC串联电路，已知R=5W，L=1H，C=（1/6）F，
1）请用Matlab绘制出该系统的单位冲激响应和单位阶跃响应的波形；
2）当输入信号
时，请画出该系统的零状态响应波形图；
图1RLC串联电路
2.绘出图1所示系统的幅频响应和相频响应特性曲线，并讨论随着R的变化，当电阻R分别为4W、2W、0.8W、0.4W时，幅频响应的变化规律。
3.给出图1所示系统的零极点分布图，并讨论随着R的变化，当电阻R分别为4W、2W、0.8W、0.4W时，系统的稳定性的变化。
离散信号与系统分析
一、离散信号的表示
1.单位取样序列
functiony=impDT（n）
y=（n==0）;
2.单位阶跃序列
functiony=uDT（n）
y=（n>=0）;
二、离散时间系统的时域分析
1.系统的零状态响应
y=filter（b,a,x）;%x输入序列，a,b为差分方程的系数向量，y输出序列与x长度相同
2.单位取样响应
y=impz（b,a,N）;%N为取样响应的样值个数
y=stepz（b,a,N）;%N为阶跃响应的样值个数
例：
a=[3-42];
b=[12];
n=0:
30;
x=impDT（n）;
h=filter（b,a,x）;%方法1，单位取样序列的零状态响应
stem（n,h,'fill'）,gridon;
title（'impulseresponse'）;
figure
（2）;
impz（b,a,30）,gridon;%方法2
3.离散信号的卷积
y=conv（x,h）%求零状态响应
例：
三、z变换
1.z=ztrans（x）;%x,z为符号表达式
x=iztrans（z）
2.部分分式展开求z反变换
[R,P,K]=residuez（b,a）
3.系统的零极点
1）roots（）%求解多项式的根得到零极点
2）[z,p,k]=tf2zp（b,a）%直接求系统的零极点
3）zplane（b,a）;%画系统的零极点图
四、离散系统的频率响应
[H,w]=freqz（b,a,N）；%w为[0,pi]范围内的N个频率等分点
[H,w]=freqz（b,a,N,’whole’）；%w范围内为[0,2pi]
五、DTFT和DFT
1.序列傅里叶变换DTFT
序列傅里叶变换的Matlab实现：
n=n1:
n2;
M=input（‘putinthenumberM=’）;
k=0:
2*M-1;%观察两个周期
X=x*（exp（-j*2*pi/M））.^（n’*k）；%序列的傅里叶变换
2.离散傅里叶变换（DFT）
DFT变换的Matlab实现：
n=0:
M-1;
k=0:
N-1;
WN=exp（-j*2*pi/N）;
kn=n'*k;
WNkn=WN.^kn;
X=x*WNkn;
3.快速傅里叶变换（FFT）
离散傅里叶变换的快速算法
fft（x）：
利用快速算法计算x的M点DFT，其中M是x的长度。
fft（x,N）：
利用快速算法计算x的N点DFT，其中N是用户指定的长度。
分两种情况：
①若x的长度M>N，则将x截短为N点序列，再作N点DFT；
②若x的长度Mifft（X）：
利用快速算法计算X的M点IDFT，其中M是X的长度。
ifft（X,N）：
利用快速算法计算X的N点IDFT，其中N是用户指定的长度。
同样分两种情况，同fft（x,N）。
4.离散傅里叶级数
定义
，周期序列的傅里叶级数（DFS）变换对为：
Matlab实现：
WN=exp（-j*2*pi/N）;
kn=n'*k;
WNkn=WN.^kn;
X=x*WNkn;
例：
设
，要求用MATLAB实现：
（1）计算
的傅里叶变换
，并绘出其幅度谱；
（2）分别计算
的4点DFT和8点DFT，并绘出其幅度谱。
并说明它们和
的关系。
n=0:
3;
x=[1,1,1,1];
k=0:
200;
W=pi/100*k;
X=x*（exp（-j*pi/100））.^（n'*k）;%计算DTFT
magX=abs（X）;
subplot（3,1,1）;
plot（W/pi,magX）;
axis（[0,2,0,5]）;
x=[1,1,1,1];
X=fft（x）;
magX=abs（X）;
k=0:
3;
W=pi/2*k;
subplot（3,1,2）;
stem（W/pi,magX）;
axis（[0,2,0,5]）;
x=[1,1,1,1,0,0,0,0];
X=fft（x）;%计算FFT
magX=abs（X）;
k=0:
7;
W=pi/4*k;
subplot（3,1,3）
stem（W/pi,magX）
axis（[0,2,0,5]）
实验二离散LTI系统的时频域分析
1.某离散LSI系统的差分方程表示式为
满足初始条件
，
，求系统的单位响应，单位阶跃响应，用filter子函数求系统输入为
时的零输入、零状态及全响应。
提示：
通过解差分方程，可以得到全响应为
，使用filter子函数对系统差分方程进行求解，同时将求解结果与理论计算的结果进行比较。
2.已知离散系统的系统函数为
求该系统的零极点及零极点分布图，并判断系统的因果性和稳定性。
3.观察系统零极点的位置对幅频响应的影响。
已知一阶离散系统的系统函数为
，
（1）假设系统的零点在原点，极点分别取0.2、0.5、0.8，比较它们的幅频响应曲线，
（2）假设系统的极点在原点，零点分别取0.2、0.5、0.8，比较它们的幅频响应曲线，从中总结零极点位置对幅频响应的影响。
IIR滤波器设计
一、模拟滤波器设计
1、特沃斯模拟滤波器设计：
①[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）
其中，参数Wp和Ws分别是通带边界频率和阻带边界频率，Wp和Ws的单位是rad/s。
Rp和Rs分别为通带最大衰减和阻带最小衰减（dB）。
返回的参数N和Wn分别为滤波器的阶数和3dB截止频率。
对于带通和带阻滤波器，Wp和Ws都是二维向量，向量的第一个元素对应低端的边界频率，第二个元素对应高端的边界频率。
②
[B,A]=butter（N,Wn,'s'）其中，N和Wn分别为滤波器的阶数和3dB截止频率。
利用此函数可以获得低通和带通滤波器系统函数的分子多项式（B）和分母多项式（A）的系数。
[B,A]=butter（N,Wn,'high','s'）可以获得高通滤波器系统函数的分子多项式（B）和分母多项式（A）的系数。
[B,A]=butter（N,Wn,'stop','s'）可以获得带阻滤波器系统函数的分子多项式（B）和分母多项式（A）的系数。
[z,p,k]=buttap（N）:
设计一个N阶的归一化的巴特沃斯原型低通模拟滤波器，返回滤波器的零点、极点和增益，此时z为空。
利用freqs函数计算模拟滤波器的频率响应：
H=freqs（B,A,w）其中，B和A分别表示滤波器系统函数的分子多项式和分母多项式的系
数。
该函数返回矢量w指定的那些频率点上的频率响应，w的单位是rad/s。
不带输出变量的freqs函数，将绘制出幅频和相频曲线。
例：
设计一个巴特沃斯模拟低通滤波器，要求通带截止频率
，通带最大衰减
，阻带截止频率
，阻带最小衰减
。
要求绘出滤波器的幅频特性曲线。
参考程序如下：
fp=5000;
wp=2*pi*fp;
fs=12000;
ws=2*pi*fs;
rp=2;
rs=30;
[N,Wn]=buttord（wp,ws,rp,rs,'s'）;
[b,a]=butter（N,Wn,'s'）;
freqs（b,a）
H=20*log10（abs（freqs（b,a,w）））;%dB表示增益
2、计切比雪夫原型低通滤波器的函数，其调用格式如下：
[N,Wn]=cheb1ord（Wp,Ws,Rp,Rs,'s'）：
其中，参数Wp和Ws分别是通带边界频率和阻带边界频率，Wp和Ws的单位是rad/s。
Rp和Rs分别为通带最大衰减和阻带最小衰减（dB）。
返回切比雪夫I型滤波器的阶数N和通带截止频率Wn。
对于带通和带阻滤波器，Wp和Ws都是二维向量，向量的第一个元素对应低端的边界频率，第二个元素对应高端的边界频率。
[N,Wn]=cheb2ord（Wp,Ws,Rp,Rs,'s'）：
参数同cheb1ord，返回切比雪夫II型滤波器的阶数N和通带截止频率Wn。
[b,a]=cheby1（N,R,Wn）：
其中，N和Wn分别为滤波器的阶数和通带截止频率,R为纹波参数
。
利用此函数可以获得低通滤波器系统函数的分子多项式（b）和分母多项式（a）的系数。
Cheby2与cheby1调用格式相同。
[B,A]=cheby1（N,R,Wn,'high'）：
可以获得高通滤波器系统函数的分子多项式（b）和分母多项式（a）的系数。
Cheby2与cheby1调用格式相同。
[B,A]=cheby1（N,R,Wn,'stop'）:
可以获得带阻滤波器系统函数的分子多项式（b）和分母多项式（a）的系数。
Cheby2与cheby1调用格式相同。
[z,p,k]=cheb1ap（N,Rp）:
返回一个N阶的归一化的切比雪夫I型低通模拟滤波器的零点、极点和增益。
切比雪夫I型低通模拟滤波器在阻带是最平坦的。
[z,p,k]=cheb2ap（N,Rs）:
返回一个N阶的归一化的切比雪夫II型低通模拟滤波器的零点、极点和增益。
切比雪夫II型低通模拟滤波器在通带是最平坦的。
例：
设计切比雪夫滤波器，设计指标要求同例7-1。
参考程序如下：
fp=5000;
wp=2*pi*fp;
fs=12000;
ws=2*pi*fs;
rp=2;
rs=30;
[N,Wn]=cheb1ord（wp,ws,rp,rs,'s'）;
[z,p,k]=cheb1ap（N,rp）;
b=k*real（poly（z））;%poly（）用来求零点向量对应的特征多项式；
a=real（poly（p））;%求极点向量对应的特征多项式
[H,omega]=freqs（b,a）;
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;
subplot（221）,plot（omega*Wn/（2*pi）,abs（H））;grid
axis（[0,15000,0,1.1]）;ylabel（'幅度'）;xlabel（'f（Hz）'）;
subplot（222）,plot（omega*Wn/（2*pi）,angle（H））;grid
axis（[0,15000,-4,4]）;ylabel（'相位'）;xlabel（'f（Hz）'）;
subplot（212）,plot（omega*Wn/（2*pi）,dbH）;grid
axis（[0,14000,-40,2]）;ylabel（'幅度（dB）'）;xlabel（'f（Hz）'）;
3、椭圆原型低通滤波器的函数，其调用格式如下：
[N,Wn]=ellipord（Wp,Ws,Rp,Rs,'s'）：
其中，参数Wp和Ws分别是通带边界频率和阻带边界频率，Wp和Ws的单位是rad/s。
Rp和Rs分别为通带最大衰减和阻带最小衰减（dB）。
返回椭圆滤波器的阶数N和通带截止频率Wn。
对于带通和带阻滤波器，Wp和Ws都是二维向量，向量的第一个元素对应低端的边界频率，第二个元素对应高端的边界频率。
[b,a]=ellip（N,Rp,Rs,Wn）：
其中，N和Wn分别为滤波器的阶数和通带截止频率Wn,R为通带的纹波参数。
利用此函数可以获得低通滤波器系统函数的分子多项式（b）和分母多项式（a）的系数。
[B,A]=ellip（N,Rp,Rs,Wn,'high'）：
可以获得高通滤波器系统函数的分子多项式（b）和分母多项式（a）的系数。
[B,A]=ellip（N,
升级会员 