FIR数字滤波器设计与实现

1、吉林大学珠海学院FIR数字滤波器设计与实现一.摘要：数字滤波器是一种具有频率选择性的离散线性系统，在信号数字处理中有着广泛的应用。其中FIR滤波器是一种常用的滤波器，它在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性，在语音分析、图像处理、雷达监测等对信号相位要求高的领域有着广泛的应用，能实现IIR滤波器不能实现的许多功能。二.关键词：FIR窗函数系统函数MATLAB三.内容提要：数字滤波器的功能就是把输入序列通过一定的运算变换成输出序列，因此数字滤波器的结构系统中就必须包括一定数量和性能的运算器件和运算单元，而运算器件和运算单元的配置必须由数字滤波器的结构特点和性能特点来决定，

2、因此在进行FIR数字滤波器的设计之前，有必要介绍和总结FIR数字滤波器的基本结构和相关特性（包括频响曲线（幅度和相位），单位冲激响应等），在介绍完其基本结构和相关特性后，就进行FIR数字滤波器的设计和实现。（一）FIR滤波器的基本结构在讨论任何一种滤波器时，都要着重分析其系统函数，FIR滤波器的系统函数为：N1H（z）h（n）zn0从该系统函数可看出，FIR滤波器有以下特点：n01）系统的单位冲激响应h（n）在有限个n值处不为零；2）系统函数H（z）在|z|>0处收敛，极点全部在z=0处（稳定系统）；3）结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含

3、有反馈的递归部分。1 .FIR滤波器实现的基本结构有：1） 横截型（卷积型、直接型）a.一般FIR滤波器的横截型（直接型、卷积型）结构：AT-1若给定差分方程为：双冷=£起口）工-。则可以直接由差分方程得出FIR滤波器结构如M-0下图所示：这就是FIR滤波器的横截型结构，又称直接型或卷积型结构。b.线性相位FIR滤波器的横截型结构页脚内容1吉林大学珠海学院若h(n)呈现对称特性，即此FIR滤波器具有线性相位，则可以简化成横截型结构，下面分情况讨论：N为奇数时线性相位FIR滤波器实现结构如图所示：N为偶数时线性相位FIR滤波器实现结构如图所示我们知道IIR滤波器的优点是可利用模拟滤波器

4、设计的结果，缺点是相位是非线性的，若需要线性相位，则要用全通网络进行校正，比较麻烦，而FIR滤波器的优点是可以方便地实现线性相位。2)、级联型将H(z)分解为若干个实系数一阶或二阶因子相乘:实现结构如下图所示：xkI斓0单中仔2L千外延时器，2L+1竹7法什2L=M令忡*器必k分行号)及苛构4可以":*级增的，点:z14每济墓本节揶制，对考此便市拄制滤波器的传输零点02L系数比直接型多，所需的乘法运算多。3)频率取样型若FIR滤波器的冲激响应为有限长(N点)序列h(n),则有如图所示的关系:因此，对h(n)可以利用DFT得到H(k),然后利用内插公式：来表示系统函数，这就为FIR滤波

5、器提供了另外一种结构：频率抽样结构，1HG这种结构由两部分级联而成：分析系统函数二一二二其中级联的第一部分为：丁-这是一个梳状滤波器，它滤掉了频率分二名及其各次谐波。N级联的第二部分为N个一阶网络并联而成，第k个一阶网络为：2L=M个或对%+%+1塔|：法南二/乂¥力嬴器特点：可以分别控制每个子系统的零点页脚内容2吉林大学珠海学院这是一个谐振频率国二红的无损耗谐振器。这个谐振器的极点正好与梳状滤波器的一个零点(i=k)相抵消，从而使这个频率上的频率响应等于H(k)o这样，N个谐振器的N个极点就和梳状滤波器的N个零点相抵消，从而在N个频率抽样点上的频率响应就分别等于N个H(k)值。有上

