数字信号处理第7章有限脉冲响应数字滤波器的设计资料

1、第7章 有限脉冲响应数字滤波器的设计 7.1 线性相位FIR数字滤波器的条件和特点 7.2 利用窗函数法设计FIR滤波器 7.3 利用频率采样法设计FIR滤波器 7.4 利用等波纹最佳逼近法设计FIR数字滤波器 7.5 IIR和FIR数字滤波器的比较 用N表示FIR滤波器单位脉冲响应h(n)的长度，其系统函数H(z)为 H(z)是z1的N1次多项式，它在z平面上有N1个零点，在原点z=0处有一个N1重极点。因此，H(z)永远稳定。稳定和线性相位特性是FIR滤波器最突出的优点。FIR滤波器的设计方法和IIR滤波器的设计方法有很大差别。FIR滤波器设计任务是选择有限长度的h(n)，使频率响应函数H

2、(ej)满足技术指标要求。7.1 线性相位FIR数字滤波器的条件和特点1 线性相位FIR数字滤波器对于长度为N的h(n)，频率响应函数为（7.1.1）（7.1.2）式中，Hg()称为幅度特性; ()称为相位特性。注意，这里Hg()不同于|H(ej)|，Hg()为的实函数，可能取负值，而|H(ej)|总是正值。线性相位FIR滤波器是指()是的线性函数，即为常数-第一类线性相位如果()满足下式：是起始相位- 第二类线性相位。 严格地说，此时()不具有线性相位特性，但以上两种情况都满足群时延是一个常数，即也称这种情况为线性相位。2. 线性相位FIR的时域约束条件线性相位FIR滤波器的时域约束条件是指

3、满足线性相位时，对h(n)的约束条件。1) 第一类线性相位对h(n)的约束条件 第一类线性相位FIR数字滤波器的相位函数()=，由式（7.1.1）和（7.1.2）得到:（7.1.5） 由式（7.1.5）得到:(7.1.6)将（7.1.6）式中两式相除得到：即 移项并用三角公式化简得到:函数h(n)sin(n)关于求和区间的中心(N1)/2奇对称，是一组解。因为sin(n)关于n=奇对称，如果取=(N1)/2，则要求h(n)关于(N1)/2偶对称，所以要求和h(n)满足如下条件: Attention: n为自变量（7.1.8） 如果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第一类

4、线性相位特性（严格线性相位特性），则h(n)应当关于n=(N1)/2点偶对称。 当N确定时，FIR数字滤波器的相位特性是一个确知的线性函数，即()=(N1)/2。N为奇数和偶数时, h(n)的对称情况分别如表7.1.1中的情况1和情况2所示。表7.1.1 线性相位FIR数字滤波器的时域和频域特性一览 2) 第二类线性相位对h(n)的约束条件第二类线性相位FIR数字滤波器的相位函数()=/2， 经过同样的推导过程可得到: 函数h(n)cos(n)关于求和区间的中心(N1)/2奇对称，是一组解，因为cos(n)关于n=偶对称，所以要求和h(n)满足如下条件：（7.1.10） 由以上推导结论可知，如

5、果要求单位脉冲响应为h(n)、长度为N的FIR数字滤波器具有第二类线性相位特性，则h(n)应当关于n=(N1)/2点奇对称。N为奇数和偶数时h(n)的对称情况分别如表7.1.1中情况3和情况4所示。2 线性相位FIR滤波器幅度特性Hg()的特点 幅度特性的特点就是线性相位FIR滤波器的频域约束条件。为了推导方便，引入两个参数符号： 表示取不大于(N1)/2的最大整数。显然，仅当N为奇数时，M=(N1)/2 。 情况1： h(n)=h(Nn1), N为奇数。将时域约束条件h(n)=h(Nn1)和()=得到:所以因为cos(n-)关于=0, , 2三点偶对称，可以看出，Hg()关于=0, , 2三

6、点偶对称。因此情况1可以实现各种（低通、高通、带通、带阻）滤波器。对于N=13的低通情况，Hg()的一种例图如表7.1.1中情况1所示。Attention: 为自变量情况2： h(n)=h(Nn1), N为偶数。仿照情况1的推导方法得到:（7.1.12）式中，。因为是偶数，所以当时而且cos(n)关于过零点奇对称，关于=0和2偶对称。所以Hg()=0，Hg()关于=奇对称，关于=0和2偶对称。因此，情况2不能实现高通和带阻滤波器。对N=12 的低通情况，Hg()如表7.1.1中情况2所示。情况3： h(n)=h(Nn1)，N为奇数。将时域约束条件h(n)=h(Nn1)和()=/2）， 并考虑，

7、得到:式中，N是奇数，=(N1)/2是整数。所以，当=0，, 2时，sin(n)=0, 而且sin(n)关于过零点奇对称。因此Hg()关于=0, , 2三点奇对称。由此可见，情况3只能实现带通滤波器。对N=13的带通滤波器举例，Hg()如表7.1.1中情况3所示。 情况4： h(n)=h(Nn1), N为偶数。用情况3的推导过程可以得到:式中，N是偶数，=(N1)/2=N/21/2。所以，当=0, 2时，sin(n)=0；当=时，sin(n)=(1)(nN/2), 为峰值点。而且sin(n)关于过零点=0和2两点奇对称，关于峰值点=偶对称。因此Hg()关于=0和2两点奇对称，关于=偶对称。由此

