数字信号处理第三版第七章.ppt

第7章 有限脉冲响应数字滤波器的设计,7.1 线性相位FIR数字滤波器的条件和特点 7.2 利用窗函数法设计FIR滤波器 7.3 利用频率采样法设计FIR滤波器 7.4 利用等波纹最佳逼近法设计FIR滤波器 7.5 IIR和FIR数字滤波器的比较 7.6 几种特殊类型滤波器简介 7.7 滤波器分析设计工具FDATool,7.1 线性相位FIR数字滤波器的条件和特点,本节主要介绍FIR滤波器具有线性相位的条件、 FIR滤波器的幅度特性、零点分布特点和网络结构的特点。 1. 线性相位条件 对于长度为N的h(n)，传输函数为：,式中，Hg()称为幅度特性，()称为相位特性。这里Hg()不同于|H(ej)|，Hg()为的实函数，可能取负值，而|H(ej)|为的正实函数。,H(ej)线性相位： ()是的线性函数，即： ()= -, 为常数 (7.1.3) 若()满足下式： ()=0 -, 0是起始相位 (7.1.4) 也称这种情况为线性相位。 以上两种情况都满足群时延是一个常数，即,一般称：满足(7.1.3)式是第一类线性相位； 满足(7.1.4)式为第二类线性相位。,2. 线性相位FIR的时域约束条件 线性相位FIR滤波器的时域约束条件是指满足线性相位时，对h(n)的约束条件。 1) 第一类线性相位对h(n)的约束条件 相位函数()=，由式（7.1.1）和（7.1.2）得到:,（7.1.5）,(7.1.6) 将（7.1.6）式中两式相除得到：,即 移项并用三角公式化简得到: （7.1.7） 函数h(n)sin(n)关于求和区间的中心(N1)/2奇对称，是满足（7.1.7）式的一组解。 因为sin(n)关于n=奇对称，如果取=(N1)/2，则要求h(n)关于(N1)/2偶对称，所以要求和h(n)满足如下条件:,（7.1.8）,表7.1.1 线性相位FIR数字滤波器的时域和频域特性一览,2) 第二类线性相位对h(n)的约束条件 相位函数()=/2，由式（7.1.1）和（7.1.2），可得到: （7.1.9） 函数h(n)cos(n)关于求和区间的中心(N1)/2奇对称，是满足式（7.1.9）的一组解， 因为cos(n)关于n=偶对称，所以要求和h(n)满足如下条件：,（7.1.10）,2 线性相位FIR滤波器幅度特性Hg()的特点 实质上，幅度特性的特点就是线性相位FIR滤波器的频域约束条件。 引入两个参数符号： ,情况1： h(n)=h(Nn1), N为奇数。 将时域约束条件h(n)=h(Nn1)和()=代入式（7.1.1）和（7.1.2），得到:,所以 （7.1.11） 因为cos(n-)关于=0, , 2三点偶对称，所以由式（7.1.11）可以看出，Hg()关于=0, , 2三点偶对称。因此情况1可以实现各种（低通、高通、带通、带阻）滤波器。,情况2： h(n)=h(Nn1), N为偶数。 仿照情况1的推导方法得到:,（7.1.12）,式中， 。因为是偶数，所以当 时,而且cos(n)关于过零点奇对称，关于=0和2偶对称。所以Hg()=0，Hg()关于=奇对称，关于=0和2偶对称。因此，情况2不能实现高通和带阻滤波器。 情况3： h(n)=h(Nn1)，N为奇数。,将时域约束条件h(n)=h(Nn1)和()= /2代入式（7.1.1）和（7.1.2）， 并考虑 ，得到:,当=0，, 2时，sin(n)=0, 且sin(n)关于过零点奇对称,Hg()关于=0, , 2三点奇对称,情况3只能实现带通滤波器,情况4： h(n)=h(Nn1), N为偶数。 用情况3的推导过程可以得到: （7.1.13） N是偶数，=(N1)/2=N/21/2。所以，当=0, 2时，sin(n)=0；当=时，sin(n)=(1)nN/2, 为峰值点。而且sin(n)关于过零点=0和2两点奇对称，关于峰值点=偶对称。因此Hg()关于=0和2两点奇对称，关于=偶对称。由此可见，情况4不能实现低通和带阻滤波器。,3. 线性相位FIR滤波器零点分布特点 线性相位的系统函数满足：,(7.1.21),“+”和“-”分别对应第一类和第二类线性相位。,一般情况： 线性相位FIR零点分布的特点：互为倒数的共轭对,图7.1.1 线性相位FIR数字滤波器的零点分布,当然，也有一些特殊情况，如图中z1、z2和z4情况。,7.2 利用窗函数法设计FIR滤波器 7.2.1 窗函数法设计原理 设计思想：从时域出发，设计h(n)，逼近hd(n) 。 设希望设计的滤波器传输函数为Hd(ej)，hd(n)是与其对应的单位脉冲响应，因此 ：,一旦Hd(ej)给定，就可得到hd(n) 。 一般情况下， Hd(ej)逐段恒定，在边界频率处有不连续点，因而hd(n)是无限时宽的，且是非因果的，为得到一有限长滤波器h(n)，最直接的方法是截断 hd(n)，即用一个窗口函数RN(n)对hd(n)进行加窗处理： h(n)= hd(n) RN(n)。 并保证截取的一段关于n=（N1）/2偶对称。,相应的单位取样响应hd(n)为：,(7.2.1),它是一个无限长的非因果序列。波形如下图所示。,设理想低通滤波器的传输函数Hd(ej) 为：,(7.2.2),实际设计的滤波器的单位脉冲响应为h(n)，长度为N，其系统函数为H(z)，,用一个有限长的序列h(n)去代替hd(n)，肯定会引起误差,图7.2.1 理想低通的单位脉冲响应及矩形窗,吉布斯（Gibbs）效应：用一个有限长的序列h(n)去代替hd(n) 会引起误差，表现在频域就是吉布斯效应。该效应引起通带内和阻带内的波动性，尤其使阻带的衰减不足，从而不能满足技术上的要求。 吉布斯效应是由于将hd(n)直接截断引起的，因此，也称为截断效应。,图7.2.2 吉普斯效应,Hd(ej)是一个以2为周期的函数，可以展为傅里叶级数，即 傅里叶级数的系数为hd(n)，当然就是Hd(ej)对应的单位脉冲响应。设计FIR滤波器就是根据要求找到N个傅里叶级数系数h(n)，n=1, 2, , N1，以N项傅氏级数去近似代替无限项傅氏级数，这样在一些频率不连续点附近会引起较大误差，这种误差就是前面说的截断效应。,因此，从这一角度来说，窗函数法也称为傅氏级数法。显然，选取傅氏级数的项数愈多，引起的误差就愈小，但项数增多即h(n)长度增加，也使成本和滤波计算量加大，应在满足技术要求的条件下，尽量减小h(n)的长度。,（7.2.4）,根据傅里叶变换的时域卷积定理，得到（7.2.3）式的傅里叶变换：,（7.2.5）,第一过零点内为主瓣,旁瓣,将Hd(ej)写成Hd(ej)=Hdg()ej， 则按照（7.2.1）式，理想低通滤波器的幅度特性函数为,将Hd(ej)和WR(ej)代入（7.2.4）式，得到：,将H(ej)写成H(ej)=Hg()ej ，则,（7.2.6）,加窗后的滤波器的幅度特性等于理想低通滤波器的幅度特性Hdg()与 矩形窗幅度特性WRg()的卷积。,图7.2.3 矩形窗加窗效应,当=c时，如(c)所示，当c 2/N时，积分近似为WRg()一半波形的积分，对Hg(0)归一化后的值近似为1/2,当=0时，Hg(0)等于图(a)与(b)两波形乘积的积分，相当于对WRg()在c之间一段波形的积分，当c2/N时，近似为之间波形的积分,当=c2/N时，如(d)所示，WR()主瓣完全在区间c, c之内，而最大的一个负旁瓣移到区间c, c之外，因此Hg(c2/N)有一个最大的正峰。,当=c+2/N时，如(e)所示，WRg()主瓣完全移到积分区间外边，由于最大的一个负旁瓣完全在区间c, c内，因此Hg(c+2/N)形成最大的负峰。,Hg()最大的正峰与最大的负峰对应的频率相距4/N,对hd(n)加矩形窗处理后，Hg()与原理想低通Hdg()的差别有以下两点：,在理想特性不连续点=c附近形成过渡带。过渡带的宽度近似等于WRg()主瓣宽度4/N。 通带内产生了波纹，最大的峰值在c2/N处。阻带内产生了余振，最大的负峰在c+2/N处。通带与阻带中波纹的情况与窗函数的幅度谱有关， WRg()旁瓣幅度的大小直接影响Hg()波纹幅度的大小。,这两点就是对hd(n)截断后，在频域的反映，即吉布斯效应。,如何减少吉布斯效应的影响，设计一个满足要求的FIR滤波器呢？ 直观上，增加矩形窗口的宽度（即加大N）可以减少吉布斯效应的影响。N 时， 在主瓣附近， WRg()近似为：,该函数的性质是： 随x (N )，主瓣幅度加高，同时旁瓣也加高，保持主瓣和旁瓣幅度相对值不变；波动的频率加快，当x 时，sinx/x趋近于 函数，因此，当N加大时， H()的波动幅度没有多大改善。 N加大带来的最大好处是过渡带变窄(过渡带:4/N)。 因此加大N并不是减少吉布斯效应的有效方法。,图7.2.4 矩形窗函数长度的影响,所以有如下结论： 调整窗口长度N可以有效地控制过渡带的宽度。 减少带内波动以及加大阻带的衰减只能从窗函数的形状上找解决方法。 为了消除吉布斯效应，取得较好的频率特性，一般采用其它类型的窗口，这些窗口函数在靠近两端的h(n)逐步衰减到零，使窗口频谱的主瓣包含更多的能量。从而降低通带内的纹波并加大阻带的衰减。 下面是几种常用的窗函数,1. 矩形窗(Rectangle Window),2. 三角形窗(Bartlett Window),3. 汉宁(Hanning)窗-升余弦窗,4. 哈明(Hamming)窗改进的升余弦窗,5. 布莱克曼(Blackman)窗,6. 凯塞贝塞尔窗(Kaiser-Basel Window),7.2.2. 几种常用的窗函数,（7.2.7）,为了描述方便，定义窗函数的几个参数: 旁瓣峰值n窗函数的幅频函数|Wg()|的最大旁瓣的最大值相对主瓣最大值的衰减值（dB）； 过渡带宽度Bg用该窗函数设计的FIR数字滤波器（FIRDF）的过渡带宽度； 阻带最小衰减s用该窗函数设计的FIRDF的阻带最小衰减。,1 矩形窗（Rectangle Window） wR(n)=RN(n) 其幅度函数为,图7.2.4 矩形窗函数长度的影响,2 三角形窗（Bartlett Window）,（7.2.8）,其频谱函数为,其幅度函数为 （7.2.10）,图7.2.5 三角窗的四种波形,主瓣能量更集中，但过渡带宽了很多,3 汉宁（Hanning）窗升余弦窗,（7.2.11）,当N1时， N1N,汉宁窗的幅度函数WHng()由三部分相加，旁瓣互相对消，使能量更集中在主瓣中。汉宁窗的四种波形如图7.2.6所示，参数为: n=31 dB; Bg=8/N; s=44 dB。,图7.2.6 汉宁窗的四种波形,4 哈明（Hamming）窗改进的升余弦窗 （7.2.12）其频谱函数WHm(ej)为,其幅度函数WHmg()为 当N时，其可近似表示为 ,图7.2.7 哈明窗的四种波形,这种改进的升余弦窗，能量更加集中在主瓣，但其主瓣宽度和汉宁窗的相同,（7.2.13）,5 布莱克曼（Blackman）窗 其频谱函数为,其幅度函数为,幅度函数由五部分组成，它们都是移位不同，且幅度也不同的WRg()函数，使旁瓣再进一步抵消。旁瓣峰值幅度进一步增加，其幅度谱主瓣宽度是矩形窗的3倍,图7.2.8 布莱克曼窗的四种波形,旁瓣进一步减小，但过渡带变宽,6 凯塞贝塞尔窗（Kaiser-Basel Window） 以上五种窗函数都称为参数固定窗函数，每种窗函数的旁瓣幅度都是固定的。凯塞贝塞尔窗是一种参数可调的窗函数，是一种最优窗函数。 （7.2.15） 式中 ,I0()是零阶第一类修正贝塞尔函数，可用下面级数计算： 一般I0()取1525项，便可以满足精度要求。 参数可以控制窗的形状。一般加大，主瓣加宽，旁瓣幅度减小，典型数据为4 9。当 =5.44时，窗函数接近哈明窗。 =7.865时，窗函数接近布莱克曼窗。在设计指标给定时，可以调整值，使滤波器阶数最低，所以其性能最优。凯塞（Kaiser）给出的估算和滤波器阶数N的公式如下:,表7.2.1 凯塞窗参数对滤波器的性能影,表7.2.2 6种窗函数的基本参数,7.2.3 用窗函数法设计FIR滤波器的步骤 根据对过渡带及阻带衰减的指标要求，选择窗函数的类型，并估计窗口长度N。,(2) 构造希望逼近的频率响应函数Hd(ej)，即,一般取,(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)。 对（7.2.19）式给出的线性相位理想低通滤波器作为Hd(ej)，由（7.2.2）式求出单位脉冲响应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。,（2） 构造Hd(ej): 式中 ,（3） 求出hd(n):,将=12代入得,（4） 加窗:,全通滤波器,理想低通滤波器的单位脉冲响应,例7.2.2 对模拟信号进行低通滤波处理，要求通带0f 1.5kHz内衰减小于1 dB， 阻带2.5 kHzf 上衰减大于40 dB。希望对模拟信号采样后用线性相位FIR数字滤波器实现上述滤波，采样频率Fs=10 kHz。用窗函数法设计满足要求的FIR数字低通滤波器，求出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，所以由式（7.2.2）和式（7.2.3）, 得到: 式中，w(n)是长度为24（ =3.395）的凯塞窗函数。,7.3 利用频率采样法设计FIR滤波器 1 用频率采样法设计滤波器的基本思想 设希望逼近的滤波器的频响函数Hd(ej) ，对它=0到2之间等间隔采样N点，得到Hd(k)：,再对Hd(k)进行N点IDFT，得到h(n)： （7.3.2） h(n)作为所设计的FIR滤波器的单位脉冲响应，其系统函数H(z)为 （7.3.3）,根据频率域采样理论，利用频率域采样值恢复原信号Z变换式，得到H(z)的内插表示形式： （7.3.4） 此式就是直接利用频率采样值Hd(k)形成滤波器的系统函数，（7.3.3）式适合FIR直接型网络结构，（7.3.4）式适合频率采样结构。,下面讨论两个问题：一个是为了设计线性相位FIR滤波器，频域采样序列Hd(k)应满足的条件；另一个是逼近误差问题及其改进措施。,2 设计线性相位滤波器时对Hd(k)的约束条件 FIR滤波器具有线性相位的条件是h(n)为实序列，且满足h(n)=h(Nn1)，其频响函数应满足的条件是：,（7.3.5）,（7.3.6）,在=02区间上N个等间隔的采样频点为 将=k代入（7.3.5）（7.3.8）式中，并写成k的函数：,k = 0，1，2，N 1,（7.3.9）,（7.3.10）,Hg(k)关于N/2点偶对称,Hg(k)关于N/2点奇对称，且Hg(N/2)=0, 设用理想低通作为希望逼近的滤波器Hd(ej)，截止频率为c，采样点数为N，Hg(k)和(k)用下列公式计算：,（7.3.13）,N=偶数时，,N=奇数时，,kc是通带内最后一个采样点的序号，所以kc的值取不大于cN/(2)的最大整数。对于高通和带阻滤波器，N只能取奇数。,3 逼近误差及其改进措施 时域角度分析： 待逼近的滤波器为Hd(ej)，对应的单位脉冲响应为hd(n)，在频域02范围等间隔采样N点，利用IDFT得到的h(n)应是hd(n)以N为周期的周期延拓的主值区序列，即 ,由于时域混叠及截断，使h(n)与hd(n)有偏差,频域角度分析：,H(ej)=FTh(n)的内插表示形式：,式中,在采样频点，(2k/N)=1，因此采样点处 与H(k)相等，逼近误差为0。在采样点之间， 由N项H(k)(2k/N)之和形成。,图7.3.1 频域幅度采样序列Hg(k)及其内插波形Hg(),dg()的误差与Hdg()特性的平滑程度有关，Hdg()特性愈平滑的区域，误差愈小；特性曲线间断点处，误差最大。,间断点变成倾斜下降的过渡带曲线，过渡带宽度近似为2/N。通带和阻带内产生震荡波纹，且间断点附近振荡幅度最大，使阻带衰减减小，往往不能满足技术要求。,增加N可以使过渡带变窄，但是通带最大衰减和阻带最小衰减随N的增大并无明显改善。且N太大，会增加滤波器的阶数，即增加了运算量和成本。,提高阻带衰减的具体方法是在频响间断点附近区间内插一个或几个过渡采样点，使不连续点变成缓慢过渡带，这样，虽然加大了过渡带，但阻带中相邻内插函数的旁瓣正负对消，明显增大了阻带衰减。,表7.3.1 过渡带采样点的个数m与滤波器阻带最小衰减 s 的经验数据,随着m， s明显,4 频率采样法设计步骤 （1） 根据阻带最小衰减 s选择过渡带采样点的个数m。 （2） 确定过渡带宽度Bt，估算频域采样点数（即滤波器长度）N。如果给定过渡带宽度Bt，则要求(m+1)2/NBt ，滤波器长度N必须满足如下估算公式:,（7.3.15）