基于窗函数法的FIR 数字低通滤波器设计
摘 要
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR ）滤波器和有限冲激响应（FIR ）滤波器。与IIR 滤波器相比，FIR 的实现是非递归的，总是稳定的；更重要的是，FIR 滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。
滤波器的设计是信号处理的核心问题之一。根据FIR 滤波器的原理，提出了FIR 滤波器的窗函数设计法，给出了在MATLAB 环境下，用窗函数法设计FIR 滤波器的过程和设计实例。通过利用不同的窗函数方法设计FIR 滤波器，对所设计的滤波器进行分析比较，得出各种方法设计的滤波器的优缺点及其不同的使用场合，从而可以在设计滤波器时能够正确的选择FIR 数字滤波器的窗函数的选取及设计方法。
关键词：FIR 滤波器，MATLAB ，窗函数
目 录
摘 要 ................................................................................................................ I 1 概述 ............................................................................................................... 1 1.1 FIR滤波器简介 . ....................................................................................... 1 1.2 窗函数设计法 . ......................................................................................... 1 2 设计原理 ....................................................................................................... 3 2.1 基本原理 .................................................................................................. 3 2.2 典型的窗函数 . ......................................................................................... 4 3 几种数字低通滤波器的窗函数设计 . .......................................................... 7 3.1 采用矩形窗设计FIR 数字低通滤波器 ................................................. 7 3.2 采用汉明窗设计FIR 数字低通滤波器 ................................................. 7 3.3 采用布莱克曼窗设计FIR 数字低通滤波器 ......................................... 9 参考文献 ......................................................................................................... 10 附 录 ............................................................................................................. 12
1 概述
1.1 FIR滤波器简介
FIR 数字滤波器设计最简单的方法是窗函数法，通常也称为傅立叶级数法。它是在时域进行的，因而必须由理想滤波器的频率响应H d (e jw ) 推导出其单位冲激响应
h d (n ) ，在设计一个FIR 数字滤波器的单位冲激响应h (n ) 去逼近h d (n ) 。根据冲激响
应的时域特性，数字滤波器可分为无限长冲激响应（IIR ）和有限长冲激响应滤波器（FIR ），FIR 的突出优点是：系统总是稳定的、易于实现线性相位、允许设计多通带（或多阻带）滤波器，但与IIR 相比，在满足同样阻带衰减的情况下需要的阶数较高，滤波器的阶数越高，占用的运算时间越多，因此在满足指标要求的情况下应尽量减少滤波器的阶数。
1.2 窗函数设计法
窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲应响应序列的设计方法，通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标，在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。 用窗函数法设计FIR 滤波器的步骤如下：
（1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N （或阶数M=N-1），窗函数类型可根据最小阻带衰减As 独立选择，因为窗口长度N 对最小阻带衰减As 没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N ，设待求滤波器的过渡带宽为Δω，它与窗口长度N 近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N ，在N 和窗函数类型确定后，即可调用MATLAB 中的窗函数求出窗函数W （n ）。
（2）根据待求滤波器的理想频率响应求出理想单位脉冲响应h d （n ），如果给出
待求滤波器频率应为H d （e jw ），则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：
h d (n ) =
12π
⎰π
-
π
H d (e jw ) e jw dw （1）
在一般情况下，h d （n ）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2π采样N 点，采用离散傅里叶反变换（IDFT ）即可求出。
（3）计算滤波器的单位脉冲响应h （n ），它是理想单位脉冲响应和窗函数的乘积。
(4)算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz 子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。
使用窗函数法设计时要满足以下两个条件：
1窗谱主瓣尽可能地窄，以获得较陡的过渡带； ○
2尽量减少窗谱的最大旁瓣的相对幅度，也就是使能量尽量集中于主瓣，减小峰肩○
和纹波，进行增加阻带的衰减。 窗函数的选择原则是：
[1]具有较低的旁瓣幅度，尤其是第一旁瓣的幅度； [2]旁瓣的幅度下降的速率要快，以利于增加阻带的衰减； [3]主瓣的宽度要窄，这样可以得到比较窄的过渡带。
通常上述的几点难以同时满足。实际中设计FIR 数字滤波器往往要求是线性相位的，因此要求w (n ) 满足线性相位的条件，即要求w(n)满足：
w (n ) =w (N -1-n ) （2）
所以，窗函数不仅有截短的作用，而且能够起到平滑的作用，在很多领域得到了应用。
2 设计原理
2.1 基本原理
设计低通FIR 数字滤波器，寻求一系统函数H (z ) ，使其频率响应H (e j ω) 逼近滤波器要求的理想频率响应H d (e j ω) ，其对应的单位脉冲响应h d (n )
H d e j ω
N -1
2
()
-j ωα⎧≤ωc ⎪e =⎨⎪⎩0，ωc
其中α=
1
h d (n )=
2π
⎰π
-
π
H d e
()
j ω
1
e d ω=
2π
j ω
⎰ω
-
ωc
e -j ωαe j ωαd ω=
c
sin [ωc (n -a )]
πn -a
如果所希望的滤波器的理想的频率响应函数为H d (e j ω)，则其对应的单位脉冲响应为
1
h d (n )=
2π
ωω
⎰πH (e )e d ω （4） π
j
j
-
d
窗函数设计法的基本原理是用有限长单位脉冲响应序列h (n )逼近h d (n )。由于
h d (n )往往是无限长序列，而且是非因果的，所以用窗函数ω(n )将h d (n )截断，并进行
加权处理，得到：
h (n )=h d (n )ω(n )
（5）
h (n )就作为实际设计的FIR 数字滤波器的单位脉冲响应序列，其频率响应函数
H e j ω为
()
H e
()=∑h (n )e
j ω
n =0
N -1
j ωn
（6）
式中，N 为所选窗函数ω(n )的长度。
2.2 典型的窗函数
（1）矩形窗(Rectangle Window)
w (n ) =R N (n ) （7）
其频率响应和幅度响应分别为：
sin(N ω/2) -j ω
W (e ) =e
sin(ω/2)
j ω
N -12
，W R (ω) =
sin(N ω/2)
（8）
sin(ω/2)
（2）三角形窗(Bartlett Window)
⎧2n ⎪N -1, w (n ) =⎨
2n ⎪2-,
N -1⎩
N -1
2 （9）
N -1
j ω
N -12
2sin(N ω/4) 2-j ω
]e 其频率响应为：W (e ) =[
N sin(ω/2)
（3）汉宁(Hanning)窗，又称升余弦窗
12n π
w (n ) =[1-)]R N (n ) （10）
2N -1
其频率响应和幅度响应分别为：
-j (2π2π
W (e ) ={0. 5W R (ω) +0. 25[W R (ω-) +W R (ω+)]}e
N -1N -1
=W (ω) e -j ωa
j ω
N -1
) ω2
W (ω) =0. 5W R (ω) +0. 25[W R (ω-
2π2π
) +W R (ω+)]N -1N -1
（4）汉明(Hamming)窗，又称改进的升余弦窗
2n π
w (n ) =[0. 54-0. 46)]R N (n ) （11）
N -1
其幅度响应为：W (ω) =0. 54W R (ω) +0. 23[W R (ω-
（5）布莱克曼(Blankman)窗，又称二阶升余弦窗
2π2π
) +W R (ω+)] N -1N -1
2n π4n π
w (n ) =[0. 42-0. 5) +0. 08)]R N (n ) （12）
N -1N -1
其幅度响应 ：
W (ω) =0. 42W R (ω) +0. 25[W R (ω-
2π2π
) +W R (ω+)]N -1N -1
4π4π
+0. 04[W R (ω-) +W R (ω+)]
N -1N -1
（6）凯塞(Kaiser)窗
I 0(β-[1-2n /(N -1)]2)
w (n ) =, 0≤n ≤N -1 （13）
I 0(β)
其中：β是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，β越大，过渡带越宽，阻带越小衰减也越大。I 0(·) 是第一类修正零阶贝塞尔函数。
若阻带最小衰减表示为A s =-20log 10δs ，β的确定可采用下述经验公式：
⎧0
⎪
β=⎨0. 5842(A s -21) 0. 4+0. 07886(A s -21)
⎪0. 1102(A -8. 7)
s ⎩
A s ≤21
2150
若滤波器通带和阻带波纹相等即δp=δs 时，滤波器节数可通过下式确定：
N =
A s -7. 95
+1
14. 36∆F
式中：∆F =
∆ωωs -ωp
= 2π2π
我们知道，用窗函数法设计的滤波器性能取决于窗函数ω(n )的类型及窗口长度N 的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见下表1。
表2-1 各种窗函数的基本参数
这样选定窗函数类型和长度N 之后，求出单位脉冲响应h (n )=h d (n )ω(n )，并按照式（6）求出H (e j ω)。H (e j ω)是否满足要求，要进行演算。一般在h (n )尾部加零使长度满足2的整数次幂，以便用FFT 计算H (e j ω)。如果要观察细节，补零点数增多即可。如果H (e j ω)不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。
如果要求线性相位特性，则h (n )还必须满足
h (n )=±h (N -1-n ) (14)
根据上式中的正、负号和长度N 的奇偶性又将线性相位FIR 滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择h (n )=h (N -1-n )这一类，而不能选择h (n )=-h (N -1-n )这一类。
我们在设计滤波器时，希望窗谱主瓣尽可能窄，以获得较陡的过渡带，同时尽可能减小最大旁瓣的相对幅度，我们发现这是不可能的，因为这本身就是一个矛盾体，所以在设计滤波器时只是根据实际情况来选择合适的窗函数。
3 几种数字低通滤波器的窗函数设计
3.1 采用矩形窗设计FIR 数字低通滤波器
100
-10-20-30-400
0.5
1
1.522.5
矩形窗滤波器幅頻响应
3
3.5
0-10
-20-30-40
00.10.2
0.30.40.50.60.7矩形窗滤波器相頻响应
0.80.91
图3-1 采用矩形窗设计FIR 数字低通滤波器的仿真图
3.2 采用汉明窗设计FIR 数字低通滤波器
0-20
-40-60-80
00.51
1.522.5
汉明窗滤波器幅頻响应
33.5
0-20
-40-60-80
00.10.2
0.30.40.50.60.7汉明窗滤波器相频响应
0.80.91
图3-2 采用汉明窗设计FIR 数字低通滤波器的仿真图
3.3 采用布莱克曼窗设计FIR 数字低通滤波器
-50
-100
-150
00.5
11.522.5布莱克曼窗滤波器幅頻响应
33.5
0-20
-40-60-80
00.10.20.3
0.40.50.60.7布莱克曼窗相频响应
0.80.91
图3-3 采用布莱克曼窗设计FIR 数字低通滤波器的仿真图
9
0.015
0.01
0.005
冲激响应序列
0.20.150.1
0.050
510冲激响应傅里叶变换
15
图3-4采用布莱克曼窗冲激响应仿真图
1510x 10
4
50
2
4
6
810
含噪语音信号
12
14
16x 10
4
3002001000
2
4
6810滤波后语言信号
12
14
16x 10
4
图3-5采用布莱克曼窗对含噪语音信号及滤波后的语音信号仿真图
10
参考文献
[1] 从玉良. 数字信号处理原理及其MATLAB 实现[M].北京：电子工业出版社.2009.7
[2] 胡广书. 数字信号处理理论、算法与实现[M].北京：清华大学出版社.2003,8 [3] 万永革. 数字信号处理的MATLAB 实现[M].北京: 科学出版社.2007. [4] 薛山.MATLAB 基础教程[M].北京: 清华大学出版社.2011.
[5] 陈怀琛. 数字信号处理教程—MATLAB 释义与实现[M].北京: 电子工业出版社.2002
[6] 缪家鼎，徐文娟． 光电技术[M]． 杭州：浙江大学出版社，1994：112-156．
[7] 李晓东，张庆红，叶瑾琳．气候学研究的若干理论问题[J]．北京大学学报：
自然科学版，1999，35(1)：101-106．
[8] 郑开青．通讯系统模拟及软件[D]．北京：清华大学无线电系，1987． [9]
Online
Computer
Library
Center,
Inc.
History
of
OCLC[EB/OL]．[2000-01-08]． http: //www. oclc. org/about/ history/default. htm.
11
附 录
矩形窗设计程序 passrad=0.4*pi; w2=boxcar(16);
n=1:1:16;
hd=sin(passrad*(n-8))./(pi*(n-8)); hd(8)=passrad/pi; h2=hd.*rot90(w2);
title('designed by Hanning window'); [mag2,rad]=freqz(h2); subplot(2,1,1);
plot(rad,20*log10(abs(mag2))); grid on ;
[h2,w2]=freqz(h2,1,100,2); subplot(2,1,2);
plot(w2,unwrap(angle(h2))); grid on ;
汉明窗设计程序 passrad=0.4*pi; w2=hamming(32); n=1:1:32;
hd=sin(passrad*(n-16))./(pi*(n-16)); hd(16)=passrad/pi; h2=hd.*rot90(w2);
title('designed by Hanning window'); [mag2,rad]=freqz(h2); subplot(2,1,1);
12
plot(rad,20*log10(abs(mag2))); grid on ;
[h2,w2]=freqz(h2,1,100,2); subplot(2,1,2);
plot(w2,unwrap(angle(h2))); grid on ;
布莱克曼窗函数设计程序 passrad=0.4*pi; w2=blackman(64); n=1:1:64;
hd=sin(passrad*(n-32))./(pi*(n-32)); hd(32)=passrad/pi; h2=hd.*rot90(w2);
title('designed by Hanning window'); [mag2,rad]=freqz(h2); subplot(2,1,1);
plot(rad,20*log10(abs(mag2))); grid on ;
[h2,w2]=freqz(h2,1,100,2); subplot(2,1,2);
plot(w2,unwrap(angle(h2))); grid on ;
13
用布莱克曼语音信号进行滤波，并得滤波前后信号的时域波形及频谱程序 clear;clc;
[x,fs]=wavread('D:\SHE.wav'); sound(x,fs); if x/2==0 x=[x' 0]; end
fnoise =10000;T=length(x)/fs; t=T/length(x):T/length(x):T; for j=1:length(x)
noise(j)=sin(2*pi*fnoise*t(j)); end x=x+noise'; fp=1000; fst=1200;
delta_w=2*pi*(fst-fp)/fs; wc=pi*(fst+fp)/fs; N=31; if N/2==0 N=N+1; end tau=(N-1)/2; for n=1:N
h(n)=sin(wc*(n-tau))/(pi*(n-tau)); end
h((N-1)/2)=(fst+fp)/fs; wn=(blackman(N)); h1=h(n).*wn';
subplot(2,1,1);stem(h1); F=abs(fft(h1));
14
subplot(2,1,2); plot(F(1:(N-1)/2)); X=abs(fft(x)); figure;subplot(2,1,1); plot(X(1:(length(x)+1)/2)); y=conv(x,h1); Y=abs(fft(y)); subplot(2,1,2);
plot(Y(1:(length(x)+1)/2)); sound(x,fs); sound(y,fs);
15