6、叙的理论分析基础可以得到FIR滤波器的频率抽样结构。FIR滤波器的频率抽样结构如图所示：频率抽样结构的特点是它的系数H(k)就是滤波器在可二军处的响应，因此控制滤波器的频率响应很方便。频率抽样结构存在问题的问题是：在有限长情况下，系数量化后极点不能和零点抵消，使FIR系统不稳定。解决方法：在r圆上进行(r<1但近似等于1)取样，即用rz1代替z1,使极点和相应的零点移到单位圆内。(a)当N为偶数时的频率取样型结构如图所示。(b)当N为奇数时频率抽样型结构如图所示。4)快速卷积结构若FIR滤波器的单位冲激响应h(n)是一个N1点有限长序列，输入x(n)是一个N2点有限长序列，那么输出y(n

7、)是x(n)与h(n)的线性卷积，它是一个L=N1+N2-1点的有限长序列。我们知道，将x(n)补上LN2个零值点，将h(n)补上LN1个零值点，然后进行L点圆周卷积，就可以代替原x(n)与h(n)的线性卷积。而圆周卷积可以用DFTffiIDFT的方法来计算，这样我们得到FIR滤波器的快速卷积结构：这里DFT和IDFT都将采用快速傅里叶变换算法，当N1和N2足够长时，比直接计算线性卷积要快得多。2 .线性相位FIR滤波器的特点从以上的讨论中可以看出，我们最感兴趣的是具有线性特性的FIR滤波器，因此在设计FIR滤波器时，需要着重研究线性相位FIR滤波器的特点和性质，在上述已经介绍了线性相位FIR

8、滤波器的横截型结构，现在介绍它的频响特性。FIR滤波器的单位冲激响应h(n)是有限长的(0&n&N-1),其Z变换为：其傅立叶变换为：页脚内容3吉林大学珠海学院其中H()是幅度函数，是一个纯实数，可正可负，9()是相位函数。可以证明，线性相位FIR滤波器的冲激响应满足对称条件：h(n)=±h(N-1-n)和H(z)z(N1)H(z1)(1)、线性相位FIR滤波器的幅度函数和相位函数：(a)当h(n)是偶对称时，其幅度函数和相位函数分别为：特点：幅度函数H()包括正负值，相位函数是严格线性相位，滤波器有(N-1)/2个抽样周期的延时，它等于单位抽样响应h(n)长度N的一

9、半。(b)当h(n)是奇对称时，其幅度函数和相位函数分别为：特点：相位函数是严格线性相位，但在零频率(必=0)处有冗/2的相移。仍有(N-1)个抽样周期的延时。因此当h(n)为奇对称时，FIR滤波器将是一个具有准确相位的正交变换网络。(2)、FIR滤波器的线性相位特性FIR滤波器的线性相位特性如图所示。(3)、任何一种线性相位FIR滤波器的群延时都为：(4)FIR滤波器幅度函数的特点分四种情况分别讨论H()的特点：(a)当h(n)偶对称，N为奇数时：幅度函数的特点：H()对必=0,宿，之龙呈偶对称。(b)当h(n)偶对称，N为偶数时：幅度函数的特点：当a二"时，巴值)二°,

