用窗函数法设计FIR数字低通滤波器要点

用窗函数法设计FIR数字低通滤波器要点

河北科技大学

课程设计报告

学生姓名： 学 号： 专业班级： 课程名称： 学年学期

指导教师：

20 年 月

课程设计成绩评定表

学生姓名 专业班级 设计题目 指 导 教 师 评 语 指导教师： 年 月 日 学 号 起止时间 成绩

目 录

1. 窗函数设计低通滤波器

1.1设计目的……………………………………………………………………1 1.2设计原理推导与计算………………………………………………………1 1.3设计内容与要求……………………………………………………………2 1.4设计源程序与运行结果……………………………………………………3 1.5思考题…………………………………………………………………… 10 1.6心得体会………………………………………………………………… 14

参考文献……………………………………………………………………… 15

1.窗函数设计低通滤波器

1.1设计目的

1. 熟悉设计线性相位数字滤波器的一般步骤。

2. 掌握用窗函数法设计FIR数字滤波器的原理和方法。

3. 熟悉各种窗函数的作用以及各种窗函数对滤波器特性的影响。 4. 学会根据指标要求选择合适的窗函数。

1.2设计原理推导与计算

如果所希望的滤波器的理想的频率响应函数为Hdej，则其对应的单位脉冲响应为

hdn12Hdejejd （4.1）

窗函数设计法的基本原理是设计设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数Hej，即

Hde1 hdn2jj，cN1e,其中

20，cHdej1ed2jcejejdcsincna na用有限长单位脉冲响应序列hn逼近hdn。由于hdn往往是无限长序列，而且是非因果的，所以用窗函数n将hdn截断，并进行加权处理，得到： hnhdnn （4.2）

hn就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响

应函数Hej为

Hejhnen0N1jn （4.3）

式中，N为所选窗函数n的长度。

用窗函数法设计的滤波器性能取决于窗函数n的类型及窗口长度N的取

- 1 -

值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表（一）。

旁瓣峰值幅度/dB 过渡带宽 阻带最小衰减/dB 过渡带带宽窗函数 2N 0.9 矩形窗 -13 4π/N 8π/N 8π/N 8π/N 12π/N 10π/N -12 三角形窗 -25 -25 2.1 汉宁窗 -31 -44 3.1 哈明窗 -41 -53 3.3 不莱克曼窗 凯塞窗(=7.865) -57 -74 5.5 -57 -80 5 表（一） 各种窗函数的基本参数

这样选定窗函数类型和长度N之后，求出单位脉冲响应hnhdn•n，并按照式（4.3）求出Hej。Hej是否满足要求，如果Hej不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。

1.3设计内容与要求

（一）设计要求：

1. 学会计算滤波器各项性能指标及如何来满足给定的指标要求。 2. 用MATLAB语言编程实现给定指标要求的滤波器的设计。

- 2 -

3. 熟悉MATLAB语言，编写程序。 4. 设计低通FIR滤波器的指标： 通带最大波动

Rp0.25dB,，

p0.2

阻带最小衰减 As50dB,（二）、设计内容：

，s0.3

1.熟悉各种窗函数,在MATLAB命令窗下浏览各种窗函数，绘出（或打印）所看到的窗函数图。

2.编写计算理想低通滤波器单位抽样响应hd(n)的m函数文件ideal.m。 3. 编写计算N阶差分方程所描述系统频响函数H(ej)的m函数文件fr.m。 4.根据指标要求选择窗函数的形状与长度N。（至少选择两种符合要求的窗函数及其对应的长度）。

H(ej)，打印在频率区间[O，π]上的幅频响应特性曲线H(ej)~，

幅度用分贝表示。

6.验证所设计的滤波器是否满足指标要求。

1.4设计的源程序及运行结果：

1、利用MATLAB窗口观察各种窗函数： %巴特利特窗 w=bartlett(20); subplot(3,2,1); plot(w);

stem(w,'y');%'y'表示黄色 %stem表示以离散图输出 title('巴特利特床窗'); xlabel('n');%横坐标为n ylabel('w(n)');%纵坐标为w(n)

%布莱克曼窗 w=blackman(20);

- 3 -

subplot(3,2,2); plot(w);

