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

单击此处编辑母版标题样式,*,*,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第一讲 信号及其运算的,MATLAB,表示,1.1,连续信号,1,、指数信号,形式：,y=A*exp,（,a*t,）,2,、正弦信号,形式：,y=A*cos,（,w,0,*t+,f,）,或：,y=A*sin,（,w,0,*t+,f,）,1,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),grid on,2,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),grid on,axis(0,4,-0.5,1.5),以,t-2*T=0,即,t=2*T,为,对称中心,3,周期性矩形脉冲信号,MATLAB,语言调用形式,:,y=square(t,duty),产生一个周期为,2,幅值为,1,的周期性方波,duty,表示占空比,(duty cycle),t=-0.0625:0.0001:0.0625;,y=square(2*pi*30*t,75);,plot(t,y),grid on,axis(-0.1,0.1,-1.5,1.5),占空比,75%,4,5,、三角脉冲信号,MATLAB,语言调用形式,:,y=tripuls(t,width,skew),产生一个最大幅值为,1,宽度为,width,、斜度为,skew,的三角波信号,.,t=-3:0.0001:3;,y=tripuls(t,4,0.5);,plot(t,y),grid on,axis(-3,3,-1.5,1.5),5,周期性三角脉冲信号,MATLAB,语言调用形式,:,y=sawtooth(t,width),产生一个周期为,2,最大幅值为,1,最小幅值为,-1,周期性,三角波信号,.,t=-5*pi:0.0001:5*pi;,y=sawtooth(t,0.5);,plot(t,y),grid on,axis(-16,16,-1.5,1.5),6,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,.)+.,7,t=0:0.0001:1;d=0:1/3:1;,y=pulstran(t,d,rectpuls,0.1);,figure(1),plot(t,y),grid on,axis(0,1,-0.1,1.5),yy=pulstran(t,d,tripuls,0.1,-1);figure(2),plot(t,yy),grid on,axis(0,1,-0.1,1.5),8,1.2,离散信号,1,、指数序列,a.k,2,、正弦序列,离散的指数序列与正弦序列用,MATLAB,的表示与连续信号类似,只不过用,stem,函数而不是用,plot,函数来画出序列的波形,.,9,3,、单位冲激序列,单位冲激序列,k=-50:50;,y=zeros(1,50),1,zeros(1,50);,stem(k,y),function y,k=impseq(k0,k1,k2);,k=k1:k2;,y=k-k0=0,写成,MATLAB,函数文件形式若,k-k0,在,k1,k,k2,范围内,10,4,、单位阶跃序列,单位阶跃序列,k=-50:50;,y=zeros(1,50),ones(1,51);,stem(k,y),axis(-60,60,-0.5,1.5),写成函数文件形式若,uk-k0,在,k1,k,k2,function y,k=stepseq(k0,k1,k2);,k=k1:k2;,y=k-k0=0,11,1.3,信号运算的,MATLAB,实现,2,、离散序列的差分与求和,1,、信号的尺度变换、反转和平移,差分,f(k)=f(k)-f(k-1),用,diff,函数,求和用,sum,函数,12,3,、连续信号的微分与积分,微分,可,用,diff,函数,积分用,quad,函数,调用形式,:,quad(function_name,a,b),例,:,对于三角波，画出,df(t)/dt,和,t,-,f(t)dt,的波形,function y=funtrip(t),y=tripuls(t,4,0.5);,13,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);grid on,xlabel(t),14,for x=1:length(t),y2(x)=quad(funtrip,-3,t(x);,end,figure(2),plot(t,y2);,title(Intergral of df);grid on,xlabel(t),15,第二讲 时域分析的,MATLAB,表示,2.1,连续时间系统的零状态响应,LTI,连续时间系统以常系数微分方程描述，系统的零状态响应可通过求解初始状态为零的微分方程得到。,形式：,y=lsim,（,sys,f,t,）,sys,是,LTI,系统模型，,f,是系统输入信号向量，是计算系统响应的抽样点向量。,LTI,系统模型,sys,要通过,MATLAB,中的,tf,函数来获得,:,sys=tf(b,a),其中,b,、,a,分别为微分方程右端和左端各项的系数向量。,16,若系统的微分方程为：,系统的输入信号为,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);grid on,17,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,分别为微分方程右端和左端各项的系数向量。,18,例：,在一水平面上物体与弹簧组成的力学系统中，弹簧一端固定，物体的质量,m=1kg,，弹簧的弹性系数,k,s,=100N/m,，物体与地面的摩擦系数,f,d,=2Ns/m,，在物体的水平方向施加一外力,f(t),，外力是强度为,10,的冲击信号，若系统的微分方程为：,求物体的位移,y(t),。,19,ts=0;te=5;dt=0.01;,sys=tf(10,1 2 100);,t=ts:dt:te;,y=impulse(sys,t);,plot(t,y);,xlabel(t(sec);,ylabel(y(t);,grid on,20,2.3,离散时间系统的零状态响应,LTI,离散时间系统常用线性常系数差分方程描述：,已知差分方程的,n,个初始状态及输入,fk,，,可由下式迭代算出系统的输出：,fk,和,yk,分别表示系统的输入和输出,n,是差分方程的阶数。,21,在零初始状态下，,MATLAB,调用的形式为：,y=filter,（,b,，,a,，,f,）,例：,受噪声干扰的信号为,fk=sk+dk,，其中,sk=(2k)0.9,k,是原始信号，,dk,是噪声，已知,M,点滑动平均系统的输入和输出关系为：,试用,MATLAB,编程实现用点滑动平均系统对受噪声干扰的信号去噪。,22,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(Time k);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(Time k);legend(sk,yk);,23,24,2.4,离散时间系统的冲击响应和阶跃响应,1,、,单位冲击响应,形式：,h=impz,（,b,a,k,）,其中,b,、,a,分别为差分方程右端和左端的系数向量，,k,是输出序列的取值范围，,h,是系统的单位冲击响应,2,、,单位阶跃响应,形式：,h=stepz,（,b,a,k,）,25,例：,用,impz,函数求下列离散时间系统的单位冲击响应,hk,，并与理论值,hk=-(-1),k,+2(-2),-k,，,k,0,进行比较：,yk+3yk-1+2yk-2=fk,k=0:10;,a=1 3 2;b=1;,h=impz(b,a,k);,subplot(2,1,1),stem(k,h);title(,单位冲击相应的近似值,);grid on;,hk=-(-1).k+2*(-2).k;,subplot(2,1,2),stem(k,hk);title(,单位冲击相应的理论值,);grid on;,26,27,2.5,卷积的计算,两个离散序列卷积和：,c=conv(a,b),28,第三讲 频域分析的,MATLAB,表示,3.1,周期信号的频域分析,1,、周期性离散时间信号,离散傅里叶级数,(DFS),和逆离散傅里叶级数（,IDFS,）表示为：,其中,W,N,=e,-j2,/N,29,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,重建的,fk,figure(3)stem(m,real(fr);title(,重建的,fk);xlabel(k);,30,31,其中,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,可以计算非周期性连续时间信号的频谱（高阶算法）。,32,例：,利用,MATLAB,采用数值方法近似计算,三角波信号 的 频谱。,function y=sfl(t,w),%,对于不同的，计算出傅里叶变换中被积函数的值,y=(abs(t)=1).*(1-abs(t).*exp(-j*w*t);,33,w=linspace(-6*pi,6*pi,512);,N=length(w);F=zeros(1,N);,for k=1:N,F(k)=quad8(sfl,-1,1,w(k);,end,figure(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(,计算误差,),34,35,2,、非周期性离散时间信号的频谱,当序列的,DTFT,（离散时间傅里叶变换）可写成,e,j,的有理多项式时，,MATLAB,的函数,freqz,可用来计算其,DTFT,的值。设,DTFT,的有理多项式为：,函数,freqz,的一般调用形式：,h=freqz(b,a,w),W,为频率抽样点向量,36,例：,利用,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);,37,38,1,、连续时间系统的频率特性,3.3,系统的频率特性分析,当频率响应,H(j,w,),是,j,w,的有理多项式时，即：,或写为,H(s),的形式：,39,MATLAB,提供,freqs,函数计算连续时间系统的频率响应,:,H=freqs(b,a,w),H,w=freqs(b,a,N),其中，,W,为频率抽样点向量,MATLAB,自动选择,N,个频率点,计算其频率响应，,N,的缺省值为,200,40,例：,三阶归一化的低通滤波器的频率响应为：,利用,MATLAB,画出该系统的幅频响应,H(j,w,),和相频响应,j,(,w,),。,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)|);grid on,41,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);grid on,42,2,、离散时间系统的频率特性,当频率响应,H(e,j,),是,e,j,的有理多项式时，即：,或写为,H(z),的形式，即：,43,MATLAB,提供,freqz,函数计算离散时间系统的频率响应,:,H=,freqz(b,a,w,),H,w,=,freqz(b,a,N,),其中，,W,为频率抽样点向量,MATLAB,在单位圆的上半周上等间距分布的,N,个频率点计算其频率响应，,N,的缺省值为,512,44,MATLAB,音频支持的函数,1,、,sound,重放矢量成为声音，,格式：,sound(y,fs,bits),2,、,soundsc,声音自动定标和重放矢量，,格式：,soundsc(y,fs,bits,slim),slim=slow,shigh,3,、,wavplay,使用视窗音频的输出装置重放声音,格式：,wavplay,(y,fs),45,4,、,wavread,读*,.wav,声音文件,x,fs,bits=wavread(file),5,、,wavrecord,使用视窗音频的输入装置记录声音,格式：,wavrecord,(n,fs,ch),6,、,wavwrite,写*,.wav,声音文件,格式：,wavwrite,(y,fs,nbits,wavefile),46,程序设计实例,1、,音频信号的分析与合成,采用,MATLAB,分析,WAV,文件。,步骤一：选择一个,WAV,文件作为分析的对象。（,*.wav）,步骤二：读,WAV,文件数据并画时域图形。,步骤三：进行,FFT,变换,并画频域图形,。,步骤四：进行该声波主要频谱的分析,。,步骤五：根据该声音的频谱，反演时域图形,。（有,失真,）,常用函数：,x，,fs,，bits=wavread(filename),d=FFT(w，l),sound(w，,fs,，bits),步骤六：进行,付立叶逆变换,IFFT,并画频域图形,。,47,音频信号的分析与合成程序,返回,%=,clear;,if,nargin,1;action=initialized;end;,%,显示打开文件的对话框,fname,pname,=,uigetfile,(SOUND1.wav,Open Wave File);,file=,pname,fname,;,x,fs,bits=,wavread,(file);%,读入声音文件,（*.,wav）,sound(x,fs,bits);%,数据通过声卡转换为声音,%=选择声道,%,y=x(:);%,双声道信号波形数据,y=x(:,1);%,左声道信号波形数据,%,y=x(:,2);%,右声道信号波形数据,48,%=,disp(,按任意键显示左声道声音波形！);,disp,(-);,pause,plot(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,(,按任意键显示左声道声音频域的幅值！);,disp,(-);,49,disp,(-);,pause,plot(abs(k);,set(,gca,FontName,Arial,FontSize,10);,title(,左声道声音频域的幅值);,xlabel,(,采样点数);,m1,i1=max(abs(k);%,找出频域最大值,F1=i1/32768*fs;F1t=num2str(F1)Hz;,text(1500,450,F1t);,m2,i2=max(abs(k(4000:5000);%,找出频域次大值,F2=(4000+i2)/32768*,fs,;,F2t=num2str(F2)Hz;text(3000,40,F2t);,m3,i3=max(abs(k(10000:12000);%,找出频域次次大值,F3=(10000+i3)/32768*,fs,;,F3t=num2str(F3)Hz;text(8000,20,F3t);,%=,50,%=,t=(0:0.0001:1);,y1=(sin(2*pi*788*t)+sin(2*pi*3174*t)*(22.65/490)+sin(2*pi*6936*t)*(2.996/490)*0.18;,disp,(,按任意键显示左声道简单合成声音波形！);,disp,(-);,pause,plot(t,y1);,set(,gca,FontName,Arial,FontSize,10);,title(,左声道简单合成声音波形);,disp,(,按任意键听简单合成左声道信号三次！);,disp,(-);,pause,sound(y1,fs,bits);%,数据通过声卡转换为声音,pause,sound(y1,fs,bits);,pause,sound(y1,fs,bits);,k1=,fft,(y1,32768);,disp,(,按任意键显示左声道简单合成声音的频域幅值！);,disp,(-);,51,