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

第一章FIR数字滤波器的基本原理及设计方法有限长单位脉冲响应数字滤波器（FIRDF，FiniteImpulseResponseDigitalFilter）IIRDF主要对幅频特性进行逼FIRDFIIRDF高的多，但是同时考虑幅频特性指标和线FIRDF。FIRDFMATLAB工具箱函数的介绍。FIR数字滤波器的设计方法有窗函数法、频率采样法和基于firls函数和remez函数的最优化方法。MATLAB语言中的数字信号处理工具箱，提供了一些滤波器的函数，使FIR滤波器的运算更加方便和快捷。在MATLAB中提供的fir1(),FIRFIRfmm分别为滤波器的期望幅频响应的频率相Firls()remez()IIIFIR滤波器，III型的区别是偶函数还是奇函数。freqz()用于求数字滤波器的频率1.1 窗函数法设计FIR设我们所要设计的FIR滤波器的传输函数是Hd位脉冲响应，因此

(ej ),hd

是与其对应的单H (ej)d

hdn

(n)ej

（1-1）h(n) 1d

H

(ej)ejnd

（1-2)如果我们能够在Hd

(ej)已知的情况下，求出hd

(n)，经过Z变换可得到滤波器的系统函数。通常情况下理想数字滤波器的单位脉冲相应hd

(n)是无限长的，且是非因果序列。获得有限脉冲响应滤波器的一种可能方法是对hd

(n)截取一段h(n)来近似代替hd差。

(n)，可是这样会改变原来的滤波器指标，出现吉布斯效应误窗函数法就是用被称为窗函数的有限加权序列w(n)来修正式(1)的傅里叶基数以求得要求的有限脉冲响应序列h(n)，即dh(n)h(n)w(n) （1-3）dw(n)是有限长序列，当n<0或n>N-1时，w(n)=0。这种方法的重点在于选择某种合适的窗函数。要求窗函数主瓣宽度尽可能衰减大。下面介绍几种常用的窗函数：1.矩形窗(RectangleWindow)其频率函数为：

w (n)RR

(n) （3-4）Rsin(N/2) jNRW(ej) sin(/2) e

（3-5）三角形窗(BartlettWindow) 2n , 0n1(N1) N1 2w (n)

2n 1

（3-6）Br 2 , (N1)nN1 N1 2其频率函数为：2

W(ej)

ej 2

（3-7)Br N2)汉宁(Hanning)窗，又称升余弦窗w (n)11cos()

(n) （3-8）Hn 2N1N利用傅里叶变换得到频率函数为： NW(e)

(N

)W

(N

)e 2 R 2W N2Hn

1 R 1

（3-9）当N 1时，N1N，所以窗函数的幅度函数为W

0.25 )W)

(3-10)N N Hn

R R 汉明(Hamming)窗，又称改进的升余弦窗w

0.46cos( )

(n)

（3-11）Hm

N1N其幅度函数为：W

(

)

(

2)NN

（3-12）Hm R R R 布莱克曼(Blankman)窗，又称二阶升余弦窗 w (n)0.420.5cos(N

)0.08cos(N

)R

(n) （3-13)Bl 1其幅度函数为：

1 NW

(

2)W

(

2)N1 N1Bl 0.04

(

R4)

(

R)

