FIR数字滤波器的基本原理及设计方法.doc

第一章 FIR数字滤波器的基本原理及设计方法有限长单位脉冲响应数字滤波器（FIRDF，Finite Impulse Response Digital Filter）的最大优点是可以实现线性相位滤波。而IIRDF主要对幅频特性进行逼近，相频特性会存在不同程度非线性。我们知道，无失真传输与滤波处理的条件是，在信号的有效频谱范围内系统幅频响应为常数，相频响应具有线性相位。在数字通信和图像处理与传输等应用场合都要求滤波器具有线性相位特性。另外FIRDF是全零点滤波器，硬件和软件实现结构简单，不用考虑稳定性问题。所以，FIRDF是一种很重要的滤波器，在数字信号处理领域得到广泛应用。当幅频特性指标相同时，FIRDF的阶数比IIRDF高的多，但是同时考虑幅频特性指标和线性相位要求时，IIRDF要附加复杂的相位校正网络，而且难以实现严格线性相位特性。所以，在要求线性相位滤波的应用场合，一般都用FIRDF。FIRDF的设计方法主要有两类：第一类是基于逼近理想滤波器特性的方法，包括窗函数法，频率采样法和等波纹最佳逼近法。第二类是最优设计法，我们主要讨论第一类设计法，侧重与滤波器的设计方法和相应的MATLAB工具箱函数的介绍。 FIR数字滤波器的设计方法有窗函数法、频率采样法和基于firls函数和remez函数的最优化方法。MATLAB语言中的数字信号处理工具箱，提供了一些滤波器的函数，使FIR滤波器的运算更加方便和快捷。在MATLAB中提供的滤波函数有fir1(),此函数以经典的方法实现加窗线性相位FIR数字滤波器设计，可以设计出低通、高通、带通和带阻滤波器；fir2函数设计的FIR滤波器，其滤波的频率特性由矢量f和m决定，f和m分别为滤波器的期望幅频响应的频率相量和幅值相量。Firls()和remez()的基本格式用于设计I型和II型线性相位FIR滤波器，I型和II型的区别是偶函数还是奇函数。freqz()用于求数字滤波器的频率响应。并且提供了各种窗函数的函数，比如，hamming()是海明窗函数，hanning()是汉宁窗函数，kaiser()是凯泽窗函数，使在设计的过程中，不用自己重新设计窗函数。1.1 窗函数法设计FIR数字滤波器设我们所要设计的FIR滤波器的传输函数是(e),(n)是与其对应的单位脉冲响应，因此 （1-1） （1-2)如果我们能够在已知的情况下，求出，经过Z变换可得到滤波器的系统函数。通常情况下理想数字滤波器的单位脉冲相应是无限长的，且是非因果序列。获得有限脉冲响应滤波器的一种可能方法是对截取一段来近似代替，可是这样会改变原来的滤波器指标，出现吉布斯效应误差。窗函数法就是用被称为窗函数的有限加权序列w(n)来修正式(1)的傅里叶基数以求得要求的有限脉冲响应序列，即 （1-3）w(n)是有限长序列，当nN-1时，w(n)=0。这种方法的重点在于选择某种合适的窗函数。要求窗函数主瓣宽度尽可能窄，以获得最小的过渡带；旁瓣相对值尽可能小，以使得通带波纹小，并且阻带衰减大。下面介绍几种常用的窗函数：1.矩形窗(Rectangle Window) （3-4）其频率函数为： （3-5）2.三角形窗(Bartlett Window) （3-6）其频率函数为： （3-7)3.汉宁(Hanning)窗，又称升余弦窗 （3-8）利用傅里叶变换得到频率函数为： （3-9）当时，所以窗函数的幅度函数为 (3-10)4.汉明(Hamming)窗，又称改进的升余弦窗 （3-11）其幅度函数为： （3-12）5.布莱克曼(Blankman)窗，又称二阶升余弦窗 （3-13)其幅度函数为： (3-14）6.凯泽(Kaiser)窗 （3-15）其中：是一个可自由选择的参数，I0( x)是第一类修正零阶贝塞尔函数10.上述窗函数的基本参数如下表窗函数旁瓣峰值幅度/db过渡带宽阻带最小衰减/db矩形窗-134/N-21三角形窗-268/N-25汉宁窗-318/N-44汉明窗-408/N-53布莱克曼窗-5712/N-74凯泽窗-5710/N-80窗函数法设计滤波器的步骤：1）根据技术要求确定待求滤波器的单位取样响应。2）根据对过渡带和阻带衰减的要求，选择窗函数的形式，并估计窗口长度N。3）计算滤波器的单位取样响应h(n)： （3-16)式中,是前面所选择好的窗函数。4）检验技术指标是否满足要求。根据下式计算： (3-17)如果不满足要求，根据具体情况重复步骤(2)(3)(4)步，直到满足要求为止。本文以一个FIR滤波器的设计为例说明如何使用MATLAB设计数字滤波器设计实例：用窗函数法设计线性相位FIR低通数字滤波器，要求通带截止频率Wp=0.4*,阻带截止频率Ws=0.5*, 通带衰减不大于3db，阻带衰减不小于40db。程序如下：Wp=0.4*pi;Ws=0.5*pi;Wdel=Ws-Wp;N=ceil(8*pi/Wdel);Wn=(0.4+0.5)*pi/2;window=hanning(N+1);b=fir1(N,Wn/pi,window);freqz(b,1,512)程序执行后得幅频和相频如下图所示： 图1.11.2 频率采样法设计FIR数字滤波器1.对理想滤波器的系统函数Hd(z)进行频率采样以得到系统的理想频响Hd(ejw)的等间隔采样值H(k)。H(k)实际上是所要求的滤波器的单位采样响应(h(n)的离散傅里叶变换（DFT），如下试： (3-18) (3-19)为了减小H（k）的通带边缘由于抽样点的变化而引起的起伏振荡，可以增加过渡点，加宽过渡带以减小通带的起伏。每一个抽样值产生一个与sin()/sin()成正比，并位移（2k）/N的频率响应，而H（k）与内插函数的线性组合就是FIR滤波器的频率响应，增加一点过渡可以使阻带衰减提高到-4454dB,二点过渡衰减-6575,三点过渡衰减-8595dB.如果不能使过渡带太宽，同时要求增大阻带衰减，可以增加取样点数N，但这样会增加计算量、延时和误差。频率取样型FIR滤波器设计步骤：(1)给定理想滤波器频率响应Hd(ejw)。(2)根据过渡带宽和阻带衰减确定过渡点数和h(n)的长度N。 (3-20)(3)由IFFT计算IDFT得到： (3-21)设计实例：率采样法设计一个带通滤波器，满足：低阻带边缘：w1s=0.2*;低通带边缘：w1p=0.35*;高通带边缘：w2p=0.65*;高阻带边缘：w2s=0.8*。设计过渡带中的频率样本值为t1和t2，取t1=0.109021,t2=0.59417456。设计程序如下：M=40;al=(M-1)/2;l=0:M-1;t1=0.109021;t2=0.59417456;Hrs=zeros(1,5),t1,t2,ones(1,7),t2,t1,zeros(1,9),t1,t2,ones(1,7),t2,t1,zeros(1,4);k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1;angh=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angh);h=real(ifft(H,M);freqz(h,1,512,1000)实验得幅频相频特性如下图所示： 图1.21.3等波纹最优化方法设计FIR数字滤波器在数字信号处理中, 利用数字滤波器可改变信号中所含频率分量的相对比例或滤除某些频率分量, 使其达到所需要的效果. 其中数字FIR 滤波器由于具有精确的线性相位, 且系统稳定, 所以广泛应用于通信、 数字图象处理、 语音信号处理、 自适应处理、 雷达/ 声纳系统等方面. 目前, FIR 滤波器 1, 2的设计主要有窗函数设计法和频率采样设计法. 但是, 这2 种方法都不易精确控制通带边界频率Wp与阻带边界频率Ws ,所以, 在实际应用中具有一定的局限性.而以最大误差最小化准则支持的切比雪夫逼近法是一种优异的设计方法, 易于精确控制Wp与Ws.与窗函数和频率采样法比较，由于这种设计法使最大误差均匀化，所以设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使阶数最低。等波纹最佳逼近法设计的数学证明复杂，已超出本科生的数学基础。所以我们略去等波纹最佳逼近法复杂的数学推导，只介绍基本思想和实现线性相位FIRDF的等波纹最佳逼近设计的MATLAB信号处理工具箱函数remez和remezord。Remez函数采用数值分析中的remez多重迭代算法求解等波纹最佳逼近问题，求得满足等波纹最佳逼近准则的FIRDF的单位脉冲相应h(n),由于切比雪夫和雷米兹(remez)对解决该问题做出了贡献，所以又称之为切比雪夫逼近法或者雷米兹逼近法。第二章 设计步骤等波纹滤波器的最优化设计方法主要有2 种,第 1 种是离散最小二乘法. 它的思路是在给定的一些离散点上,使实际的幅频特性和理想幅频特性之间的误差的平方和为最小. 第 2 种是最优化等波纹设计法,也称为雷米兹法或切比雪夫逼近法. 该类型滤波器幅频特性在通带和阻带上的误差峰值是均匀分布的,其误差具有等波纹特性,因而把波纹的幅度控制到最小, 或在同等指标下减小它的阶次. 第 1 种方法是连续最小的平方法的推广, 容易理解,但它的指标与滤波器没有直接关联,误差平方小的滤波器不能保证没有窄而大的波纹出现,像吉布斯效应 那样. 第 2 种方法直接控制通带波动和阻带衰减, 最具针对性,是滤波器的最优化设计方法. 因此,采用 MAT LAB信号处理工具箱提供的函数,运用最优化等波纹设计法实现数字 FIR 滤波器的设计和仿真.。完整的最优化等波纹滤波器设计, 除了切比雪夫等波纹逼近公式外,还要考虑:( 1)滤波器采样点数n 的确定.( 2)极值数目的确定. 最优化等波纹滤波器的误差函数在给定的频率上有( L + 2)或( L + 3)个极值,L 为多项式的阶数. 对于某些Wp, Ws 的组合,可能得到有( L + 3)个极值的滤波器. 此时, L = f loor ( ( n-1) / 2) .( 3)建立频率修正的算法. 在程序中自动进行反复的迭代修正,直到达到要求的精度为止.交替定理能保证切比雪夫逼近问题的解存在并且惟一, 但它并没有说如何得到该解, 既不知 n(或L ) ,也不知极值的频率和波纹系数W. ParksMcClellan提供了利用Remez 交换算法导出的迭代算法. 它假设已知滤波器的n,并且设1= 2 , 就可得到解,其中, 为最大波纹系数, 2 为阻带波动系数.和n 是相关的, n 越大,越小. 滤波器的技术指标中给出了通带波动系数 1和2 , W p , W s ,因此需要设定n. 凯泽提出逼近 n 的公式: Parks 和McClellan 算法首先猜设( L + 2)个极值频率W i ,估计极值频率上的最大误差及其 Wi , 由这些新的频率点拟合出一个新的L 阶多项式, 重复以上过程, 直到找出最优集Wi和全局最大波纹系数1为止. 迭代结果得到多项式系数d(n) ,最后算出滤波器脉冲响应 h(n) .由于n是近似的, 1可能不等于2 . 若12 ,则增加n;若12 ,则减小n,然后再次用Remez 算法确定一个新的.重复此过程, 直至12,这样就可得到等波纹滤波器.性能分析比较分别采用窗口法、频率取样法和优化设计法设计双带滤波器，其指标要求为： 1）第一带通滤波器 fp1=4 kHz, fp2 =5 kHz， fr1=3.6 kHz， fr2=5.4 kHz；2）第二带通滤波器 fp3=6 kHz，fp4=7 kHz，fr3=5.6 kHz，fr4=7.4 Hz；3）通带衰耗：Ap40 dB，模拟取样频率 fs=20 kHz。 MATLAB 模拟结果如下。 图1.3 图1.4 图1.5比较以上三种滤波器的设计方法，在同样的阶数下，优化设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。频率取样法设计过程简单，但阻带衰减明显，若适当选取过渡带样点值，会取得较窗口法略好的衰耗特性。窗口设计法在阶数较低时，阻带特性不满足设计要求，只有当滤波器阶数较高时，使用哈明窗和凯塞窗基本可以达到阻带衰耗要求。FIR 滤波器以其稳定和容易实现严格的线性相位，使信号处理后不产生相位畸变，而在实际中获得广泛的应用。第三章MATLAB的实现3.1MATLAB函数介绍在介绍MATLAB工具箱函数remezord和remez之前，先介绍等波纹滤波器的技术指标及其描述参数，等波纹滤波器的幅频特性函数曲线如图7.4.2所示。图7.4.2给出了等波纹滤波器技术指标的两种描述参数。图7.4.2(a)用衰减的dB数描述，即p=2，p=2dB, s=1120, s=20dB.这是工程实际中常用的指标描述方法。但是，用等波纹最佳逼近设计法求滤波器阶数N和误差加权函数W()时，要求给出滤波器通带和阻带的振荡波纹幅度1和2.图7.4.2(b)给出了用通带和阻带的振荡波纹幅度1和2描述的技术指标。显然，两种描述参数之间可以换算。如果设计指标以p和s给出，为了调用MATLAB工具箱函数进行设计，就必须由p和s换算出通带和阻带的振荡波纹幅度1和2。对比图7.4.2(a)和图7.4.2(b)得出关系式：p=-20lg（1-1）（1+1）=20lg（1+1）（1-1） (7.4.2)s=-20lg2（1+1） -20lg2 (7.4.3)由式(7.4.2)和(7.4.3)得到1=(10p/20-1)/(10p/20+1) (7.4.4)2=10(-s)/20 (7.4.5)按照式和(7.4.4)和(7.4.5)计算得到图7.4.1（b）中的参数：1=0.1146，2=0.1.实际中1和2一般很小，这里是为了观察等波纹特性及参数1和2的含义，特意取较大值。1 Remez函数采用remez算法实现线性相位FIR数字滤波器的等波纹最佳逼近设计。其调用格式为：hn=remez(M，f,m,w)返回单位脉冲响应响亮hn。调用参数含义如下：M：FIRDF阶数，hn长度为N=M+1。f和m：给出希望逼近的幅度特性。f为边界频率向量，0f1，要求为f为单调增向量（即f(k)f（k+1），而且从0开始，以1结束，1对应数字频率=（或模拟频率Fs2，Fs表示时域采样频率）。m是与f对应的幅度向量，m与f长度相等，m（k）表示频点f(k)的频响应值。如果用命令Plot（f,m）画出幅频响应曲线，则k为奇数时，频段f(k),f(k+1)上的幅频响应就是期望逼近的幅频响应值，频段f（k+1），f（k+2）为无关区。简言之，Plot（f，m）命令画出的幅频响应就是期望逼近的幅频响应值，频段f（k+1），f（k+2）为无关区。简言之，Plot(f,m)命令画出的幅频响应曲线中，起始频段为第一段，奇数频段为逼近区，偶数频段为无关区。W/g10除了设计选频FIRDF，remez函数还可以设计两种特殊滤波器，调用格式如下。设计希尔伯特变换器的调用格式为：hn=remez（M，f,m,w,hilbert）hn具有奇对称特性：hn(n)= hn（M+1-n）。在边界频率向量f给定的通带上希望幅度逼近常数1.一般f=a,b，m=1,1,且f中第一个边界频率a不能为0，必须满足0ab1.这是因为希尔伯特变换器属于第二类线性相位滤波器，且滤波器长度N为奇数，所以只能实现带通滤波器。设计数字微分器调用格式为：hn=remez（M,f,m,w,defferentiator）下一节将介绍希尔伯特变换器和数字微分器的设计和应用。Remez函数的调用参数（M，f,m,w），一般都是调用remezord函数来计算。2.remezord功能：根据逼近指标估算等波纹最佳逼近FIRDF的最低阶数m、误差加权向量M和归一化边界频率向量f，从而是滤波器在满足指标的前提下造价最低。其返回参数作为remez函数的调用参数。其调用格式为： M,fo,mo,w=remezord（f,m,rip,Fs） 参数说明如下。 返回参数作为remez函数的调用参数，设计的滤波器可以满足由参数f,m,rip和Fs买哦数的指标。f与remez中类似，这里f可以是模拟频率（Hz）或归一化数字频率，但必须以0开始，以Fs2（用归一化频率时对应1）结束，而且其中省略了0和Fs2两个频点。Fs为采样频率，缺省时默认Fs=2Hz。但是这里f的长度（包括省略的0和Fs2两个频点）是m的两倍，即m中的每个元素表示f给定的一个逼近频段上希望逼近的幅度值。例如，对图7.4.4，f=14，516；m=1，0。 图3.1 注意：省略Fs时，f中必须为归一化频率。有事估算的阶数M略小，使设计结果达不到指标要求，这时要取M+1或M+2（必须注意对滤波器长度N=M+1的奇偶性要求）。所以必须检验设计结果。如果无关区（过渡带）太窄，或截止频率接近零频率和FS2时，设计结果不正确。 rip表示f和m描述的各逼近频段允许的波纹振幅（幅频响应最大偏差），f的长度是rip的两倍。 一般以N,fo,mo,w=remezord（f,m,rip,Fs）返回的参数作为remez的调用参数，计算单位脉冲响应：hn=remez（N，fo,mo,w），对比前面介绍的remez调用参数，可清楚地看出remezord返回参数N，fo,mo和w的含义。 综上所述，调用remez和remezord函数设计线性相位FIRDF，关键是根据设计指标求出remezord函数的调用参数f,m,rip和Fs，其中Fs一般是题目给定的，或根据实际信号处理要求（按照采样定理）确定。 下面给出由给定的各种滤波器设计指标确定remezord调用参数f,m,和rip的公式，编程时直接诶套用即可。 低通滤波器设计指标逼近通带：0，p，通带最大衰减：pdB；逼近阻带：s，阻带最小衰减：sdB。remezord调用参数：f=p，s，m= 1，0，rip= 1，2 （7.4.6）其中，f向量省去了起点频率0和终点频率1，1和2分别为通带和阻带波纹幅度，由式（7.4.4）和式（7.4.5）计算得到，下面相同。 高通滤波器设计指标逼近通带：p，通带最大衰减：pdB；逼近阻带：0，s，阻带最小衰减：sdB。remezord调用参数： f= s，p，m= 0，1，rip= 2，1 (7.4.7) 带通滤波器设计指标逼近通带：p1，pu，通带最大衰减：pdB；逼近阻带：0，s1，su,，阻带最小衰减：sdB。remezord调用参数： f= s1，p1，pu，us，m= 0，1,0，rip= 2，1, 2 (7.4.8) 带阻滤波器设计指标逼近阻带：s1，su，阻带最大衰减：sdB；逼近阻带：0，p1，pu,，通带最小衰减：pdB。remezord调用参数： f= p1，s1，su，pu，m= 1，0,1，rip= 1, 2, 1 (7.4.9)3.2设计举例1利用等波纹最佳逼近法设计FIR带通滤波器，其指标如下：阻带下截止频率s=0.2，通带下截止频率p=0.35，通带上截止频率p=0.65, 阻带上截止频率s=0.8,通带最大衰减p=1dB,阻带最小衰减s=60dB。解：调用remezord和remez函数求解。由调用格式知道，首先要根据设计指标确定remezord函数的调用参数，在直接编写程序调用remezord和remez函数设计得到h(n)。将设计指标带入式（7.4.8）得到remezord函数的调用参数。本例设计程序为ep741.m。%ep741.m:例7.4.1 用remez函数设计带通滤波器f=0.2,0.35,0.65,0.8; %省去了0，1m=0,1,0; rp=1;rs=60;%由式(7.4.4)和式(7.4.5)求通带和阻带波纹幅度dat1和dat2Dat1=(10(rp20)-1) (10(rp20)+1);Dat2=10(-rs20);Rip=dat2,dat1,dat2;M,fo,mo,w=remezord(f,m,rip);M=M+2;%remezord估算的阶数偏小，加2才满足要求Hn=remez(M,fo,mo,w);绘图部分略去程序运行结果：M=28.即h(n)的长度N=29。h(n)及其耗损函数曲线如图7.4.5所示。而用窗函数设计的FIR滤波器的阶数为80，由此例可见，等波纹最佳逼近法可以使滤波器的阶数大大降低。图3.2图3.3 2 FIR希尔伯特变换器3 例 利用等波纹最近逼近法设计14阶FIR希尔伯特变换器，要求通带为0.1,0.9.调用remez函数的设计程序为ep744.m. %ep744.m:例7.4.4 调用remez函数设计FIR希