stem(w,'b');%'b'表示蓝色 title('布莱克曼窗'); xlabel('n'); ylabel('w(n)'); %矩形窗 w=boxcar(20); subplot(3,2,3); plot(w); stem(w,'r'); title('矩形窗');

xlabel('n'); ylabel('w(n)'); %海明窗 w=hamming(20); plot(w);

stem(w,'m');%'m'表示紫色 title('海明窗'); xlabel('n'); ylabel('w(n)'); %汉宁窗 w=hanning(20); subplot(3,2,5); plot(w);

stem(w,'g');%'g'表示绿色

巴特利特床窗11title('汉宁窗'); xlabel('n'); ylabel('w(n)'); %凯泽窗 beta=5.6533; w=kaiser(20,beta); subplot(3,2,6); plot(w);

stem(w,'k');%'k'表示黑色 title('凯泽窗,beta=5.6533');

xlabel('n'); ylabel('w(n)');

布莱克曼窗w(n)0.50w(n)0510n矩形窗15200.500510n海明窗152011w(n)0.500510n汉宁窗1520w(n)0.5001015n凯泽窗,beta=5.653352011w(n)0.50w(n)0510n15200.500510n1520 常用窗函数的图形

2、理想低通滤波器单位抽样响应hd(n)的m函数文件ideal.m。

- 4 -

function hd=ideal(wc,M) %理想低通滤波器计算

%hd为0到M-1之间的理想脉冲响应 %wc为截止频率 %M为理想滤波器的长度 alpha=(M-1)/2; n=0:M-1; m=n-alpha+eps; hd=sin(wc*m)./(pi*m);

3、N阶差分方程所描述的系统频响函数的m函数文件fr.m。function[db,mag,pha,gfd,w]=fr(b,a) %求解系统响应 %db为相位振幅（db) %mag为绝对振幅 %pha为相位响应 %grd为群延时 %w为频率样本点矢量

%b为Ha(z)分析多项式系数(对FIR而言，b=h) %a为Hz(z)分母多项式系数(对FIR而言，a=1) [H,w]=freqz(b,a,1000,'whole'); H=(H(1:501))'; w=(w(1:501))'; mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=angle(H); gfd=grpdelay(b,a,w); 4、实际低通滤波器FIR： %用海明窗设计低通滤波器 wp=0.2*pi;

- 5 -

ws=0.3*pi; tr_width=ws-wp;

disp(['海明窗设计低通滤波器参数:']); M=ceil(6.6*pi/tr_width)+1;

disp(['滤波器的长度为',num2str(M)]); n=0:M-1;

wc=(ws+wp)/2; %理想LPF的截止频率 hd=ideal(wc,M); w_ham=(hamming(M))'; h=hd.*w_ham;

[db,mag,pha,gfd,w]=fr(h,[1]); delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1))); %求出实际通带波动 disp(['实际带通波动为',num2str(Rp)]);

As=-round(max(db(ws/delta_w+1:1:501))); %求出最小阻带衰减 disp(['最小阻带衰减为-',num2str(As)],’db’); %绘图

subplot(1,1,1) subplot(2,6,1) stem(n,hd);

title('理想冲击响应'); axis([0 M-1 -0.1 0.3]); ylabel('hd(n)'); subplot(2,6,2) stem(n,w_ham); title('海明窗'); axis([0 M-1 0 1.1]); ylabel('w(n)'); subplot(2,6,7)

- 6 -

stem(n,h);

title('实际冲激响应'); axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)'); subplot(2,6,8) plot(w/pi,db);

title('幅度响应（db）'); axis([0 1 -100 10]); grid;

xlabel('以pi为单位的频率'); ylabel('分贝数');

图（1）海明窗设计的FIR

海明窗设计低通滤波器参数: 滤波器的长度为67 实际带通波动为0.03936

- 7 -

最小阻带衰减为-52db

%用布莱克曼窗设计低通滤波器

wp=0.2*pi; ws=0.3*pi; tr_width=ws-wp;

disp(['布莱克曼窗设计低通滤波器的参数:']); M=ceil(11.0*pi/tr_width)+1;

disp(['滤波器的长度为',num2str(M)]); n=0:M-1;

