数字信号处理结课论文

基于MATLAB的数字滤波器设计摘要数字滤波器的实现是数字信号处理中的重要组成部分，设计过程较为复杂，牵涉到模型逼近、指标选择、计算机仿真、性能分析及可行性分析等一系列的工作，本文从设计原理以及数学软件matlab出发阐述数字滤波器的设计原理与方法。关键词：MATLAB，数字信号处理，数字滤波器第1页1绪论数字滤波器是数字信号处理的重要应用，21世纪数字滤波器及其相关技术广泛的应用于通信、电子、自动控制等领域，是一种有效的抑制噪音、提取有用信号的方法。它本身可以用硬件实现，也可以通过软件来实现，还可以通过专用的DSP处理器配合相应的软件，即软硬结合的方式来实现。数字滤波器可以分为有限冲击响应（FIR）和无限冲激响应(IIR)两种。通过MATLAB程序，实现输入相应技术指标及滤波器模型，输出相应数字滤波器的参数的功能。2无限长单位脉冲响应滤波器（IIR）的设计根据模拟滤波器设计数字滤波器，就是通过已知的模拟滤波器系统的系统函数H(s)来设计数字滤波器的系统函数H(z)，主要是通过脉冲响应不变法，或双线性变换法完成S平面到Z平面的转换。通过典型的模拟滤波器(诸如：巴特沃斯滤波器、切比雪夫滤波器等)可以实现一定参数要求的数字滤波器。根据已有的数字滤波器设计不同参数或者不同频带通断类型的数字滤波器。例如已知数字低通滤波器的模型，通过变量代换得到不同截止频率的数字低通滤波器，或通过已知低通滤波器的模型设计高通、高阻、带通、带阻滤波器，这里主要是通过来完成相应的变量代换来实现滤波器类型的变换和参数的变换。3有限长单位脉冲响应滤波器(FIR)的设计IIR滤波器可用于较少的阶数达到所要求的幅度特性，且实现时所需的运算次数及存储单元都很少，十分适合于对于相位特性没有严格要求的场合，如果对相位特性有要求，这时选用FIR滤波器较好。3.1窗函数法从时域出发，把理想的窗口函数hd(n)截取成有限长的，以此h(n)来逼近理想的窗口函数hd(n)，从而频率响应H(jw)也近似于理想的频率响应Hd(jw)。主要窗函数有：矩形窗，汉宁(Hanning)窗，汉明(Hamming)窗，布莱克曼(Blackman)窗，凯塞(Kaiser)窗。3.2频率采样法第2页从频率出发，对理想的频率响应Hd(ejw)加以等间隔采样，以此Hd(k)实际FIR滤波器的频率特性的离散样本H(k)。4利用MATLAB具体实现数字滤波器4.1MATLAB中IIR数字滤波器设计的设计函数信号处理工具箱提供的IIR经典设计方法是基于经典的低通滤波器到具有相同性能指标的数字滤波器的变换。基本原理就是先根据滤波器的技术指标设计出相应的模拟滤波器，然后再将设计好的模拟滤波器变换为满足指标的数字滤波器。利用表1的设计函数，可以很容易地产生任何阶数的高通、低通、带通、带阻数字滤波器。表1滤波器设计函数滤波器类型设计函数Bessel（用于模拟滤波器）[b,a]=besself(n,Wn,options)[z,p,k]=besself(n,Wn,options)[A,B,C,D]=besself(n,Wn,options)Butterworth[b,a]=butter(n,Wn,options)[z,p,k]=butter(n,Wn,options)[A,B,C,D]=butter(n,Wn,options)ChebyshevTypeI[b,a]=cheby1(n,Wn,options)[z,p,k]=cheby1(n,Wn,options)第3页[A,B,C,D]=cheby1(n,Wn,options)ChebyshevTypeII[b,a]=cheby2(n,Wn,options)[z,p,k]=cheby2(n,Wn,options)[A,B,C,D]=cheby2(n,Wn,options)椭圆滤波器[b,a]=ellip(n,Wn,options)[z,p,k]=ellip(n,Wn,options)[A,B,C,D]=ellip(n,Wn,options)表中，n为滤波器阶次，Wn为滤波器的归一化截止频率；函数默认为低通或带通滤波器。b,a分别为滤波器传递函数的分子和分母的系数向量；z,p,k分别为滤波器的零点、极点和增益。options为滤波器类型参数，high为高通滤波器，截止频率为Wn；stop为带阻滤波器，截止频率为Wn=[W1,W2]。4.2脉冲响应不变法设计数字滤波器调用格式：[bz,az]=impinvar(b,a,Fs)，再给定模拟滤波器参数b,a和取样频率Fs的前提下，计算数字滤波器的参数。两者的冲激响应不变，即模拟滤波器的冲激响应按Fs取样后等同于数字滤波器的冲激响应。4.3利用双线性变换法设计数字ButterWorth滤波器调用格式：[bz,az]=bilinear[b,a,Fs]，根据给定的分子b,分母系数a和取样频率Fs，根据双线性变换将模拟滤波器变换成离散滤波器，具有分子系数向量bz和分母向量az。4.4IIR设计实例利用脉冲响应不变法设计一个巴特沃斯数字低通滤波器，要求：ωp=0.25π，Rp=1dB，第4页ωs=0.4π，As=15dB，滤波器采样频率Fs=2000Hz。程序如下：wp=0.25*pi;ws=0.4*pi;Rp=1;As=15;ripple=10^(-Rp/20);%滤波器的通带衰减对应的幅度值Attn=10^(-As/20);%滤波器的阻带衰减对应的幅度值%转换为模拟滤波器的技术指标Fs=2000;T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs;%模拟原型滤波器的计算[n,Omgc]=buttord(Omgp,Omgs,Rp,As,’s’)%计算阶数n和截止频率[z0,p0,k0]=buttap(n);%设计归一化的巴特沃斯模拟滤波器原型ba1=k0*real(poly(z0));%求原型滤波器的系数baa1=real(poly(p0));%求原型滤波器的系数a[ba,aa]=lp2lp(ba1,aa1,Omgc);%变换为模拟低通滤波器%用脉冲响应不变法计算数字滤波器的系数[bd,ad]=impinvar(ba,aa,Fs)[C,B,A]=dir2par(bd,ad)[H,w]=freqz(bd,ad);dbH=20*log10((abs(H)+eps)/max(abs(H));subplot(2,2,1);plot(w/pi,abs(H));ylable(‘|H|’);title(‘幅度响应’);axis([0,1,0,1.1]);set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);set(gca,’YTickMode’,’manual’,’YTick’,[0,Attn,ripple,1]);gridsubplot(2,2,2);plot(w/pi,dbH);title(‘幅度响应(dB)’)；ylabel(‘dB’);xlable(‘频率’);axis([0,1,-40,5]);set(gca,’XTickMode’,’manual’,’XTick’,[0,0.25,0.4,1]);set(gca,’YTickMode’,’manual’,’YTick’,[-50,-15,-1,0]);gridsubplot(2,2,4);zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title(‘零极点图’);自定义dir2par函数：functionI=cplxcomp(p1,p2)I=[];forj=1:length(p2)fori=1:length(p1)if(abs(p1(i)-p2(j))0.0001)I=[I,i];endendendI=I';function[C,B,A]=dir2par(num,den)M=length(num);第5页N=length(den);[r1,p1,C]=residuez(num,den);p=cplxpair(p1,10000000*eps);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);ifK*2==N;fori=1:2:N-2;Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix(i+1)/2,:)=real(Brow);A(fix(i+1)/2,:)=real(Arow);end[Brow,Arow]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(Brow),0];A(K,:)=[real(Arow),0];elsefori=1:2:N-1;Brow=r(i:1:i+1,:);Arow=p(i:1:i+1,:);[Brow,Arow]=residuez(Brow,Arow,[]);B(fix(i+1)/2,:)=real(Brow);A(fix(i+1)/2,:)=real(Arow);endend程序运行结果如下：n=5Omgc=103.2016bd=0.00720.03620.07250.07250.03620.0072ad=1.0000-1.94341.9680-1.07020.3166-0.0392sos=1.00000.995601.0000-0.319301.00002.00721.00721.0000-0.69840.20531.00001.99720.99731.0000-0.92570.5976g=0.0072频率特性如图1所示：第6页00.250.4100.17780.89131|H|幅度响应00.250.41-101相位响应00.250.41-15-10幅度响应(dB)dB频率()-101-1-0.500.51RealPartImaginaryPart零极点图图1频率特性曲线4.5MATLAB中FIR数字滤波器的设计的函数在MATLAB中产生窗函数十分简单：根据长度n产生一个相应类型窗函数的类型：（1）矩形窗(RectangleWindow)：调用格式：w=boxcar(n)。（2）三角窗(TriangularWindow)：调用格式：w=triang(n)。（3）汉宁窗(HanningWindow)：调用格式：w=hanning(n)。（4）哈明窗(HammingWindow)：调用格式：w=hamming(n)。（5）布莱克曼窗(BlackmanWindow)：调用格式：w=blackman(n)。（6）凯塞窗(KaiserWindow)：调用格式：w=kaiser(n,beta)。基于窗函数的FIR滤波器设计利用MATLAB提供的函数firl来实现。调用格式：firl(n,Wn,’ftype’,Window)，n为阶数，Wn是截止频率（如果输入是形如[W1W2]的矢量时，本函数将设计带通滤波器，其通带为W1ωW2），ftype第7页是滤波器的类型(低通省略该参数，高通high，带阻stop)，Window是窗函数。格式中：n为FIR滤波器的阶数，对于高通、带阻滤波器，n取偶数。Wn为滤波器的截止频率，0~1；对于带通、带阻滤波器，Wn=[W1,W2]，且W1W2,对于多带滤波器Wn=[W1,W2,W3,W4],频率分段0WW1，W1WW2,W2WW3,...‘ftype’为滤波器类型，默认为低通滤波器或带通滤波器。其选择如下：‘high’为高通滤波器，’stop’为带阻滤波器，’DC-1’使多带的第一频带为通带，’DC-0’使多带的第一频带为阻带。Window为窗函数，列向量，其长度为n+1。4.6FIR设计实例用矩形窗设计一个FIR数字低通滤波器，要求：N=64，截止频率wc=0.4，描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。程序如下：wc=0.4*pi;N=64;n=0:N-1;hd=ideal_