FIR数字滤波器设计与实现

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

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

3、基本结构有：1） 横截型（卷积型、直接型）a.一般FIR滤波器的横截型(直接型、卷积型)结构：若给定差分方程为： 。 则可以直接由差分方程得出FIR滤波器结构如下图所示：这就是FIR滤波器的横截型结构，又称直接型或卷积型结构。b.线性相位FIR滤波器的横截型结构 若h(n)呈现对称特性，即此FIR滤波器具有线性相位，则可以简化成横截型结构，下面分情况讨论：N为奇数时线性相位FIR滤波器实现结构如图所示：N为偶数时线性相位FIR滤波器实现结构如图所示我们知道IIR滤波器的优点是可利用模拟滤波器设计的结果，缺点是相位是非线性的，若需要线性相位，则要用全通网络进行校正，比较麻烦，而FIR滤波器的优点

4、是可以方便地实现线性相位。2)、级联型将H(z)分解为若干个实系数一阶或二阶因子相乘：实现结构如下图所示：该结构图中有2L=M个延迟器，2L+1=M+1个乘法器，2L=M个加法器2L=M个延迟器，2L+1=M+1个乘法器，2L=M个加法器特点：可以分别控制每个子系统的零点。分析H(z)及结构图可以得出级联型的特点：每个基本节控制一对零点，便于控制滤波器的传输零点。系数比直接型多，所需的乘法运算多。3）频率取样型 若FIR滤波器的冲激响应为有限长(N点)序列h(n)，则有如图所示的关系： 因此，对h(n)可以利用DFT得到H(k)，然后利用内插公式：来表示系统函数，这就为FIR滤波器提供了另外一

5、种结构：频率抽样结构， 这种结构由两部分级联而成：分析系统函数 其中级联的第一部分为： 这是一个梳状滤波器，它滤掉了频率 及其各次谐波。 级联的第二部分为N个一阶网络并联而成，第k个一阶网络为： 它在单位圆上有一个极点： 这是一个谐振频率 的无损耗谐振器。这个谐振器的极点正好与梳状滤波器的一个零点(i=k)相抵消，从而使这个频率上的频率响应等于H(k)。这样，N个谐振器的N个极点就和梳状滤波器的N个零点相抵消，从而在N个频率抽样点上的频率响应就分别等于N个H(k)值。有上叙的理论分析基础可以得到FIR滤波器的频率抽样结构。FIR滤波器的频率抽样结构如图所示：频率抽样结构的特点是它的系数H(k)

6、就是滤波器在 处的响应，因此控制滤波器的频率响应很方便。频率抽样结构存在问题的问题是：在有限长情况下，系数量化后极点不能和零点抵消，使FIR系统不稳定。解决方法：在r圆上进行(r1但近似等于1)取样，即用r代替，使极点和相应的零点移到单位圆内。（a）当N为偶数时的频率取样型结构如图所示。（b）当N为奇数时频率抽样型结构如图所示。4）快速卷积结构若FIR滤波器的单位冲激响应h(n)是一个N1点有限长序列，输入x(n)是一个N2点有限长序列，那么输出y(n)是x(n)与h(n)的线性卷积，它是一个LN1+N2-1点的有限长序列。 我们知道，将x(n)补上LN2个零值点，将h(n)补上LN1个零值点

7、，然后进行L点圆周卷积，就可以代替原x(n)与h(n)的线性卷积。而圆周卷积可以用DFT和IDFT的方法来计算，这样我们得到FIR滤波器的快速卷积结构：这里DFT和IDFT都将采用快速傅里叶变换算法，当N1和N2足够长时，比直接计算线性卷积要快得多。2.线性相位FIR滤波器的特点从以上的讨论中可以看出，我们最感兴趣的是具有线性特性的FIR滤波器，因此在设计FIR滤波器时，需要着重研究线性相位FIR滤波器的特点和性质，在上述已经介绍了线性相位FIR滤波器的横截型结构，现在介绍它的频响特性。 FIR滤波器的单位冲激响应h(n)是有限长的（0nN-1）,其Z变换为： 其傅立叶变换为： 其中H()是幅

8、度函数，是一个纯实数，可正可负， ()是相位函数。可以证明，线性相位FIR滤波器的冲激响应满足对称条件： h(n) =h(N-1-n)和（1）、线性相位FIR滤波器的幅度函数和相位函数：（a）当h(n)是偶对称时，其幅度函数和相位函数分别为： 特点： 幅度函数H()包括正负值，相位函数是严格线性相位，滤波器有(N-1)/2个抽样周期的延时，它等于单位抽样响应h(n)长度N的一半。（b）当h(n)是奇对称时，其幅度函数和相位函数分别为： 特点： 相位函数是严格线性相位，但在零频率(0)处有/2的相移。仍有(N-1)个抽样周期的延时。因此当h(n)为奇对称时，FIR滤波器将是一个具有准确相位的正交