%理想LPF的截止频率 wc=(ws+wp)/2; hd=ideal(wc,M); w_bla=(blackman(M))'; h=hd.*w_bla;

[db,mag,pha,gfd,w]=fr(h,[1]); delta_w=2*pi/1000;

Rp=-(min(db(1:1:wp/delta_w+1))); %求出实际通带波动 disp(['实际带通波动为',num2str(Rp)]);

As=-round(max(db(ws/delta_w+1:1:501))); %求出最小阻带衰减 disp(['最小阻带衰减-',num2str(As)],’db’); %绘图

subplot(2,6,3) stem(n,hd);

title('理想冲击响应'); axis([0 M-1 -0.1 0.3]); ylabel('hd(n)'); subplot(2,6,4) stem(n,w_bla); title('布莱克曼窗');

- 8 -

axis([0 M-1 0 1.1]); ylabel('w(n)'); subplot(2,6,9) stem(n,h);

title('实际冲激响应'); axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)'); subplot(2,6,10) plot(w/pi,db);

title('幅度响应（db）'); axis([0 1 -100 10]); grid;

xlabel('以pi为单位的频率'); ylabel('分贝数');

图（2）布莱克曼窗设计的FIR

- 9 -

布莱克曼窗设计低通滤波器的参数: 滤波器的长度为111 实际带通波动为0.0033304 最小阻带衰减为-73db 5、技术指标比较：

（1）海明窗设计低通滤波器参数: 滤波器的长度为67 实际带通波动为0.03936 最小阻带衰减为-52db

（2）布莱克曼窗设计低通滤波器的参数: 滤波器的长度为111 实际带通波动为0.0033304 最小阻带衰减为-73db

在相同的技术指标下用布莱克曼窗设计的低通滤波器实际带通波动实际带通波动最小，最小阻带衰减，滤波器的长度最大；海明窗和凯泽窗最小阻带衰减差不多，滤波器的长度页差不多，但是海明窗实际波动小于凯泽窗；所以用布莱克曼窗用设计的FIR最逼近理想单位冲击响应。这三个窗设计的低通滤波器都符合要求。

1.5思考题：

1. 设计线性相位数字滤波器的一般步骤：

（1）、给定所要求的频率响应函数Hd(ejw)以及技术指标阻带衰减和; （2）、求hd=IDTFT[Hd(ejw)];

由过渡带带宽及带阻最小衰减的要求，利用六种常见的窗函数基本参数的比较表或凯泽窗设计FIR滤波器的经验公式，选择窗函数的形状及N的大小（一般N要通过几次试探而最后确定）；

（3）求得所设计的FIR滤波器的单位抽样响应：

h(n)=hd(n)w(n),N=0,1,2,3…,N-1

（4）、求Hd(ejw)=DIFT[h(n)],校验是否满足设计要求，如果不满足，则重

- 10 -

新设计。

2. 窗函数有哪些指标要求?对给定指标要求的低通滤波器，理论计算所需窗函数的长度N。

答:窗函数的指标要求：主瓣宽度，旁瓣峰值。海明窗设计低通滤波器的

长度为67，凯泽窗设计低通滤波器的长度为60，布莱克曼窗设计低通滤波器的长度为111。

3. 用窗函数法设计FIR滤波器，滤波器的过渡带宽度和阻带衰减与哪些因素有关？

答：过渡带宽度与窗函数的形状和窗的宽度有关；阻带衰减只有窗函数的

形状决定，不受N的影响。

4、计算理想带通滤波器单位抽样响应hd(n)时取N为奇数和N为偶数有没有区别？你计算时所用的方法是仅适合于N为奇数或偶数的一种还是两种都可以用？

答： 以海明窗为例说明：

理想冲击响应N=670.30.2hd(n)0.10-0.10102030理想冲击响应N=684050600.30.2hd(n)0.10-0.10102030理想冲击响应N=0.30.2405060hd(n)0.10-0.10102030405060 滤波器的长度为 67 实际带通波动为 0.03936

- 11 -

最小阻带衰减为 52 滤波器的长度为 68 实际带通波动为 0.0324 最小阻带衰减为 53 滤波器的长度为 实际带通波动为 0.068677 最小阻带衰减为 45

