MATLAB在信号处理中的应用

MATLAB在信号处理中的应用4.1信号及其表示4.1.1连续时间信号的表示连续时间信号：时间变化连续。如y=x(t)离散时间信号(序列)：时间离散，如x(nT)=x(t)|t=nT.4.1.2工具箱中的信号产生数函数名功能函数名功能sawtooth产生锯齿波或三角波信号pulstran产生冲激串square产生方波信号rectpule产生非周期的方波信号sinc产生sinc函数波形tripuls产生非周期的三角波信号chirp产生调频余弦信号diric产生Dirichlet或周期sinc函数gauspuls产生高斯正弦脉冲信号gmonopuls产生高斯单脉冲信号vco电压控制振荡器例：产生f=50Hz的锯齿波、三角波Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=sawtooth(2*pi*50*t,0);x2=sawtooth(2*pi*50*t,1);x3=sawtooth(2*pi*50*t,0.5);subplot(311);plot(t,x1);subplot(312);plot(t,x2);subplot(313);plot(t,x3)1、sawtooth函数产生锯齿波或三角波用法：　　x=sawtooth(t)将产生周期为2π的锯齿波。以0~2π这个周期内为例，当t=0时，x=-1，当t=2π时，x=1。　　x=sawtooth(t,width)　　width是0到1之间的标量。在0到2π×width区间内，x的值从-1线性变化到1；在2π×width~2π区间内，x的值又从1线性变化到-1。sawtooth(t,1)和sawtooth(t)是等价的。2、square函数产生矩形波例：产生50Hz占空比分别为20％和50％的矩形波。Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=square(2*pi*50*t,20);x2=square(2*pi*50*t,50);subplot(211);plot(t,x1);subplot(212);plot(t,x2);3、sinc函数产生sinc波形或sin(πt)/(πt)波形例：t=linspace(-10,+10,200);x=sinc(t);plot(t,x);4.1.3离散时间信号的表示在MATLAB中，离散时间信号x(n)的表示：需用一个向量x表示序列幅值，用另一个等长的定位时间变量n，才能完整地表示一个序列。[例4-10]绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。MATLAB源程序为：n=-3:5;%定位时间变量x=[0,0,-1,1,2,1,-1,0,0];stem(n,x);grid;%绘制棒状图line([-3,5],[0,0]);%画x轴线xlabel('n');ylabel('x[n]')运行结果如图4.11所示。图4.11离散时间信号图形line([起点横坐标，终点横坐标],[起点纵坐标，终点纵坐标])4.1.4几种常用离散时间信号的表示1．单位脉冲序列直接实现：x=zeros(1,N);x(1,n0)=1;2．单位阶跃序列直接实现：n=[ns:nf];x=[(n-n0)>=0];x=zeros(1,5);x(1,3)=1;>>n=-2:2;>>stem(n,x)n=[1:4];x=[(n-2)>=0];stem(n,x)3．实指数序列直接实现：n=[ns:nf];x=a.^n;4．复指数序列直接实现：n=[ns:nf];x=exp((sigema+jw)*n);5．正(余)弦序列直接实现：n=[ns:nf];x=cos(w*n+sita);4.2信号的基本运算4.2.1信号的相加与相乘y(n)=x1(n)+x2(n)y(n)=x1(n)×x2(n)MATLAB实现：y=x1+x2;y=x1.*x24.2.2序列移位与周期延拓运算序列移位：y(n)=x(n-m)。MATLAB实现：y=x;ny=nx-m序列周期延拓：y(n)=x((n))M，MATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1)4.2.3序列翻褶与序列累加运算序列翻褶：y(n)=x(-n)。MATLAB可实现:y=fliplr(x)序列累加的数学描述为：MATLAB实现：y=cumsum(x)4.2.4两序列的卷积运算两序列卷积运算：MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。4.2.5两序列的相关运算两序列相关运算：。MATLAB实现：y=xcorr(x1,x2)。已知离散信号x(n)和h(n)，求y(n)=x(n)*h(n)，并用图形表示。Nh=20;Nx=10;m=5;%设定Nx，Nh和位移值mn=0:Nh-1;h1=(0.9).^n;%产生h1(n)h2=h1;nx=0:Nx-1;x1=ones(1,Nx);%产生x1(n)x2=zeros(1,Nx+m);fork=m+1:m+Nx%产生x2(n)=x1(n-m)x2(k)=x1(k-m);end%产生x2(n)y1=conv(x1,h1);%计算y1(n)=x1(n)*h1(n)y2=conv(x2,h2);%计算y2(n)=x2(n)*h2(n)subplot(3,2,1)stem(nx,x1,'.')axis([03001.2]),title(‘x1(n)’)%绘图……（以下省略）4.4线性时不变系统4.4.1系统的描述系统传递函数　　在Matlab中，传递函数用分子、分母两个多项式的系数表示，系数为降幂排列。线性时不变（LTI）系统的常用表示方法1、传递函数表示法　　在Matlab中，传递函数用分子、分母两个多项式的系数表示，系数为降幂排列。分子（Numerator）：B=[b(1)b(2)…b(m+1)]分母（Denominator）：A=[a(1)a(2)…b(n+1)]例：num=[10.21];den=[10.51];2、零－极点模型表示法　　在Matlab中，增益系数、零点向量、极点向量用三个列向量表示。增益系数（Gain）：k零点向量（Zero）：z=[z1z2…zn]'极点向量（Pole）：p=[p1p2…pn]'sys=zpk(z,p,k)%获得零-极点模型表达式>>k=3;>>z=[234]';>>p=[567]';>>sys=zpk(z,p,k,0.1)Zero/pole/gain:3(z-2)(z-3)(z-4)-------------------(z-5)(z-6)(z-7)Samplingtime:0.1线性系统模型的变换函数4.4.2系统模型的转换函数函数名功能说明函数名功能说明ss2tf状态空间模型转换为传递函数模型zp2tf零－极点增益模型转换为传递函数模型ss2zp状态空间模型转换为零－极点增益模型zp2ss零－极点增益模型转换为状态空间模型ss2sos状态空间模型转换为二次分式模型zp2sos零－极点增益模型转换为二次分式模型tf2ss传递函数模型转换为状态空间模型sos2tf二次分式模型转换为传递函数模型tf2zp传递函数模型转换为零－极点增益模型sos2zp二次分式模型转换为零－极点增益模型tf2sos传递函数模型转换为二次分式模型sos2ss二次分式模型转换为状态空间模型[例4-18]求离散时间系统的零、极点向量和增益系数。在命令窗口输入：>>num=[2,3];den=[1,0.4,1];>>[num,den]=eqtflength(num,den);　%使长度相等>>[z,p,k]=tf2zp(num,den)屏幕显示为z=　0-1.5000p=-0.2000+0.9798i　-0.2000-0.9798ik=24.5线性时不变系统的响应4.5.1　线性时不变系统的时域响应1．连续LTI系统的响应2．离散LTI系统的响应用MATLAB中的卷积函数conv()来实现。用MATLAB中的卷积函数conv()来实现。4.5.2LTI系统的单位冲激响应1.求连续LTI系统的单位冲激响应函数impulse()格式：[Y,T]=impulse(sys)　或impulse(sys)功能：返回系统的响应Y和时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。2.求离散系统的单位冲激响应函数dimpulse()格式：[y,x]=dimpulse(num,den)功能：返回项式传递函数的单位冲激响应y向量和时间状态历史记录x向量。1、impulse函数求连续系统的单位冲击响应。b=[23];a=[10.41];sys=tf(b,a)Transferfunction:2s+3---------------s^2+0.4s+1>>[yt]=impulse(sys);>>plot(t,y)>>*2、impz函数求离散系统（数字滤波器）的单位冲击响应。（注：Matlab7.0不再支持dimpulse函数）4.5.3时域响应的其它函数1.求连续LTI系统的零输入响应函数initial()格式：[y,t,x]=initial(a,b,c,d,x0)功能：计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x为状态记录，t为仿真所用的采样时间向量。2.求离散系统的零输入响应函数dinitial()格式：[y,x,n]=dinitial(a,b,c,d,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。n为仿真所用的点数。3.求连续系统的单位阶跃响应函数step()格式：[Y,T]=step(sys)　功能：返回系统的单位阶跃响应Y和仿真所用的时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数(TF)、零极增益模型(ZPK)或状态空间模型(SS)。4.求离散系统的单位阶跃响应函数dstep()格式：[y,x]=dstep(num,den)功能：返回多项式传递函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。4.6线性时不变系统的频率响应1．求模拟滤波器Ha(s)的频率响应函数freqs()格式：H＝freqs(B,A,W)功能：计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(jΩ)，结果存于H向量中。[例4-31]已知某模拟滤波器的系统函数求该模拟滤波器的频率响应。MATLAB源程序如下。B=1;A=[12.61313.41422.61311];W=0:0.1:2*pi*5;freqs(B,A,W)图4.30模拟滤波器的频率响应[例4-32]已知某滤波器的系统函数为求该滤波器的频率响应。MATLAB源程序为：B=[10000000–1];A=1;freqz(B,A)该程序运行所绘出的幅频与相频性曲线如图4.31所示。图4.31滤波器幅度和相位曲线2．求数字滤波器H(z)的频率响应函数freqz()格式：H=freqz(B,A,W)功能：计算由向量W(rad)指定的数字频率点上(通常指[0,π]范围)在H(z)的频率响应H(ejw)。3．滤波函数filter格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。[例4-33]设系统差分方程为MATLAB源程序为：B=1;A=[1,-0.8];N=0:31;x=0.8.^n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)该程序运行所得结果如图4.32所示。，求该系统对信号的响应。图4.32系统对信号的响应4.7傅里叶(Fourier)变换4.7.1连续时间、连续频率－傅里叶变换4.7.2连续时间、离散频率－傅里叶级数正变换：逆变换：正变换：逆变换：4.7.3时间离散、连续频率－序列傅里叶变换4.7.4离散时间、离散频率－离散傅里叶级数4.7.5离散时间、离散频率－离散傅里叶变换(DFT)正变换：逆变换：正变换：逆变换：正变换：逆变换：1．一维快速正傅里叶变换函数fft格式：X=fft(x,N)功能：采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用混合算法。2．一维快速逆傅里叶变换函数ifft格式：x=ifft(X,N)功能：采用FFT算法计算序列向量X的N点IDFT变换。[例4-36]用快速傅里叶变换FFT计算下面两个序列的卷积。,并测试直接卷积和快速卷积的时间。图4.35快速卷积框图MATLAB程序(部分)：%线性卷积xn=sin(0.4*[1:15]);　%对序列x(n)赋值,M=15hn=0.9.^(1:20);　　%对序列h(n)赋值,N=20yn=conv(xn,hn);%直接调用函数conv计算卷积%园周卷积L=pow2(nextpow2(M+N-1));Xk=fft(xn,L);Hk=fft(hn,L);Yk=Xk.*Hk;yn=ifft(Yk,L);图4.36x(n),h(n)及其线性卷积波形4.8IIR数字滤波器的方法1.数字滤波器的频率响应函数幅度响应：相位响应：图4.37理想低通、高通、带通、带阻数字滤波器幅度特性2.滤波器的技术指标幅度响应指标、相位响应指标图4.38数字低通滤波器的幅度特性通带要求：阻带要求：通带最大衰减：阻带最小衰减：4.8.1冲激响应不变法2.MATLAB信号处理工箱中的专用函数impinvar()：格式：[BZ,AZ]=impinvar(B,A,Fs)功能：把具有[B,A]模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的传递函数模型[BZ,AZ]。采样频率Fs的默认值为Fs=1。1.冲激响应不变法设计IIR数字滤波器的基本原理：[例4-37]MATLAB源程序如下：num=[1];%模拟滤波器系统函数的分子den=[1,sqrt(5),2,sqrt(2),1];%模拟滤波器系统函数的分母[num1,den1]=impinvar(num,den)%求数字低通滤波器的系统函数程序的执行结果如下：num1=-0.00000.09420.21580.0311den1=1.0000-2.00321.9982-0.76120.1069MATLAB信号处理工具箱中的专用双线性变换函数bilinear()格式：[numd,dend]＝bilinear(num,den,Fs)功能：把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。4.8.2双线性变换法双线性变换利用频率变换关系：[例4-38]MATLAB源程序如下：num=1;%模拟滤波器系统函数的分子den=[1,sqrt(3),sqrt(2),1];%模拟滤波器系统函数的分母[num1,den1]=bilinear(num,den,1)%求数字滤波器的传递函数运算的结果如下：　num1=0.05330.15990.15990.0533den1=1.0000-1.33820.9193-0.15464.8.3IIR数字滤波器的频率变换设计法1.IIR数字滤波器的频率变换设计法的基本原理根据滤波器设计要求，设计模拟原型低通滤波器，然后进行频率变换，将其转换为相应的模拟滤波器(高通、带通等)，最后利用冲激响应不变法或双线性变换法，将模拟滤波器数字化成相应的数字滤波器。图4.39IIR数字滤波器MATLAB设计步骤流程图1．MATLAB的典型设计利用在MATLAB设计IIR数字滤波器可分以下几步来实现(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标；(2)根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc；(3)利用最小阶数N产生模拟低通滤波原型；(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器；(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。[例4-39]设计一个数字信号处理系统，它的采样率为Fs＝100Hz，希望在该系统中设计一个Butterworth型高通数字滤波器，使其通带中允许的最小衰减为0.5dB，阻带内的最小衰减为40dB，通带上限临界频率为30Hz，阻带下限临界频率为40Hz。MATLAB源程序设计如下：%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;[N,Wc]=buttord(wp,ws,rp,rs,'s');　%选择滤波器的最小阶数[Z,P,K]=buttap(N);　%创建Butterworth低通滤波器原型[A,B,C,D]=zp2ss(Z,P,K);%零－极点增益模型转换为状态空间模型[AT,BT,CT,DT]=lp2hp(A,B,C,D,Wc);　%实现低通向高通的转变[num1,den1]=ss2tf(AT,BT,CT,DT);%状态空间模型转换为传递函数模型%运用双线性变换法把模拟滤波器转换成数字滤波器[num2,den2]=bilinear(num1,den1,100);　[H,W]=freqz(num2,den2);　%求频率响应plot(W*Fs/(2*pi),abs(H));grid;%绘出频率响应曲线xlabel('频率/Hz');ylabel('幅值')程序运行结果如图4.40所示。2．MATLAB的直接设计图4.39IIR数字滤波器MATLAB设计步骤流程图[例4-41]试设计一个带阻IIR数字滤波器，其具体的要求是：通带的截止频率：wp1＝650Hz、wp2＝850Hz；阻带的截止频率：ws1＝700Hz、ws2＝800Hz；通带内的最大衰减为rp＝0.1dB；阻带内的最小衰减为rs＝50dB；采样频率为Fs＝2000Hz。MATLAB源程序设计如下：wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp=[wp1,wp2]/(Fs/2);ws=[ws1,ws2]/(Fs/2);　%利用Nyquist频率频率归一化[N,wc]=ellipord(wp,ws,rp,rs,'z');%求滤波器阶数[num,den]=ellip(N,rp,rs,wc,'stop');　　　　　%求滤波器传递函数[H,W]=freqz(num,den);　　　　　　　　　%绘出频率响应曲线plot(W*Fs/(2*pi),abs(H));grid;xlabel('频率/Hz');ylabel('幅值')该程序运行后的幅频响应曲线如图4.42所示。4.9FIR数字滤波器设计格式：w=boxcar(M)功能：返回M点矩形窗序列。MATLAB信号处理工具箱中的窗函数法设计FIR数字滤波器的专用命令fir1()。格式：B＝fir1(N,wc)功能：设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器，返回的向量B为滤波器的系数(单位冲激响应序列)，其长度为N+1。4.9.1窗函数设计法窗函数设计的基本原理：h(n)=w(n)hd(n)w(n)为窗函数，hd(n)理想数字滤波器的单位冲激响应。在MATLAB信号处理工具箱中为用户提供了Boxcar(矩形)、Bartlet(巴特利特)、Hanning(汉宁)等窗函数，这些窗函数的调用格式相同。FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称h(n)=h(N-n-1)或h(n)=-h(N-n-1)FIR数字滤波器具有线性相位：或[例4-43]用矩形窗设计线性相位FIR低通滤波器。该滤波器的通带截止频率wc=pi/4，单位脉冲响h(n)的长度M=21。并绘出h(n)及其幅度响应特性曲线。MATLAB源程序为：M=21;wc=pi/4;　　　　　　　%理想低通滤波器参数n=0:M-1;r=(M-1)/2;nr=n-r+eps*((n-r)==0);hdn=sin(wc*nr)/pi./nr;　%计算理想低通单位脉冲响应hd(n)ifrem(M,2)~=0,hdn(r+1)=wc/pi;end　%M为奇数时，处理n=r点的0/0型wn1=boxcar(M);%矩形窗hn1=hdn.*wn1';%加窗subplot(2,1,1);stem(n,hn1,'.');line([0,20],[0,0]);xlabel('n'),ylabel('h(n)'),title('矩形窗设计的h(n)');hw1=fft(hn1,512);w1=2*[0:511]/512;　　　%求频谱subplot(2,1,2),plot(w1,20*log10(abs(hw1)))xlabel('w/pi'),ylabel('幅度(dB)');title('幅度特性(dB)');程序运行结果如图4.44所示。4.9.2频率抽样法1.频率抽样法的基本原理对所期望的滤波器的频率响应Hd(ejw)进行等间隔采样获得H(k)，利用h(n)=IDFT[H(k)]求得FIR的单位冲激响应。2.MATLAB信号处理工具箱中的频率抽样法专用函数命令fir2()格式：B＝fir2(N,F,A)功能：设计一个N阶的FIR数字滤波器，其频率响应由向量F和A指定，滤波器的系数（单位冲激响应）返回在向量B中，长度为N+1。向量F和A分别指定滤波器的采样点的频率及其幅值，F中的频率必须在0.0到1.0之间，1.0对应于采样频率的一半。它们必须按递增的顺序从0.0开始到1.0为结束。　　[例4-47]试用频率抽样法设计一个FIR低通滤波器，该滤波器的截止频率为0.5pi，频率抽样点数为33。MATLAB源程序为：N=32;F=[0:1/32:1];%设置抽样点的频率，抽样频率必须含0和1。A=[ones(1,16),zeros(1,N-15)];%设置抽样点相应的幅值B=fir2(N,F,A);freqz(B);%绘制滤波器的幅相频曲线figure(2);stem(B,'.');%绘制单位冲激响应的实部line([0,35],[0,0]);xlabel('n');ylabel('h(n)');图4.49滤波器的频率响应和单位冲激响应序列4.9.3MATLAB的其它相关函数1．最小二乘逼近法设计线性相位FIR滤波器函数fircls()2．有限制条件的最小二乘逼近法设计低通和高通FIR数字滤波器函数fircls1()3．最小二乘逼近法设计线性相位FIR数字滤波器函数firls()4．升余弦FIR滤波器设计函数firrcos()5．Parks-McClellan优化等纹波FIR滤波器设计函数remez()