1、有限长单位冲激响应数字滤波器的特点有限长单位冲激响应数字滤波器的特点:有限长单位冲激响应(FIR)可以做成具有严格的线性相位,同时又可以具有任意的幅度特性。FIR滤波器的单位抽样响应是有限长的,因而FIR滤波器一定是稳定的。只要经过一定的延时,任何非因果有限长序列都能变成因果的有限长序列,总能用因果系统来实现。FIR滤波器由于单位冲激响应是有限长的,因而可以用快速傅里叶变换(FFT)算法来实现过滤信号,从而可大大提高运算效率。但是,要取得很好的衰减特性,FIR滤波器的阶次比HR滤波器的要高。本章主要讨论线性相位滤波器的设计。本章主要讨论线性相位滤波器的设计。7.1 线性相位线性相位FIR滤波器
2、的特点滤波器的特点7.1.1 线性相位条件线性相位条件如果一个线性移不变系统的频率响应有如下形式:(7.1)则其具有线性相位。这里 是一个实数。因而,线性相位系统有一个恒定的群延时 (7.2)()()()|()|jjjjH eHeH ee 在实际应用中,有两类准确的线性相位,分别要求满足 (7.3)(7.4)FIR滤波器具有式(7.3)的线性相位的充分必要条件是:单位抽样响应 关于群延时 偶对称,即满足(7.5)(7.6)()()()h n12N()(1)01h nh NnnN 满足式(7.5)和式(7.6)的偶对称条件的FIR滤波器分别称为I型线性相位滤波器和型线性相位滤波器。FIR滤波器具
3、有式(7.4)的线性相位的充分必要条件是:单位抽样响应 关于群延时 奇对称,即满足(7.7)(7.8)(7.9)()h n12N2()(1)01h nh NnnN 把满足式(7.7)、(7.8)和式(7.9)的奇对称条件的FIR滤波器分别称为型线性相位滤波器和型线性相位滤波器。1I型线性相位滤波器型线性相位滤波器 7.1.2 线性相位滤波器频率响应的特点线性相位滤波器频率响应的特点由于偶对称性,一个I型线性相位滤波器的频率响应可表示为(7.10)其中(1)/2(1)/20()()cos()Njj NnH eea kk11()2()1,2,.,22NNa khkk1(0)()2Nah幅度函数为
4、(7.11)相位函数为 (7.12)(1)/20()()cos()NnHa kk()(1)2N I型线性相位滤波器的幅度函数和相位函数的特点:幅度函数和相位函数的特点:幅度函数幅度函数对 偶对称,同时对 也呈偶对称;相位函数相位函数为准确的线性相位。12N0,22型线性相位滤波器型线性相位滤波器 一个型线性相位滤波器,由于N是偶数,所以,的对称中心在半整数点 。其频率响应可以表示为:(7.13)其中()h n12N/2(1)/20()1()cos()2Njj NnH eeb kk()2()1,2,.,22NNb khkk幅度函数为 (7.14)相位函数为 (7.15)()(1)2N/20()1
5、()cos()2NnHb kk型线性相位滤波器的幅度函数和相位函数的特点:幅度函数和相位函数的特点:幅度函数的特点:幅度函数的特点:(1)当 时,=0,也就是说 在 处必然有一个零点;(2)对 呈奇对称,对 呈偶对称。相位函数的特点:同相位函数的特点:同I I型线性相位滤波器。型线性相位滤波器。()H()Hz1z ()H0,23型线性相位滤波器型线性相位滤波器 由于型线性相位滤波器关于 奇对称,且 为整数,所以,其频率响应可以表示为 (7.16)其中 12N(1)/2(1)/21()()sin()Njj NnH ejec kk11()2()1,2,.,22NNc khkk幅度函数为 (7.17
6、)相位函数为 (7.18)(1)/21()()sin()NnHc kk()(1)22N 型线性相位滤波器的幅度函数和相位函数的特点幅度函数和相位函数的特点:幅度函数的特点:幅度函数的特点:(1)当 时,=0,也就是说 在 处都为零点;(2)对 均呈奇对称。相位函数的特点:相位函数的特点:既是准确的线性相位,又包括 的相移,所以又称 移相器,或称正交变换网络。()H()Hz()H0,21z 0,2/2904型线性相位滤波器型线性相位滤波器 型线性相位滤波器关于 奇对称,且N为偶数,所以为非整数。其频率响应可以表示为 (7.19)其中 12N/2(1)/21()1()sin()2Njj NnH e
7、jed kk()2()1,2,.,22NNd khkk幅度函数为 (7.20)相位函数为 (7.21)()(1)22N/21()1()sin()2NnHd kk型线性相位滤波器的幅度函数和相位函数的特点幅度函数和相位函数的特点:幅度函数的特点幅度函数的特点:(1)在 处必为零,也就是说 在 处为零点;(2)在 处呈奇对称,在 处呈偶对称 相位函数相位函数的特点:同型线性相位滤波器。()Hz()H0,2()H1z 0,27.1.3 零点位置零点位置对于I型或型线性相位滤波器,意味着()(1)h nh Nn(1)1()()NH zzH z对于型或型线性相位滤波器,意味着()(1)h nh Nn (
8、1)1()()NH zzH z 在上述两种情况下,如果 在 处等于零,则在处也一定等于零。所以 的零点呈倒数对出现。另外,若 是实值的,则复零点呈共轭倒数对出现,或者说是共轭镜像的。()H z0zz01/zz()H z()h n一个线性相位滤波器零点的约束条件一个线性相位滤波器零点的约束条件 线性相位滤波器的级联结构实现 5N 122N在此情况下,111112212211()(1)(1)(1)(1)1 12(cos)12(cos)iiiijjjjiiiiiiiiiiiH zz rez rezezerrrzr zrzzr(1)零点 既不在实轴上,也不在单位圆上,零点是两组互为倒数的共轭对,其基本
9、因子为 izz(7.22),1,0ijiiiizrer在此情况下,3N 112N(2)零点 在单位圆上,但不在实轴上,即 ,此时零点的共轭值就是他的倒数,其基本因子为 izz(7.23)1112()(1)(1)12(cos)iijjiiH zz ez ezz 1,0,iiir在此情况下,3N 112N(3)零点 在实轴上,但不在单位圆上,即 ,此时零点是实数,他没有复共轭部分,只有倒数,倒数也在实轴上,其基本因子为 izz(7.24)1,0 iir或 111211()(1)(1)1()iiiiiH zrzzrzzrr 式中“”号相当于 ,零点在负实轴上,“”相当于零点在正实轴上。i0i(4)零
10、点 既在单位圆上,但在实轴上,即 ,此时零点只有两种情况,即 ,这时零点既是自己的复共轭,又是倒数,其基本因子为 izz(7.25)1,0iir 或1()1iH zz 式中“”号相当于 ,零点在负实轴上,“”相当于零点在正实轴上。i0i在此情况下,即有半个抽样的延时。2N 11/22N线性相位滤波器只能由以上这几种因子的组合而成。11zz 或7.2 窗函数设计法窗函数设计法 7.2.1 设计方法设计方法 给出所要求的理想低滤波器频率响应()jdeH设计一个FIR滤波器频率响应 10()()Njj nneeHh n逼近()jdeH窗函数设计法设计是在时域进行(7.26)()jdeH由 的傅里叶反
11、变换导出 sin()11()()22ccjj nj ncddeeenh nHddn()jdeH()dh n()h n由于 是矩形频率特性,故 一定是无限长的序列,而且是非因果的,而要设计的是FIR滤波器,必然 是有限长的。所以要用有限长的 来逼近无限长的 ,最有效的方法是截断 ,即用一个有限长度的窗函数序列 来截取 ,并将截短后的 移位,得()h n()dh n()dh n()n()dh n()dh n11()()()22dNNh nnh n(7.27)窗函数序列的形状及长度的选择很关键。取 ,即取矩形窗 例例7.1 设计一低通滤波器,所希望的频率响应截止频率 在 之间为1,在 之间为0,分别
12、取N=11,21,41,观察其频谱响应的特点。()jdeH00.250.25解:解:(1)/2 00.25()0 0.25j NjdeeH1 01()0 nNn其它 由式(7.27)1sin0.25()12()()12()2dNnNh nh nNn当N=11时,求得(0)(10)0.045,(1)(9)0,(2)(8)0.075(3)(7)0.1592,(4)(6)0.2251,(5)0.25hhhhhhhhhhh 显然 ,满足对称关系。152N根据序列 ,分别求得N=11,21,41时的幅频特性()h n|()|jeH由图可以看出,当N取的过小时,通频带过窄,且阻带内波纹较大,过渡带较宽,当
13、N增大时,与 的近似程度越来越好。但当N增大时,通带内出现了波纹,而且随着N的继续增大,这些波纹并不消失,只是最大的尖峰处越来越接近于间断点,这种现象称作吉布斯现象。()jeH()jdeH吉布斯现象的产生是由于对 突然截短的结果。()dh n为了减少吉布斯现象,应选取旁瓣较小的窗函数。1矩形窗矩形窗 7.2.2 各种窗函数各种窗函数 窗函数为()()NnRn(7.28)幅度函数为(7.29)sin()2()|()|sin()2jRRNWW e主瓣宽度 ,过渡带宽 。4/2 2/NN0.9 2/N2汉宁(汉宁(Hanning)窗(又称升余弦窗)窗(又称升余弦窗)窗函数为 2()0.50.5cos
14、()()1NnnRnN(7.30)幅度函数为(7.31)22()0.5()0.25()()11RRRnnWWWWNN主瓣宽度 ,过渡带宽 。4 2/8/NN3.1 2/N3海明(海明(Hamming)窗(又称改进的升余弦窗)窗(又称改进的升余弦窗)窗函数为(7.32)幅度函数为(7.33)主瓣宽度 ,过渡带宽 。4.4.凯泽(凯泽(Kaiser)窗)窗 窗函数为(7.34)2()0.540.46cos()()1NnnRnN22()0.54()0.23()()11RRRnnWWWWNN4 2/8/NN3.3 2/N20021(1)1(),01()nINnnNI其中 为第一类变形零阶贝塞尔函数,是
15、一个可自由选择的参数,改变 值就可对主瓣宽度与旁瓣衰减进行选择,一般选择 。过渡带宽 。0I495 2/N 窗函数窗谱性能指标加窗后滤波器性能指标旁瓣峰值(dB)主瓣宽度()过渡带宽()阻带最小衰减(dB)矩形窗汉宁窗海明窗凯泽窗-13-31-41-57244 0.93.13.35-21-44-53-80最小阻带衰减只由窗形决定,不受N的影响,而过渡带宽则随N的增加而减小。表表7.1 几种窗函数的基本参数比较几种窗函数的基本参数比较2/N2/N1高通数字滤波器的设计高通数字滤波器的设计 7.2.3 其他各型其他各型FIR滤波器的设计方法滤波器的设计方法 令(7.35)(1)/2 ()0 0j
16、NjcdceeH则 11()()2211()22ccNNjnjndeeh ndd求得 11sinsin22()1()2cdNNnnh nNn(7.36)2高通数字滤波器的设计高通数字滤波器的设计 令(7.37)则 求得 (7.38)(1)/2 ()0 j NjlhdeeH其它11()()2211()22lhhlNNjnjndeeh ndd11sinsin22()1()2hldNNnnh nNn3带阻数字滤波器的设计带阻数字滤波器的设计 令(7.39)则 求得 (7.40)(1)/2|,|()0 j NjlhdeeH其它 111()()()222111()222hllhNNNjnjnjndeee
17、h nddd111sinsinsin222()1()2lhdNNNnnnh nNn比较式(7.36)、(7.38)、(7.40)可知,一个高通滤波器相当于用一个全通滤波器减去一个低通滤波器;一个带通滤波器相当于两个低通滤波器相减,其中一个截止频率在 ,另一个在 ;一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器,低通滤波器的截止频率 ,高通滤波器在 。hlhl选取一个满意的窗函数,令()()(),0,1,.,1dh nn h nnN(7.41)则 即为要设计的滤波器的单位抽样响应。()h n按上述方法设计的滤波器,由于满足了 的对称关系,因此都具有线性相位。()(1)h nh Nn 7.3
18、 频率抽样设计法频率抽样设计法 频率抽样法是从频域出发,把给定的理想频率响应 加以等间隔抽样,即()jdeH2()()jddnkeHHk令 2()()(),0,1,.,1jddnkeH kHkHkN(7.42)由DFT定义,得2101()(),0,1,.,1NjnkNdkh nHk ekNN(7.43)可求得滤波器的系统函数(7.44)21110002110012011()()()1 ()11 ()1NNNjnknnNdnnkNNjnknNdknNNdjnkkNH zh n zHk ezNHkezNzHkNez该系统的频率响应为12011()()()1jNNjdjjnkz ekjNeH eH
19、zHkNee(7.45)经过推导,有(7.46)由式(7.46)可知,是由内插函数(7.47)1(1)/2(1)/0sin(2/)/2()()sin(2/)/2Njj Nj NkNdkNk NH eeHk eNk N()jeH(1)2sin2()sin2NjNeN的插值所决定的,即 102()()()NjdkH eHkkN 由内插公式(7.47)可知,在各频率抽样点上,滤波器的实际频率响应严格地和理想频率响应值相等。但是在抽样点之间的频率响应则是由N个离散值 作为权重和插值函数 线性组合的结果。显然抽样点N取得越大,近似程度越好,N的选取要视在通带和阻带内的技术要求而定。()dHk()(7.4
20、7)102()()()NjdkH eHkkN 的指定原则的指定原则()dHk(3)由 求出的 应具有线性相位。(1)在通带内可令|=1,阻带内|=0,且在通带内赋给 一相位函数;()dHk()dHk()dHk(2)指定的 应保证由式(7.43)求出的 是实序列;()dHk()h n()h n()jeH 的指定的指定()dHk由式(7.46)知,若保证(1)/()j NkNdHk e实数则 就具有线性相位,。()jeH1()2N 并考虑|=1,等效地指定()dHk(1)/(),0,1,.,1j NkNdHkekN(7.48)根据DFT的性质可知,为保证 是实序列,应满足下列对称关系 ()dHk*
21、()()()dddHkHkHNk(7.49)由于(1)()/(1)(1)/(1)*()()j NN kNj Nj NkNj NddeeeeHNkHk(7.50)当N为偶数时,;当N为奇数时,。这样当N为偶数时,若按式(7.48)对 赋值,就不能满足式(7.49)的对应关系。由此,按如下原则对 赋值。(1)1j Ne(1)1j Ne()dHk()dHkN为偶数时(1)/(1)/0,1,.,/2 1()0 /2 /2 1,.,1j NkNdj NkNeekNHkkNkNNN为奇数时(1)/()0,1,.,1j NkNdeHkkN(7.51)(7.52)用频率抽样法设计用频率抽样法设计FIR数字滤波
22、器的步骤:数字滤波器的步骤:(1)根据所设计的滤波器的通带与阻带的要求,根据N为偶数还是奇数,按式(7.51)、(7.52)指定 ,在阻带内,=0;()dHk()dHk(2)由指定的 构成所设计的滤波器的转移函数,也可由式(7.44)求得频率响应 。()dHk()jeH例例7.2 用频率抽样法设计一个低通滤波器,其截止频率是抽样频率的1/10,取N=20。解:解:此处N为偶数,且在通带内对 抽样时,仅得两个点,由式(7.51),有()jeH(0)1dH19/20(1)jdeH*19/20(19)(20 1)(1)jdddeHHH在其它点处,()0dHk 将 代入式(7.43),得 序列如下()
23、dHk()h n(0)(19)0.04877hh(1)(18)0.0391hh(2)(17)0.0207hh(3)(16)0.0046hh(4)(15)0.03436hh(5)(14)0.0656hh(6)(13)0.0954hh(7)(12)0.12071hh(8)(11)0.1391hh(9)(10)0.14877hh7.4 应用应用MATLABMATLAB设计设计FIRFIR数字滤波器数字滤波器 1窗函数窗函数(1)bartlett.m (三角窗)(三角窗)(2)blackman.m (布莱克曼窗)(布莱克曼窗)(3)boxcar.m(矩形窗)(矩形窗)(4)hamming.m(海明窗)
24、(海明窗)(5)hanning.m(汉宁窗)(汉宁窗)(6)triang.m(三角窗)(三角窗)(7)chebwin.m(切比雪夫窗)(切比雪夫窗)(8)kaiser.m(凯泽窗)(凯泽窗)7.4.1 与本章内容有关的与本章内容有关的MATLABMATLAB文件文件 (1)fir1.m本文件采用窗函数法设计FIR数字滤波器,其调用格式是1)b=fir1(N,Wn)2)b=fir1(N,Wn,high)3)b=fir1(N,Wn,stop)2FIR数字滤波器的文件数字滤波器的文件 式中N为滤波器的阶次,因此滤波器的长度为N+1;Wn是通带截止频率,其值在01之间,1对应抽样频率的一半;b是设计好
25、的滤波器系数。对于格式(1),若Wn是一标量,则可用来设计低通滤波器;若Wn是 的向量,则用来设计带通滤波器;若Wn是 的向量,则可用来设计带滤波器,此时,格式将变为:b=fir1(N,Wn,DC-1)或b=fir1(N,Wn,DC-0)格式(2)用来设计高通滤波器;格式(3)用来设计带阻滤波器。21L1(2)fir2.m本文件采用窗函数法设计具有任意幅频特性的FIR数字滤波器。其调用格式是b=fir1(N,F,M)其中F是频率向量,其值在01之间,M是与F相对应的所希望的幅频响应。不指定窗函数的类型,则自动选择汉明窗。(3)remez.m本文件用来设计采用切比雪夫最佳一致逼近FIR数字滤波器
26、。同时,还可以用来设计希尔伯特变换器和差分器。其调用格式是1)b=remez(N,F,A)2)b=remez(N,F,A,W)3)b=remez(N,F,A,W,hilbert)4)b=remez(N,F,A,W,differentiator)其中,N是给定的滤波器的阶次;b是设计的滤波器的系数,其长度为N+1;F是频率向量,其值在01之间;A是对应F的各频段上的理想幅频响应;W是各频段上的加权向量。注意:若注意:若b的长度为偶数,设计高通和带阻滤波器时有可能出现错误,的长度为偶数,设计高通和带阻滤波器时有可能出现错误,因此最好保证因此最好保证b的长度为奇数,即的长度为奇数,即N应为偶数。应为
27、偶数。(4)remexord.m 本文件采用切比雪夫一致逼近设计FIR数字滤波器时所需要的滤波器阶次。其调用格式是N,Fo,Ao,W=remexord(F,A,DEV,Fs)式中,F、A的含义同文件(3),是通带和阻带上的偏差;该文件输出的是符合要求的滤波器阶次N、频率向量Fo、幅度向量Ao和加权向量W。若设计者事先不能确定自己要设计的滤波器的阶次,那么,调用remexord后,就可利用这一族参数再调用remez,即b=remez(N,Fo,Ao,W),从而设计出所需要的滤波器。因此,通常remez和remexord结合使用。说明:说明:remexord给出的阶次给出的阶次N有可能偏低,这时适
28、当增加有可能偏低,这时适当增加N即可;另外,即可;另外,若若N为奇数,就可令其加为奇数,就可令其加1,使其变为偶数,这样,使其变为偶数,这样b的长度为奇数。的长度为奇数。(5)sgolay.m本文件用来设计Savitzky-Golay平滑滤波器。其调用格式是b=sgolay(k,f)式中k是多项式的阶次,f是拟合的双边点数。要求 kf,且f为奇数。(6)firls.m本文件用最小平方法设计线性相位FIR数字滤波器。可设计任意给定的理想幅频特性。(7)fircls.m用带约束的最小平方法设计线性相位FIR数字滤波器。可设计任意给定的理想幅频特性。(8)fircls1.m用带约束最小平方法设计线性
29、相位FIR低通和高通滤波器。可设计任意给定的理想幅频特性。(9)firrcos.m用来设计低通线性相位FIR数字滤波器,其过渡带为余弦函数形状。7.4.2 应用应用MATLABMATLAB设计设计FIRFIR数字滤波器数字滤波器 例例7.3 令N=10,分别用矩形窗和海明窗重复例7.1。解:解:根据要求编制MATLAB程序如下:clear all;N=10;b1=fir1(N,0.25,boxcar(N+1);%用矩形窗作为冲激响应的窗函数b2=fir1(N,0.25,hamming(N+1);%用Hamming窗作为冲激响应的窗函数%M=128;h1=freqz(b1,1,M);h2=fre
30、qz(b2,1,M);%分别求两个滤波器的频率响应;t=0:10;subplot(221)stem(t,b2,.);hold on;plot(t,zeros(1,11);grid;f=0:0.5/M:0.5-0.5/M;M1=M/4;for k=1:M1 hd(k)=1;hd(k+M1)=0;hd(k+2*M1)=0;hd(k+3*M1)=0;endsubplot(222)plot(f,abs(h1),b-,f,abs(h2),g-,f,hd,-);grid;运行结果如图7.5所示。图图7.5 运行结果运行结果 例例7.4设计一多带滤波器,要求理想幅频响应在归一化频率0.20.3,0.60.8
31、之间为1,其余均为0。解:解:程序如下:clear all;f=0 0.19 0.2 0.3 0.31 0.59 0.6 0.8 0.81 1;%给定频率轴分点;m=0 0 1 1 0 0 1 1 0 0;%给定在这些频率分点上理想的幅频响应N1=30;N2=90;%取两种不同的滤波器长度;b1=fir2(N1,f,m);b2=fir2(N2,f,m);%得到两个滤波器;subplot(311);stem(b1,.);grid;subplot(312);stem(b2,.);grid;M=128;h1,w=freqz(b1,1,M,1);h2,w=freqz(b2,1,M,1);subplot
32、(313);plot(w,abs(h1),b-,w,abs(h2),g-);grid;图7.6 滤波器单位抽样响应及其幅频响应曲线例例7.5 利用切比雪夫最佳一致逼近法设计一低通滤波器,要求通带边缘频率 ,阻带边缘频率 。解:解:clear all;f=0.6.7 1;%给定频率轴分点;A=1 1 0 0;%给定在这些频率分点上理想的幅频响应;weigh=1 10;%给定在这些频率分点上的加权;6.0p7.0sb=remez(32,f,A,weigh);%设计出切比雪夫最佳一致逼近滤波器;%h,w=freqz(b,1,256,1);h=abs(h);h=20*log10(h);figure(1
33、)stem(b,.);grid;figure(2)plot(w,h);grid;例例7.6 利用切比雪夫最佳一致逼近法设计一个多阻带陷波器,数字系统的抽样频率为500Hz,去掉工频信号(50Hz)及二次、三次谐波的干扰。解:解:clear all;%用切比雪夫最佳一致逼近设计线性相位多带FIR滤波器;f=0.14.18.22.26.34.38.42.46.54.58.62.66 1;A=1 1 0 0 1 1 0 0 1 1 0 0 1 1;weigh=8 1 8 1 8 1 8;b=remez(64,f,A,weigh);%h,w=freqz(b,1,256,1);hr=abs(h);h=abs(h);h=20*log10(h);figure(1)stem(b,.);grid;figure(2)plot(w,h);grid;