8、可见，情况4不能实现低通和带阻滤波器。对N=12的高通滤波器举例，Hg()如表7.1.1中情况4所示。3. 线性相位FIR数字滤波器的零点分布特点将h(n)=h(N1n)代入上式, 得到:由（7.1.14）式可以看出，如z=zi是H(z)的零点，其倒数也必然是其零点；又因为h(n)是实序列，H(z)的零点必定共轭成对，因此也是其零点。这样，线性相位FIR滤波器零点必定是互为倒数的共轭对，确定其中一个，另外三个零点也就确定了, 如图7.1.1中。当然，也有一些特殊情况，如图7.1.1中z1、z2和z4情况。图7.1.1 线性相位FIR数字滤波器的零点分布7.2 利用窗函数法设计FIR滤波器7.2

9、.1 窗函数法设计原理设希望逼近的滤波器频率响应函数为Hd(ej)，其单位脉冲响应是hd(n)。如果能够由已知的Hd(ej)求出hd(n)，经过Z变换可得到滤波器的系统函数。但通常以理想滤波器作为Hd(ej)，其幅度特性逐段恒定，在边界频率处有不连续点，因而hd(n)是无限时宽的，且是非因果序列。例如，线性相位理想低通滤波器的频率响应函数Hd(ej)为 其单位脉冲响应hd(n)为 理想低通滤波器的单位脉冲响应hd(n)是无限长，且是非因果序列。hd(n)的波形如图7.2.1（a）所示。为了构造一个长度为N的第一类线性相位FIR滤波器，只有将hd(n)截取一段，并保证截取的一段关于n=（N1）/

10、2偶对称。设截取的一段用h(n)表示，即式中, RN(n)是一个矩形序列，长度为N，波形如图7.2.1（b）所示。由该图可知，当取值为（N1）/2时，截取的一段h(n) 关于n=（N1）/2偶对称，保证所设计的滤波器具有线性相位。我们实际设计的滤波器的单位脉冲响应为h(n)，长度为N，其系统函数为H(z)，这样用一个有限长的序列h(n)去代替hd(n)，肯定会引起误差，表现在频域就是通常所说的吉布斯（Gibbs）效应。该效应引起过渡带加宽以及通带和阻带内的波动，尤其使阻带的衰减小，从而满足不了技术上的要求。这种吉布斯效应是由于将hd(n)直接截断引起的，因此，也称为截断效应。图7.2.1 窗函

11、数设计法的时域波形（矩形窗，N=30）图7.2.2 吉普斯效应 Hd(ej)是一个以2为周期的函数，可以展为傅里叶级数，即傅里叶级数的系数为hd(n)，为Hd(ej)对应的单位脉冲响应。设计FIR滤波器就是根据要求找到N个傅里叶级数系数h(n)，n=1, 2, , N1，以N项傅氏级数去近似代替无限项傅氏级数，这样在一些频率不连续点附近会引起较大误差，这种误差就是前面说的截断效应。窗函数法又称为傅氏级数法显然，选取傅氏级数的项数愈多，引起的误差就愈小，但项数增多即h(n)长度增加，也使成本和滤波计算量加大，应在满足技术要求的条件下，尽量减小h(n)的长度。RN(n)（矩形序列）就是起对无限长序

12、列的截断作用，RN(n)可看做一个窗口，h(n)则是从窗口看到的一段hd(n)序列，所以称h(n)=hd(n)RN(n)为用矩形窗对hd(n)进行加窗处理。 用w(n)表示窗函数，用下标表示窗函数类型，矩形窗记为wR(n)。用N表示窗函数长度。（7.2.4）根据傅里叶变换的时域卷积定理，h(n)=hd(n)WR(n)的傅里叶变换：式中，Hd(ej)和WR(ej)分别是hd(n)和WN(n)的傅里叶变换，即（7.2.5）WRg()称为矩形窗的幅度函数，将图中2/N, 2/N区间上的一段波形称为WRg()的主瓣，其余较小的波动称为旁瓣。将Hd(ej)写成Hd(ej)=Hdg()ej， 则按照（7.

13、2.1）式，理想低通滤波器的幅度特性函数（如图7.2.3（a）所示）为式中图7.2.3 矩形窗加窗效应将Hd(ej)和WR(ej)代入（7.2.4）式，得到：将H(ej)写成H(ej)=Hg()ej ，则加窗后的滤波器的幅度特性等于理想低通滤波器的幅度特性Hdg()与矩形窗幅度特性WRg()的卷积。Summary： 1， 在理想特性不连续点=c附近形成过渡带。过渡带的宽度近似等于WRg()主瓣宽度4/N。2， 通带内产生了波纹，最大的峰值在c2/N处。阻带内产生了余振，最大的负峰在c+2/N处。通带与阻带中波纹的情况与窗函数的幅度谱有关， WRg()旁瓣幅度的大小直接影响Hg()波纹幅度的大小

14、。 吉布斯效应直接影响滤波器的性能。通带内的波纹影响滤波器通带的平稳性，阻带内的波纹影响阻带内的衰减，可能使最小衰减不满足技术指标要求。 一般滤波器都要求过渡带愈窄愈好。直观上，好像增加矩形窗的长度，即加大N，就可以减少吉布斯效应的影响。只要分析一下N加大时WRg()的变化，就可以看到这一结论不是完全正确。我们讨论在主瓣附近的情况。在主瓣附近，按照式（7.2.5），WRg()可近似为该函数的性质是随x加大（N加大），主瓣幅度加高，同时旁瓣也加高，保持主瓣和旁瓣幅度相对值不变；另一方面，N加大时，WRg()的主瓣和旁瓣宽度变窄，波动的频率加快。三种不同长度的矩形窗函数的幅度特性WRg()曲线如图