(3-14）N1 W N1 R R凯泽(Kaiser)窗I(1[12nI(1[12n/(N1)]2)0

0nN1 （3-15）k I()0其中：β是一个可自由选择的参数，I0(x)是第一类修正零阶贝塞尔函数[10].上述窗函数的基本参数如下表窗函数旁瓣峰值幅度/db过渡带宽阻带最小衰减/db矩形窗-134/N-21三角形窗-268/N-25汉宁窗-318/N-44汉明窗-408/N-53布莱克曼窗-5712/N-74凯泽窗-5710/N-80窗函数法设计滤波器的步骤：根据技术要求确定待求滤波器的单位取样响应h(n。dN。h(n)：h(n)h(n)(n) （3-16)d式中,(n)是前面所选择好的窗函数。检验技术指标是否满足要求。根据下式计算：H(e

j)N1h(n)ejnn0

(3-17)H(e不满足要求，根据具体情况重复步骤(2)(3)(4)为止。本文以一个FIR滤波器的设计为例说明如何使用MATLAB设计数字滤波器FIRWp=0.4*,Ws=0.5*,3db40db。程序如下：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)程序执行后得幅频和相频如下图所示：M

500-50-100-1500 0.1 0.2 0.3 0.4 0.5 0.6 0.7 NormalizedFrequency( rad/sample)

0.9 1edeaP

0-1000-2000-3000-40000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 NormalizedFrequency( rad/sample)

0.9 1图1.1频率采样法设计FIR数字滤波器对理想滤波器的系统函数 Hd(z)进行频率采样以得到系统的理想频响Hd(ejw)的等间隔采样值H(k)H(k)实际上是所要求的滤波器的单位采样响应(h(n))的离散傅里叶变换（DF，如下试：N1

1

jwN1N

j

sin(Nw)2H(ejw)

Nk0N

H(k)(wk ) eN N

2k0

k)2 N

(3-18)h(n)

1N1H(k)Nk0

2j knN (3-19)为了减小H（k）的通带边缘由于抽样点的变化而引起的起伏振荡，可以增加过渡点，加宽过渡带以减小通带的起伏。sin(Nw)/sin(w)成正比，并位移（2k）/N的频2 2H（k）FIR点过渡可以使阻带衰减提高到-44~54dB,二点过渡衰减-65~75,三点过渡衰减-85~95dB.如果不能使过渡带太宽，同时要求增大阻带衰减，可以增加取样点数N，但这样会增加计算量、延时和误差。FIR滤波器设计步骤：(1)Hd(ejw)(2)根据过渡带宽和阻带衰减确定过渡点数和h(n)的长度N。 H(k)H (z) d ze2 jk/N(3)由IFFT计算IDFT得到：H(n)

1N1 j2knH(k)H(k)ek0低通带边缘：w1p=0.35*;高通带边缘：w2p=0.65*;高阻带边缘：w2s=0.8*。t1t2t1=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)实验得幅频相频特性如下图所示：(uaM

2000-200-4000 50 100 150 200 250 Frequency(Hz)

350 400 450 500e(hP

10000-1000-20000 50 100 150 200 250 Frequency(Hz)

350 400 450 500图1.2等波纹最优化方法设计FIR数字滤波器在数字信号处理中,利用数字滤波器可改变信号中所含频率分量的相对比例或滤除某些频率分量,使其达到所需要的效果.其中数字FIR滤波器由于具有精确的线性相位,且系统稳定,所以广泛应用于通信、数字图象处理、语音信号处理、自适应处理、雷达/声纳系统等方面.目前,FIR滤波器[1,2]的设计主要有窗函数设计法和频率采样设计法.但是,这2种方法都不易精确控制通带边界频率Wp与阻带边界频率Ws,所以,在实际应用中具有一定的局限性.而以最大误差最小化准则支持的切比雪夫逼近法是一种优异的设计方法,易于精确控制Wp与Ws.最低。等波纹最佳逼近法设计的数学证明复杂已超出本科生的数学基础所以们略去等波纹最佳逼近法复杂的数学推导，只介绍基本思想和实现线性相位FIRDF的等波纹最佳逼近设计的 MATLAB信号处理工具箱函数remezremezordRemezremezFIRDFh(n),由于切比雪夫和雷米兹(remez)对解决该问题做出了贡献，所以又称之为切比雪夫逼近法或者雷米兹逼近法。第二章设计步骤等波纹滤波器的最优化设计方法主要有2种,第1.它的思路是在给定的一些离散点上,使实际的幅频特性和理想幅频特性之间的误差的平方和为最小.第2种是最优化等波纹设计法,也称为雷米兹法或切比雪夫逼近法.该类型滤波器幅频特性在通带和阻带上的误差峰值是均匀分布的,差具有等波纹特性,因而把波纹的幅度控制到最小,或在同等指标下减小它的阶次.第1种方法是连续最小的平方法的推广,容易理解,但它的指标与滤波器没有直接关联,误差平方小的滤波器不能保证没有窄而大的波纹出现,像吉布斯效应那样.第2种方法直接控制通带波动和阻带衰减,最具针对性,是滤波器的最优化设计方法.因此,采用MATLAB信号处理工具箱提供的函数,纹设计法实现数字FIR滤波器的设计和仿真.。完整的最优化等波纹滤波器设计,除了切比雪夫等波纹逼近公式外,还要考虑:(1)滤波器采样点数n的确定.(2)极值数目的确定.最优化等波纹滤波器的误差函数在给定的频率上有(L+2)或(L+3)个极值,L为多项式的阶数.Wp,Ws的组合得到有(L+个极值的滤波器.此时,L=floor((n-1)/2).(3)建立频率修正的算法.在程序中自动进行反复的迭代修正,直到达到要求的精度为止.交替定理能保证切比雪夫逼近问题的解存在并且惟一,并没有说如何得到该解,既不知n(L),也不知极值的频率δ和波纹系W.ParksMcClellanRemez.n,并且设δ1=δ2,就可得到解,其中,δ系数,δ2为阻带波动系数.δ和n是相关的,n越大,δ越小.滤波器的技术指标中给出了通带波动系数δ1和δ2,Wp,Ws,因此需要设定n.凯泽提出逼近n的公式:ParksMcClellan算法首先猜设(L+2)个极值频率{Wi},最大误差及其Wi},由这些新的频率点拟合出一个新的L阶多项式,过程,直到找出最优集和全局最大波纹系数δ1为止.d(n),最后算出滤波器脉冲响应h(n).n是近似的,δ1可能不等于δ2若δ1>δ2则增加n;若δ1<δ2则减小n,然后再次用Remezδ.重复此过程,直至δ1≈δ2,这样就可得到等波纹滤波器.性能分析比较分别采用窗口法、频率取样法和优化设计法设计双带滤波器，其指标要求为：1）第一带通滤波器fp1=4kHz,fp2=5kHz， fr1=3.6kHz，fr2=5.42）第二带通滤波器fp3=6kHz，fp4=7kHz，fr3=5.6kHz，fr4=7.4Hz；3）通带衰耗：Ap<3rad，阻带衰耗：Ar>40dBfs=20kHzMATLAB模拟结果如下。图FIR获得广泛的应用。第三章MATLAB的实现MATLAB函数介绍MATLABremezordremez7.4.2所示。图7.4.2给出了等波纹滤波器技术指标的两种描述参数。7.4.2(a)dB数描述，即ωp=π／2，αp=2dB,ωs=11π／20,s=20dB.这是工程实际中常用的指标描述方法。但是，用等波纹最佳逼近设计法NW(ω)纹幅度δ1和δ(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=-20lg[δ2／（1+δ1）]≈-20lgδ2 (7.4.3)由式(7.4.2)和(7.4.3)得到δ1=(10^αp/20-1)/(10^αp/20+1) (7.4.4)δ2=10^(-αs)/20 (7.4.5)按照式和(7.4.4和(7.4.5)7.4.（bδ1=0.114δ2=0.1.实际中δ1和δ2δ1和δ2特意取较大值。1Remez函数remezFIR用格式为：hn=remez(M，f,m,w)hnM：FIRDF阶数，hnN=M+1。f和m：给出希望逼近的幅度特性。f为边界频率向量，0≤f≤1，要求为f为单调增向量（即f(k)＜f（k+1），01结束，1对应数字频率ω=π（Fs／2，Fs时域采样频率。mff表示频点f(k)的频响应值。如果用命令画出幅频响应曲线，则k[f(k),f(k+1)]上的幅频响应就是期望逼近的幅频响应值，频段[f（k+，（k+）]简言之，Plot（f，m）命令画出的幅频响应就是期望逼近的幅频响应值，频段[fk+（k+]Plot(f,m始频段为第一段，奇数频段为逼近区，偶数频段为无关区。W/g10下。设计希尔伯特变换器的调用格式为：hn=remez（M，f,m,w,’hilbert’）hn具有奇对称特性：hn(n)=﹣hn（M+1-f给定的通带1.f=[a,b]，m=[1,1],fa0＜a＜b＜1.这是因为希尔伯特变换器属于第二类线性相位滤波器，且N为奇数，所以只能实现带通滤波器。设计数字微分器调用格式为：hn=remez（M,f,m,w,’defferentiator’）下一节将介绍希尔伯特变换器和数字微分器的设计和应用。Remez函数的调用参数（M，f,m,remezord2.remezordFIRDFm、误差加权Mfremez函数的调用参数。其调用格式为：[M,fo,mo,w]=remezord（f,m,rip,Fs）参数说明如下。返回参数作为remezf,m,ripFsfremez中类似，这里f可以是模拟频率0Fs／2（1）结束，而且其0Fs／2FsFs=2Hz。但是这里f（0Fs／2两个频点mmf给定的一个逼近频段上希望逼近的幅度值。例如，对图16]；m=[1，0]。图3.1Fs时，fM略小，使设计结果达不到指标要求，这时要取M+1或M+2（必须注意对滤波器长度N=M+1的奇偶性要求（过渡带FS／2时，设计结果不正确。ripfm描述的各逼近频段允许的波纹振幅（幅频响应最大偏差rip的两倍。一般以[N,fo,mo,w]=remezord（f,m,rip,Fs）remez数，计算单位脉冲响应：hn=reme（，fo,mo,wremezremezordN，fo,mow的含义。remezremezordFIRDFremezordf,m,ripFsFs的，或根据实际信号处理要求（按照采样定理）确定。下面给出由给定的各种滤波器设计指标确定remezord调用参数f,m,和rip的公式，编程时直接诶套用即可。低通滤波器设计指标最小衰减：αsdB。remezord调用参数：f=[ωp／π，ωs／π]，m=1，0]，rip=δ1，δ2] 0和δ2分别为通带和阻带波纹幅度，由式（7.4.4）和式（7.4.5）计算得到，下面相同。高通滤波器设计指标最小衰减：αsdB。remezord调用参数：f=(7.4.7)带通滤波器设计指标逼近通带：[ωp1，ωpu]，通带最大衰减：αpdB；逼近阻带：[0，ωs1]，[ωsu,π]，阻带最小衰减：αsdB。remezord调用参数：f=[δ2，δ1,δ2] (7.4.8)带阻滤波器设计指标逼近阻带：[ωs1，ωsu]，阻带最大衰减：αsdB；逼近阻带：[0，ωp1]，[ωpu,π]，通带最小衰减：αpdB。remezord调用参数：f=[δ1,δ2,δ1] (7.4.9)设计举例1FIR率ω₂p=0.65π,阻带上截止频率ω₂s=0.8π,通带最大衰减αp=1dB,阻带最小衰减αs=60dB。remezordremezremezordremezordremezh(n)。将设计指标带入式remezordep741.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)dat1Dat1=(10∧(rp／20)-1)／(10∧(rp／20)+1);Dat2=10∧(-rs／20);Rip=[dat2,dat1,dat2];[M,fo,mo,w]=remezord(f,m,rip);M=M+2;%remezord估算的阶数偏小，加2才满足要求Hn=remez(M,fo,mo,w);绘图部分略去M=28h(nN=2h(n)7.4.5所示。而用窗函数设计的FIR80佳逼近法可以使滤波器的阶数大大降低。图3.2图