MATLAB在信号处理中的应用中.ppt

1,第一讲信号及其运算的MATLAB表示,1.1连续信号,1、指数信号形式：y=A*exp（a*t）,2、正弦信号形式：y=A*cos（w0*t+f）或：y=A*sin（w0*t+f）,2,3、抽样函数用sinc函数表示，其定义为：sinc(t)=sin(t)/(t)MATLAB语言调用形式:y=sinc(t)t=-3*pi:pi/100:3*pi;y=sinc(t/pi);plot(t,y),gridon,3,4、矩形脉冲信号MATLAB语言调用形式:y=rectpuls(t,width)产生一个幅值为1,宽度为width、相对于t=0点左右对称的矩形波信号t=0:0.001:4;T=1;y=rectpuls(t-2*T,2*T);plot(t,y),gridonaxis(0,4,-0.5,1.5)以t-2*T=0,即t=2*T为对称中心,4,周期性矩形脉冲信号MATLAB语言调用形式:y=square(t,duty)产生一个周期为2,幅值为1的周期性方波,duty表示占空比(dutycycle)t=-0.0625:0.0001:0.0625;y=square(2*pi*30*t,75);plot(t,y),gridonaxis(-0.1,0.1,-1.5,1.5)占空比75%,5,5、三角脉冲信号MATLAB语言调用形式:y=tripuls(t,width,skew)产生一个最大幅值为1,宽度为width、斜度为skew的三角波信号.t=-3:0.0001:3;y=tripuls(t,4,0.5);plot(t,y),gridonaxis(-3,3,-1.5,1.5),6,周期性三角脉冲信号MATLAB语言调用形式:y=sawtooth(t,width)产生一个周期为2,最大幅值为1,最小幅值为-1周期性三角波信号.t=-5*pi:0.0001:5*pi;y=sawtooth(t,0.5);plot(t,y),gridonaxis(-16,16,-1.5,1.5),7,6、一般周期性脉冲信号MATLAB语言调用形式:y=pulstran(t,d,func)依据名为func的连续函数并以之为一个周期,从而产生一串周期性的连续函数。用于指定周期性的偏移量，这个func函数会被计算length(d)次,最后值相当于:y=func(t-d(1)+func(t-d(2)+.,更一般的调用形式为:y=pulstran(t,d,func,p1,p2)。其中p1,p2为需要传送给func函数的额外输入参数值,最后值相当于:y=func(t-d(1),p1,p2,.)+func(t-d(2),p1,p2,.)+.,8,t=0:0.0001:1;d=0:1/3:1;y=pulstran(t,d,rectpuls,0.1);figure(1),plot(t,y),gridonaxis(0,1,-0.1,1.5)yy=pulstran(t,d,tripuls,0.1,-1);figure(2),plot(t,yy),gridonaxis(0,1,-0.1,1.5),9,1.2离散信号,1、指数序列a.k,2、正弦序列,离散的指数序列与正弦序列用MATLAB的表示与连续信号类似,只不过用stem函数而不是用plot函数来画出序列的波形.,10,3、单位冲激序列,单位冲激序列,k=-50:50;y=zeros(1,50),1,zeros(1,50);stem(k,y),functiony,k=impseq(k0,k1,k2);k=k1:k2;y=k-k0=0,写成MATLAB函数文件形式若k-k0在k1kk2范围内,11,4、单位阶跃序列,单位阶跃序列,k=-50:50;y=zeros(1,50),ones(1,51);stem(k,y),axis(-60,60,-0.5,1.5),写成函数文件形式若uk-k0在k1kk2,functiony,k=stepseq(k0,k1,k2);k=k1:k2;y=k-k0=0,12,1.3信号运算的MATLAB实现,2、离散序列的差分与求和,1、信号的尺度变换、反转和平移,差分f(k)=f(k)-f(k-1)用diff函数求和用sum函数,13,3、连续信号的微分与积分,微分可用diff函数积分用quad函数,调用形式:quad(function_name,a,b),例:对于三角波，画出df(t)/dt和t-f(t)dt的波形,functiony=funtrip(t)y=tripuls(t,4,0.5);,14,h=0.001,t=-3:h:3;y1=diff(funtrip(t)*1/h;figure(1),plot(t(1:length(t)-1),y1);title(df/dt);gridon,xlabel(t),15,forx=1:length(t)y2(x)=quad(funtrip,-3,t(x);endfigure(2),plot(t,y2);title(Intergralofdf);gridon,xlabel(t),16,第二讲时域分析的MATLAB表示,2.1连续时间系统的零状态响应,LTI连续时间系统以常系数微分方程描述，系统的零状态响应可通过求解初始状态为零的微分方程得到。形式：y=lsim（sys,f,t）sys是LTI系统模型，f是系统输入信号向量，是计算系统响应的抽样点向量。LTI系统模型sys要通过MATLAB中的tf函数来获得:sys=tf(b,a),其中b、a分别为微分方程右端和左端各项的系数向量。,17,若系统的微分方程为：,系统的输入信号为f(t)=10sin(t)，计算物体的位移y(t)。,ts=0;te=5;dt=0.01;sys=tf(1,1,2,100);t=ts:dt:te;f=10*sin(2*pi*t);y=lsim(sys,f,t);plot(t,y);xlabel(t(sec);ylabel(y(t);gridon,18,2.2连续时间系统的冲击响应和阶跃响应,1、冲击响应形式：y=impulse（sys,t）sys是LTI系统模型,sys=tf(b,a),其中b、a分别为微分方程右端和左端各项的系数向量。,2、阶跃响应形式：y=step(sys,t)sys是LTI系统模型,sys=tf(b,a),其中b、a分别为微分方程右端和左端各项的系数向量。,19,例：在一水平面上物体与弹簧组成的力学系统中，弹簧一端固定，物体的质量m=1kg，弹簧的弹性系数ks=100N/m，物体与地面的摩擦系数fd=2Ns/m，在物体的水平方向施加一外力f(t)，外力是强度为10的冲击信号，若系统的微分方程为：,求物体的位移y(t)。,20,ts=0;te=5;dt=0.01;sys=tf(10,12100);t=ts:dt:te;y=impulse(sys,t);plot(t,y);xlabel(t(sec);ylabel(y(t);gridon,21,2.3离散时间系统的零状态响应,LTI离散时间系统常用线性常系数差分方程描述：,已知差分方程的n个初始状态及输入fk，可由下式迭代算出系统的输出：,fk和yk分别表示系统的输入和输出,n是差分方程的阶数。,22,在零初始状态下，MATLAB调用的形式为：y=filter（b，a，f）,23,R=51;d=rand(1,51)-0.5;k=0:R-1;s=2*k.*(0.9.k);f=s+d;subplot(1,2,1)plot(k,d,r-.,k,s,b:,k,f,k-);xlabel(Timek);legend(dk,sk,fk);M=5;b=ones(M,1)/M;a=1;y=filter(b,a,f);subplot(1,2,2)plot(k,s,b:,k,y,r-);xlabel(Timek);legend(sk,yk);,24,25,2.4离散时间系统的冲击响应和阶跃响应,1、单位冲击响应形式：h=impz（b,a,k）其中b、a分别为差分方程右端和左端的系数向量，k是输出序列的取值范围，h是系统的单位冲击响应,2、单位阶跃响应形式：h=stepz（b,a,k）,26,例：用impz函数求下列离散时间系统的单位冲击响应hk，并与理论值hk=-(-1)k+2(-2)-k，k0进行比较：yk+3yk-1+2yk-2=fk,k=0:10;a=132;b=1;h=impz(b,a,k);subplot(2,1,1)stem(k,h);title(单位冲击相应的近似值);gridon;hk=-(-1).k+2*(-2).k;subplot(2,1,2)stem(k,hk);title(单位冲击相应的理论值);gridon;,27,28,2.5卷积的计算,两个离散序列卷积和：c=conv(a,b),29,第三讲频域分析的MATLAB表示,3.1周期信号的频域分析,1、周期性离散时间信号离散傅里叶级数(DFS)和逆离散傅里叶级数（IDFS）表示为：,其中WN=e-j2/N,30,MATLAB表示为：F=fft(f)f=ifft(F),例：利用MATLAB计算一周期矩形波序列的DFS系数。,N=32;M=4;%定义周期矩形波序列的参数%产生该周期矩形波序列f=ones(1,M+1),zeros(1,N-2*M-1),ones(1,M);F=fft(f);%计算DFS系数m=0:N-1;figure(1),stem(m,real(F);title(Fm的实部);xlabel(m);figure(2),stem(m,imag(F);title(Fm的虚部);xlabel(m);fr=ifft(F);%计算由DFS系数Fm重建的fkfigure(3)stem(m,real(fr);title(重建的fk);xlabel(k);,31,32,其中func是被积函数的函数文件名，a、b分别为定积分的下上限，Tol表示指定允许的相对或绝对积分误差，非零的Trace表示以被积函数的点绘图形式来跟踪函数生成的返回值，若Tol和Trace均赋予“”，则两者均自动使用缺省值；”p1,p2,.”表示被积函数所需的多个额外输入参数,y=quad8(func,a,b)y=quad8(func,a,b,Tol,Trace,p1,p2,.),3.2非周期信号的频域分析,1、非周期性连续时间信号的频谱quad8和quadl可以计算非周期性连续时间信号的频谱（高阶算法）。,33,例：利用MATLAB采用数值方法近似计算三角波信号的频谱。,functiony=sfl(t,w)%对于不同的，计算出傅里叶变换中被积函数的值y=(abs(t)=1).*(1-abs(t).*exp(-j*w*t);,34,w=linspace(-6*pi,6*pi,512);N=length(w);F=zeros(1,N);fork=1:NF(k)=quad8(sfl,-1,1,w(k);endfigure(1),plot(w,real(F);xlabel(omega);ylabel(F(jomega);title(三角波信号近似频谱)；figure(2),plot(w,real(F)-sinc(w/2/pi).2);%sinc函数值即为理论计算结果xlabel(omega);ylabel(error);title(计算误差),35,36,2、非周期性离散时间信号的频谱当序列的DTFT（离散时间傅里叶变换）可写成ej的有理多项式时，MATLAB的函数freqz可用来计算其DTFT的值。设DTFT的有理多项式为：,函数freqz的一般调用形式：h=freqz(b,a,w),W为频率抽样点向量,37,例：利用MATLAB画出a=0.9时的幅度频谱，,b=1;a1=1,-0.9;a2=1,0.9;w=linspace(0,2*pi,512);h1=freqz(b,a1,w);h2=freqz(b,a2,w);plot(w/pi,abs(h1),w/pi,abs(h2),k:)xlabel(omega/pi);legend(alpha=0.9,alpha=-0.9);,38,39,1、连续时间系统的频率特性,3.3系统的频率特性分析,当频率响应H(jw)是jw的有理多项式时，即：,或写为H(s)的形式：,40,MATLAB提供freqs函数计算连续时间系统的频率响应:H=freqs(b,a,w)H,w=freqs(b,a,N)其中，W为频率抽样点向量,MATLAB自动选择N个频率点计算其频率响应，N的缺省值为200,41,w=linspace(0,6*pi,200);b=1;a=1,2,2,1;H=freqs(b,a,w);subplot(2,1,1)plot(w,abs(H)set(gca,xtick,0,1,2,3,4,5);set(gca,ytick,0,0.4,0.707,1);xlabel(omega(rad/s);ylabel(|H(jomega)|);gridon,42,subplot(2,1,2)plot(w,angle(H)set(gca,xtick,0,1,2,3,4,5);set(gca,ytick,-4:2:4);xlabel(omega(rad/s);ylabel(phi(rad);gridon,43,2、离散时间系统的频率特性,当频率响应H(ej)是ej的有理多项式时，即：,或写为H(z)的形式，即：,44,MATLAB提供freqz函数计算离散时间系统的频率响应:H=freqz(b,a,w)H,w=freqz(b,a,N)其中，W为频率抽样点向量,MATLAB在单位圆的上半周上等间距分布的N个频率点计算其频率响应，N的缺省值为512,45,MATLAB音频支持的函数,1、sound重放矢量成为声音，格式：sound(y,fs,bits)2、soundsc声音自动定标和重放矢量，格式：soundsc(y,fs,bits,slim),slim=slow,shigh3、wavplay使用视窗音频的输出装置重放声音格式：wavplay(y,fs),46,4、wavread读*.wav声音文件x,fs,bits=wavread(file)5、wavrecord使用视窗音频的输入装置记录声音格式：wavrecord(n,fs,ch)6、wavwrite写*.wav声音文件格式：wavwrite(y,fs,nbits,wavefile),47,程序设计实例,1、音频信号的分析与合成,采用MATLAB分析WAV文件。,步骤一：选择一个WAV文件作为分析的对象。（*.wav）,步骤二：读WAV文件数据并画时域图形。,步骤三：进行FFT变换并画频域图形。,步骤四：进行该声波主要频谱的分析。,步骤五：根据该声音的频谱，反演时域图形。（有失真）,常用函数：,x，fs，bits=wavread(filename),d=FFT(w，l),sound(w，fs，bits),步骤六：进行付立叶逆变换IFFT并画频域图形。,48,音频信号的分析与合成程序,返回,%=clear;ifnargin1;action=initialized;end;%显示打开文件的对话框fname,pname=uigetfile(SOUND1.wav,OpenWaveFile);file=pname,fname;x,fs,bits=wavread(file);%读入声音文件（*.wav）sound(x,fs,bits);%数据通过声卡转换为声音%=选择声道%y=x(:);%双声道信号波形数据y=x(:,1);%左声道信号波形数据%y=x(:,2);%右声道信号波形数据,49,%=disp(按任意键显示左声道声音波形！);disp(-);pauseplot(y);set(gca,FontName,Arial,FontSize,10);title(左声道声音波形);xlabel(采样点数);disp(按任意键听左声道信号三次！);disp(-);pause,sound(y,fs,bits);%数据通过声卡转换为声音pause,sound(y,fs,bits);pause,sound(y,fs,bits);k=fft(y,32768);disp(按任意键显示左声道声音频域的幅值！);d