15、7.2.4(a)、(b)、(c)所示。用这三种窗函数设计的FIR滤波器的幅度特性Hg()曲线如图7.2.4(d)、(e)、(f)所示。因此，当N加大时，Hg()的波动幅度没有多大改善，带内最大肩峰比H(0)高8.95%，阻带最大负峰值为H(0)的8.95%，使阻带最小衰减只有21 dB。加大N只能使Hg()过渡带变窄（过渡带近似为主瓣宽度4/N）。因此加大N, 并不是减小吉布斯效应的有效方法。图7.2.4 矩形窗函数长度的影响调整窗口长度N只能有效地控制过渡带的宽度，而要减少带内波动以及增大阻带衰减，只能从窗函数的形状上找解决问题的方法。构造新的窗函数形状，使其谱函数的主瓣包含更多的能量，相应

16、旁瓣幅度更小。旁瓣的减小可使通带、阻带波动减小，从而加大阻带衰减。但这样总是以加宽过渡带为代价的。7.2.2 典型窗函数介绍几种常用窗函数的时域表达式、时域波形、幅度特性函数（衰减用dB计量）曲线，以及用各种窗函数设计的FIR数字滤波器的单位脉冲响应和损耗函数曲线。为了叙述简单，把这组波形图简称为“四种波形”。下面均以低通为例，Hd(ej)取理想低通，c=/2，窗函数长度N=31。1 矩形窗（Rectangle Window）wR(n)=RN(n)前面已分析过，按照（7.2.5）式，其幅度函数为（7.2.7） 为了描述方便，定义窗函数的几个参数: 旁瓣峰值n窗函数的幅频函数|Wg()|的最大旁

17、瓣的最大值相对主瓣最大值的衰减值（dB）；过渡带宽度Bg用该窗函数设计的FIR数字滤波器（FIRDF）的过渡带宽度；阻带最小衰减s用该窗函数设计的FIRDF的阻带最小衰减。图7.2.4所示的矩形窗的参数为: n=13 dB；Bg=4/N; s=21 dB。2 三角形窗（Bartlett Window）（7.2.8） 其频谱函数为其幅度函数为 （7.2.10）三角窗的四种波形如图7.2.5所示，参数为: n=25 dB; Bg=8/N; s=25 dB。图7.2.5 三角窗的四种波形3 汉宁（Hanning）窗升余弦窗（7.2.11） 当N1时， N1N汉宁窗的幅度函数WHng()由三部分相加，

18、旁瓣互相对消，使能量更集中在主瓣中。汉宁窗的四种波形如图7.2.6所示，参数为: n=31 dB; Bg=8/N; s=44 dB。图7.2.6 汉宁窗的四种波形4 哈明（Hamming）窗改进的升余弦窗 （7.2.12）其频谱函数WHm(ej)为其幅度函数WHmg()为当N时，其可近似表示为这种改进的升余弦窗，能量更加集中在主瓣中，主瓣的能量约占99.96%，瓣峰值幅度为40 dB，但其主瓣宽度和汉宁窗的相同，仍为8/N。可见哈明窗是一种高效窗函数，所以MATLAB窗函数设计函数的默认窗函数就是哈明窗。哈明窗的四种波形如图7.2.7所示，参数为: n=41 dB; Bg=8/N; s=53

19、dB。图7.2.7 哈明窗的四种波形（7.2.13）5 布莱克曼（Blackman）窗其频谱函数为其幅度函数为这样其幅度函数由五部分组成，它们都是移位不同，且幅度也不同的WRg()函数，使旁瓣再进一步抵消。旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍。布莱克曼窗的四种波形如图7.2.8所示，参数为: n=57 dB; B=12/N; s=74 dB。（7.2.14）图7.2.8 布莱克曼窗的四种波形 6 凯塞贝塞尔窗（Kaiser-Basel Window）以上五种窗函数都称为参数固定窗函数，每种窗函数的旁瓣幅度都是固定的。凯塞贝塞尔窗是一种参数可调的窗函数，是一种最优窗函数。 （7

20、.2.15）式中I0()是零阶第一类修正贝塞尔函数，可用下面级数计算：一般I0()取1525项，便可以满足精度要求。 参数可以控制窗的形状。一般加大，主瓣加宽，旁瓣幅度减小，典型数据为4 9。当 =5.44时，窗函数接近哈明窗。 =7.865时，窗函数接近布莱克曼窗。在设计指标给定时，可以调整值，使滤波器阶数最低，所以其性能最优。凯塞（Kaiser）给出的估算和滤波器阶数N的公式如下:（7.2.17） 式中，Bt=|sp|, 是数字滤波器过渡带宽度。应当注意，因为式（7.2.17）为阶数估算，所以必须对设计结果进行检验。另外，凯塞窗函数没有独立控制通带波纹幅度，实际中通带波纹幅度近似等于阻带波

21、纹幅度。凯塞窗的幅度函数为（7.2.16） （7.2.18）对的8种典型值，将凯塞窗函数的性能列于表7.2.1中，供设计者参考。由表可见, 当 =5.568时, 各项指标都好于哈明窗。6种典型窗函数基本参数归纳在表7.2.2中，可供设计时参考。表7.2.1 凯塞窗参数对滤波器的性能影 表7.2.2 6种窗函数的基本参数 表中过渡带宽和阻带最小衰减是用对应的窗函数设计的FIR数字滤波器的频率响应指标。随着数字信号处理的不断发展，学者们提出的窗函数已多达几十种，除了上述6种窗函数外，比较有名的还有Chebyshev窗、Gaussian窗。wn=boxcar(N)列向量wn中返回长度为N的矩形窗函数

22、w(n)wn=bartlett(N) 列向量wn中返回长度为N的三角窗函数w(n)wn=hanning(N) 列向量wn中返回长度为N的汉宁窗函数w(n)wn=hamming(N) 列向量wn中返回长度为N的哈明窗函数w(n)wn=blackman(N) 列向量wn中返回长度为N的布莱克曼窗函数w(n)wn=kaiser(N, beta) 列向量wn中返回长度为N的凯塞贝塞尔窗函数w(n)7.2.3 用窗函数法设计FIR滤波器的步骤用窗函数法设计FIR滤波器的步骤如下：(1) 根据对过渡带及阻带衰减的指标要求，选择窗函数的类型，并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰

