数字信号||用窗函数法设计FIR数字滤波器（5）

实验五  用窗函数法设计FIR数字滤波器

一、实验目的

(1)加深对窗函数法设计FIR数字滤波器的基本原理的理解。
(2)学习用MATLAB语言的窗函数法编写设计FIR数字滤波器的程序。
(3)了解MATLAB有关窗函数法设计的常用子函数。

二、实验涉及的MATLAB子函数

1.boxcar
功能：矩形窗。
格式：w＝boxcar(n)
说明：boxcar(n)函数可产生一长度为n的矩形窗函数。

2.triang

功能：三角窗。
格式：w＝triang(n)
说明：triang(n)函数可得到n点的三角窗函数。三角窗系数为：当n为奇数

当n为偶数时，

1. bartlett
   功能：bartlett(巴特利特)窗。
   格式：w＝Bartlett(n)
   说明：bartlett(n)可得到n点的bartlett窗函数。bartlett窗函数系数为
2. hamming
   功能：hamming(哈明)窗。
   格式：w＝hamming(n)
   说明：hamming(n)可产生n点的hamming窗。hamming窗函数系数为

1. hanning
   功能：hanning(汉宁)窗。
   格式：w＝hanning(n)
   说明：hanning(n)可产生n点的hanning窗。hanning窗函数系数为

1. blackman
   功能：blackman(布莱克曼)窗。
   格式：w＝blackman(n)
   说明：blackman(n)可产生n点的blackman窗。blackman窗函数系数为
   
   与等长度的hamming和hanning窗相比，blackman窗的
2. chebwin
   功能：chebyshev(切比雪夫)窗。
   格式：w＝chebwin(n，r)
   说明：chebwin(n，r)可产生n点的chebyshev窗函数，其傅里叶变换后的旁瓣波纹低于主瓣r dB。注意：当n为偶数时，窗函数的长度为n＋1。

8.kaiser
功能：kaiser(凯塞)窗。
格式：w＝kaiser(n，beta)
说明：kaiser(n，beta)可产生n点的kaiser窗函数，其中，beta为影响窗函数旁瓣的b参数，其最小的旁瓣抑制a与b之间的关系为

9.fir1
功能：基于窗函数的FIR数字滤波器设计——标准频率响应，以经典方法实现加窗线性相位FIR滤波器设计，可设计出标准的低通、带通、高通和带阻滤波器。
格式：※b＝fir1(n，Wn)；设计截止频率为Wn的hamming(哈明)加窗线性相位滤波器，滤波器系数包含在b中。当0≤Wn≤1(Wn＝1相应于0.5fs)时，可得到n阶低通FIR滤波器。
　　当Wn＝［W1W2］时，fir1函数可得到带通滤波器，※b＝fir1(n，Wn，¢ftype¢)；可设计高通和带阻滤波器，由ftype决定：
·当ftype＝high时，设计高通FIR滤波器；
·当ftype＝stop时，设计带阻FIR滤波器。
　　在设计高通和带阻滤波器时，fir1函数总是使用偶对称N为奇数(即第一类线性相位FIR滤波器)的结构，因此当输入的阶次为偶数时，fir1函数会自动加1。

※b＝fir1(n，Wn，Window)；利用列矢量Window中指定的窗函数进行滤波器设计，Window长度为n＋1。如果不指定Window参数，则fir1函数采用hamming窗。
※b＝fir1(n，Wn，¢ftype¢，Window)；可利用ftype和Window参数设计各种加窗的滤波器。
　　由fir1函数设计的FIR滤波器的群延迟为n/2。

三、实验原理

1.运用窗函数法设计FIR数字滤波器
FIR数字滤波器的系统函数为

这个公式也可以看成是离散LSI系统的系统函数：

分母a0为1，其余ak全都为0时的一个特例。由于极点全部集中在0点，稳定和线性相位特性是FIR滤波器突出的优点，因此在实际中具有更广泛的使用价值。
FIR滤波器的设计任务是选择有限长度的h(n)，使传输函数H(ejw)满足技术要求。主要设计方法有窗函数法、频率采样法和切比雪夫等波纹逼近法等。本实验主要介绍用窗函数法设计FIR数字滤波器。用窗函数法设计FIR数字滤波器的基本步骤是：
(1)根据过渡带和阻带衰减设计指标选择窗函数的类型，估算滤波器的阶数N。
(2)由数字滤波器的理想频率响应H(ejw)求出其单位冲激响应hd(n)。
对于理想的数字低通滤波器频率响应，有下列的子程序可以实现(该程序名为ideal－lp.m)：
functionhd＝ideal－lp(wc，N)
%hd＝点0到N－1之间的理想脉冲响应
%wc＝截止频率(弧度)

%N＝理想滤波器的长度
tao＝(N－1)/2；
n＝［0：(N－1)］；
m＝n－tao＋eps；  %加一个小数以避免0作除数
hd＝sin(wc*m)./(pi*m)；
　　其它选频滤波器则可以由低通频响特性合成。如一个通带在wc1～wc2之间的带通滤波器，在给定N值的条件下，可以用下列程序实现：
hd＝ideal－lp(wc2，N)－ideal－lp(wc1，N)；

(3)计算数字滤波器的单位冲激响应h(n)＝w(n)hd(n)。
(4)检查设计出的滤波器是否满足技术指标。
　　如果不满足技术指标，则需要重新选择或调整窗函数的类型，估算滤波器的阶数N。再重复前面的四个步骤，直到满足指标为止。
　　常用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、切比雪夫窗、布莱克曼窗、凯塞窗等，MATLAB均有相应的子函数可以调用。
　　另外，MATLAB信号处理工具箱还提供了fir1子函数，可以用于窗函数法设计FIR滤波器。

