第4章5-7-数字滤波器的原理和设计方法

第四章数字滤波器的原理和设计方法TheoryandDesignMethodofDigitalFilter4.5IIR数字滤波器的频率变换前面的一些例子已经说明了如何应用冲激不变法和双线性变换法，根据具有选频特性的低通模拟系统函数来设计无限冲激响应数字滤波器。下图画出了四种常见的选频滤波器的理想频率响应，(a)、(b)、(c)和(d)分别表示低通、高通、带通和带阻滤波器的理想频率响应。设计这样一些选频滤波器的传统方法有两种。方法1：首先设计一个模拟原型低通滤波器，然后通过频率变换把它变换成所需要的模拟高通、带通或带阻滤波器，最后再使用冲激不变法或双线性变换法变换成相应的数字高通、带通和带阻滤波器。方法2：先设计了一个模拟原型低通滤波器，然后采用冲激响应不变法或双线性变换法将它转换成数字原型低通滤波器，最后通过频率变换把数字原型低通滤波器变换成所需要的数字高通、带通或带阻滤波器。方法1的缺点是，由于产生混叠失真，因此不能用冲激不变法来变换成高通或带阻滤波器。因此本节只讨论方法2。在方法2中，从模拟低通滤波器到数字低通滤波器的转换，前面已经讨论过了，因此下面只讨论数字低通滤波器到数字高通、带通和带阻滤波器的转换问题。我们曾经用双线性变换将模拟系统函数变换成数字系统函数，对于低通、高通、带通和带阻型选频滤波器来说，也可以利用一种和双线性变换非常相象的有理变换从低通数字滤波器变换出来。设Hl(v)是数字原型低通滤波器的系统函数，Hd(z)是所要求的滤波器的系统函数。数字域的频率变换，就是要找出一个变量代换的公式，使得所要求的系统函数这里用v-1是因为系统函数的标准形式,一般写成z-1的形式，换到v平面即是v-1。频率变换中的变量代换公式必须满足下列条件：(1)F(z-1)必须是z-1的有理函数；(2)v平面的单位圆内部映射到z平面的单位圆内部。从这些条件出发，我们可推导出频率变换的实用公式。设v平面单位圆是v=ejθ，z平面单位圆是z=ejω，则对比等式两边，有如果把变量代换的有理函数F(z-1)看成是一个系统函数，那么该系统的幅频特性曲线在任何ω处恒为1，这样的函数就是全通函数。任何全通函数都可表示为其中αk是F(z-1)的极点。为了满足稳定性的要求，必须有|αk|<1。这样，通过选择适当的N值和αk值，可以得出各种各样的映射。1)低通→低通的z平面变换已知数字低通原型H(v)的截止频率是θp，如果要把它变成截止频率为ωp的另一个低通H(z)。频率点的对应关系为：v平面 z平面θ=0 ω=0θ=π ω=πθ=θp ω=ωp当ω由0→π时，θ由0→π，变化量为1个π，所以F(z-1)是一阶全通，具有下列形式：将代入上式得：以代入，解出α，得到2)低通→高通的z平面变换3)低通→带通的z平面变换4)低通→带阻的z平面变换4.6FIR数字滤波器的设计方法IIR数字滤波器的优点是可以利用模拟滤波器设计的结果，而模拟滤波器的设计有大量图表可查，方便简单。但是它也有明显的缺点，就是相位的非线性，若须线性相位，则要采用全通网络进行相位校正。在图象处理以及数据传输中，都要求信道具有线性相位特性。而有限长单位冲激响应(FIR)数字滤波器就可以做成具有严格的线性相位，而同时可以具有任意的幅度特性。此外，FIR滤波器的单位抽样响应是有限长的，极点都位于原点，因而滤波器一定是稳定的。再有，只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而总能用因果系统来实现。最后，FIR滤波器由于单位冲激响应是有限长的，因而可以用快速傅里叶变换(FFT)算法来实现过滤信号，这可大大提高运算效率。IIR滤波器设计中的各种变换法，对FIR滤波器设计是不适用的，这是因为那里是利用有理分式的系统函数，而FIR滤波器的系统函数只是z-1的多项式。FIR滤波器的主要缺点是，必须用很长冲激响应的FIR滤波器才能很好地逼近锐截止滤波器，这意味着需要很大的运算量。另一个缺点是，线性相位FIR滤波器的时延不一定总是样本间隔的整数倍，在某些信号处理应用中，这种非整数时延会带来一些不希望有的问题。从以上讨论看出，我们最感兴趣的是具有线性相位的FIR滤波器。对非线性相位的FIR滤波器，一般可以用IIR滤波器来代替，因而这里不去讨论它。FIR数字滤波器的设计方法与IIR数字滤波器的设计方法很不一样，它不能利用模拟滤波器的设计方法。FIR数字滤波器的设计方法主要有窗函数法、频率取样法和等波纹逼近法等3种，本章主要介绍窗函数法，也简要地介绍频率取样法。4.6.1窗函数法这种方法也称为傅里叶级数法。设计方法：一般是先给定所要求的理想的滤波器频率响应Hd(ejω)，要求设计一个FIR滤波器频率响应来逼近Hd(ejω)。但是设计是在时域进行的，因而先由Hd(ejω)的傅里叶反变换导出hd(n)。由于理想滤波器的频率响应Hd(ejω)具有矩形频率特性，故hd(n)一定是无限长的序列，且是非因果的，而我们要设计的是FIR滤波器，其h(n)必然是有限长的，所以要用有限长的h(n)来逼近无限长的hd(n)，最有效的方法是截断hd(n)，或者说用一个有限长度的窗口函数序列ω(n)来截取hd(n)，因而窗口函数序列的形状及长度的选择就很关键。我们以一个截止频率为ωc的线性相位的理想矩形幅度特性的低通滤波器为例来加以讨论，设时延为α，即对应的冲激响应为：显然，hd(n)是以α为中心的无限长非因果序列，如右图(a)所示。现在需要寻找一个有限长序列h(n)来逼近hd(n)，h(n)应满足FIR滤波器的基本条件，即它是偶对称或奇对称的，以满足线性相位的要求，它还应当是因果的。这样，有和可以把h(n)看作是hd(n)与一矩形序列ωR(n)(如上图b所示)相乘的结果，即其中相乘的结果h(n)如上图c所示。ωR(n)称为矩形窗函数。窗函数不一定是矩形窗函数，也可以是其它窗函数，因此一般将h(n)表示成根据傅里叶变换的卷积性质，h(n)的频谱函数可表示为(4.79)即FIR数字滤波器的频谱函数是理想低通滤波器的频谱函数与窗函数频谱的卷积。采用不同的窗函数，对应的H(ejω)有不同的形状。矩形窗ωR(n)的频谱为其中矩形窗ωR(n)的频谱的图形如下图所示。ω从-2π/N到-2π/N之间的WR(ω)称为窗函数的主瓣，主瓣两侧呈衰减振荡的部分称为旁瓣。理想低通滤波器的频率响应可表示为：其幅度响应Hd(ω)为将Hd(ejω)和WR(ejω)分别代入4.79得到因此FIR数字滤波器的幅度响应为(4.81)上式表明，由理想低通滤波器的冲激响应加窗得到的FIR滤波器，它的幅度响应等于理想低通滤波器的幅度响应与窗函数频谱的幅度响应的周期卷积，如下图所示。只要看几个特殊的频率点，就可以看出H(ω)的一般情况。特别要注意卷积过程给H(ω)造成的起伏现象。1、先来看ω=0时零频率处的响应值H(0)。根据式(4.81)，H(0)等于图4.49中(a)与(b)两个函数乘积的积分，即WR(θ)在θ=-ωc到θ=+ωc这一段的面积，当ωc>>2π/N时(这个条件一般能满足)，H(0)实际上就很近似于WR(θ)全部θ从-π到+π的面积。2、再看ω=ωc时的卷积值，这时Hd(θ)正好与WR(ω-θ)的一半重叠，如图4.49(d)所示，因此卷积值正好是零频响应H(0)的一半，即H(ωc)/H(0)＝0.5，如图4.49(f)所示。3、当ω在通带截止频率ωc以内，即ω=ωc-2π/N时，WR(ω-θ)的整个主瓣都在Hd(θ)的通带内，如图4.49(d)所示，因此卷积结果有最大值，这时频率响应出现正肩峰。4、对于ω=ωc-2π/N、WR(ω-θ)的主瓣全部在Hd(θ)的通带外，如图4.49(e)所示，在通带内旁瓣负的面积大于正的面积，因此卷积值达到最大负值，H(ω)在这里出现负肩峰。5、当ω进一步增大时，卷积值也将随着WR(ω-θ)的旁瓣在通带内的面积的变化而变化，这样就造成H(ω)以零值为中心的上下起伏波动。6、当ω由ω-2π/N向通带内减小时，WR(ω-θ)的右旁瓣进入Hd(ω)的通带，这时，卷积值H(ω)在WR(ω-θ)的主瓣和左右旁瓣的共同作用下将以H(0)为中心上下波动。从以上分析及图4.49(f)可以看出，理想低通滤波器经加窗处理后，主要受到加窗处理两方面的影响。第一，使滤波器的理想频率特性在不连续点处边沿加宽，出现过渡带，这主要是由窗函数频谱的主瓣引起的，过渡带的宽度取决于窗函数主瓣的宽度，矩形窗对应的过渡带的宽度Δω=4π/N。一般来说，过渡带的宽度与N成反比；第二，滤波器在通带和阻带内产生波纹，这种现象称为吉布斯(Gibbs)现象，主要是由窗函数的频谱的旁瓣造成的。根据以上讨论，可知在一般情况下，对窗函数的要求有二：①旁瓣高度尽可能小，即尽可能让能量集中于主瓣，以减少通带和阻带中的波纹；②主瓣宽度尽量窄，以获得尽可能陡的过渡带。但是，这两个要求是互相矛盾的，不可能同时满足。具体来说，降低旁瓣高度必然会使主瓣变宽；反之，压窄主瓣宽度，不可避免地会使旁瓣变高。往往是增加主瓣宽度以换取对旁瓣的抑制。以矩形窗为例，它的频谱为可见，改变N，只能改变窗谱的主瓣宽度、改变ω坐标的比例以及改变WR(ω)的绝对值大小，但是不能改变主瓣与旁瓣的相对比例(当然N太小时，会影响旁瓣的相对值)，这个相对比例是由sinx/x决定的，或者说只由窗函数的形状来决定的。因而，当截取长度N增加时，只会减小过渡带宽(4π/N)，而不会改变肩峰的相对值。用矩形窗截取无限长序列hd(n)来得到有限长序列h(n)，由于突然将h(n)截短，因而破坏了序列hd(n)的均匀收敛性，这意味着人为地强迫hd(n)收敛。不均匀收敛性在频谱中是以吉布斯现象反映出来的。矩形窗所形成的FIR滤波器的频率响应的波纹幅度很大，最大肩峰值达8.95％，如图4.49(f)所示。为了减小波纹幅度，一方面可以加大窗的长度N，但效果并不显著；另一方面可采用不同的窗函数来改善不均匀收敛性。图4.50所示的是几种常用的窗函数：它们的定义式和频谱函数分述如下：1、矩形窗2、Bartlett窗（三角形窗）3、汉宁(Hanning)窗(升余弦窗)或利用傅里叶变换的调制特性，即利用和考虑到RN(n)的傅里叶变换为则得当N>>1时，N-1≈N，得到窗谱的幅度函数为因此可以认为汉宁窗的频谱由图4.51所示的3部分组成，3部分频谱相加的结果使旁瓣大大抵消，而使能量有效地集中在主瓣内，代价是使主瓣的宽度加大了一倍，即为8π/N。4、哈明(Hamming)窗，又称改进的升余弦窗把升余弦窗加以改进，可以得到旁瓣更小的效果，窗形式为：其频率响应的幅度函数为结果可将99.963％的能量集中在窗谱的主瓣内，与汉宁窗相比，主瓣宽度相同为8π/N，但旁瓣幅度更小，旁瓣峰值小于主瓣峰值的1%。5、布莱克曼(Blackman)窗为了更进一步抑制旁瓣，可再加上余弦的二次谐波分量，得到Blackman窗：其频谱的幅度函数为此时主瓣宽度为矩形窗谱主瓣宽度的三倍，即为12π/N。图4.52描绘的是N=51时上列5种窗函数的频谱函数图形，图中以相对衰减A=201g|W(ω)/W(0)|dB为纵坐标。从图中可以看出，这5种窗函数的旁瓣衰减依次增大，主瓣宽度依次加宽。图4.53所示的是用这5种窗函数设计的低通FIR数字滤波器的频率响应特性。窗函数的长度N＝51，理想低通滤波器的截止频ωc=π/2。从图中可看出，用矩形窗设计的滤波器的过渡带最窄，但阻带衰减指标最差，仅有-21dB左右。而用布莱克曼窗设计的阻带衰减指标最好，可达-74dB，但过渡带最宽，约为矩形窗的3倍。6、凯泽(Kaiser)窗这是一种适应性较强的窗，其窗函数的表示式为其中，I0(x)是第一类修正零阶贝塞尔函数。它可用以下的级数来计算：在实际应用中，级数取15～25项就可以达到足够的精度。凯泽窗是一族窗函数。β是可调参数，调节β值可以改变主瓣的宽度和旁瓣的幅度，β的典型值在4<β<9范围内。凯泽窗的曲线示于图4.55中。当n=(N-1)/2(即中点)时，ω(n)=1，当n从中点向两边变化时，ω(n)逐渐减小，β越大，ω(n)变化愈快。当n=0和n=N-1时，ω(0)=ω(N-1)=1/I0(β)。在上图中，β=5.44的曲线接近于哈明窗，β=8.5的曲线与Blackman窗相近，而β=0的曲线就是矩形窗。参数β选得越大，则ω(n)窗越窄，而频谱的旁瓣越小，但主瓣宽度也相应增加。因而改变β值就可对主瓣宽度与旁瓣衰减进行选择。可以看出，阻带最小衰减只由窗形状决定，不受N的影响，而过渡带的宽度则随窗宽N的增加而减小。窗函数设计FIR数字滤波器的步骤：(1)给出希望设计的滤波器的频率响应函数Hd(ejω)；(2)根据允许的过渡带宽度及阻带衰减，初步选定窗函数和N值；(3)计算以下积分，求出hd(n)；或(4)将hd(n)与窗函数相乘得FIR数字滤波器的冲激响应h(n)；(5)计算FIR数字滤波器的频率响应，并验证是否达到所要求的指标或从而得到H(ω)和φ(ω)。窗函数法计算中的主要问题其一，需要预先确定窗函数的形式和窗序列的点数N，以满足给定的频率响应指标，从而很难准确控制滤波器通带的边缘；其二，若Hd(ejω)不能用简单函数表示，则计算式(4.97)的积分非常困难。第一个问题只有通过多次设计来解决。如图4.56所示，理想低通滤波器的截止频率为ωc，由于窗函数主瓣的作用而产生过渡带，出现了通带截止频率ω1和阻带截止频率ω2。在ω1和ω2处的衰减是否满足通带和阻带的要求，也就是ω1和ω2是否就是所需要的通带和阻带的截止频率，这是不一定的。为了得到满意的结果，不得不假设不同的ωc，进行多次设计。例如，假设要求设计的FIR数字滤波器在ω=0.25π处有-3dB的衰减，设计时先选择理想低通滤波器的ωc=0.25π，得到的FIR数字滤波器在ω=0.25π处衰减为-6dB左右，不符合要求；然后改选ωc=0.30π进行设计，得到的FIR数字滤波器在ω=0.2π处的衰减恰好为-3dB左右。第二个问题的解决办法是用求和来代替积分，以便在计算机上计算，也就是要计算离散傅里叶反变换，一般都采用FFT来计算。由式(4.97)知将积分限分成M段，也就是令抽样频率为则有频域的抽样，造成时域序列的周期延拓，延拓周期是M，即实际上，由于hd(n)随n增加衰减很快，一般只要M足够大，即M>>N，在窗口范围内能很好地逼近hd(n)。窗函数法的优点是简单，有闭合形式的公式可循，因而很实用。其缺点是通带、阻带的截止频率不易控制。例4.7设计一个满足下式要求的FIR线性相位低通数字滤波器