10、在z=-1处有一个零点，区(切对毋="是奇对称；如果滤波器在的=癖处幅度不为零(如高通滤波器)，则不能用这种滤波器。(c)当h(n)奇对称，N为奇数时幅度函数的特点：页脚内容4吉林大学珠海学院H()在必=0,比，2加处都为零，也就是H(z)在Z二±1处为零；H()对必=0,比，2笈都成奇对称。(d)当h(n)奇对称，N为偶数时：幅度函数的特点：H()在必=0,2充处为零,即H(z)在z=1处为零点；H(co)对3=0,2结呈奇对称，对3=既呈偶对称。(5)、零点位置：线性相位FIR滤波器的系统函数有以下关系：可见，若zZi是H(z)的零点，则z1/Zi也一定是H(z)的零点

11、。又由于当h(n)是实数时，H(z)的零点必成共腕对出现，所以zzi及z1/乙也一定是H(z)的零点。因而线性相位FIR滤波器的零点必是互为倒数的共腕对。其有四种可能性：(1)巴既不在实轴上，也不在单位园上，则是互为倒数的两组共腕对。(2)与不在实轴上，但是在单位园上，则共腕对的倒数是它们本身，故只有一组共腕对。(3)/在实轴上而不在单位园上，只有倒数部分，无复共腕部分。(4) 4既在实轴上又在单位园上，有两种可能，z=1或z=1。(二)FIR数字滤波器的设计在介绍和总结完FIR数字滤波器的基本结构和相关特性(包括频响曲线(幅度和相位)，单位冲激响应等)后，就是FIR数字滤波器的设计和实现，F

12、IR数字滤波器的设计步骤有：1 .技术要求(预期性能指标)页脚内容5吉林大学珠海学院技术要求是由实际用途决定的，它可由理想滤波器的系统函数Hd(z)、脉冲响应h(n)和差分方程描述。2 .逼近(近似)在数字滤波器性能分析的基础上，利用已经学过的概念和数学知识提供数字滤波器的表述，它是理想滤波器的一种近似。3 .实现上面一步的结果是一个滤波器的表述，它既可能是一个系统函数、也可能为差分方程，或者是单位脉冲响应h(n),依据这个结果进行数字滤波器结构的实际和软硬件的实现。下面重点介绍目前最主要的三种FIR数字滤波器的设计方法：(1) .窗函数设计法(时间窗口法)这种方法也称为傅立叶级数法。其设计是

13、在时域进行的，先用傅氏反变换求出理想数字滤波器的单位抽样响应hd(n),然后时域移位并加时间窗w(n)对其截断，从而求得FIR滤波器的单位抽样响应h(n)在设计过程中，将无限长序列变为有限长是通过时域加矩形窗乘积实现数据的截断的。时域乘积对应了频域卷积，从而对频响特征发生的改变。常见的窗函数有：矩形窗、三角形(Bartlertt)窗、汉宁(Hanning)窗。海明(Hamming)窗、布拉克曼(Blackman)窗、凯泽(kaiser)窗等，下面介绍几种常用的窗函数：矩形窗例：用矩形窗设计低通数字滤波器程序及运行结果如下：omegac=0.37;N=81;m=(N-1)/2;n=0:2*m+1

14、0;h=omegac/pi*sinc(omegac*(n-m)/pi);w=ones(1,N)zeros(1,length(n)-N);hd=h.*w;omega=-pi:2*pi/300:pi;Hd=freqz(hd,1,omega);页脚内容6吉林大学珠海学院plot(omega,abs(Hd);汉宁窗(升余弦窗).Rif舞.StM:州*'1-匚。M(m=S5/5)0.25小口+mM%利用傅氏变换的移位特性，汉宁窗频谱的幅度函数W()可用矩形窗的幅度函数表示为：三部分矩形窗频谱相加，使旁瓣互相抵消，能量集中在主瓣，旁瓣大大减小，主瓣宽度增加1倍。汉明窗(改进的开余弦窗)对汉宁窗的改

15、进，在主瓣宽度(对应第一零点的宽度)相同的情况下，旁瓣进一步减小，可使99.96%的能量集中在主瓣内。布莱克曼窗(三阶开余弦窗)增加一二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，增加N可减少过渡带。频谱的幅度函数为：凯塞窗以上四种窗函数，都是以增加主瓣宽度为代价来降低旁瓣。凯塞窗则可自由选择主瓣宽度和旁瓣衰减。10(x)是零阶贝塞尔函数，参数B可自由选择，决定主瓣宽度与旁瓣衰减。B越大，w(n)窗越窄,其频谱的主瓣变宽，旁瓣变小。一般取4VB<9B=5.44接近汉明B=8.5接近布莱克曼B=0为矩形除了上述窗口外，还有所谓的“组构窗口”，即由一些简单的窗函数组合构成，如常见

16、的三角形窗(Bartlett窗)它是由二个长度为(N/2)的矩形窗进行线性卷积而成的。几种窗函数的比较四种常用的窗口函数矩形窗；汉宁窗；汉明窗；布莱克曼窗图窗口函数频谱(N=51,A=20lg|W()/W(0)|)页脚内容7吉林大学珠海学院从(a)一(d),旁弁的衰减逐步增加，主弁相应加宽(用对数可放大小旁瓣)。图在同一指标下用四种窗口设计的LP滤波器频率特性(N=51,wc=0.8tt)可见，矩形窗设计的过渡带最窄，但阻带最小衰减也最差，仅-21dB;布莱克曼窗设计的阻带最小衰减最好，达-74dB,但过渡带最宽，约为矩形窗设计的三倍。几种窗口函数的具体性能比较见下表。窗函数主瓣宽度旁瓣峰值衰

17、减(dB)阻带最小衰减(dB)矩形4兀/N-13-21汉宁8兀/N-31-44汉明8兀/N-41-53布莱克曼12-/N-57-74设计要点：如果希望得到的滤波器的理想频率响应为H(ejw)那么FIR滤波器的设计就在于寻找一个传递函数去逼近H(ejw),在这种逼近中有两种直接的方法，一是从时域入手，即刚才讨论的窗函数设计法。另一种从频域入手，就是即将讨论的频率取样法窗函数设计法是从单位脉冲响应序列着手，使h(n)逼近理想的单位脉冲响应序列hd(n)Mn)。我们知道hd(n)可以从理想频响Hd(ejw)通过傅氏反变换获得但一般来说，理想频响都为分段恒定，在边界频率处有突变点，所以，这样得到的理想

18、单位脉冲响应hd(n)往往都是无限长序列，而且是非因果的。但FIR的h(n)是有限长的，问题是怎样用一个有限长的序列去近似无限长的hd(n)o最直接简单的办法是直接截取其一段来代替。这种截取可以形象地想象为h(n)就好象是通过一个“窗口”所看到一段hd(n),因此，h(n)也可表达为hd(n)和一个“窗函数”的乘积，即页脚内容8吉林大学珠海学院h(n)=w(n)hd(n)在这里窗口函数就是矩形脉冲函数RnS),当然以后我们还可看到，为了改善设计滤波器的特性，它还可以有其它的形式，相当于在矩形窗内对hd(n)作一定的加权处理。我们从一个截止频率为c的线性相位理想低通滤波器为例来讨论FIR的设计问