　由于第一类线性相位滤波器(类型Ⅰ)能进行低通、高通、带通、带阻滤波器的设计，因此，本实验所有滤波器均采用第一类线性相位滤波器。

2.各种窗函数特性的比较
例24-1 在同一图形坐标上显示矩形窗、三角形窗、汉宁窗、哈明窗、布莱克曼窗、凯塞窗的特性曲线。
　　解 程序如下：
N＝64；beta＝7.865；n＝1：N；   %输入N、凯塞窗需要的b值
wbo＝boxcar(N)；%矩形窗
wtr＝triang(N)；%三角形窗
whn＝hanning(N)；%汉宁窗

whm＝hamming(N)；%哈明窗
wbl＝blackman(N)；%布莱克曼窗
wka＝kaiser(N，beta)；%凯塞窗
plot(n¢，［wbo，wtr，whn，whm，wbl，wka］)； %在同一界面上作图
axis(［0，N，0，1.1］)；
legend(¢矩形¢，¢三角形¢，¢汉宁¢，¢哈明¢，¢布莱克曼¢，¢凯塞¢) %线型标注

N=64;beta=7.865;n=1:N;   %输入N、凯塞窗需要的b值

wbo=boxcar(N);%矩形窗

wtr=triang(N);%三角形窗

whn=hanning(N);%汉宁窗

whm=hamming(N);%哈明窗

wbl=blackman(N);%布莱克曼窗

wka=kaiser(N,beta);%凯塞窗

plot(n',[wbo,wtr,whn,whm,wbl,wka]); %在同一界面上作图

axis([0,N,0,1.1]);

legend('矩形','三角形','汉宁','哈明','布莱克曼','凯塞') %线型标注

　　程序结果如图24-1所示，MATLAB将自动用不同颜色标出各条曲线。因黑白印刷无法分辨，故改用不同线型表示。

图24-1 常用窗函数形状比较

　　为了便于滤波器设计，表24-1给出了六种窗函数的特性参数。

四、实验任务

(1)阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。
(2)认真阅读实验原理，明确本次实验任务，读懂例题程序，了解实验方法。
(3)列写调试通过的实验程序，打印或描绘实验程序产生的曲线图形。

3.用窗函数法设计FIR数字低通滤波器
例24-2 用矩形窗设计一个FIR数字低通滤波器，要求：N＝64，截止频率为wc＝0.4p，描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。
　　解 程序如下：
wc＝0.4*pi；          %输入设计指标
N＝64；n＝0：N－1；
hd＝ideal－lp(wc，N)；%建立理想低通滤波器
windows＝(boxcar(N))¢；%使用矩形窗，并将列向量变为行向量

b＝hd.*windows；%求FIR系统函数系数
　　［H，w］＝freqz(b，1)；%求解频率特性
dbH＝20*log10((abs(H)＋eps)/max(abs(H)))； %化为分贝值
%作图
subplot(2，2，1)，stem(n，hd)；
axis(［0，N，1.1*min(hd)，1.1*max(hd)］)；title(¢理想脉冲响应¢)；
xlabel(¢n¢)；ylabel(¢hd(n)¢)；
subplot(2，2，2)，stem(n，windows)；
axis(［0，N，0，1.1］)；title(¢窗函数特性¢)；

xlabel(¢n¢)；ylabel(¢wd(n)¢)；
subplot(2，2，3)，stem(n，b)；
axis(［0，N，1.1*min(b)，1.1*max(b)］)；title(¢实际脉冲响应¢)；
xlabel(¢n¢)；ylabel(¢h(n)¢)；

subplot(2，2，4)，plot(w/pi，dbH)；
axis(［0，1，－80，10］)；title(¢幅度频率响应¢)；
xlabel(¢频率(单位：＼pi)¢)；ylabel(¢H(e^{j＼omega})¢)；
set(gca，¢XTickMode¢，¢manual¢，¢XTick¢，［0，wc/pi，1］)；
set(gca，¢YTickMode¢，¢manual¢，¢YTick¢，［－50，－20，－3，0］)；grid

wc=0.4*pi;          %输入设计指标

N=64;n=0:N-1;

hd=ideal-lp(wc, N);%建立理想低通滤波器

windows=(boxcar(N))';%使用矩形窗，并将列向量为行向量

b=hd.*windows;%求FIR系统函数系数

[H,w]=freqz(b,1);%求解频率特性

dbH=20*log10((abs(H)+eps)/max(abs(H))); %化为分贝值%作图

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

axis([0, N, 1.1*min(hd), 1.1*max(hd)]);

title('理想脉冲响应');

xlabel('n');ylabel('hd(n)');

subplot(2,2,2),stem(n, windows);

axis([0,N,0,1.1]);title('窗函数特性');

xlabel('n');ylabel('wd(n)');

subplot(2,2,3),stem(n,b);

axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');

xlabel('n');ylabel('h(n)');

subplot(2,2,4),plot(w/pi,dbH);

axis([0,1,-80,10]);title('幅度频率响应');

xlabel('频率(单位：＼pi)');ylabel('H(e^{j\omega})');

set(gca,'XTickMode','manual','XTick',[0,wc/pi,1]);

set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid

程序运行结果如图24-2所示。

图24-2 例24-2设计的数字低通滤波器特性