23、减满足要求的情况下，尽量选择主瓣窄的窗函数。然后根据过渡带宽度估计窗口长度N。待求滤波器的过渡带宽度Bt近似等于窗函数主瓣宽度，且近似与窗口长度N成反比，NA/Bt，A取决于窗口类型，例如，矩形窗的A=4，哈明窗的A=8等，参数A的近似和精确取值参考表7.2.2。 (2) 构造希望逼近的频率响应函数Hd(ej)，即 对所谓的“标准窗函数法”，就是选择Hd(ej)为线性相位理想滤波器（理想低通、理想高通、理想带通、理想带阻）。以低通滤波器为例，Hdg()应满足： （7.2.19） 由图7.2.2知道，理想滤波器的截止频率c近似位于最终设计的FIRDF的过渡带的中心频率点，幅度函数衰减一半（约6

24、dB）。所以如果设计指标给定通带边界频率和阻带边界频率p和s， 一般取（7.2.20）(3) 计算hd(n)。 如果给出待求滤波器的频响函数为Hd(ej)，那么单位脉冲响应用下式求出：（7.2.21）如果Hd(ej)较复杂，或者不能用封闭公式表示，则不能用上式求出hd(n)。我们可以对Hd(ej)从=0到=2采样M点，采样值为 ，k=0 ,1, 2, , M1，进行M点IDFT(IFFT)，得到： （7.2.22）根据频域采样理论，hdM(n)与hd(n)应满足如下关系：因此，如果M选得较大，可以保证在窗口内hdM(n)有效逼近hd(n)。对线性相位理想低通滤波器作为Hd(ej)，其单位脉冲响

