1、1第第6章章 有限脉冲响应数字滤波器的设计有限脉冲响应数字滤波器的设计 IIR的设计SpecificationsDesired IIRButterworth,Chebyshev,Ellipse脉冲响应不变法阶跃响应不变法双线性变换法LP to HP,BP,BS2回顾IIR优点：利用模拟滤波器成熟的理论和设计图表，保留了模拟滤波器优良的幅度特性缺点：只考虑了幅度 特性，没考虑相位特性，一般具有非线性相位特性如果要得到线性相位特性，需要增加相位校正网络，使滤波器变得复杂，增加了成本，而且难以保证严格的线性相位特性引入FIR优点：在满足幅度特性技术要求的同时，很容易实现严格的线性相位特性3FIR滤波
2、器概述系统函数10()()NnnH zh n zH(z)在z平面有N-1个零点，在原点处有一个N-1重极点，因此H(z)永远稳定。稳定和线性相位特性是FIR滤波器的最突出优点。FIR滤波器设计方法不同于IIR。IIR借助模拟滤波器，FIR直接选择有限长度的h(n)，使得频响函数满足要求。三种主要设计方法：窗函数法、频率采样法和Chebyshev等波纹法46.1 线性相位线性相位FIR数字滤波器回顾数字滤波器回顾6.2 利用窗函数法设计利用窗函数法设计FIR滤波器滤波器6.3 利用频率采样法设计利用频率采样法设计FIR滤波器滤波器6.4 利用切比雪夫逼近法设计利用切比雪夫逼近法设计FIR滤波器滤
3、波器6.5 IIR和和FIR数字滤波器的比较数字滤波器的比较6.6 FDATool56.1 线性相位线性相位FIR数字滤波器回顾数字滤波器回顾 考虑长度为N的h(n)，系统函数为：1()0()()()Njj njgnH eh n eHe 什么是线性相位FIR？10()()NnnH zh n z 频率响应函数为：jzeHg()称为幅度特性，()称为相位特性。注意，Hg()不同于|H(ej)|，Hg()为的实函数，可能取负值，而|H(ej)|总是正值。(6.1.1)(6.1.2)6 H(ej)线性相位是指：()是的线性函数，即 为常数 (6.1.3)或()满足下式：,0是起始相位 (6.1.4)严
4、格地说，(7.1.4)中()不具有线性相位，但以上两种情况都满足群时延是一个常数，即()dd 0 第一类线性相位第二类线性相位7第一类线性相位：h(n)是实序列且对(N-1)/2偶对称，即线性相位条件：()(1)h nh Nn第二类线性相位：h(n)是实序列且对(N-1)/2奇对称，即()(1)h nh Nn 注意：充分条件(6.1.5)(6.1.6)87.2 利用窗函数法设计利用窗函数法设计FIR滤波器滤波器 理想低通、高通、带通、带阻数字滤波器幅度特性)(ejH)(ejH)(ejH)(ejH0低通0高通0带通0带阻222222229理想滤波器为Hd(z)()()jj nddnHeh n e
5、基本思路hd(n)是与其对应的单位脉冲响应。频率响应为()()nddnHzh n z求和上下限为什么要取无穷？对是连续的，需要无穷个这样的函数，才能逼近分段连续的理想滤波器的频率响应j ne10Unfortunately,hd(n)序列不可能无限长，我们只能选取其中的N个，组成实际的序列h(n)，比如：()()()dNh nh n wn基本思路wN(n)是长度为N的窗口。得到：注意：通常选择构造长度为N的第一类线性相位FIR滤波器，即要求h(n)关于n=(N-1)/2偶对称。可以选择多种形状的窗口10()()NnnH zh n z11图6.2.1 理想低通的单位脉冲响应及矩形窗如果选择矩形窗口
6、RN(n)窗函数法的一些性质 在时域上，截取了理想滤波器的部分单位单位响应序列，作为实际滤波器的单位脉冲序列 用有限长的h(n)代替无限长的hd(n)，必然引起误差，即Gibbs效应，使得在频域上 过渡带加宽 通带、阻带波动 N越大，误差越小，但系统越复杂、成本越高 设计的基本目标 构造窗函数，减少Gibbs效应 满足线性相位特性12窗函数法导致的频域误差 h(n)对应的频率响应是什么样子？1310()()NnnH zh n z可以这样计算：10()()Njj nnH eh n e也可以这样计算：根据傅立叶变换的频域卷积定理，有()()()dNh nh n wn()()11()()()()()
7、22jjjjjdNdNH eHeWeHeRed 该用哪种方法呢？14wN(n)=RN(n)11002221(1)22221()()1sin(/2)()sin(/2)j NNNjj nj nNNjnnNNNjjjjNjaNjjjeReRn eeeeeeNeReeee矩形窗sin(/2)1(),sin(/2)2NNNRRN()称为矩形窗的幅度函数15RN()的特性RN()在-2/N,2/N区间的一段波形称为主瓣，其余较小的波动称为旁瓣。主瓣旁瓣16设Hd(z)为希望逼近的线性相位低通滤波器，其频率响应为：1,()0,cdcH得到实际滤波器H(z)的频率响应为：1()()()2dNHHRd其中，其幅
8、度特性Hd()为()()jjddHeHe()1()()()21()()2jj ajadNj adNH eHeRedeHRd 最终得到实际滤波器H(z)的幅度特性H()为：线性相位Hd()和RN()的卷积17Hd()和RN()的卷积形成的波形：=0，H(0)为(a)、(b)两波形的积分，即RN()在-c,c之间积分。当c 2/N，近似于RN()在-,之间积分，归一化该积分为1。=c，当c 2/N，近似为RN()一半波形的积分，近似值为1/2=c-2/N，主瓣在-c,c之内，最大负旁瓣在-c,c之外，Hd()有一个最大正峰=c+2/N，主瓣在-c,c之外，最大负旁瓣在-c,c之内，Hd()有一个最
9、大负峰最大正负峰距离4/N18 通过以上分析可知，对hd(n)加矩形窗处理后，H()和原理想低通Hd()差别有以下三点：(1)在理想特性不连续点=c附近形成过渡带。过渡带的宽度，近似等于RN()主瓣宽度，即4/N。(2)通带内增加了波动，最大的峰值在c-2/N处。阻带内产生了余振，最大的负峰在c+2/N处。(3)通带和阻带中波纹的情况与窗函数的幅度函数有关，RN()旁瓣幅度的大小直接影响了H()波纹幅度的大小。选择合适的窗函数会改变旁瓣幅度，因此可以改变H()波纹幅度19 增加矩形窗长度（即加大N）的效果增加N并不能减小波纹幅度 在主瓣附近，RN()可近似为sin(/2)sin()/2NNxR
10、NxN增大时：主瓣幅度增加，旁瓣也增加，相对值不变主瓣、旁瓣距离变小，波动频率加快20 几种典型窗函数四种波形图：窗函数时域波形、幅度特性曲线理想低通滤波器加窗后的单位脉冲响应、幅度特性曲线 定义几个参数：旁瓣峰值n最大旁瓣的最大值相对于主瓣最大值的衰减值(dB)过渡带宽度Bg该窗函数设计的FIR数字滤波器的过渡带宽度阻带最小衰减s该窗函数设计的FIR数字滤波器的阻带最小衰减21 1.矩形窗(Rectangle Window)sin(/2)()/2NNR其频率响应为()()jjaNNReRe幅度特性为指标：n=-13dB;Bg=4/N;s=-21dB22 2.三角形窗(Bartlett Win
11、dow)21,0(1)12()212,(1)112BrnnNNwnnNnNN(6.2.8)其频率响应为 1()22sin()24()sin(/2)NjjBNWeeN(6.2.9)2sin()24()sin(/2)BgNWN幅度特性为23指标：n=-25dB;Bg=8/N;s=-25dB24 3.汉宁(Hanning)窗升余弦窗2()0.51 cos()()1HnNnwnRnN当N1时，N-1N，幅度特性为22()0.5()0.25()()HngNNNWRRRNN12()()NjjHnHngWeWe其频率响应为 旁瓣对消，能量集中在主瓣25图6.2.3 汉宁窗的幅度特性26指标：n=-31dB;
12、Bg=8/N;s=-44dB27 4.哈明(Hamming)窗改进的升余弦窗2()0.540.46cos()()1HmNnnRnN其频响函数WHm(ej)为22()()11()0.54()0.23()0.23()22()0.54()0.23()0.23()11jjjjNNHmNNNHmgNNNWeReReReWRRRNN其幅度函数WHmg()为当N1时，可近似表示为22()0.54()0.23()0.23()HmgNNNWRRRNN能量更集中在主瓣，约占99.96%。是Matlab默认窗函数28指标：n=-41dB;Bg=8/N;s=-53dB29 5.布莱克曼(Blackman)窗24()0
13、.420.5cos0.08cos()11BlNnnnRnNN其频率响应函数为22()()1122()()11()0.42()0.25()()0.04()()jjjjNNBlNNNjjNNNNWeReReReReRe其幅度函数为lg22()0.42()0.25()()11440.04()()11BNNNNNWRRRNNRRNN旁瓣进一步抵消，但过渡带增宽30指标：n=-57dB;Bg=12/N;s=-74dB31常用的窗函数的时域波形32 图6.2.5 常用窗函数的幅度特性(a)矩形窗；(b)巴特利特窗(三角形窗)；(c)汉宁窗；(d)哈明窗；(e)布莱克曼窗 33 图6.2.6 理想低通加窗后
14、的幅度特性(N=51,c=0.5)(a)矩形窗；(b)巴特利特窗(三角形窗)；(c)汉宁窗；(d)哈明窗；(e)布莱克曼窗34 6.凯塞贝塞尔窗(Kaiser-Basel Window)002201()(),01()21(1)11()1()!2kkkIw nnNInNxIxk 式中 I0(x)是零阶第一类修正贝塞尔函数，可用下面级数计算：前五种为参数固定窗函数，旁瓣幅度是固定的。凯塞贝塞尔窗参数可调，是一种最优窗函数。35 一般I0(x)取1525项，便可以满足精度要求。参数可以控制窗的形状。一般加大，主瓣加宽，旁瓣幅度减小，典型数据为49。当=5.44时，窗函数接近哈明窗。=7.865时，窗
15、函数接近布莱克曼窗。凯塞窗的幅度函数为(1)/21()(0)2()cosNkkknWww nn36凯塞窗参数对滤波器的性能影响 37六种窗函数的基本参数 38 (1)根据技术要求确定待求滤波器的单位取样响应hd(n)。如果给出待求滤波器的频响为Hd(ej)，那么单位取样响应用下式求出：1()()2jjddh nHee d22101()()()MjkjknMMMddkhnIDFT HkHeeM根据频率采样定理，hM(n)与hd(n)应满足如下关系：()()()MdMrhnh nrM Rn用窗函数设计FIR滤波器的步骤39 例如，理想低通滤波器的单位取样响应hd(n)为：(2)根据对过渡带及阻带衰
16、减的要求，选择窗函数的形式，并估计窗口长度N。设待求滤波器的过渡带宽度近似等于窗函数主瓣宽度。(3)计算滤波器的单位取样响应h(n)，h(n)=hd(n)w(n)sin()()()cdnh nn (4)验算设计出的滤波器频率响应是否满足指标：10()()Njj nnH eh n e40 例6.2.1 用矩形窗、汉宁窗和布莱克曼窗设计FIR低通滤波器，设N=11,c=0.2rad。解 用理想低通作为逼近滤波器，有sin()(),010()1(1)52sin(0.2(5)(),010(5)cddnh nnnNnh nnn41 用汉宁窗设计：()()(),0102()0.5(1 cos)10dHnH
17、nh nh n wnnnwn用布莱克曼窗设计：11()()()22()(0.420.5cos0.08cos)()1010dBlBlh nh n wnnnwnRn42图6.2.7 例6.2.1的低通幅度特性Using Matlab例：用Hanning窗设计高通FIRDF，要求p=/2，s=/4，ap=1dB，as=40dB。43%HP_Windowwp=.5;ws=.25;Bt=wp-ws;N0=ceil(6.2/Bt);N=N0+mod(N0+1,2);wc=(wp+ws)/2;hn=FIRl(N-1,wc,high,hanning(N);figure(1);stem(0:(length(hn
18、)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);查表可知，Hanning窗满足40dB要求。4400.10.20.30.40.50.60.70.80.91-120-100-80-60-40-20020Frequency(pi)Magnitude(dB)N=25,wc=0.375pi;hn=-0.0004 -0.0006 0.0028 0.0071 -0.0000 -0.0185 -0.0210 0.016
19、5 0.0624 0.0355 -0.1061 -0.2898 0.6249 -0.2898 -0.1061 0.0355 0.0624 0.0165 -0.0210 -0.0185 -0.0000 0.0071 0.0028 -0.0006 -0.00040510152025-0.3-0.2-0.100.10.20.30.40.50.60.7例：对模拟信号进行低通滤波处理。要求p=1.5kHz,s=2.5kHz,ap=1dB,as=40dB。采样频率Fs=10kHz。45%LP_Windowwp=2*1.5/10;ws=2*2.5/10;as=40;Bt=ws-wp;alpha=0.5842
20、*(as-21)+0.07886*(as-21);M=ceil(as-8)/2.285/Bt);wc=(wp+ws)/2;hn=FIRl(M,wc,kaiser(M+1,alpha);figure(1);stem(0:(length(hn)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);用凯塞窗46N=24;hn=0.0039 0.0041 -0.0062 -0.0147 0.0000 0.0286 0.0
21、242 -0.0332 -0.0755 0.0000 0.1966 0.3724 0.3724 0.1966 0.0000 -0.0755 -0.0332 0.0242 0.0286 0.0000 -0.0147 -0.0062 0.0041 0.00390510152025-0.1-0.0500.050.10.150.20.250.30.350.400.10.20.30.40.50.60.70.80.91-80-70-60-50-40-30-20-100Frequency(pi)Magnitude(dB)例：设计带阻FIRDF，要求pl=.2,sl=.35,su=.65,pu=.8，ap=1
22、dB，as=60dB。47%BS_Windowwp=.2,.8;ws=.35,.65;Bt=ws(1)-wp(1);M=ceil(12/Bt)-1;wp=(ws(1)+wp(1)/2,(ws(2)+wp(2)/2;hn=FIRl(M,wp,stop,blackman(M+1);figure(1);stem(0:(length(hn)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);48N=80;0102030
23、405060708090-0.100.10.20.30.40.50.600.10.20.30.40.50.60.70.80.91-100-80-60-40-20020Frequency(pi)Magnitude(dB)496.3 利用频率采样法设计利用频率采样法设计FIR滤波器滤波器 设待设计的滤波器的传输函数用Hd(ej)表示，对它在=0到2之间等间隔采样N点，得到Hd(k)，2210()(),0,1,2,11()IDFT()(),0,1,2,1jddkDNjknNddkHkHekNh nHkHk ekNN再对N点Hd(k)进行IDFT，得到h(n)，50 式中，h(n)作为所设计的滤波器的
24、单位脉冲响应，其系统函数H(z)为 10()()NnnH zh n z12011()()1N NdjkkNzHkH zNez51 1.用频率采样法设计线性相位滤波器的条件 FIR滤波器具有线性相位的条件是h(n)是实序列，且满足h(n)=h(N-n-1)，在此基础上我们已推导出其传输函数应满足的条件是：()()()1()2()(2),odd()(2),evenjjdgggggHeHeNHHNHHN (6.3.5)(6.3.6)(6.3.7)(6.3.4)52 在=02之间等间隔采样N点，2,0,1,2,1kk kNN 将=k代入(7.3.4)(7.3.7)式中，并写成k的函数：()()()1
25、21()2jkdgHkHk eNNkkkNN(6.3.8)(6.3.9)()(),odd()(),evenggggHkHNk kHkHNk k(6.3.10)(6.3.11)53 设用理想低通作为希望设计的滤波器，截止频率为c，采样点数N，Hg(k)和(k)用下面公式计算：N=奇数时，()()1,0,1,2,()0,1,2,11(),0,1,2,1ggcgcccHkHNkkkHkkkkNkNkk kNN(6.3.12)注意：Hg(N)=Hg(0)ckc是通带内最后一个采样点序号54 N=偶数时，()1,0,1,2,()0,1,2,1()1,1,2,1(),0,1,2,1gcgcccgcHkkk
26、HkkkkNkHNkkkNkk kNN (6.3.13)55 2.逼近误差及其改进措施 如果待设计的滤波器为Hd(ej),对应的单位脉冲响应为hd(n)，1()()2jj nddh nHeed 则由频率域采样定理知道，在频域02之间等间隔采样N点，利用IDFT得到的h(n)应是hd(n)以N为周期，周期性延拓乘的主值区序列，即()()()dNrh nh nrN Rn 由于理想滤波器对应的hd(n)序列无限长，因此有限的周期N必然带来时域混叠。N越大，混叠越小。56 由采样定理表明，频率域等间隔采样H(k)，经过IDFT得到h(n)，其Z变换H(z)和H(k)的关系为10122()()()1 s
27、in(/2)()sin(/2)NjkNjH eH kkNNeN12011()()1NNjkNzH kH zNez 频率响应可以写成内插的形式：57在采样点上，逼近误差为0在采样点之间，有误差在Hdg()平滑的区域，误差较小在Hdg()突变的区域，误差较大在频响间断点附近区间内插几个过渡采样点，是不连续点变成过渡带，可增大阻带衰减，但会增加过渡带宽度阻带衰减和过渡带宽度是一对矛盾58 图6.3.1 理想低通滤波器增加过渡点 过渡带过渡带采样点采样点数数m123as4454dB6575dB8595dB59 例6.3.1 利用频率采样法设计线性相位低通滤波器，要求截止频率c=/3，as=40dB，过
28、渡带宽度Bt/16。解：as=40dB时需要增加一个过渡带采样点，即m1。21tBmN 计算N：21tNmBUsing Matlab60%LP_FreqSampleT=input(T=);wp=1/3;Bt=1/16;m=1;N=ceil(m+1)*2/Bt)+1;N=N+mod(N+1,2);F=0,wp,wp+2/N,wp+4/N,1;A=1,1,T,0,0;hn=fir2(N-1,F,A,boxcar(N);figure(1);stem(0:(length(hn)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*
29、log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);61T=.38as=-43.4411dB010203040506070-0.1-0.0500.050.10.150.20.250.30.350.400.10.20.30.40.50.60.70.80.91-100-80-60-40-20020Frequency(pi)Magnitude(dB)00.10.20.30.40.50.60.70.80.91-120-100-80-60-40-20020Frequency(pi)Magnitude(dB)62T=.5as=-29.6896dB
30、00.10.20.30.40.50.60.70.80.91-120-100-80-60-40-20020Frequency(pi)Magnitude(dB)T=.6as=-25.0690dB过渡带插值T不同，as不同窗函数和频率采样法的特点优点：简单方便，易于实现缺点：边界频率不易精确控制 窗函数法使得通带、阻带波纹幅度相等，频率采样法只能控制阻带波纹幅度，它们都不能分别控制通带和阻带的幅度 所设计的滤波器在阻带边界附近衰减最小，距离阻带边界越远，衰减越大。如果阻带边界的衰减刚好满足要求，那么阻带其它频段有很大富裕量，性价比较低63646.4 切比雪夫最佳逼近法切比雪夫最佳逼近法 设希望设计的
31、滤波器幅度特性为Hd()，实际设计的滤波器幅度特性为Hg()，其加权误差E()用下式表示：E()=W()Hd()-Hg()什么是最佳逼近？设计具有线性相位的FIR滤波器，使E()最优其单位脉冲响应序列h(n)必须满足一定条件，例如：h(n)=h(n-N-1)且N为奇数6512()()NjjgH eeH0()()()()cosMdnEWHa nnM=(N-1)/2 minmax(),passbandstopbanda nAEA1(1)20()()cosNgnHa nn最佳逼近问题：选择M+1个系数a(n),使加权误差E()的最大值为最小W()可以在不同频段选择不同权重66 该定理指出最佳一致逼近
32、的充要条件是E()在A上至少呈现M+2个“交错”(alternations)，使得10121()()max(),iiAMEEEA Alternation Theorem属于多项式逼近问题，证明已超出要求67对于低通滤波器0,psA 性能指标11221120lg120lg20lg1psaa 2012020210110110ppsaaaUsing Matlab设计带阻滤波器，通带0,0.2pi、0.8pi,pi，ap=1dB；阻带0.35pi,0.65pi，as=60dB。%BS_Equiripf=0.2,0.35,0.65,0.8;m=1,0,1;ap=1;as=60;delta1=(10(ap
33、/20)-1)/(10(ap/20)+1);delta2=10(-as/20);rip=delta1,delta2,delta1;M,fo,mo,w=remezord(f,m,rip);hn=remez(M,fo,mo,w);figure(1);stem(0:(length(hn)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);69M=28,N=29,fo=0,.2,.35,.65,.8,1,mo=1,1,
34、0,0,1,1,w=1,57.5011,1051015202530-0.100.10.20.30.40.50.600.10.20.30.40.50.60.70.80.91-100-80-60-40-20020Frequency(pi)Magnitude(dB)设计数字低通滤波器,p=1.5kHz,ap=1dB;s=2.5kHz,as=40dB;采样频率10kHz%LP_EquiripFs=10000;f=1500,2500;m=1,0;ap=1;as=40;delta1=(10(ap/20)-1)/(10(ap/20)+1);delta2=10(-as/20);rip=delta1,delta
35、2;M,fo,mo,w=remezord(f,m,rip,Fs);M=M+1;hn=remez(M,fo,mo,w);figure(1);stem(0:(length(hn)-1),hn);fk=0:1/200:1;hk=freqz(hn,1,fk*pi);figure(2);plot(fk,20*log10(abs(hk);xlabel(Frequency(pi);ylabel(Magnitude(dB);71M=15,N=16,w=1,5.7501051015-0.1-0.0500.050.10.150.20.250.30.350.400.10.20.30.40.50.60.70.80.9
36、1-350-300-250-200-150-100-50050Frequency(pi)Magnitude(dB)726.5 IIR和和FIR数字滤波器的比较数字滤波器的比较 性能：IIR滤波器传输函数的极点可位于单位圆内的任何地方，因此可用较低的阶数获得高的选择性，所用的存贮单元少，所以经济而效率高。但是这个高效率是以相位的非线性为代价的对于同样的性能指标，FIR滤波器所要求的阶数一般比IIR滤波器高5-10倍，成本高，信号时延大，但是为线性相位结构：IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定；另外，由于运算过程中的舍入处理，有限字长效应可能会引起寄生振荡FIR滤波器采用非递归结构，在理论和实际中都不存在稳定性问题，运算误差引起的噪声较小；此外，FIR可以采用FFT算法实现，运算速度可以大大提高73设计工具：IIR滤波器可以借助于模拟滤波器的成果，因此一般都有有效的封闭形式的设计公式可供准确计算，计算工作量比较小，对计算工具的要求不高。FIR滤波器计算通带、阻带衰减等无显式表达式，其边界频率不易精确控制；FIR的设计依赖计算工具灵活性：IIR滤波器虽然设计简单，但主要用于低通、高通、带通和带阻等滤波器，脱离不了典型滤波器的频响特性FIR滤波器则要灵活很多，易于适应某些特殊的应用746.6 FDATool Matlab：FDATool
