窗函数法设计滤波器

1、实验六 用窗函数法设计fir数字滤波器一. 实验目的(1) 掌握用窗函数法设计fir数字滤波器的原理与方法。(2) 熟悉线性相位fir数字滤波器的特性。(3) 了解各种窗函数对滤波特性的影响。二. 实验内容和要求(1) 复习用窗函数法设计fir数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2) 用升余弦窗设计一线性相位低通fir数字滤波器，截止频率。窗口长度n =15，33。要求在两种窗口长度情况下，分别求出，打印出相应的幅频特性和相频特性曲线，观察3db带宽和20db带宽。总结窗口长度n 对滤波器特性的影响。设计低通fir数字滤波器时，一般以理想低通滤波特性为逼近函数，即其中(3) ，

2、用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性曲线，观察3db带宽和20db带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。三. 实验方法、步骤及结果测试如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为 （4.1）窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到：（4.2）就作为实际设计的fir数字滤波器的单位脉冲响应序列，其频率响应函数为（4.3）式中，n为所选窗函数的长度。我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度n的取值。设计过程中，要根据对阻带最小

3、衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度n 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。表1 各种窗函数的基本参数窗函数旁瓣峰值幅度/db过渡带宽阻带最小衰减/db矩形窗-134/n-12三角形窗-258/n-25汉宁窗-318/n-44哈明窗-418/n-53不莱克曼窗-5712/n-74凯塞窗(=7.865)-5710/n-80这样选定窗函数类型和长度n之后，求出单位脉冲响应，并按照式（4.3）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用fft计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函数类型和长度n

4、 ，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足根据上式中的正、负号和长度n的奇偶性又将线性相位fir滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。主程序框图如图4.1所示。其中幅度特性要求用db表示。开始读入窗口长度n计算hd(n)调用窗函数子程序求w(n)调用子程序（函数）计算h(k)=dfth(n)调用绘图子程序（函数）绘制h(k)幅度相位曲线结束图6-1 主程序框图计算h(n)= hd(n) w(n)设画图时，用打印幅度特性。第k点对应的频率。为使曲线包络更接近的幅度特性曲线，dft变换区间要选大些

5、。例如窗口长度n=33时，可通过在末尾补零的方法，使长度变为64，再进行64点dft，则可以得到更精确的幅度衰减特性曲线。下面给出matlab主程序：%实验四，用窗函数法设计fir数字滤波器b=1;close all;i=0;while(b); temp=menu(选择窗函数长度n,n=10,n=15,n=20,n=25,n=30,n=33,n=35,n=40,n=45,n=50,n=55,n=60,n=64); menu1=10,15,20,25,30,33,35,40,45,50,55,60,64; n=menu1(temp); temp=menu(选择逼近理想低通滤波器截止频率wc,wc

6、=pi/4,wc=pi/2,wc=3*pi/4,wc=pi,wc=0.5,wc=1.0,wc=1.5,wc=2.0,wc=2.5,wc=3.0); menu2=pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3; w=menu2(temp); n=0:(n-1); hd=ideal(w,n); %得到理想低通滤波器 k=menu(请选择窗口类型:,boxcar,hamming,hanning,blackman); if k=1 b=boxcar(n); string=boxcar,n=,num2str(n); else if k=2 b=hamming(n); stri

7、ng=hamming,n=,num2str(n); else if k=3 b=hanning(n); string=hanning,n=,num2str(n); else if k=4 b=blackman(n); string=blackman,n=,num2str(n); end end end end h=hd.*(b); %得到fir数字滤波器 h,m=freqz(h,1,1024,whole); %求其频率响应 mag=abs(h); %得到幅值 db=20*log10(mag+eps)/max(mag); pha=angle(h); %得到相位 i=i+1; figure(i)

