【MATLAB课件】第七章matlab在信号处理中的应用

1、1MATLAB语言2本章目标本章目标 熟悉熟悉MATLABMATLAB软件平台、工具箱、高效的数值软件平台、工具箱、高效的数值计算及符号计算功能；计算及符号计算功能； 熟悉熟悉MATLABMATLAB软件的信号处理编程方法和结果软件的信号处理编程方法和结果的可视化；的可视化； 了解数字信号处理的计算机仿真方法；了解数字信号处理的计算机仿真方法； 进一步加深对信号与系统的基本原理、方法进一步加深对信号与系统的基本原理、方法及应用的理解。及应用的理解。MATLAB语言3主要内容主要内容 7.1 7.1 基本信号的表示及可视化基本信号的表示及可视化 7.2 7.2 连续信号的时域运算与时域变换连续信

2、号的时域运算与时域变换 7.3 7.3 线性系统的时域分析线性系统的时域分析 7.4 7.4 连续时间信号的频域分析连续时间信号的频域分析 7.5 7.5 连续系统的复频域分析连续系统的复频域分析 7.6 7.6 信号采样与重构信号采样与重构MATLAB语言47.17.1基本信号的表示及可视化基本信号的表示及可视化 利用利用MatlabMatlab软件信号处理工具箱软件信号处理工具箱(Signal Processing (Signal Processing Toolbox)Toolbox)中的专用函数产生信号并绘出波形。中的专用函数产生信号并绘出波形。产生正弦波产生正弦波t=0:0.01:3t

3、=0:0.01:3* *pi;pi;y=sin(2y=sin(2* *t);t);plot(t,y)plot(t,y)产生矩形脉冲信号产生矩形脉冲信号t=-3:0.01:3;t=-3:0.01:3;y=rectpuls(t-1,2);%y=rectpuls(t-1,2);%以以t=1t=1为中心，宽度为为中心，宽度为2 2plot(t,y)plot(t,y)MATLAB语言5产生周期方波产生周期方波t=0:0.01:1;t=0:0.01:1;y=square(4y=square(4* *pipi* *t);t);plot(t,yplot(t,y) )产生周期锯齿波产生周期锯齿波t=0:0.00

4、1:2.5;t=0:0.001:2.5;y=sawtooth(2y=sawtooth(2* *pipi* *3030* *t);t);plot(t,yplot(t,y) )axis(0 0.2 -1 1)axis(0 0.2 -1 1)产生指数函数波形产生指数函数波形x=linspace(0,1);x=linspace(0,1);%x=0:0.01:1;%x=0:0.01:1;y=exp(-x);y=exp(-x);plot(x,yplot(x,y) )MATLAB语言6产生抽样信号产生抽样信号t=-10:1.5:10; t=-10:1.5:10; % %向量向量t t时间范围时间范围t=t1