0.98≤|H(ejω)|≤1.02，0≤|ω|≤0.18π|H(ejω)|≤0.003，0.22π≤|ω|≤π(1)选择合适的窗函数。(2)求滤波器的阶数。(3)求理想低通滤波器的截止频率ωc和滤波器的时延α。(4)求滤波器的单位取样响应。解：阻带衰减为

20log(0.003)=-50.5dB查P146，表4.2，选择哈明窗可满足阻带衰减的要求。(2)题目要求的过渡带宽度为

Δω=0.22π-0.18π=0.04π或Δf=Δω/2π=0.02由表4.2知道，哈明窗的过渡带宽度Δω与滤波器的阶数N有下列关系：Δω=8π/N，从而N=8π/Δω=8π/0.04π=200(3)ωc=(ωP+ωT)/2=(0.18π+0.22π)/2=0.20π滤波器的时延α=(N-1)/2=199/2=99.5，取α=100。(4)滤波器的单位取样响应因滤波器的频率特性为H(ejω)=|H(ejω)|e-jωα所以滤波器的单位取样响应为4.6.2频率取样法窗函数法是从时域出发，把理想的hd(n)用一定形状的窗函数截取成有限长的h(n)，以此h(n)来近似理想的hd(n)，这样得到的频率响应H(ejω)逼近于我们所要求的理想的频率响应Hd(ejω)。频率抽样法则是从频域出发，把给定的理想频率响应Hd(ejω)加以等间隔抽样，所以称之为频率取样法。然后以此Hd(k)作为实际FIR数字滤波器的频率特性的抽样值H(k)，即令知道H(k)后，由DFT定义，可以用领域的这N个抽样值H(k)来唯一确定有限长序列h(n)，而由频域采样的内插公式(3.50)知道，利用这N个频域抽样值H(k)同样可求得FIR滤波器的系统函数H(z)及频率响应H(ejω)。这个H(z)或H(ejω)将逼近H(d)或Hd(ejω)。根据书本P.78频率取样，可以得到内插公式：(4.105)FIR滤波器的频率响应也可用H(k)表示：(4.106)其中φ(ω)为内插函数，将φ(ω)代入4.106，化简后可得式(4.105)和式(4.106)是用频率取样法设计FIR数字滤波器的理论基础。从内插公式我们看到，在各频率抽样点上，滤波器的实际频率响应是严格地和理想频率响应数值相等，即但是在抽样点之间的频响则是由各抽样点的加权内插函数的延伸叠加而形成，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状，理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小，如下图(b)梯形理想频率特性所示。反之，如果抽样点之间的理想频率特性变化越陡，则内插值与理想值之误差就越大，因而在理想频率特性的不连续点附近，就会产生肩峰和波纹，如图(a)矩形理想频率特性所示。一般来说，在过渡带的频率响应特性变化较陡，在通带或阻带内的变化比较平缓，因此我们重点考虑过渡带抽样的优化设计问题。为了提高逼近质量，使逼近误差更小，也就是减小在通带边缘由于抽样点的陡然变化而引起的起伏振荡(这种起伏振荡就使得阻带内最小衰减很小)，和窗函数法的平滑截断一样，这里是使理想频率响应的不连续点的边缘，加上一些过渡的抽样点，这些点上抽样的最佳值由计算机算出。这样的结果，就增加了过渡带，减小了频带边缘的突变，也就减小了起伏振荡，因而增大了阻带最小衰减。这些抽样点上的取值不同，效果也就不同，因为由频率抽样公式看出，每一个频率抽样值，都要产生一个与常数sin(ωN/2)/sin(ω/2)成正比并且在频率上位移2πk/N的频率响应，而FIR滤波器的频率响应就是H(k)与内插函数φ(ω-2πk/N)的线性组合。如果精心设计过渡带的抽样值，就有可能使它的相邻频带(通带、阻带)的波纹得以减小，从而设计出较好的滤波器。一般过渡带取一、二、三点抽样值即可得到满意结果。在低通设计中，不加过渡抽样点时，阻带最小衰减为-20dB，一点过渡抽样的最优设计，阻带最小衰减可提高到-44dB到-54dB左右，二点过渡抽样的最优设计可达-65dB到-75dB左右，而加三点过渡抽样的最优设计则可达-85dB到-95dB左右。频率取样型FIR滤波器设计步骤给定理想滤波器频率响应。根据过渡带宽和阻带衰减确定过渡点数和h(n)长度N。。由IFFT计算IDFT得到。频率取样方法的优点是，可以直接在频域中用选择过渡取样值的方法，来得到良好的设计效果，适合于最优化设计。这种方法的缺点是，频率控制点的位置受到频率轴上的N个取样点的限制，因而滤波器的截止频率不易控制。如果要自由地选择截止频率，就必须增加取样点数N，这样做是不经济的。例利用频率取样法，设计一个低通FIR数字滤波器，其理想频率特性是矩形，已知ωc=0.5π，抽样点数为奇数N=33，要求滤波器具有线性相位。解：根据指标，可画出频率取样后的H(k)序列，如下图。由于|H(k)|是对称于ω=π的，我们只对0≤ω≤π即0≤k≤16的区间感兴趣。故可将π≤ω≤2π即17≤k≤32的图形略去不画。截止频率ωc=0.5π满足16π/33≤ωc≤17π/33，按频率取样方式设计，N=33，则将以上值代入H(ejω)，得到按此式计算的结果如下图所示。由图看出过渡带宽为2π/33，而最小阻带衰减则约为-20dB。这一衰减在大多数情况下是不令人满意的。为了改善频率特性，以满足指标要求，可在通带和阻带交界处安排一个或几个不等于零也不等于l的抽样值。本例中用优化算法算出在k=9处，|H(9)|=0.5，则得如下图(a)所示结果。这相当于加宽过渡带，其宽度为2×2π/N=4π/33，算出阻带最小衰减约为-40dB左右，如下图(b)所示。如果要进一步增加阻带衰减，可再添上第二个不等于1也不等于零的取样，这样过渡带又加宽了，如果不允许再增大过渡带宽，而又希望增大阻带衰减，还可增加抽样点数N，例如同样是