8、subplot(2,2,1); n=0:n-1; stem(n,h,.); axis(0,n-1,-0.1,0.3); hold on; n=0:n-1; x=zeros(n); plot(n,x,-); xlabel(n); ylabel(h(n); title(实际低通滤波器的h(n); text(0.3*n),0.27,string); hold off; subplot(2,2,2); plot(m/pi,db); axis(0,1,-100,0); xlabel(w/pi); ylabel(db); title(衰减特性(db); grid; subplot(2,2,3); plot

9、(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,-); title(相频特性); xlabel(频率(rad); ylabel(相位(rad); axis(0,3.15,-4,4); subplot(2,2,4); plot(m,mag); title(频率特性); xlabel(频率w(rad); ylabel(幅值); axis(0,3.15,0,1.5); text(0.9,1.2,string); b=menu(do you want to continue ?,yes,no); if b=2 b=0; endendtemp=menu(clo

10、se all figure ?,yes,no);if temp=1 close allend程序运行结果：运行程序，根据实验内容要求和程序提示选择你要进行的实验参数。三个实验参数选定后，程序运行输出用所选窗函数设计的实际fir低通数字滤波器的单位脉冲响应h(n)、幅频衰减特性(20lg|h(ejw)|)、相频特性及幅频特性|h(ejw)|的波形，h(n)和|h(ejw)|图中标出了所选窗函数类型及其长度n值。对四种窗函数（n=15和n=33）的程序运行结果如图4-2到图4-9所示，由图可以看出用各种窗函数设计的fir滤波器的阻带最小衰减及过渡带均与教材中一致。在通带内均为严格相位特性。图6-2

11、 矩形窗（n=15）图6-3 矩形窗（n=33）图6-4 哈明窗（n=15）图6-5 哈明窗（n=33）图6-6 汉宁窗（n=15）图6-7 汉宁窗（n=33）图6-8 布莱克曼窗（n=15）图6-9 布莱克曼窗（n=33）四. 实验报告要求(1) 简述实验原理及目的。(2) 按照实验步骤以及要求，比较各种情况下的滤波性能，说明窗口长度n和窗函数类型对滤波特性的影响。(3) 总结用窗函数法设计fir滤波器的主要特点。(4) 简要回答思考题。五. 思考题(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器，写出设计步骤。(2) 如果要求用窗函数法设计带

12、通滤波器，而且给定上、下边带截止频率为和，试求理想带通的单位脉冲响应27%实验中的子函数：产生理想低通滤波器单位脉冲响应hd(n)function hd=ideal(w,n);alpha=(n-1)/2;n=0:(n-1);m=n-alpha+eps;hd=sin(w*m)./(pi*m);相关练习：例1 用凯塞窗设计一fir低通滤波器，低通边界频率 ，阻带边界频率 ，阻带衰减 不小于50db。 解 首先由过渡带宽和阻带衰减 来决定凯塞窗的n和 图4.1给出了以上设计的频率特性，(a) 为n=30直接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的matlab程序为： wn=kai

13、ser(30,4.55); nn=0:1:29; alfa=(30-1)/2; hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa); h=hd.*wn; h1,w1=freqz(h,1); plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10); grid;xlabel(归一化频率/p) ylabel(幅度/db) 例2 利用雷米兹交替算法，设计一个线性相位低通fir数字滤波器，其指标为：通带边界频率fc=800hz，阻带边界fr=1000hz，通带波动 阻带最小衰减at=40db，采样频率fs=4000hz。 解 在matlab中

14、可以用remezord 和remez两个函数设计，其结果如图4.2，matlab程序如下： fedge=800 1000; mval=1 0; dev=0.0559 0.01; fs=4000; n,fpts,mag,wt=remezord(fedge,mval,dev,fs); b=remez(n,fpts,mag,wt); h,w=freqz(b,1,256); plot(w*2000/pi,20*log10(abs(h); grid; xlabel(频率/hz) ylabel(幅度/db) 函数remezord中的数组fedge为通带和阻带边界频率，数组mval是两个边界处的幅值，而数组

15、dev是通带和阻带的波动，fs是采样频率单位为hz。7.7 设计fir数字滤波器的窗函数法窗函数法的基本思想是用一具有有限长度样值响应、并具有线性相位的系统函数逼近理想滤波器的系统函数。以低通滤波器为例，设系统函数在奈奎斯特区间为（7.7-1）频率响应如图7.7-1所示。对取idtft，有（7.7-2） 图7.7-1理想低通滤波器的频率响应式（7.7-2）给出的单位样值响应说明理想低通滤波器是非因果系统，实际中无法实现。如果把对称地截断为有限长度，设长度为n，为奇数，截取后的样值响应为，当足够大时，在某种程度上可近似表示，为（7.7-3）其中。式（7.7-3）仍然为偶对称函数，而且是非因果的，

16、如果将其右移点，则有（7.7-4） 样值响应的截断与移位如图7.7-2所示。(a)(b)图7.7-2截断并移位的样值响应以上通过对理想滤波器的无限长样值响应截断为有限长并移位得到了物理可实现滤波器的样值响应，那么的幅频响应是否为理想幅频响应的逼近，其相位响应是否与频率成线性关系，下面对其进行分析。根据dtft变换公式，、和的频率响应分别为（7.7-5）（7.7-6）（7.7-7）式（7.7-6）相当于对做级数截断，当足够大时，逼近于。由式（7.7-4）和傅里叶变换的移位性质可得（7.7-8）由于为实偶对称函数，根据傅里叶变换的对称性，的变换也为实偶对称函数，因此，有（7.7-9）（7.7-10

17、） 可见，的幅频响应逼近于理想滤波器的幅频响应，而相频响应为的分段线性函数，在改变符号的频率处跳变。通过以上分析看出，样值响应为偶对称实序列、且为奇数的这种fir滤波器不仅在幅度上可以满足滤波要求，而且能够取得分段线性相位特性。当为奇对称的实序列时，根据傅里叶变换的对称性，的实部为零，这时（7.7-9）（7.7-10） 则的相位也为频率的分段线性函数。当的长度为偶数时，若仍满足对称关系（7.7-11a）或（7.7-11b）也可证明的相位为频率的分段线性函数。设取，为21和时，由式（7.7-4），分别为 把代入式（7.7-7），用计算机求得的幅频响应如图7.7-3所示。由图可见：在处幅度总为0.

18、5；越大，过渡带越窄，幅频响应越接近于理想响应；通带与阻带波纹在接近过渡带处较大，在两侧总有比较大的肩峰出现，这是因级数的截断造成的，这种现象称为傅里叶级数的吉伯斯（gibbs）现象。理论分析表明，两肩峰的幅度近似为8.9，与的取值无关，因此，滤波器的阻带衰减为。(a) n = 21 (b) n = 51图7.7-3 为21和时幅频响应上述fir滤波器的样值响应还可表示为（7.7-12）其中为长度为的矩形窗函数，即（7.7-13）由图7.7-3，矩形窗的逼近性能不是很好，为了加大阻带衰减也即减小肩峰的影响，可以考虑采用其他形状的窗函数，常用的有1.三角窗（bartlett窗）（7.7-14） 2.汉宁（hanning）窗（升余弦窗）（7.7-15） 3.哈明（hamming）窗（改进升余弦窗） （7.7-16）4. 布莱克曼（blackman）窗 （7.7-17）用以上几种窗函数设计的滤波器的指标如表7.7-1所示。设计时一般根据阻带衰减指标选取符合要求的窗函数，再根据过渡带宽度确定样值响应长度。当仿真结果不符合要求时，调整重新设计。 表7.7-1 常用