窗函数法设计FIR数字低通滤波器w

1、基于窗函数法的FIR数字低通滤波器设计基于窗函数法的FIR数字低通滤波器设计摘 要数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。与IIR滤波器相比，FIR的实现是非递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。滤波器的设计是信号处理的核心问题之一。根据FIR滤波器的原理，提出了FIR滤波器的窗函数设计法，

2、给出了在MATLAB环境下，用窗函数法设计FIR滤波器的过程和设计实例。通过利用不同的窗函数方法设计FIR滤波器，对所设计的滤波器进行分析比较，得出各种方法设计的滤波器的优缺点及其不同的使用场合，从而可以在设计滤波器时能够正确的选择FIR数字滤波器的窗函数的选取及设计方法。关键词：FIR滤波器，MATLAB，窗函数I目 录摘 要I1 概述11.1 FIR滤波器简介11.2 窗函数设计法12 设计原理32.1 基本原理32.2 典型的窗函数43 几种数字低通滤波器的窗函数设计73.1 采用矩形窗设计FIR数字低通滤波器73.2 采用汉明窗设计FIR数字低通滤波器73.3 采用布莱克曼窗设计FIR

3、数字低通滤波器9参考文献10附 录121 概述1.1 FIR滤波器简介 FIR数字滤波器设计最简单的方法是窗函数法，通常也称为傅立叶级数法。它是在时域进行的，因而必须由理想滤波器的频率响应推导出其单位冲激响应，在设计一个FIR数字滤波器的单位冲激响应去逼近。根据冲激响应的时域特性，数字滤波器可分为无限长冲激响应（IIR）和有限长冲激响应滤波器（FIR），FIR的突出优点是：系统总是稳定的、易于实现线性相位、允许设计多通带（或多阻带）滤波器，但与IIR相比，在满足同样阻带衰减的情况下需要的阶数较高，滤波器的阶数越高，占用的运算时间越多，因此在满足指标要求的情况下应尽量减少滤波器的阶数。1.2 窗

4、函数设计法窗函数设计法是一种通过截短和计权的方法使无限长非因果序列成为有限长脉冲应响应序列的设计方法，通常在设计滤波器之前，应该先根据具体的工程应用确定滤波器的技术指标，在大多数实际应用中，数字滤波器常常被用来实现选频操作，所以指标的形式一般为在频域中以分贝值给出的相对幅度响应和相位响应。 用窗函数法设计FIR滤波器的步骤如下： （1）根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过

5、渡带宽为，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数W（n）。 （2）根据待求滤波器的理想频率响应求出理想单位脉冲响应hd（n），如果给出待求滤波器频率应为Hd（ejw），则理想的单位脉冲响应可以用下面的傅里叶反变换式求出： （1）在一般情况下，hd（n）是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样N点，采用离散傅里叶反变换（IDFT）即可求出。 （3）计算滤波器的单位脉冲响

6、应h（n），它是理想单位脉冲响应和窗函数的乘积。 (4)算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。使用窗函数法设计时要满足以下两个条件：窗谱主瓣尽可能地窄，以获得较陡的过渡带；尽量减少窗谱的最大旁瓣的相对幅度，也就是使能量尽量集中于主瓣，减小峰肩和纹波，进行增加阻带的衰减。窗函数的选择原则是：1具有较低的旁瓣幅度，尤其是第一旁瓣的幅度；2旁瓣的幅度下降的速率要快，以利于增加阻带的衰减；3主瓣的宽度要窄，这样可以得到比较窄的过渡带。通常上述的几点难以同时满足。实际中设计FIR数字滤波器

7、往往要求是线性相位的，因此要求满足线性相位的条件，即要求w(n)满足： （2）所以，窗函数不仅有截短的作用，而且能够起到平滑的作用，在很多领域得到了应用。2 设计原理2.1 基本原理设计低通FIR数字滤波器，寻求一系统函数，使其频率响应逼近滤波器要求的理想频率响应，其对应的单位脉冲响应 （3）其中 如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为 （4）窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到： （5）就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为 （6）式中，N

8、为所选窗函数的长度。2.2 典型的窗函数（1）矩形窗(Rectangle Window) （7） 其频率响应和幅度响应分别为：， （8）（2）三角形窗(Bartlett Window) （9） 其频率响应为： （3）汉宁(Hanning)窗，又称升余弦窗 （10）其频率响应和幅度响应分别为：（4）汉明(Hamming)窗，又称改进的升余弦窗 （11）其幅度响应为： （5）布莱克曼(Blankman)窗，又称二阶升余弦窗 （12）其幅度响应 ： （6）凯塞(Kaiser)窗 （13） 其中：是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，越大，过渡带越宽，阻带越小衰减也越大。

9、I0(·)是第一类修正零阶贝塞尔函数。 若阻带最小衰减表示为，的确定可采用下述经验公式： 若滤波器通带和阻带波纹相等即p=s时，滤波器节数可通过下式确定： 式中：我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见下表1。表2-1 各种窗函数的基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134/N-12汉明窗-418/N-53布莱克曼窗-5712/N-74这样选定窗函数类型和长度N之后，求出单位脉冲响应，并

10、按照式（6）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足 (14)根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。我们在设计滤波器时，希望窗谱主瓣尽可能窄，以获得较陡的过渡带，同时尽可能减小最大旁瓣的相对幅度，我们发现这是不可能的，因为这本身就是一个矛盾体，所以在设计滤波器时只是根据实际

11、情况来选择合适的窗函数。3 几种数字低通滤波器的窗函数设计3.1 采用矩形窗设计FIR数字低通滤波器图3-1 采用矩形窗设计FIR数字低通滤波器的仿真图3.2 采用汉明窗设计FIR数字低通滤波器图3-2 采用汉明窗设计FIR数字低通滤波器的仿真图3.3 采用布莱克曼窗设计FIR数字低通滤波器图3-3 采用布莱克曼窗设计FIR数字低通滤波器的仿真图图3-4采用布莱克曼窗冲激响应仿真图图3-5采用布莱克曼窗对含噪语音信号及滤波后的语音信号仿真图参考文献1 从玉良.数字信号处理原理及其MATLAB实现M.北京：电子工业出版社.2009.72 胡广书.数字信号处理理论、算法与实现M.北京：清华大学出版

12、社.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. Hist

13、ory of OCLCEB/OL2000-01-08 http: /www. oclc. org/about/ history/default. htm. 附 录矩形窗设计程序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);gr

14、id 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);plot(rad,20*log10(abs(mag2);g

15、rid 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(m

16、ag2);grid on; h2,w2=freqz(h2,1,100,2);subplot(2,1,2);plot(w2,unwrap(angle(h2);grid on;用布莱克曼语音信号进行滤波，并得滤波前后信号的时域波形及频谱程序clear;clc;x,fs=wavread('D:SHE.wav');sound(x,fs);if x/2=0 x=x' 0;endfnoise =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);endx=x+noise'fp=1000;fst=1200;delta_w=2*pi*(fst-fp)/fs;wc=pi*(fst+fp)/fs;N=31;if N/2=0N=N+1;endtau=(N-1)/2;for n=1:Nh(n)=sin(wc*(n-tau)/(pi*(n-tau);endh(N-1)/2)=(fst+fp)