9、变换网络。(2)、 FIR滤波器的线性相位特性 FIR滤波器的线性相位特性如图所示。(3)、 任何一种线性相位FIR滤波器的群延时都为： （4）FIR滤波器幅度函数的特点 分四种情况分别讨论H（）的特点：（a）当h(n)偶对称，N为奇数时： 幅度函数的特点：H（）对0， ， 呈偶对称。（b）当h(n)偶对称，N为偶数时： 幅度函数的 特点： 当 时， ，在z-1处有一个零点， 对 是 奇对称； 如果滤波器在 处幅度不为零（如高通滤波器）， 则不能用这种滤波器。 （c）当h(n)奇对称，N为奇数时幅度函数的特点：H（）在0， ，处都为零，也就是H(z)在 处为零；H（）对0， ， 都成奇对称。

10、（d）当h(n)奇对称，N为偶数时：幅度函数的特点： H（）在0， 处为零，即H(z)在z=1处为零点； H（）对0， 呈奇对称，对 呈偶对称。(5)、零点位置:线性相位FIR滤波器的系统函数有以下关系： 可见，若是H(z)的零点，则也一定是H(z)的零点。又由于当h(n)是实数时，H(z)的零点必成共轭对出现，所以及也一定是H(z)的零点。因而线性相位FIR滤波器的零点必是互为倒数的共轭对。 其有四种可能性： （1） 既不在实轴上，也不在单位园上，则是互为倒数的两组共轭对。 （2） 不在实轴上，但是在单位园上，则共轭对的倒数是它们本身，故只有一组共轭对。（3） 在实轴上而不在单位园上，只有倒

11、数部分，无复共轭部分。 （4） 既在实轴上又在单位园上，有两种可能，z1或z1。（二）FIR数字滤波器的设计在介绍和总结完FIR数字滤波器的基本结构和相关特性（包括频响曲线（幅度和相位），单位冲激响应等）后，就是FIR数字滤波器的设计和实现，FIR数字滤波器的设计步骤有：1.技术要求（预期性能指标） 技术要求是由实际用途决定的，它可由理想滤波器的系统函数Hd(z)、脉冲响应h(n)和差分方程描述。2.逼近（近似） 在数字滤波器性能分析的基础上，利用已经学过的概念和数学知识提供数字滤波器的表述，它是理想滤波器的一种近似。3.实现 上面一步的结果是一个滤波器的表述，它既可能是一个系统函数、也可能为