19、题，设低通滤波器的时延为a,即理想特性的hd(n弁DHd(w)这是一个以a为中心的偶对称的无限长非因果序列，如果截取一段n=0N-1的hd(n)作为h(n),则为要保证所得到是线性相位FIR滤波器，延时a应为h(n)长度N的一半，a=(N-1)/2其中W(n尸Rn”)(为了表示一致)设W(ejw)为该窗口函数的频谱：页脚内容9吉林大学珠海学院用幅度函数和相函数来表示，则有。其线性相位部分ej则是表示延时一半长度a=(N-1)/2,对频响起作用的是它的幅度函数图矩形窗函数及其幅度函数理想频响也可以写成幅度函数和相位函数的表示形式Hd(ejw)=Hd(w)ejwn其中幅度函数为二个信号时域乘积对应

20、于频域卷积，所以有如果也以幅度函数H()和相位函数来表示H(ejw),H(ejw)=H(co)ejwn则实际FIR滤波器的幅度函数H()为正好是理想滤波器幅度函数与窗函数幅度函数的卷积。=0的响应，因一般有c>>2冗/N,恸)在-wc,wc内近似包含全频域的值，所以，H(0)=H(0);二c时，一半重叠，H(c)/H(0)=0.5;二c-2冗/N时，第一旁瓣(负数)在通带外，出现正肩峰；=c+2tt/N时，第一旁瓣(负数)在通带内，出现负肩峰。由最后的频响函数图我们可以看到，页脚内容10吉林大学珠海学院在通带截止频率的两旁=入土期的地方，H()出现最大正负肩峰值；在这两个肩峰之间形