25、应hd(n)：为保证线性相位特性， =(N1)/2。(4) 加窗得到设计结果：h(n)=hd(n)w(n)。【例7.2.1】 用窗函数法设计线性相位高通FIRDF，要求通带截止频率p=/2 rad，阻带截止频率s=/4 rad，通带最大衰减 p=1 dB，阻带最小衰减 s=40 dB。解 (1） 选择窗函数w(n)，计算窗函数长度N。已知阻带最小衰减 s=40 dB，由表（7.2.2）可知汉宁窗和哈明窗均满足要求，我们选择汉宁窗。本例中过渡带宽度Btps=/4， 汉宁窗的精确过渡带宽度Bt=6.2/N，所以要求Bt=6.2/N/4，解之得N24.8。对高通滤波器N必须取奇数，取N=25。由式（

26、7.2.11）, 有（2） 构造Hd(ej):式中（3） 求出hd(n):将=12代入得 (n12)对应全通滤波器，是截止频率为3/8的理想低通滤波器的单位脉冲响应，二者之差就是理想高通滤波器的单位脉冲响应。这就是求理想高通滤波器的单位脉冲响应的另一个公式。 （4） 加窗:7.2.4 窗函数法的MATLAB设计函数简介实际设计时一般用MATLAB工具箱函数。可调用工具箱函数fir1实现窗函数法设计步骤（2）（4）的解题过程。(1) fir1是用窗函数法设计线性相位FIR数字滤波器的工具箱函数，以实现线性相位FIR数字滤波器的标准窗函数法设计。这里的所谓“标准”，是指在设计低通、高通、带通和带阻

27、FIR滤波器时，Hd(ej)分别表示相应的线性相位理想低通、高通、带通和带阻滤波器的频率响应函数。因而将所设计的滤波器的频率响应称为标准频率响应。Fir1的调用格式及功能： hn=fir1(M, wc)，返回6 dB截止频率为wc的M阶（单位脉冲响应h(n)长度N=M+1）FIR低通（wc为标量）滤波器系数向量hn，默认选用哈明窗。滤波器单位脉冲响应h(n)与向量hn的关系为h(n)=hn(n+1) n=0, 1, 2, , M而且满足线性相位条件: h(n)=h(N1n)。其中wc为对归一化的数字频率，0wc1。 当wc=wcl, wcu时，得到的是带通滤波器，其6 dB通带为wclwcu。

28、 hn=fir1(M, wc, ftype)，可设计高通和带阻FIR滤波器。当ftype=high时，设计高通FIR滤波器；当ftype=stop，且wc=wcl, wcu时，设计带阻FIR滤波器。 应当注意，在设计高通和带阻FIR滤波器时，阶数M只能取偶数（h(n)长度N=M+1为奇数）。不过，当用户将M设置为奇数时，fir1会自动对M加1。 hn=fir1(M, wc, window)，可以指定窗函数向量window。如果缺省window参数，则fir1默认为哈明窗。例如:hn=fir1(M, wc, bartlett(M+1)，使用Bartlett窗设计；hn =fir1(M, wc,

29、blackman(M+1)，使用blackman窗设计；hn=fir1(M, wc, ftype, window)，通过选择wc、ftype和window参数（含义同上），可以设计各种加窗滤波器。(2) fir2为任意形状幅度特性的窗函数法设计函数，用fir2设计时，可以指定任意形状的Hd(ej)，它实质是一种频率采样法与窗函数法的综合设计函数。主要用于设计幅度特性形状特殊的滤波器（如数字微分器和多带滤波器等）。用help命令查阅其调用格式及调用参数的含义。例7.2.1 的设计程序ep721.m如下: ep721.m: 例7.2.1 用窗函数法设计线性相位高通FIR数字滤波器wp=pi/2;

30、ws=pi/4;Bt=wp-ws; 计算过渡带宽度N0=ceil(6.2*pi/Bt); 根据表7.2.2汉宁窗计算所需h(n)长度N0，ceil(x)取大于等于x的最小整数N=N0+mod(N0+1, 2); 确保h(n)长度N是奇数wc=(wp+ws)/2/pi; 计算理想高通滤波器通带截止频率(关于归一化)hn=fir1(N-1, wc, high, hanning(N); 调用fir1计算高通FIR数字滤波器的h(n)略去绘图部分运行程序得到h(n)的25个值:h(n)= 0.0004 0.0006 0.0028 0.0071 0.0000 0.0185 0.0210 0.0165 0

31、.0624 0.0355 0.10610.2898 0.6249 0.2898 0.1061 0.0355 0.0624 0.0165 0.0210 0.0185 0.0000 0.0071 0.00280.00060.0004高通FIR数字滤波器的h(n)及损耗函数如图7.2.9所示。图7.2.9 高通FIR数字滤波器的h(n)波形及损耗函数曲线【例7.2.2】 对模拟信号进行低通滤波处理，要求通带0f 1.5kHz内衰减小于1 dB， 阻带2.5kHzf 上衰减大于40 dB。希望对模拟信号采样后用线性相位FIR数字滤波器实现上述滤波，采样频率Fs=10 kHz。用窗函数法设计满足要求的F

32、IR数字低通滤波器，求出h(n)，并画出损耗函数曲线。为了降低运算量，希望滤波器阶数尽量低。解 （1） 确定相应的数字滤波器指标: 通带截止频率为 阻带截止频率为 阻带最小衰减为 s =40dB（2） 用窗函数法设计FIR数字低通滤波器，为了降低阶数选择凯塞窗。根据式（7.2.16）计算凯塞窗的控制参数为指标要求过渡带宽度Bt=sp=0.2，根据式（7.2.17）计算滤波器阶数为取满足要求的最小整数M=23。所以h(n)长度为N=M+1=24。但是，如果用汉宁窗，h(n)长度为N=40。理想低通滤波器的通带截止频率c=(s+p)/2=0.4，所以, 得到:式中，w(n)是长度为24（ =3.3

33、95）的凯塞窗函数。实现本例设计的MATLAB程序为ep722.m。ep722.m: 例7.2.2 用凯塞窗函数设计线性相位低通FIR数字滤波器 fp=1500;fs=2500;rs=40; wp=2*pi*fp/Fs;ws=2*pi*fs/Fs; Bt=ws-wp; %计算过渡带宽度 alph=0.5842*(rs-21)0.4+0.07886*(rs-21); %根据(7.2.16)式计算kaiser窗的控制参数 N=ceil(rs-8)/2.285/Bt); %根据(7.2.17)式计算kaiser窗所需阶数N wc=(wp+ws)/2/pi; %计算理想高通滤波器通带截止频率(关于归一

34、化） hn=fir1(N,wc,kaiser(N+1,alph); %调用kaiser计算低通FIRDF的h(n) 以下绘图部分省去运行程序得到h(n)的24个值:h(n)=0.00390.00410.00620.01470.00000.02860.0242 0.03320.07550.00000.1966 0.37240.37240.1966 0.00000.07550.0332 0.0242 0.02860.00000.0147 0.00620.00410.0039低通FIR数字滤波器的h(n)波形和损耗函数曲线如图7.2.10所示。图7.2.10 低通FIR数字滤波器的h(n)波形及损耗

35、函数曲线【例7.2.3】 窗函数法设计一个线性相位FIR带阻滤波器。要求通带下截止频率lp =0.2，阻带下截止频率ls=0.35，阻通带上截止频率us=0.65，通带上截止频率up=0.8， 通带最大衰减 p=1 dB，阻带最小衰减 s=60 dB。解 本例直接调用fir1函数设计。因为阻带最小衰减 s=60 dB，所以选择布莱克曼窗，再根据过渡带宽度选择滤波器长度N，布莱克曼窗的过渡带宽度Bt=12/N，所以解之得N=80。调用参数设计程序为ep723.m，参数计算也由程序完成。ep723.m: 例7.2.3 用窗函数法设计线性相位带阻FIR数字滤波器wlp=0.2*pi;wls=0.35

36、*pi;wus=0.65*pi;wup=0.8*pi; %设计指标参数赋值B=wls-wlp; %过渡带宽度N=ceil(12*pi/B); %计算阶数N,ceil(x)为大于等于x的最小整数wp=(wls+wlp)/2/pi,(wus+wup)/2/pi; %设置理想带通截止频率hn=fir1(N,wp,stop,blackman(N+1); %带阻滤波器要求h(n)长度为奇数，所以取N+1省略绘图部分程序运行结果: N=81由于h(n)数据量太大，因而仅给出h(n)的波形及损耗函数曲线，如图7.2.11所示。图7.2.11 带阻FIR数字滤波器的h(n)波形及损耗函数曲线7.3 利用频率采

37、样法设计FIR滤波器1 用频率采样法设计滤波器的基本思想设希望逼近的滤波器的频响函数用Hd(ej)表示，对它在=0到2之间等间隔采样N点，得到Hd(k)：（7.3.1）再对Hd(k)进行N点IDFT，得到h(n)： （7.3.2）式中，h(n)作为所设计的FIR滤波器的单位脉冲响应，其系统函数H(z)为（7.3.3）另外根据频率域采样理论，得到H(z)的内插表示形式：（7.3.4）此式就是直接利用频率采样值Hd(k)形成滤波器的系统函数，（7.3.3）式适合FIR直接型网络结构，（7.3.4）式适合频率采样结构。 一个是为了设计线性相位FIR滤波器，频域采样序列Hd(k)应满足的条件；另一个是

38、逼近误差问题及其改进措施。2 设计线性相位滤波器时对Hd(k)的约束条件FIR滤波器具有线性相位的条件是h(n)为实序列，且满足h(n)=h(Nn1)，在此基础上我们已推导出其频响函数应满足的条件是：（7.3.5）（7.3.6）(7.3.7)(7.3.8)N = 奇数 N = 偶数 在=02区间上N个等间隔的采样频点为将=k代入（7.3.5）（7.3.8）式中，并写成k的函数：（7.3.12） k = 0，1，2，N 1 （7.3.9）（7.3.10）（7.3.11）N = 奇数 N = 偶数 N等于奇数时Hg(k)关于N/2点偶对称。N等于偶数时，Hg(k)关于N/2点奇对称，且Hg(N/2

39、)=0。 设用理想低通作为希望逼近的滤波器Hd(ej)，截止频率为c，采样点数为N，Hg(k)和(k)用下列公式计算：N=奇数时，（7.3.13） N=偶数时，上面公式中的kc是通带内最后一个采样点的序号，所以kc的值取不大于cN/(2)的最大整数。另外，对于高通和带阻滤波器，这里N只能取奇数。（7.3.14）3 逼近误差及其改进措施如果待逼近的滤波器为Hd(ej)，对应的单位脉冲响应为hd(n)，则由频率域采样定理知道，在频域02范围等间隔采样N点，利用IDFT得到的h(n)应是hd(n)以N为周期的周期延拓的主值区序列，即如果Hd(ej)有间断点，那么相应的单位脉冲响应hd(n)应是无限长

40、的。这样，由于时域混叠及截断，使h(n)与hd(n)有偏差。所以，频域的采样点数N愈大，时域混叠愈小，设计出的滤波器频响特性愈逼近Hd(ej)。上面是从时域分析其设计误差的来源，下面从频域分析。频域采样定理表明，频率域等间隔采样H(k)，经过IDFT得到h(n)，由（3.3.7）和（3.3.8）式，得到H(ej)=FTh(n)的内插表示形式： 式中上式表明，在采样频点，k=2k/N, k=0，1，2，N1, (2k/N)=1，因此采样点处与H(k)相等，逼近误差为0。在采样点之间，由N项H(k)(2k/N)之和形成。频域幅度采样序列Hg(k)及其内插波形Hg()如图7.3.1所示。图7.3.1

41、的上图 7.3.1(a) 中, 实线表示希望逼近的理想幅度函数Hdg(), 黑点表示幅度采样序列Hg(k); 下图7.3.1(b)中, 实线Hg(）与虚线dg()的误差与Hdg()特性的平滑程度有关，Hdg()特性愈平滑的区域，误差愈小；特性曲线间断点处，误差最大。表现形式为间断点变成倾斜下降的过渡带曲线，过渡带宽度近似为2/N。通带和阻带内产生震荡波纹，且间断点附近振荡幅度最大，使阻带衰减减小，往往不能满足技术要求。 当然，增加N可以使过渡带变窄，但是通带最大衰减和阻带最小衰减随N的增大并无明显改善。且N太大，会增加滤波器的阶数，即增加了运算量和成本。N=15和N=75两种情况下的幅度内插波

42、形Hg()如图7.3.2所示，图中的空心圆和实心圆点分别表示N=15和N=75时的频域幅度采样。运行本书程序库程序fig731-2.m，即可绘制出图7.3.1和图7.3.2，并分别输出两种采样点数（N=15和N=75）的通带最大衰减 p和阻带最小衰减 s。N=15时，通带最大衰减 p=0.8340 dB，阻带最小衰减 s=15.0788 dB； N=75时，通带最大衰减 p=1.0880 dB，阻带最小衰减 s=16.5815 dB。 所以，直接对理想滤波器的频率响应采样的“基本频率采样设计法”不能满足一般工程对阻带衰减的要求。 图7.3.1 频域幅度采样序列Hg(k)及其内插波形Hg()图7

43、.3.2 N=15和N=75的幅度内插波形Hg() 在窗函数设计法中，通过加大过渡带宽度换取阻带衰减的增加。频率采样法同样满足这一规律。提高阻带衰减的具体方法是在频响间断点附近区间内插一个或几个过渡采样点，使不连续点变成缓慢过渡带，这样，虽然加大了过渡带，但阻带中相邻内插函数的旁瓣正负对消，明显增大了阻带衰减。过渡带采样点的个数与阻带最小衰减 s的关系以及使阻带最小衰减 s最大化的每个过渡带采样值求解都要用优化算法解决。其基本思想是将过渡带采样值设为自由量，用一种优化算法（如线性规划算法）改变它们，最终使阻带最小衰减 s最大。表7.3.1 过渡带采样点的个数m与滤波器阻带最小衰减 s的经验数据

44、 4 频率采样法设计步骤综上所述，可归纳出频率采样法的设计步骤: （1） 根据阻带最小衰减 s选择过渡带采样点的个数m。 （2） 确定过渡带宽度Bt，估算频域采样点数（即滤波器长度）N。如果增加m个过渡带采样点，则过渡带宽度近似变成(m+1)2/N。当N确定时，m越大，过渡带越宽。如果给定过渡带宽度Bt，则要求(m+1)2/NBt ，滤波器长度N必须满足如下估算公式:（7.3.15） （3） 构造一个希望逼近的频率响应函数: 设计标准型片断常数特性的FIR数字滤波器时，一般构造幅度特性函数Hdg()为相应的理想频响特性，且满足表7.1.1要求的对称性。 （4） 按照（7.3.1）式进行频域采样

45、: （7.3.16） （7.3.17） 并加入过渡带采样。过渡带采样值可以设置为经验值，或用累试法确定，也可以采用优化算法估算。 （5） 对H(k)进行N点IDFT，得到第一类线性相位FIR数字滤波器的单位脉冲响应：（7.3.18） （6） 检验设计结果。如果阻带最小衰减未达到指标要求，则要改变过渡带采样值，直到满足指标要求为止。如果滤波器边界频率未达到指标要求，则要微调Hdg()的边界频率。 上述设计过程中的计算相当繁琐，所以通常借助计算机设计。MATLAB就是一种很有效的语言。【例7.3.1】 用频率采样法设计第一类线性相位低通FIR数字滤波器，要求通带截止频率p=/3，阻带最小衰减大于4

46、0 dB，过渡带宽度Bt/16。解 查表7.3.1， s=40 dB时，过渡带采样点数m=1。将m=1和Bt/16代入（7.3.15）式估算滤波器长度: N(m+1)2/Bt=64，留一点富余量，取N=65。构造Hd(ej)=Hdg()ej(N1)/2为理想低通特性，其幅度响应函数Hdg()如图7.3.3(a)中实线所示。图7.3.3 一个过渡点的设计结果(T=0.38)设计由以下程序ep731.m完成：ep732.m: 用频率采样法设计FIR低通滤波器 T=input(T= ) 输入过渡带采样值T Bt=pi/16;wp=pi/3; 过渡带宽度为pi/16，通带截止频率为pi/3 m=1;N

47、=ceil(m+1)*2*pi/Bt)+1; 按式（7.3.15）估算采样点数N Np=fix(wp/(2*pi/N); Np+1为通带0, wp上采样点数 Ns=N2*Np1; Ns为阻带wp, 2*piwp上采样点数 Hk=ones(1, Np+1), zeros(1, Ns), ones(1, Np);N为奇数，幅度采样向量偶对称H(k)=H(Nk) Hk(Np+2)=T;Hk(NNp)=T; 加一个过渡采样 thetak=pi*(N1)*(0:N1)/N; 相位采样向量(k)=(N1)k/N, 0kN1 Hdk=Hk.*exp(j*thetak); 构造频域采样向量Hd(k) hn=r

48、eal(ifft(Hdk); h(n)=IDFTH(k) Hw=fft(hn, 1024); 计算频率响应函数:DFTh(n)wk=2*pi* 0:1023/1024; Hgw=Hw.*exp(j*wk*(N1)/2); 计算幅度响应函数Hg() 计算通带最大衰减Rp和阻带最小衰减Rs Rp=max(20*log10(abs(Hgw) hgmin=min(real(Hgw);Rs=20*log10(abs(hgmin) 以下绘图部分略去运行程序，输入T=0.38，得到设计结果如图7.3.3所示，并输出通带最大衰减 p=0.4767 dB， 阻带最小衰减 s=43.4411 dB。但是，如果过渡

49、带采样值T=0.5和0.6，则得到阻带最小衰减 s=29.6896 dB和25.0690 dB。由此可见，当过渡带采样点数给定时，过渡带采样值不同，则逼近误差不同。所以，对过渡带采样值进行优化设计才是有效的方法。hn= fir2(M,F,A,window(M+1)设计一个M阶线性相位FIR数字滤波器，返回长度为N=M+1的单位脉冲响应序列向量hn。window表示窗函数名，缺省该项时默认选用Hamming窗。可供选择的窗函数有Boxcar、 Bartlett、Hann、Hamm、 Blackman、Kaiser和 Chebwin。当window=boxcar时，fir2就是纯粹的频率采样设计法

50、。希望逼近的幅度特性由边界频率向量F和相应的幅度向量A确定，plot(F, A)画出的就是希望逼近的幅度特性曲线。图7.3.4 例7.3.1希望逼近的幅度特性（T=0.38）F为对归一化的数字频率向量，0F1。而且F的元素必须是单调递增的，以0开始，以1结束，1对应于模拟频率Fs/2。对例7.3.1，F=0,wc/pi, wc/pi+2/N,wc/pi+4/N,1; A=1,1,T,0,0。其中，wc/pi+2/N为过渡带采样点频率，T为过渡带采样值。plot(F,A)画出的就是希望逼近的幅度特性曲线, 如图7.3.4所示。%调用fir2求解例7.3.1的程序ep731b.mT=input(T

51、= ) %键入过渡采样值TBt=pi/16;wp=pi/3; %过渡带宽度pi/16,通带截止频率为pi/3；m=1;%过渡点个数m=1N=ceil(m+1)*2*pi/Bt)+1 %按（7.3.15）式估算采样点数NF=0,wp/pi,wp/pi+2/N,wp/pi+4/N,1;A=1,1,T,0,0;%设置调用参数向量F和Ahn=fir2(N-1,F,A,boxcar(N);%选用矩形窗函数以下与ep731.m相同（省略）。频率采样法的缺点: 滤波器边界频率不易精确控制。 窗函数设计法总使通带和阻带波纹幅度相等，频率采样法只能依靠优化过渡带采样点的取值控制阻带波纹幅度，所以两种方法都不能分

52、别控制通带和阻带波纹幅度。但是工程上对二者的要求是不同的，希望能分别控制。 所设计的滤波器在阻带边界频率附近的衰减最小，距阻带边界频率越远，衰减越大。所以，如果在阻带边界频率附近的衰减刚好达到设计指标要求，则阻带中其他频段的衰减就有很大富余量。这就说明这两种设计法存在较大的资源浪费，或者说所设计滤波器的性能价格比低。7.4 利用等波纹最佳逼近法设计FIR数字滤波器7.4.1 等波纹最佳逼近法的基本思想用Hd()表示希望逼近的幅度特性函数，要求设计线性相位FIR数字滤波器时，Hd()必须满足线性相位约束条件。用Hg()表示实际设计的滤波器幅度特性函数。定义加权误差函数E()为（7.4.1） W(

53、)称为误差加权函数，用来控制不同频段（一般指通带和阻带）的逼近精度。等波纹最佳逼近基于切比雪夫逼近，在通带和阻带以|E()|的最大值最小化为准则，采用Remez多重交换迭代算法求解滤波器系数h(n)。 W()取值越大的频段, 逼近精度越高，开始设计时应根据逼近精度要求确定W()，在Remez多重交换迭代过程中W()是确知函数。 等波纹最佳逼近设计中，把数字频段分为“逼近（或研究）区域”和“无关区域”。逼近区域一般指通带和阻带，而无关区域一般指过渡带。设计过程中只考虑对逼近区域的最佳逼近。应当注意，无关区宽度不能为零，即Hd()不能是理想滤波特性。图7.4.1 等波纹滤波器的幅频特性函数曲线及指

54、标参数（7.4.2） （7.4.3） 由式（7.4.2）和（7.4.3）得到（7.4.4） 按照式（7.4.4）和（7.4.5）计算得到图7.4.1(b)中的参数: 1= 0.1146 ，2= 0.1。实际中, 1和2一般很小，这里为了观察等波纹特性及参数1和2的含义，特意取较大值。（7.4.5） 图7.4.2 误差加权函数W()和滤波器阶数N对逼近精度的影响步骤： （1） 根据给定的逼近指标估算滤波器阶数N和误差加权函数W()； （2） 采用remez算法得到滤波器单位脉冲响应h(n)。MATLAB工具箱函数remezord和remez就是完成以上2个设计步骤的有效函数。7.4.2 reme

55、z和remezord函数及滤波器设计指标1. remez和remezord函数1) remez采用remez算法可实现线性相位FIR数字滤波器的等波纹最佳逼近设计。其调用格式为 hn=remez(M, f, m, w)调用结果返回单位脉冲响应向量hn。remez函数的调用参数（M, f, m, w）一般通过调用remezord函数来计算。调用参数含义如下: M为 FIR数字滤波器阶数，hn长度N=M+1。 f和m给出希望逼近的幅度特性。f为边界频率向量，0f1，要求f为单调增向量（即f(k)f(k+1), k=1, 2, ），而且从0开始，以1结束，1对应数字频率=（模拟频率Fs/2，Fs表示

56、时域采样频率）。m是与f对应的幅度向量，m与f长度相等。如果用命令Plot(f, m)画出幅频响应曲线，则k为奇数时，频段f(k)，f(k+1)上的幅频响应就是期望逼近的幅频响应值，频段f(k+1)，f(k+2)为无关区。起始频段为第一段，奇数频段为逼近区，偶数频段为无关区。例如，对图7.4.2，f=0, 1/4, 5/16, 1；m=1, 1, 0, 0；plot(f, m)画出的幅度特性曲线如图7.4.3所示， w为误差加权向量，其长度为f的一半。w(i)表示对m中第i个逼近频段的误差加权值。图7.4.2(a)中, w=1, 10。缺省w时，默认w为全1（即每个逼近频段的误差加权值相同）。

57、 除了设计选频FIR数字滤波器，remez函数还可以设计两种特殊滤波器：希尔伯特变换器和数字微分器，调用格式分别为 hn=remez(M, f, m, w, hilbert) hn=remez(M, f, m, w, defferentiator)2) remezord采用remezord函数，可根据逼近指标估算等波纹最佳逼近FIR数字滤波器的最低阶数M、误差加权向量w和归一化边界频率向量f，使滤波器在满足指标的前提下造价最低。其返回参数作为remez函数的调用参数。其调用格式为 M, fo, mo, w=remezord(f, m, rip, Fs) 参数说明：f与remez中的类似，这里f

58、可以是模拟频率（单位为Hz）或归一化数字频率，但必须从0开始，到Fs/2（用归一化频率时对应1）结束，而且其中省略了0和Fs/2两个频点。Fs为采样频率，缺省时默认Fs=2Hz。但是这里f的长度（包括省略的0和Fs/2两个频点）是m的两倍，即m中的每个元素表示f给定的一个逼近频段上希望逼近的幅度值。例如，对图7.4.3，f=1/4, 5/16 ，m=1, 0。注意: 省略Fs时，f中必须为归一化频率。 有时估算的阶数M略小，使设计结果达不到指标要求，这时要取M+1或M+2（必须注意对滤波器长度N=M+1的奇偶性要求）。所以必须检验设计结果。 如果无关区（过渡带）太窄，或截止频率太接近零频率和F

59、s/2时，设计结果可能不正确。rip表示f和m描述的各逼近频段允许的波纹幅度（幅频响应最大偏差）. 一般以N, fo, mo, w=remezord(f, m, rip, Fs)返回的参数作为remez的调用参数，计算单位脉冲响应: hn=remez(N, fo, mo, w)。综上所述，调用remez和remezord函数设计线性相位FIR数字滤波器，关键是根据设计指标求出remezord函数的调用参数f、m、rip和Fs，其中Fs一般是题目给定的，或根据实际信号处理要求（按照采样定理）确定。2. 滤波器设计指标1) 低通滤波器设计指标逼近通带: 0, p，通带最大衰减: pdB；逼近阻带:

60、 s, ，阻带最小衰减: sdB。remezord调用参数:（7.4.6）其中，f向量省去了起点频率0和终点频率1，1和2分别为通带和阻带波纹幅度，由式（7.4.4）和（7.4.5）计算得到，下面相同。2) 高通滤波器设计指标逼近通带: p,，通带最大衰减: pdB；逼近阻带: 0, s，阻带最小衰减: sdB。remezord调用参数:（7.4.7）3) 带通滤波器设计指标逼近通带: pl, pu，通带最大衰减: pdB；逼近阻带: 0, sl，su, 阻带最小衰减: sdB。remezord调用参数:（7.4.8）4) 带阻滤波器设计指标逼近阻带: sl, su，阻带最大衰减: sdB；逼