12、差分方程，或者是单位脉冲响应h(n)，依据这个结果进行数字滤波器结构的实际和软硬件的实现。下面重点介绍目前最主要的三种FIR数字滤波器的设计方法：（1）窗函数设计法（时间窗口法）这种方法也称为傅立叶级数法。其设计是在时域进行的，先用傅氏反变换求出理想数字滤波器的单位抽样响应hd(n)，然后时域移位并加时间窗w(n)对其截断，从而求得FIR 滤波器的单位抽样响应h(n)；在设计过程中，将无限长序列变为有限长是通过时域加矩形窗乘积实现数据的截断的。时域乘积对应了频域卷积，从而对频响特征发生的改变。常见的窗函数有：矩形窗、三角形(Bartlertt)窗、汉宁(Hanning)窗。海明(Hamming

13、)窗、布拉克曼(Blackman)窗、凯泽(kaiser)窗等，下面介绍几种常用的窗函数： 矩形窗 例：用矩形窗设计低通数字滤波器程序及运行结果如下：omegac=0.37;N=81;m=(N-1)/2;n=0:2*m+10;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);plot(omega,abs(Hd); 汉宁窗（升余弦窗）=0.5 利用傅氏变换的移位特性，汉宁窗频谱的幅度函数W（）可用矩形窗的幅度函数

14、表示为： 三部分矩形窗频谱相加，使旁瓣互相抵消，能量集中在主瓣，旁瓣大大减小，主瓣宽度增加1倍。 汉明窗（改进的升余弦窗） 对汉宁窗的改进，在主瓣宽度（对应第一零点的宽度）相同的情况下，旁瓣进一步减小，可使99.96%的能量集中在主瓣内。 布莱克曼窗（三阶升余弦窗） 增加一二次谐波余弦分量，可进一步降低旁瓣，但主瓣宽度进一步增加，增加N可减少过渡带。频谱的幅度函数为： 凯塞窗 以上四种窗函数，都是以增加主瓣宽度为代价来降低旁瓣。凯塞窗则可自由选择主瓣宽度和旁瓣衰减。 I0(x)是零阶贝塞尔函数，参数可自由选择，决定主瓣宽度与旁瓣衰减。越大，w(n)窗越窄，其频谱的主瓣变宽，旁瓣变小。一般取 4

15、2/N，WR()在-c, c内近似包含全频域的值，所以，H(0)Hd(0)； =c时，一半重叠，H(c)/H(0)=0.5； =c-2/N时，第一旁瓣（负数）在通带外，出现正肩峰； =c+2/N时，第一旁瓣（负数）在通带内，出现负肩峰。 由最后的频响函数图我们可以看到， 在通带截止频率的两旁=c 的地方，H()出现最大正负肩峰值； 在这两个肩峰之间形成一个过渡带，过渡带的宽度等于WR()的主瓣宽度。在最大肩峰的两侧，则形成长长的余振； 看公式 其中x=N/2,所以N的改变不能改变主瓣与旁瓣的比例关系，最多只能改变WR（）的绝对值大小和起伏的密度，当N增加时，幅值拉高，频率轴变密，而最大肩峰经计

16、算可知总为8.95%，这种现象称为吉布斯（Gibbs）效应。 肩峰值的大小决定了滤波器通带内的平稳程度和阻带内的衰减，所以对滤波器的性能有很大的影响。 为了改善滤波器的这些特性，只有改变窗函数的形状，标准为： 窗谱主瓣宽度要窄，以获得较陡的过渡带； 相对于主瓣幅度，旁瓣要尽可能小，使能量尽量集中在主瓣中，这样就可以减小肩峰和余振，以提高阻带衰减和通带平稳性。 但实际上这两个标准不能兼得，一般总是通过增加主瓣宽度来换取对旁瓣的抑制。 窗口法的优点是简单，有闭合的公式可用，性能及参数都有表格资料可查，计算程序简单，较为实用。缺点是当较为复杂时，就不容易由反付里叶变换求得。边界频率因为加窗的影响而不

17、易控制。（2）频率取样法工程上，常给定频域上的技术指标，所以采用频域设计更直接。 一、基本思想： 使所设计的FIR数字滤波器的频率特性在某些离散频率点上的值准确地等于所需滤波器在这些频率点处的值，在其他频率处的特性则有较好的逼近。 二.设计方法：1）确定 2）计算h(n) 3）计算H(Z) 三、 约束条件 ： 为了设计线性相位的FIR滤波器，采样值H（k）要满足一定的约束条件。 前已指出,具有线性相位的FIR滤波器，其单位采样响应h(n)是实序列，且满足h(n)=+/-h(N-1-n)，由此得到的幅频和相频特性，就是对H(k)的约束。（见P101，表4.1）。 例如，要设计第一类线性相位FIR

18、滤波器，即N为奇数，h(n)偶对称，则 幅度函数H（）应具有偶对称性； 令 由此可得， 必须取为： 而必须满足对称性： 同样，若要设计第二种线性相位FIR滤波器，N为偶数，h(n)偶对称，相位关系同上，由于幅度特性是奇对称的， 因此，Hk也必须满足对称要求： 其他两种线性相位FIR数字滤波器的设计，同样也要满足幅度与相位的约束条件。 四、设计误差 :FIR设计步骤： 给定指标 或H（z） 关心的是，由上述设计过程得到的 与H（k）的逼近程度，以及 与H（k）的关系。 令 ，则 单位圆上的频响为： 这是一个内插公式，式中 为内插函数令 则 所以，在每个采样点上，频响 严格地与理想特性H(k)一致

19、，在采样点之间，频响由各采样点的内插函数延伸迭加而形成，因而有一定的逼近误差，误差大小与理想频率响应的曲线形状有关，理想特性平滑，则误差小；反之，误差大，在理想频率响应的不连续点会产生肩峰和波纹。N增加，则采样点变密，内插误差减小。 （3）.切比雪夫逼近法除了窗函数设计法（时间窗口法）和频率采样法，还可以用切比雪夫逼近法设计FIR滤波器，且切比雪夫逼近法是一种等波纹逼近法，在用切比雪夫逼近法设计FIR滤波器时，需要用到雷米兹（Remez）交替算法，并且需要遵从两个准则：均方误差最小准则和最大误差最小化准则。雷米兹（Remez）交替算法：能很好的解决通带截止频率 和阻带截止频率 不能精确控制的问

20、题。下图就是雷米兹（Remez）交替算法的流程图(三). 算例描述基于窗函数的FIR数字滤波器设计的方法十分简单，其主要步骤为：（1） 对滤波器的理想特性进行傅立叶逆变换获得理想滤波器的单位脉冲响应h(n)（2） 由性能指标确定窗函数W(n)和窗口长度N，由过度带宽近似于窗函数主瓣宽求得窗口长度N（3） 求得实际滤波器的单位脉冲响应h(n)（4） 检验滤波器的性能下面举例说明：例1：用窗函数法设计一个线性相位FIR低通滤波器，并满足性能指标。通带边界频率，阻带边界频率，阻带衰减不小于40通带波纹不大于3。由题意，阻带衰减不小于40，选取汉宁窗。用MATLAB编程如下：% MATLAB PROG

21、RAM eg-1% Window-based FIR filter design%Perfoment parameterwp =0.5*pi;ws=0.66*pi;%Width of transition band wdelta =ws-wp;%Length of the filterN= ceil(8*pi/wdelta)if rem(N,2)=0N=N+1;end% Length of the windowNw =N;%Cutoff Frequency of the filterwc =(wp+ws)/2;%Compute impulse response of ideal filter

22、n =0: N-1;alpha =(N-1)/2;m =n-alpha+0.00001;hd =sin(wc*m)./(pi*m);%Compute time response of the Hanning Windowwin =hanning(Nw);%Compute acture impulse response of the filterh=hd.*win;b=h; freqz(b,1,512)程序运行即得所设计FIR线形相位滤波器频率特性。如图所示：该例题运用了上叙步骤进行编程，其实在MATLAB里，信号处理工具箱提供了基于上叙原理设计标准型 FIR滤波器的工具函数FIR1，函数FIR

23、1是采用经典窗函数设计线性相位 FIR数字滤波器的函数，其调用格式为：b=fir1(n,n) b=fir1(n, n,ftype)b=fir1(n, n,window)b=fir1(n, n,ftype,windows)下面的程序就是调用工具函数FIR1对上题进行重新编程% MATLAB PROGRAM eg-1% Window-based FIR filter design%Perfoment parameterwp =0.5*pi;ws=0.66*pi;%Width of transition band wdelta =ws-wp;%Length of the filterN= ceil(

24、8*pi/wdelta)%Cutoff Frequency of the filterWn=(0.5+0.66)*pi/2;%usingFIR1Design the filterb=fir1(N,Wn/pi,hanning(N+1);freqz(b,1,512)由上可知，用函数FIR1设计的FIR数字滤波器在通带内具有很好的线性相位特性例二：设计一个24阶FIR带通滤波器，通带频率为0.350.65利用工具函数FIR1编写程序如下：% MATLAB PROGRAM eg-2% Window-based FIR filter design%Perfoment parameterwp=0.35 0

25、.65;N=24;%usingFIR1Design the filterb=fir1(2*N, wp);freqz(b,1,512) 这例也使用了函数fir1的设计方法。例三： 用CLS法设计一个25阶低通滤波器，截止频率0.3% MATLAB PROGRAM eg-6%desired FIR lower filtern=55;wo=0.3;dp=0.02; ds=0.008; b=fircls1(n,wo,dp,ds,plot);此例调用了fircls1函数。在实际应用中，remez算法是一种比较常用的设计滤波器的算法。其优点有：它靠一次次的迭代求得一组交错点组频率，避免了直接的烦琐的计算下

26、面列举使用REMEZ函数设计滤波器的例子例四：用雷米兹算法设计一个滤波器.程序及波形如下。%design of equiripple liner-phrase fir filtersN=14;f=0 0.2 0.4 0.7 0.8 1.0;m=0.1 0.1 1 1 0.1 0.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低通滤波器，低通边界频率 ，阻带边界频率 ，阻带衰减 不小于50dB。 解 ：首先

27、由过渡带宽和阻带衰减 来决定凯塞窗的N和 凯塞窗设计对应的MATLAB程序为： 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(归一化频率/p) ； ylabel(幅度/dB) 例六：利用雷米兹交替算法，设计一个线性相位低通FIR数字滤波器，其指标为：通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动 阻带最小衰减At=40dB，采样频率fs=4000Hz。 解： 在MATLAB中可以用rem