21、成一个过渡带，过渡带的宽度等于W()的主瓣宽度。在最大肩峰的两侧，则形成长长的余振；看公式町=sin(冷；：鼻状过;=r巴士其中x=Nw/2,所以N的改变不能改变主瓣与旁瓣的比例关系，最多只能改变W()的绝对值大小和起伏的密度，当N增加时，幅值拉高，频率轴变密，而最大肩峰经计算可知总为8.95%,这种现象称为吉布斯(Gibbs)效应。肩峰值的大小决定了滤波器通带内的平稳程度和阻带内的衰减，所以对滤波器的性能有很大的影响。为了改善滤波器的这些特性，只有改变窗函数的形状，标准为：窗谱主瓣宽度要窄，以获得较陡的过渡带；相对于主瓣幅度，旁瓣要尽可能小，使能量尽量集中在主瓣中，这样就可以减小肩峰和余振，

22、以提高阻带衰减和通带平稳性。但实际上这两个标准不能兼得，一般总是通过增加主瓣宽度来换取对旁瓣的抑制。窗口法的优点是简单，有闭合的公式可用，性能及参数都有表格资料可查，计算程序简单，较为实用。缺点是当Hd(ejw)较为复杂时，hd(n)就不容易由反付里叶变换求得。边界频率因为加窗的影响而不易控制。(2) .频率取样法工程上，常给定频域上的技术指标，所以采用频域设计更直接。一、基本思想：使所设计的FIR数字滤波器的频率特性在某些离散频率点上的值准确地等于所需滤波器在这些频率点处的值，在其他频率处的特性则有较好的逼近。二.设计方法：1)确定狂上&上页脚内容11吉林大学珠海学院2左4=日伏)=

23、2)计算h(n)h(n)=3)计算H(Z)三、约束条件：为了设计线性相位的FIR滤波器，采样值H(k)要满足一定的约束条件。前已指出，具有线性相位的FIR滤波器，其单位采样响应h(n)是实序列，且满足h(n)=+/-h(N-1-n),由此得到的幅频和相频特性，就是对H(k)的约束。(见P101,表4.1)。例如，要设计第一类线性相位FIR滤波器，即N为奇数，h(n)偶对称，则幅度函数H()应具有偶对称性；令.)由此可得，&必须取为：而Hk必须满足对称性：HkHnk同样，若要设计第二种线性相位FIR滤波器，N为偶数，h(n)偶对称，相位关系同上，由于幅度特性是奇对称的，因此，H也必须满足

24、对称要求：HkHnk其他两种线性相位FIR数字滤波器的设计，同样也要满足幅度与相位的约束条件。四、设计误差：FIR设计步骤：给定指标先时二耶)二琲抽)或H(z)页脚内容12吉林大学珠海学院关心的是，由上述设计过程得到的月，口与H(k)的逼近程度，以及丹/扣)与H(k)的关系令卬=或门”,则单位圆上的频响为:这是一个内插公式，式中)=尸、一1'为内插函数7Nsin32成/W)/2i=0,N1所以，在每个采样点上，频响严格地与理想特性H(k)一致，在采样点之间，频响由各采样点的内插函数延伸迭加而形成，因而有一定的逼近误差，误差大小与理想频率响应的曲线形状有关，理想特性平滑，则误差小；反之，

