2.2 FIR数字滤波器的设计
FIR滤波器设计的任务是选择有限长度的h(n),使传输函数H(ejw)满足一定的幅度特性和线性相位要求。由于FIR滤波器很容易实现严格的线性相位,所以FIR数字滤波器设计的核心思想是求出有限的脉冲响应来逼近给定的频率响应。
设计过程一般包括以下三个基本问题[7]: (1)根据实际要求确定数字滤波器性能指标;
(2)用一个因果稳定的系统函数去逼近这个理想性能指标; (3)用一个有限精度的运算去实现这个传输函数。
设计方法一般分为窗函数法、频率采样法、等波纹逼近法,本文主要以窗函数法举例。 2.2.1利用窗函数法设计FIR数字滤波器
设计FIR数字滤波器的最简单的方法是窗函数法,通常也称之为傅立叶级数法。这种方法一般是先给定所要求的理想滤波器的频率响应Hd(ejw),要求设计一个FIR数字滤波器频率响应H(ejw)=
N?1?h(n)en?0?jwn去逼近理想的频率响应Hd(ejw)。然而,窗函数法设计FIR数字滤波器是
在时域中进行的,因此必须首先由理想的频率响应Hd(ejw)的傅里叶反变换推导出对应的单位脉冲响应hd(n)
hd(n)=
12?????Hd(ejw)ejwndw
(2-2)
再设计一个FIR数字滤波器的单位取样响应h(n)去逼近hd(n)。
由于Hd(ejw)是矩形频率特性,故hd(n)一定是无限长的序列,且是非因果的。而要设计的是FIR数字滤波器,其h(n)必定是有限长的,所以要用有限长的h(n)来逼近无限长的hd(n),最简单且有效的方法就是截断hd(n)。
设计过程如下:
Hd(e) ????hd(n)????h(n)????H(e)
j?IDTFFT*w(n)DTFTj?加窗的作用是通过把理想滤波器的无限长脉冲响应hd(n)乘以窗函数w(n)来产生一个被截断的脉冲响应,即h(n)?hd(n)w(n)并且对频率响应进行平滑[15]。
窗函数主要用来减少序列因截断而产生的Gibbs效应。但当这个窗函数为矩形时,得到的FIR滤波器幅频响应会有明显的Gibbs效应,并且任意增加窗函数的长度(即FIR滤波器的抽头数)Gibbs效应也不能得到改善。为了克服这种现象,窗函数应该使设计的滤波器:
(1)频率特性的主瓣宽度应尽量窄,以获取较陡的过渡带;
(2)尽量减少频率特性的最大旁瓣的相对幅度。也就是能量尽量集中于主瓣,增大阻带的衰减。
6
但这两项要求是不能同时满足的。当选用主瓣宽度较窄时,虽然得到较陡的过渡带,但是通带和阻带的波动明显增加;当选用最小的旁瓣幅度时,虽然能得到平坦的幅度响应和较小的阻带波纹,但过渡带加宽,即主瓣会加宽。因此,实际所选用的窗函数往往是它们的折中。在保证主瓣宽度达到一定要求的前提下,适当的牺牲主瓣宽度以换取相对旁瓣的抑制[6]。
设计FIR数字滤波器常用的窗函数有: 1.矩形窗(Rectangle window)
w(n)?RN0?n?N?1?1,???0,其他
(2-3)
其频率响应为
WR(ejw)?WR(w)e?j(N?1)w2
(2-4) (2-5)
WR(w)?WR(ejwsin(wN/2)sin(w/2)
)主瓣宽度为4?/N。
2.三角形窗(Bartlett window)
N?1?2n,0?n???N?12w(n)??N?1?2-2n,?n?N?1?N?12?(2-6)
其频率响应为
jwW(e)?2N?1sin[({N?1)]wN?1?j()w2sin(Nw/4)2?j(2)w242? }e()esin(w/2)Nsin(w/2)N?1(2-7)
近似结果在N>>1时成立。此时,主瓣宽度为8?瓣却小很多。
3.汉宁窗(Hanning window) 汉宁窗又称升余弦窗。
w(n)?sin(2/N,比矩形窗主瓣宽度增加一倍,但旁
?nN?1)RN(n)?12[1?cos(2?nN?1)]RN(n)
(2-8)
其频率响应为
W(ejw)?{0.5WR(w)?0.25[WR(w?2?N?1)?WR(w?2?N?1)]}e?j(N?12)w?W(w)e?j(N?12)w
(2-9)
7
当N>>1时,N-1?N,所以,窗函数的幅度函数为
W(w)?0.5WR(w)?0.25[WR(w?2?N)?WR(w?2?N)] (2-10)
这三部分之和,使旁瓣互相抵消,能量更集中在主瓣。但是代价是主瓣宽度比矩形窗的主瓣宽度增加一倍,即为8?/N。
4.海明窗(Hamming window) 海明窗又称改进的升余弦窗。
把升余弦窗加以改进,可以得到旁瓣更小的效果,形式为
w(n)?[0.54?0.46cos(2?nN?1)]RN(n)
(2-11)
其频率响应为
W(w)?0.54WR(w)?0.23[WR(w??0.54WR(w)?0.23[WR(w?2?N2?N?1)?WR(w?2?N)]2?N?1)] (2-12)
)?WR(w?与汉宁窗相比,主瓣宽度相同,为8?量集中在窗谱的主瓣内。
5.布莱克曼窗(Blackman window) 布莱克曼窗又称二阶升余弦窗。
/N,但旁瓣又进一步压低,结果可将99.963%的能
为了进一步抑制旁瓣,对升余弦窗再加上一个二次谐波的余弦分量,变成布莱克曼窗,故又称二阶升余弦窗。
w(n)?[0.42?0.5cos(2?nN?1)?0.08cos(4?nN?1)]RN(n)
(2-13)
其频率响应为
W(w)?0.42WR(w)?0.25[WR(w??0.04[WR(w?4?N?1)?WR(w?2?N?1)]/N)?WR(w?2?N?1)]4?N?1
(2-14)
此时主瓣宽度为矩形窗谱主瓣宽度的3倍,即为12?6.凯塞窗(Kaiser window)
。
以上几种窗函数是各以一定主瓣加宽为代价,来换取某种程度的旁瓣抑制,而凯塞窗则是全面地反映主瓣与旁瓣衰竭之间的交换关系,可以在它们两者之间自由的选择它们的比重。
凯塞窗是一种适应性较强的窗,其窗函数的表达式为
w(n)?I0(?1?[1?2n/(N?1)])I0(?)8
2,0?n?N-1
(2-15)
式中,I0(x)是第一类变形零阶贝塞尔函数,?是一个可自由选择的参数,它可以同时调整主瓣宽度与旁瓣电平,?越大,则w(n)窗越窄,而频谱的旁瓣越小,但主瓣宽度也相应增加。因而改变?值就可对主瓣宽度与旁瓣衰减进行选择。一般选择4<9,这相当于旁边幅度与主瓣幅度的比值由3.1%变到0.047%(-30dB—-67dB)[9]。 2.2.2窗函数法设计线性相位FIR数字滤波器的一般步骤
窗函数法的设计步骤归纳如下:
(1)给定希望逼近的频率响应函数Hd(ejw); (2)根据式(2-2)求单位脉冲响应hd(n)
hd(n)=
12???H??d(ejw)ejwndw
如果Hd(ejw)很复杂或不能直接计算积分,则必须用求和代替积分,以便在计算机上计算,也就是要计算离散傅里叶反函数。
(3)由过渡带宽及阻带最小衰减的要求,可选定窗形状,并估计窗口长度。因过渡带?w近似与窗口长度成反比,N的窗函数。
(4)计算所设计的FIR数字滤波器的单位脉冲响应。
h(n)?hd(n)w(n)?A/?w,A决定于窗口形式,A参数选择参考表1。按照过渡带及
阻带衰减情况,选择窗函数形式。原则是在保证阻带衰减满足要求的情况下,尽量选择主瓣窄
,0?n?N?1
(5)由h(n)求FIR数字滤波器的系统函数H(z)
N?1H(z)??n?0h(n)z?n
表1 6种窗函数基本参数的比较
过度带宽?w 8?/N 8?/N 8?/N 12?/N 10?/N
4?/N
阻带最小衰减/dB
-21 -25 -44 -53 -74 -80
窗函数
矩形窗 三角形窗 汉宁窗 海明窗 布莱克曼窗
凯塞窗 (β=7.865)
旁瓣峰值幅度/dB
-13 -25 -31 -41 -57 -57
根据窗函数法设计一个FIR数字低通滤波器,技术指标如下:
9
通带截止频率ωp=0.2πrad,通带允许波动Rp=0.25dB; 阻带截止频率ωs=0.3πrad,阻带衰减Rs=50dB。
思路:查表1可知,海明窗、布莱克曼窗以及凯塞窗均可提供大于50dB的衰减,但是海明窗具有较小的过渡带从而具有较小的长度N,故采用海明窗设计该低通滤波器
基于MATLAB对FIR数字低通滤波器设计实现的编程如下: clc; clear;
Wp=0.2*pi;Rp=0.25;Ws=0.3*pi;Rs=50; wt=Ws-Wp;%计算过渡带宽
N=8*pi/wt;%按海明窗计算滤波器的长度 n=[0:1:N-1];
Wc=(Wp+Ws)/2;%计算3dB通带截止频率 hd=ideal_lp(Wc,N); w_han=(hamming(N))'; h=hd.*w_han;
[db,mag,pha,grd,w]=freqz_m(h,1); subplot(1,1,1); subplot(2,2,1);
plot(w/pi,db);title('幅频特性图'); xlabel('w/pi'); ylabel('dB');
axis([0, 1, -150, 50]); subplot(2,2,2);
plot(w/pi,pha);title('相频特性图'); xlabel('w/pi'); ylabel('pha'); subplot(2,2,3);
stem(n,h,'.');title('实际脉冲响应'); xlabel('n');ylabel('h(n)'); subplot(2,2,4);
stem(n,w_han,'.');title('海明窗'); xlabel('n');ylabel('海明窗');
基于MATLAB对FIR数字低通滤波器的仿真结果如下图2.1所示:
10