信号与系统的matlab表示.doc
- 文档编号:4871524
- 上传时间:2023-05-07
- 格式:DOC
- 页数:23
- 大小:365KB
信号与系统的matlab表示.doc
《信号与系统的matlab表示.doc》由会员分享,可在线阅读,更多相关《信号与系统的matlab表示.doc(23页珍藏版)》请在冰点文库上搜索。
连续信号与系统分析
一、典型信号的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)
23
实验一连续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的长度M ifft(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为通带的
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 信号 系统 matlab 表示