25、误差大，在理想频率响应的不连续点会产生肩峰和波纹。n增加，则采样点变密，内插误差减小。(3) .切比雪夫逼近法除了窗函数设计法(时间窗口法)和频率采样法，还可以用切比雪夫逼近法设计FIR滤波器,且切比雪夫逼近法是一种等波纹逼近法，在用切比雪夫逼近法设计FIR滤波器时，需要用到雷米兹(RemeZ交替算法，并且需要遵从两个准则：均方误差最小准则和最大误差最小化准则。雷米兹(RemeZ交替算法：能很好的解决通带截止频率声和阻带截止频率巩不能精确控制的问题。下图就是雷米兹(RemeZ交替算法的流程图(三).算例描述基于窗函数的FIR数字滤波器设计的方法十分简单，其主要步骤为：(1)对滤波器的理想特性进

26、行傅立叶逆变换获得理想滤波器的单位脉冲响应h(n)(2)由性能指标确定窗函数W(n)和窗口长度N,由过度带宽近似于窗函数主瓣宽求得窗口长度N(3) 求得实际滤波器的单位脉冲响应h(n)页脚内容13吉林大学珠海学院(4) 检验滤波器的性能下面举例说明：例1:用窗函数法设计一个线性相位FIR低通滤波器，并满足性能指标。通带边界频率p0.5,阻带边界频率ws0.66,阻带衰减不小于40dB,通带波纹不大于3dB。由题意，阻带衰减不小于40dB,选取汉宁窗。用MATLA蝙程如下：%MATLABPROGRAMeg-1%Window-basedFIRfilterdesign%Perfomentparame

27、terwp=0.5*pi;ws=0.66*pi;%Widthoftransitionbandwdelta=ws-wp;%LengthofthefilterN=ceil(8*pi/wdelta)ifrem(N,2)=0N=N+1;end%LengthofthewindowNw=N;%CutoffFrequencyofthefilterwc=(wp+ws)/2;%Computeimpulseresponseofidealfiltern=0:N-1;alpha=(N-1)/2;页脚内容14吉林大学珠海学院m=n-alpha+0.00001;hd=sin(wc*m)./(pi*m);%Computet

28、imeresponseoftheHanningWindowwin=hanning(Nw);%Computeactureimpulseresponseofthefilterh=hd.*win'b=h;freqz(b,1,512)程序运行即得所设计FIR线形相位滤波器频率特性。如图所示：该例题运用了上叙步骤进行编程，其实在MATLAIfi,信号处理工具箱提供了基于上叙原理设计标准型FIR滤波器的工具函数FIR1,函数FIR1是采用经典窗函数设计线性相位FIR数字滤波器的函数，其调用格式为：b=fir1(n,n)b=fir1(n,n,"type')b=fir1(n,“win

29、dow)b=fir1(n,n,'ftype',windows)下面的程序就是调用工具函数FIR1对上题进行重新编程%MATLABPROGRAMeg-1%Window-basedFIRfilterdesign%Perfomentparameterwp=0.5*pi;ws=0.66*pi;%Widthoftransitionbandwdelta=ws-wp;%LengthofthefilterN=ceil(8*pi/wdelta)%CutoffFrequencyofthefilterWn=(0.5+0.66)*pi/2;%using'FIR1'Designthefi

30、lter页脚内容15吉林大学珠海学院b=fir1(N,Wn/pi,hanning(N+1);freqz(b,1,512)由上可知，用函数FIR1设计的FIR数字滤波器在通带内具有很好的线性相位特性例二：设计一个24阶FIR带通滤波器，通带频率为0.350300.65利用工具函数FIR1编写程序如下：%MATLABPROGRAMeg-2%Window-basedFIRfilterdesign%Perfomentparameterwp=0.350.65;N=24;%using'FIR1'Designthefilterb=fir1(2*N,wp);freqz(b,1,512)这例也使

31、用了函数fir1的设计方法。例三：用CLS法设计一个25阶低通滤波器，截止频率0.3冗%MATLABPROGRAMeg-6%desiredFIRlowerfiltern=55;wo=0.3;dp=0.02;ds=0.008;b=fircls1(n,wo,dp,ds,'plot');此例调用了fircls1函数。在实际应用中，remez算法是一种比较常用的设计滤波器的算法。其优点有：它靠一次次的迭代求得一组交错点组频率，避免了直接的烦琐的计算下面列举使用REME商数设计滤波器的例子例四：用雷米兹算法设计一个滤波器.程序及波形如下。%designofequirippleliner-

32、phrasefirfilters页脚内容16吉林大学珠海学院N=14;f=00.20.40.70.81.0;m=0.10.1110.10.1;b=remez(N,f,m);>>h,w=freqz(b,1,256);>>mag=20*log10(abs(h);>>plot(w/pi,mag);grid;>>figure;>>k=0:N;>>stem(k,b);例五：用凯塞窗设计一FIR低通滤波器，低通边界频率吗=0.左,阻带边界频率%=0.5?1r,阻带衰减必不小于50dB。解：首先由过渡带宽和阻带衰减&来决定凯塞窗的N和凯塞窗设计对应的MATLA程序为：wn=kaiser(30,4.55);nn=0:1:29;alfa=(30-1)/2;hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa);h=hd.*wn'h1,w1=freqz(h,1);plot(w1/pi,20*log10(abs(h1);axis(0,1,-80,10);grid;xlabel('归