N取奇偶数时，低通滤波器的幅度函数是不同的，如上图所示（海明窗设计低通滤波器N取奇偶数时的图），通过比较，当滤波器的长度大于技术指标要求的长度时，选择偶数也是满足要求的，所以海明窗奇偶都满足。同理，N取奇偶数，选择其他的窗函数也满足设计低通滤波器的指标要求。

比较所选窗长N相同但窗形状不同对滤波器设计结果的影响以及选同一种窗函数但窗长N不同时对滤波器设计结果的影响，将结论写在报告中。

海明窗 布莱克曼窗 滤波器的长度为67 实际带通波动0.29523 最小阻带衰减为30 凯泽窗 滤波器的长度为67 实带通波动0.047919 最小阻带衰减为51 线性斜率beta=4.5513 以凯泽窗设计的低通滤波器的长度为60 滤波器的长度N=60，实际带通波动为0.12166 用布莱克曼窗和海明最小阻带衰减为39 窗窗设计低通滤波器 滤波器的长度为60 实际带通波动为0.47079 最小阻带衰减为26 滤波器的长度为60 实际带通波动0.053711 最小阻带衰减为51 线性斜率beta=4.5513 以海明窗设计的低通滤波器的长度为67 滤波器的长度N=67，实际带通波动0.03936 用布莱克曼窗和凯泽最小阻带衰减为52 窗设计低通滤波器 - 12 -

以布莱克曼窗设计的滤波器的长度为111 低通滤波器的长度实际带通波动为0.024 滤波器的长度为111 实际带通波动0.0033304 最小阻带衰减为73 滤波器的长度为111 实际带通波动0.034426 最小阻带衰减为58 线性斜率beta=4.5513 N=111，用凯泽窗和海最小阻带衰减为59 明窗窗设计低通滤波器 结论：

1、当以海明窗设计的低通滤波器的长度N，用布莱克曼窗和凯泽窗设计低通滤波器，即N=67时，三个窗函数设计的低通滤波器中布莱克曼窗设计的不满足要求，凯泽窗设计的带通波动大于海明窗设计的低通滤波器带通波动，最小阻带衰减小于海明窗设计的低通滤波器。凯泽窗设计的结果没有用海明窗设计结果好。

当N取67 、60、111时，用海明窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<67时不满足设计技术指标。

2、以凯泽窗设计的低通滤波器的长度N，用布莱克曼窗和海明窗窗设计低通滤波器，即N=60时，海明窗设计的低通滤波器的最小阻带衰减不满足指标要求，布莱克曼窗设计的低通滤波器的最小阻带衰减和带通波动都不满足要。

当N取67 、60、111时，用布莱克曼窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N<111时不满足设计技术指标。

3、以布莱克曼窗设计的低通滤波器的长度N，用凯泽窗和海明窗窗设计低通滤波器，即N=111时，三个都满足要求。但是用布莱克曼窗函数设计的低通滤波器的性能更好。

当N取67 、60、111时，用凯泽窗设计的低通滤波器的性能不同，通过上表比较，N也越大性能越好，满足要求，当N>60时都满足设计技术指标。

- 13 -

1.6心得体会：

通过运用MATLAB设计低通滤波器，学习了MATLAB 语言的编程以及熟悉MATLAB软件环境，会利用MATLAB编写基本程序，熟悉绘制图形命令的正确使用，以及掌握MATLAB处理数字信号的有关命令。在命令窗口运行编写好的m文件，学习M文件调用其他的M文件，文件之间的调用必须满足这两个文件在一个文件夹里。

学习生成常用窗函数的MATLAB语言函数，并通过调用生成各种窗函数MATLAB内部函数，以及编写窗函数M文件，在figure窗口观察常用窗函数的图形。

- 14 -

参考文献

1、王华、李有军、刘建存，《MATLAB 电子仿真与应用教程》，国防科技大学出版社 ，2007

2、闻新、周露、张鸿《MATLAB 科学图形构建基础与应用》，科技出版社，2002 3、程佩青，数字信号处理及应用，清华大学出版社，2009

4、王宏，MATLAB及其在信号处理中的应用，清华大学出版社，2004

- 15 -