MATLAB在信号处理中的应用中

1第一讲信号及其运算的MATLAB表示1.1连续信号1、指数信号形式：y=A*exp（a*t）2、正弦信号形式：y=A*cos（w0*t+f）或：y=A*sin（w0*t+f）23、抽样函数用sinc函数表示，其定义为：sinc(t)=sin(πt)/(πt)MATLAB语言调用形式:y=sinc(t)t=-3*pi:pi/100:3*pi;y=sinc(t/pi);plot(t,y),gridon34、矩形脉冲信号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为对称中心00.511.522.533.54-0.500.511.54周期性矩形脉冲信号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%-0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1-1.5-1-0.500.511.555、三角脉冲信号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])-3-2-10123-1.5-1-0.500.511.56周期性三角脉冲信号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])-15-10-5051015-1.5-1-0.500.511.576、一般周期性脉冲信号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,...)+...8t=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])00.10.20.30.40.50.60.70.80.9100.511.500.10.20.30.40.50.60.70.80.9100.511.591.2离散信号1、指数序列a.^k2、正弦序列离散的指数序列与正弦序列用MATLAB的表示与连续信号类似,只不过用stem函数而不是用plot函数来画出序列的波形.103、单位冲激序列单位冲激序列k=-50:50;y=[zeros(1,50),1,zeros(1,50)];stem(k,y)-50-40-30-20-100102030405000.10.20.30.40.50.60.70.80.91function[y,k]=impseq(k0,k1,k2);k=k1:k2;y=[k-k0==0]写成MATLAB函数文件形式若δ[k-k0]在k1≤k≤k2范围内,0,0,...}{...,0,0,10k0kk01][114、单位阶跃序列单位阶跃序列k=-50:50;y=[zeros(1,50),ones(1,51)];stem(k,y),axis([-60,60,-0.5,1.5])-60-40-200204060-0.500.511.5写成函数文件形式若u[k-k0]在k1≤k≤k2function[y,k]=stepseq(k0,k1,k2);k=k1:k2;y=[k-k0=0],1,1,...}{...,0,0,10k0kku01][121.3信号运算的MATLAB实现2、离散序列的差分与求和1、信号的尺度变换、反转和平移差分∇f(k)=f(k)-f(k-1)用diff函数求和用sum函数133、连续信号的微分与积分微分可用diff函数积分用quad函数,调用形式:quad(‘function_name’,a,b)例:对于三角波，画出df(t)/dt和t-∞f(t)dt的波形functiony=funtrip(t)y=tripuls(t,4,0.5);14h=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-3-2-1012300.511.522.5Integralofdftforx=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)()()()(tftydttdydttyd100222若系统的微分方程为：系统的输入信号为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)');gridon00.511.522.533.544.55-0.25-0.2-0.15-0.1-0.0500.050.10.150.2t(sec)y(t)182.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的冲击信号，若系统的微分方程为：)()()()(ttydttdydttyd10100222求物体的位移y(t)。2000.511.522.533.544.55-0.8-0.6-0.4-0.200.20.40.60.81t(sec)y(t)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)');gridon212.3离散时间系统的零状态响应LTI离散时间系统常用线性常系数差分方程描述：nimjjijkfbikya00][][nimjjijkfabikyaaky1000][)/(][)/(][已知差分方程的n个初始状态及输入f[k]，可由下式迭代算出系统的输出：f[k]和y[k]分别表示系统的输入和输出,n是差分方程的阶数。22在零初始状态下，MATLAB调用的形式为：y=filter（b，a，f）例：受噪声干扰的信号为f[k]=s[k]+d[k]，其中s[k]=(2k)0.9k是原始信号，d[k]是噪声，已知M点滑动平均系统的输入和输出关系为：101MnnkfMky][][试用MATLAB编程实现用点滑动平均系统对受噪声干扰的信号去噪。23R=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('d[k]','s[k]','f[k]');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('s[k]','y[k]');240204060-1012345678Timekd[k]s[k]f[k]0204060012345678Timeks[k]y[k]252.4离散时间系统的冲击响应和阶跃响应1、单位冲击响应形式：h=impz（b,a,k）其中b、a分别为差分方程右端和左端的系数向量，k是输出序列的取值范围，h是系统的单位冲击响应2、单位阶跃响应形式：h=stepz（b,a,k）26例：用impz函数求下列离散时间系统的单位冲击响应h[k]，并与理论值h[k]=-(-1)k+2(-2)-k，k≥0进行比较：y[k]+3y[k-1]+2y[k-2]=f[k]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;27282.5卷积的计算两个离散序列卷积和：c=conv(a,b)29第三讲频域分析的MATLAB表示3.1周期信号的频域分析1、周期性离散时间信号离散傅里叶级数(DFS)和逆离散傅里叶级数（IDFS）表示为：10101NmmkNNkmkNWmFNmFIDFSkfWkfkfDFSmF][]}[{][][]}[{][其中WN=e-j2π/N30MATLAB表示为：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('F[m]的实部');xlabel('m');figure(2),stem(m,imag(F));tit