5、:p:t2,pt=t1:p:t2,p为时间间隔为时间间隔f=f=sinc(tsinc(t); ); % %sinc(tsinc(t)=)=sin(pisin(pi* *t)/pit)/pi* *t tplot(t,fplot(t,f); %); %显示该信号的时域波形显示该信号的时域波形title(f(t)=Sa(t);title(f(t)=Sa(t);xlabel(txlabel(t)axis(-10,10,-0.4,1.1)axis(-10,10,-0.4,1.1)产生单位阶跃信号产生单位阶跃信号 定义阶跃函数定义阶跃函数 function f=function f=heaviside(t

6、heaviside(t) )f=(t0);f=(t0); 调用阶跃函数调用阶跃函数t=-1:0.01:3;t=-1:0.01:3;f=f=heaviside(theaviside(t););plot(t,fplot(t,f););axis(-1,3,-0.2,1.2);axis(-1,3,-0.2,1.2); 1000ttt MATLAB语言7产生单位冲激信号产生单位冲激信号 定义冲激函数定义冲激函数function chongji(t1,t2,t0)function chongji(t1,t2,t0)dt=0.01;t=t1:dt:t2;dt=0.01;t=t1:dt:t2;n=length

7、(t);n=length(t);x=zeros(1,n);x=zeros(1,n);x(1,(-t0-t1)/dt+1)=1/dt;%x(1,(-t0-t1)/dt+1)=1/dt;%表示在表示在t=t+t0t=t+t0时刻存在冲激时刻存在冲激stairs(t,xstairs(t,x););axis(t1,t2,0,1.2/dt);title(axis(t1,t2,0,1.2/dt);title(单位冲击信号单位冲击信号(t(t)调用调用chongji(-1,5,0)chongji(-1,5,0)；可以试着给别的可以试着给别的t1,t2,t0.t1,t2,t0. 10,0t dtttMATLA

8、B语言87.2 7.2 连续信号的时域运算与时域变换连续信号的时域运算与时域变换 1.1.加加( (减减) )、乘运算、乘运算要求二个信号序列长度相同要求二个信号序列长度相同. .t=0:0.01:2;t=0:0.01:2;f1=exp(-3f1=exp(-3* *t);t);f2=0.2f2=0.2* *sin(4sin(4* *pipi* *t);t);f3=f1+f2;f3=f1+f2;f4=f1.f4=f1.* *f2;f2;subplot(2,2,1);plot(t,f1);title(f1(t);subplot(2,2,1);plot(t,f1);title(f1(t);subpl

9、ot(2,2,2);plot(t,f2);title(f2(t);subplot(2,2,2);plot(t,f2);title(f2(t);subplot(2,2,3);plot(t,f3);title(f1+f2);subplot(2,2,3);plot(t,f3);title(f1+f2);subplot(2,2,4);plot(t,f4);title(f1subplot(2,2,4);plot(t,f4);title(f1* *f2);f2);MATLAB语言92.2.用用matlabmatlab的符号函数实现信号的反折、移位、的符号函数实现信号的反折、移位、尺度变换尺度变换. .由由

10、f(tf(t) )到到f(-at+b)(af(-at+b)(a0)0)步骤步骤: : f(t)f(tb)f(atb)f( atb) 移位尺度反折MATLAB语言10例例: :已知已知f(tf(t)=)=sin(t)/tsin(t)/t, ,试通过反折、移位、尺度变换由试通过反折、移位、尺度变换由f(tf(t) )得到得到f(-2t+3) f(-2t+3) 的波形的波形. .symssyms t; t;f=f=sin(t)/tsin(t)/t; %; %定义符号函数定义符号函数f(tf(t)=)=sin(t)/tsin(t)/tf1=subs(f,t,t+3);f1=subs(f,t,t+3);

11、 % %对对f f进行移位进行移位f2=subs(f1,t,2f2=subs(f1,t,2* *t);t); % %对对f1f1进行尺度变换进行尺度变换f3=subs(f2,t,-t);f3=subs(f2,t,-t); % %对对f2f2进行反褶进行反褶subplot(2,2,1);ezplot(f,-8,8);grid on;subplot(2,2,1);ezplot(f,-8,8);grid on;% % ezplotezplot是符号函数绘图命令是符号函数绘图命令subplot(2,2,2);ezplot(f1,-8,8);grid on;subplot(2,2,2);ezplot(f

12、1,-8,8);grid on;subplot(2,2,3);ezplot(f2,-8,8);grid on;subplot(2,2,3);ezplot(f2,-8,8);grid on;subplot(2,2,4);ezplot(f3,-8,8);grid on;subplot(2,2,4);ezplot(f3,-8,8);grid on;注注: :也可用一条指令也可用一条指令:subs(f,t,-2:subs(f,t,-2* *t+3)t+3)实现实现f(tf(t) )到到f(-2t+3)f(-2t+3)的变换的变换MATLAB语言11例：已知信号例：已知信号用用matlabmatlab求

13、求f(t+2),f(t-2),f(-t),f(2t),-f(t),f(t+2),f(t-2),f(-t),f(2t),-f(t),并绘出时域波并绘出时域波形。形。symssyms t tf=sym(t/2+1)f=sym(t/2+1)* *(heaviside(t+2)-heaviside(t-2);(heaviside(t+2)-heaviside(t-2);subplot(2,3,1);ezplot(f,-3,3);subplot(2,3,1);ezplot(f,-3,3);y1=subs(f,t,t+2);y1=subs(f,t,t+2); subplot(2,3,2),ezplot(y

14、1,-5,1); subplot(2,3,2),ezplot(y1,-5,1);y2=subs(f,t,t-2);y2=subs(f,t,t-2); subplot(2,3,3),ezplot(y2,-1,5); subplot(2,3,3),ezplot(y2,-1,5);y3=y3=subs(f,t,-tsubs(f,t,-t);); subplot(2,3,4),ezplot(y3,-3,3); subplot(2,3,4),ezplot(y3,-3,3);y4=subs(f,t,2y4=subs(f,t,2* *t);t); subplot(2,3,5),ezplot(y4,-2,2)

15、; subplot(2,3,5),ezplot(y4,-2,2);y5=-f;y5=-f; subplot(2,3,6),ezplot(y5,-3,3) subplot(2,3,6),ezplot(y5,-3,3) 1222tf tttMATLAB语言123. 3. 卷积运算卷积运算Y=Y=conv(x,hconv(x,h) )实现实现x,hx,h二个序列的卷积二个序列的卷积. .例例: :求二个信号的卷积求二个信号的卷积. .t11=0; t12=1;t11=0; t12=1;t21=0; t22=2;t21=0; t22=2;t1=t11:0.001:t12;ft1=2.t1=t11:0.

16、001:t12;ft1=2.* *rectpuls(rectpuls(t1-0.5,1);t1-0.5,1);subplot(2,2,1);plot(t1,ft1);axis(0 3 0 4)subplot(2,2,1);plot(t1,ft1);axis(0 3 0 4)t2=t21:0.001:t22;ft2=t2;t2=t21:0.001:t22;ft2=t2;subplot(2,2,2);plot(t2,ft2);axis(0 3 0 4)subplot(2,2,2);plot(t2,ft2);axis(0 3 0 4)t3=t11+t21:0.001:t12+t22;t3=t11+t

17、21:0.001:t12+t22;ft3=conv(ft1,ft2)ft3=conv(ft1,ft2)ft3=ft3ft3=ft3* *0.0010.001subplot(2,2,3);plot(t3,ft3);axis(0 3 0 4)subplot(2,2,3);plot(t3,ft3);axis(0 3 0 4) 12,01( )0,tf telse 2,020,ttftelse MATLAB语言137.37.3线性系统的时域分析线性系统的时域分析1. 1. 脉冲响应脉冲响应impulseimpulse(num,den,T)(num,den,T)y=y=impulseimpulse(nu

18、m,den,T)(num,den,T)T:T:为等间隔的时间向量为等间隔的时间向量, ,指明要计算响应的时间点指明要计算响应的时间点. .2. 2. 阶跃响应阶跃响应stepstep(num,den,T)(num,den,T)y=step(y=step(num,den,T)num,den,T)T:T:为等间隔的时间向量为等间隔的时间向量, ,指明要计算响应的时间点指明要计算响应的时间点. .3)3)对任意输入的响应对任意输入的响应lsimlsim(num,den,U,T)(num,den,U,T)y=y=lsimlsim(num,den,U,T)(num,den,U,T)U:U:任意输入信号任

19、意输入信号. . T:T:为等间隔的时间向量为等间隔的时间向量, ,指明要计算响应的时间点指明要计算响应的时间点. .MATLAB语言14例：已知描述某连续系统的微分方程为例：已知描述某连续系统的微分方程为：试用试用MatMatlablab绘出该系统冲激响应和阶跃响应；输入信号绘出该系统冲激响应和阶跃响应；输入信号为为 , ,该系统零状态响应该系统零状态响应y y（t t）。b=1 2;b=1 2;a=1 2 1;a=1 2 1;subplot(1,3,1);subplot(1,3,1);impulseimpulse(b,a); %(b,a); %冲激响应冲激响应subplot(1,3,2);

20、subplot(1,3,2);stepstep(b,a) %(b,a) %阶跃响应阶跃响应p=0.5; %p=0.5; %定义取样时间间隔定义取样时间间隔t=0:p:5; % t=0:p:5; % 定义时间范围定义时间范围x=exp(-2x=exp(-2* *t);t); % %定义输入信号定义输入信号subplot(1,3,3);subplot(1,3,3);lsim(b,a,x,tlsim(b,a,x,t);); % %对系统的输出信号进行仿真对系统的输出信号进行仿真 22yty ty tftf t 2tf tetMATLAB语言15例例: :对系统对系统 分别求脉冲响应、阶跃响应及对输入

21、分别求脉冲响应、阶跃响应及对输入u(tu(t)=)=sin(tsin(t) )的响应的响应. .num=1,1;num=1,1;den=1,1.3,0.8;den=1,1.3,0.8;T=0:0.1:3;T=0:0.1:3;y1=y1=impulseimpulse(num,den,T);(num,den,T);y2=y2=stepstep(num,den,T);(num,den,T);U=sin(T);U=sin(T);y3=y3=lsimlsim(num,den,U,T);(num,den,U,T);subplot(1,3,1);plot(T,y1);title(subplot(1,3,1)

22、;plot(T,y1);title(脉冲响应脉冲响应)subplot(1,3,2);plot(T,y2);title(subplot(1,3,2);plot(T,y2);title(阶跃响应阶跃响应)subplot(1,3,3);plot(T,y3);title(subplot(1,3,3);plot(T,y3);title(输入为输入为u=sintu=sint的响应的响应) 21( )1.30.8sHsssMATLAB语言167.4 7.4 连续时间信号的频域分析连续时间信号的频域分析1 1连续信号的傅立叶变换连续信号的傅立叶变换信号的傅立叶变换定义为：信号的傅立叶变换定义为： ()( )j

23、 tF jf t edt ( )()j tf tF jedMATLAB语言17（1 1）符号解法）符号解法调用函数：调用函数：fourierfourier、ifourierifourier例：求函数例：求函数 的傅立叶变换。的傅立叶变换。symssyms t tF=fourier(exp(-2F=fourier(exp(-2* *abs(tabs(t)F=4/(4+w2)F=4/(4+w2)例：求例：求傅立叶逆变换傅立叶逆变换syms t wsyms t wf =ifourier(1/( 1+ w2),t)f =ifourier(1/( 1+ w2),t)f =f = 1/2 1/2* *ex

24、p(-exp(-t)t)* *Heaviside(t)+1/2Heaviside(t)+1/2* *exp(t)exp(t)* *Heaviside(-t)Heaviside(-t) 2( )tf te 21()1F jMATLAB语言18例：设例：设 ，画出，画出f(tf(t) )及其幅频谱。及其幅频谱。syms t;syms t;f=1/2f=1/2* *exp(-2exp(-2* *t)t)* *sym(Heaviside(t);sym(Heaviside(t);F=F=fourierfourier(f(f););subplot(211);ezplot(f);subplot(211);e

25、zplot(f);subplot(212);ezplot(abs(F)subplot(212);ezplot(abs(F)符号解法的局限符号解法的局限： 运算结果是符号表达式，必须使用运算结果是符号表达式，必须使用ezplotezplot来作图。如果来作图。如果结果中含有冲激函数等项，那么结果中含有冲激函数等项，那么ezplotezplot无法作图。无法作图。 返回结果有可能包含一些不能直接表达的式子，甚至会返回结果有可能包含一些不能直接表达的式子，甚至会出现出现“未被定义函数或变量未被定义函数或变量”项，更加无法作图。项，更加无法作图。 很多场合，尽管连续，但却不可能表示成符号表达式。很多场

26、合，尽管连续，但却不可能表示成符号表达式。 21( )( )2tf tet( )f tMATLAB语言19（2 2）数值解法）数值解法从定义出发：从定义出发： 用用MATLABMATLAB作数值计算，它不能计算无限区间。可以根作数值计算，它不能计算无限区间。可以根据波形情况，确定积分上下限，将据波形情况，确定积分上下限，将t t分分N N等份，用求和等份，用求和代替积分代替积分。 0()( )lim()j tj nnF jf t edtf n e 12112()( ) ( ),( ),. (),.iNNj tiij tj tj tNF jf t ef tf tf teeeMATLAB语言20

27、同样，在频域内也是取一系列的样本点，然后计算出不同同样，在频域内也是取一系列的样本点，然后计算出不同w w处处的的F F值，然后就得到了傅立叶变换的数值解。而不同值，然后就得到了傅立叶变换的数值解。而不同w w处的处的F F值值都用上面同一个公式计算，因此可以使用都用上面同一个公式计算，因此可以使用MATLABMATLAB的矩阵计算功的矩阵计算功能。将能。将w w设为一个行向量，则设为一个行向量，则: : 其中，其中，F F是与是与w w等长的行向量，等长的行向量，tt是列向量，是列向量，w w是行向量，是行向量，tt* *w w是一个矩阵，其行数与是一个矩阵，其行数与t t相同，列数与相同，

28、列数与w w相同。相同。 *exp(* * )*Ffj twdtMATLAB语言21例：求门信号例：求门信号的傅立叶变换。的傅立叶变换。dtdt=0.02; =0.02; taotao=2;=2;t=-t=-tao:dt:taotao:dt:tao; ;f=Heaviside(t+1)-Heaviside(t-1);f=Heaviside(t+1)-Heaviside(t-1);w1=10w1=10* *pi;N=500;k=0:N;pi;N=500;k=0:N;w=kw=k* *w1/N;w1/N;F=fF=f* *exp(-jexp(-j* *tt* *w)w)* *dt;dt;F=F=r

29、eal(Freal(F););w=-fliplr(w),w(2:501); %w=-fliplr(w),w(2:501); %双边谱双边谱F=fliplr(F),F(2:501);F=fliplr(F),F(2:501);subplot(211),plot(t,f)subplot(211),plot(t,f)subplot(212),plot(w,F)subplot(212),plot(w,F) 211( )( )01tf tg tt MATLAB语言222.2.频响特性频响特性系统在正弦激励下稳态响应随信号频率变化的特性系统在正弦激励下稳态响应随信号频率变化的特性. .| |H(jH(j )

30、|:)|:幅频响应特性幅频响应特性. . ( ( ):):相频响应特性相频响应特性( (或相移特性或相移特性).).MatlabMatlab求系统频响特性函数求系统频响特性函数freqsfreqs的调用格式的调用格式: :h=h=freqs(num,denfreqs(num,den, , ) ) 为等间隔的角频率向量为等间隔的角频率向量, ,指明要计算响应的频率点指明要计算响应的频率点. . )()()()(jjsejHsHjHMATLAB语言23例例: :求系统求系统的频响特性的频响特性. . num=1,1; num=1,1;den=1,1.3,0.8;den=1,1.3,0.8;W=0:

31、0.1:100;W=0:0.1:100;h=freqs(num,den,W);h=freqs(num,den,W);subplot(2,1,1);plot(W,abs(h);title(subplot(2,1,1);plot(W,abs(h);title(幅频特性幅频特性)axis(0,20,0,1.5);axis(0,20,0,1.5);set(gca,xtick,0,10,20);set(gca,xtick,0,10,20);set(gca,ytick,0,1/sqrt(2),1.25);grid on;set(gca,ytick,0,1/sqrt(2),1.25);grid on;sub

32、plot(2,1,2);plot(W,angle(h);title(subplot(2,1,2);plot(W,angle(h);title(相频特性相频特性)axis(0,20,-pi/2,0.2);axis(0,20,-pi/2,0.2);set(gca,xtick,0,10,20);set(gca,xtick,0,10,20);set(gca,ytick,-pi/2,-pi/4,0);grid on;set(gca,ytick,-pi/2,-pi/4,0);grid on; 21( )1.30.8sH sssMATLAB语言247.5 7.5 连续系统的复频域分析连续系统的复频域分析1.

33、 1. 系统函数系统函数H(sH(s) )定义定义系统零状态响应的拉氏变换与激励的拉氏变换之比系统零状态响应的拉氏变换与激励的拉氏变换之比. .H(sH(s)=)=R(s)/E(sR(s)/E(s) )在在matlabmatlab中中, , 传递函数描述法是通过传递函数分子和分母关传递函数描述法是通过传递函数分子和分母关于于s s降幂排列的多项式系数来表示的降幂排列的多项式系数来表示的. .例如例如, ,某系统传递函数：某系统传递函数：则可用如下二个向量则可用如下二个向量numnum和和denden来表示来表示: :num=1,1num=1,1den=1,1.3,0.8den=1,1.3,0.

34、8 21( )1.30.8sH sssMATLAB语言252. 2. 系统零、极点分布与系统稳定性关系系统零、极点分布与系统稳定性关系 系统函数系统函数H(sH(s) )集中表现了系统的性能集中表现了系统的性能, ,研究研究H(sH(s) )在在S S平面中极平面中极点分布的位置点分布的位置, ,可很方面地判断系统稳定性可很方面地判断系统稳定性. .u稳定系统稳定系统: : H(sH(s) )全部极点落于全部极点落于S S左半平面左半平面( (不包括虚轴不包括虚轴),),即若即若满足，则系统是稳定的满足，则系统是稳定的. .u不稳定系统不稳定系统: : H(sH(s) )极点落于极点落于S S

35、右半平面右半平面, ,或在虚轴上具有二阶或在虚轴上具有二阶以上极点以上极点, ,则在足够长时间后则在足够长时间后, ,h(th(t) )仍继续增长仍继续增长, , 系统是不稳定系统是不稳定u临界稳定系统临界稳定系统: : H(sH(s) )极点落于极点落于S S平面虚轴上平面虚轴上, ,且只有一阶且只有一阶, ,则在则在足够长时间后足够长时间后, ,h(th(t) )趋于一个非零数值或形成一个等幅振荡趋于一个非零数值或形成一个等幅振荡. . H(sH(s) )零、极点零、极点可用可用matlabmatlab的多项式求根函数的多项式求根函数rootsroots求得：求得：极点极点: :p=p=r

36、oots(denroots(den) )零点零点: :z=z=roots(numroots(num) )根据根据p p和和z z，用，用plotplot命令命令即可画出系统零极点分布图即可画出系统零极点分布图, ,进而分析判进而分析判断系统稳定性。断系统稳定性。MATLAB语言26例例: : 系统函数系统函数画出系统零、极点分布图画出系统零、极点分布图, , 判断该系统稳定性。判断该系统稳定性。num=1,0,-4;num=1,0,-4;den=1,2,-3,2,1;den=1,2,-3,2,1;p=p=roots(denroots(den) )z=z=roots(numroots(num)

37、)plot(real(p),imag(pplot(real(p),imag(p),),* *);hold on;);hold on;plot(real(z),imag(z),o);gridplot(real(z),imag(z),o);grid on on 由系统零极点分布图可知由系统零极点分布图可知, ,该系统有一对极点位于该系统有一对极点位于s s右半右半平面平面, ,故系统是不稳定的故系统是不稳定的. . 24324( )2321sH sssssMATLAB语言277.6 7.6 信号采样与重构信号采样与重构Nyquist 采样定理： 对于带限信号 f t ,如果其带宽为mw , 那么以

38、大于 2mw 的频率对其进行采样，得到采样信号 sft，则该信号 f t 就可由采样信号 sft 完全确定，将采样信号通过一个低通滤波器，可以恢复出原始信号 f t。 那么，若采样间隔为sT ，应有smTw，则 ssssnnftf ttnTf nTtnT MATLAB语言28信号重构： 下面我们看如何由 sft ，来重构 f t 。 sft 经过一低通滤波器，那么 *sf tfth t 对于低通滤波器，有： 1cscwh tFHjwTSa w t 则有：( )( )()()()cscsscscsnwT wf tf tTSa w tf nT Sa w tnT 因为sin()sin ( )tc t

39、Satt，所以t( )sincSa t（ ） 这样，上面的式子可以写成：( )()sin ()ccsssnwwf tTf nTctnT MATLAB语言29例：信号的采样与重构。例：信号的采样与重构。wm=1 ; %wm=1 ; %信号带宽信号带宽wcwc=wm; %=wm; %滤波器截至频率滤波器截至频率tsts=pi/wm %=pi/wm %采样间隔，临界采样采样间隔，临界采样 wsws=2=2* *pi/pi/tsts; %; %采样角频率采样角频率n=-100:100; %n=-100:100; %时域采样点数目时域采样点数目ntsnts=n=n* *tsts; %; %时域采样点时域采样点f=f=sinc(ntssinc(nts/pi); %/pi); %采样信号采样信号dtdt=0.005;t=-15:dt:15; =