第24章--用窗函数法设计FIR数字滤波器

1、1.1 市场与市场营销1.2 我国汽车市场的开展与现状复习思考题实验24 用窗函数法设计FIR数字滤波器一、实验目的(1)加深对窗函数法设计FIR数字滤波器的根本原理的理解。(2)学习用MATLAB语言的窗函数法编写设计FIR数字滤波器的程序。(3)了解MATLAB有关窗函数法设计的常用子函数。二、实验涉及的MATLAB子函数1.boxcar功能：矩形窗。格式：wboxcar(n)说明：boxcar(n)函数可产生一长度为n的矩形窗函数。2.triang功能：三角窗。格式：wtriang(n)说明：triang(n)函数可得到n点的三角窗函数。三角窗系数为：当n为奇数时，当n为偶数时，3.ba

2、rtlett功能：bartlett(巴特利特)窗。格式：wBartlett(n)说明：bartlett(n)可得到n点的bartlett窗函数。bartlett窗函数系数为4.hamming功能：hamming(哈明)窗。格式：whamming(n)说明：hamming(n)可产生n点的hamming窗。hamming窗函数系数为5.hanning功能：hanning(汉宁)窗。格式：whanning(n)说明：hanning(n)可产生n点的hanning窗。hanning窗函数系数为6.blackman功能：blackman(布莱克曼)窗。格式：wblackman(n)说明：blackma

3、n(n)可产生n点的blackman窗。blackman窗函数系数为与等长度的hamming和hanning窗相比，blackman窗的主瓣稍宽，旁瓣稍低。7.chebwin功能：chebyshev(切比雪夫)窗。格式：wchebwin(n，r)说明：chebwin(n，r)可产生n点的chebyshev窗函数，其傅里叶变换后的旁瓣波纹低于主瓣r dB。注意：当n为偶数时，窗函数的长度为n1。8.kaiser功能：kaiser(凯塞)窗。格式：wkaiser(n，beta)说明：kaiser(n，beta)可产生n点的kaiser窗函数，其中，beta为影响窗函数旁瓣的b参数，其最小的旁瓣抑制

4、a与b之间的关系为增加b可使主瓣变宽、旁瓣的幅度降低。9.fir1功能：基于窗函数的FIR数字滤波器设计标准频率响应，以经典方法实现加窗线性相位FIR滤波器设计，可设计出标准的低通、带通、高通和带阻滤波器。格式：bfir1(n，Wn)；设计截止频率为Wn的hamming(哈明)加窗线性相位滤波器，滤波器系数包含在b中。当0Wn1(Wn1相应于0.5fs)时，可得到n阶低通FIR滤波器。当WnW1W2时，fir1函数可得到带通滤波器，其通带为w1ww2。bfir1(n，Wn，ftype)；可设计高通和带阻滤波器，由ftype决定：当ftypehigh时，设计高通FIR滤波器；当ftypestop

5、时，设计带阻FIR滤波器。在设计高通和带阻滤波器时，fir1函数总是使用偶对称N为奇数(即第一类线性相位FIR滤波器)的结构，因此当输入的阶次为偶数时，fir1函数会自动加1。bfir1(n，Wn，Window)；利用列矢量Window中指定的窗函数进行滤波器设计，Window长度为n1。如果不指定Window参数，那么fir1函数采用hamming窗。bfir1(n，Wn，ftype，Window)；可利用ftype和Window参数设计各种加窗的滤波器。由fir1函数设计的FIR滤波器的群延迟为n/2。三、实验原理1.运用窗函数法设计FIR数字滤波器FIR数字滤波器的系统函数为这个公式也可

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

7、的子程序可以实现(该程序名为ideallp.m)：functionhdideallp(wc，N)%hd点0到N1之间的理想脉冲响应%wc截止频率(弧度)%N理想滤波器的长度tao(N1)/2；n0：(N1)；mntaoeps； %加一个小数以防止0作除数hdsin(wc*m)./(pi*m)；其它选频滤波器那么可以由低通频响特性合成。如一个通带在wc1wc2之间的带通滤波器，在给定N值的条件下，可以用以下程序实现：hdideallp(wc2，N)ideallp(wc1，N)；(3)计算数字滤波器的单位冲激响应h(n)w(n)hd(n)。(4)检查设计出的滤波器是否满足技术指标。如果不满足技术指

8、标，那么需要重新选择或调整窗函数的类型，估算滤波器的阶数N。再重复前面的四个步骤，直到满足指标为止。常用的窗函数有矩形窗、三角形窗、汉宁窗、哈明窗、切比雪夫窗、布莱克曼窗、凯塞窗等，MATLAB均有相应的子函数可以调用。另外，MATLAB信号处理工具箱还提供了fir1子函数，可以用于窗函数法设计FIR滤波器。由于第一类线性相位滤波器(类型)能进行低通、高通、带通、带阻滤波器的设计，因此，本实验所有滤波器均采用第一类线性相位滤波器。2.各种窗函数特性的比较例24-1 在同一图形坐标上显示矩形窗、三角形窗、汉宁窗、哈明窗、布莱克曼窗、凯塞窗的特性曲线。解 程序如下：N64；beta7.865；n1

9、：N； %输入N、凯塞窗需要的b值wboboxcar(N)；%矩形窗wtrtriang(N)；%三角形窗whnhanning(N)；%汉宁窗whmhamming(N)；%哈明窗wblblackman(N)；%布莱克曼窗wkakaiser(N，beta)；%凯塞窗plot(n，wbo，wtr，whn，whm，wbl，wka)； %在同一界面上作图axis(0，N，0，1.1)；legend(矩形，三角形，汉宁，哈明，布莱克曼，凯塞) %线型标注程序结果如图24-1所示，MATLAB将自动用不同颜色标出各条曲线。因黑白印刷无法分辨，故改用不同线型表示。图24-1 常用窗函数形状比较为了便于滤波器设

10、计，表24-1给出了六种窗函数的特性参数。表24-1 六种窗函数的特性参数表3.用窗函数法设计FIR数字低通滤波器例24-2 用矩形窗设计一个FIR数字低通滤波器，要求：N64，截止频率为wc0.4p，描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。解 程序如下：wc0.4*pi； %输入设计指标N64；n0：N1；hdideallp(wc，N)；%建立理想低通滤波器windows(boxcar(N)；%使用矩形窗，并将列向量变为行向量bhd.*windows；%求FIR系统函数系数H，wfreqz(b，1)；%求解频率特性dbH20*log10(abs(H)eps)/max(a

11、bs(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

12、)，plot(w/pi，dbH)；axis(0，1，80，10)；title(幅度频率响应)；xlabel(频率(单位：pi)；ylabel(H(ejomega)；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设计的数字低通滤波器特性例24-3 选择适宜的窗函数设计一个FIR数字低通滤波器，要求：通带截止频率为wp0.3p，Rp0.05 dB；阻带截止频率为ws0.45p，As50 dB。描绘该滤波器的脉冲响应、窗

13、函数及滤波器的幅频响应曲线和相频响应曲线。解 查表24-1，选择哈明窗。程序如下：wp0.3*pi；ws0.45*pi； %输入设计指标deltawwswp；%计算过渡带的宽度N0ceil(6.6*pi/deltaw)； %按表24-1所示哈明窗数据，求滤波器长度N0NN0mod(N01，2)%为实现FIR类型I偶对称滤波器，应确保N为奇数windows(hamming(N)；%使用哈明窗，并将列向量变为行向量wc(wswp)/2；%截止频率取通阻带频率的平均值hdideallp(wc，N)；%建立理想低通滤波器bhd.*windows；%求FIR系统函数系数db，mag，pha，grd，wf

14、reqzm(b，1)； %求解频率特性n0：N1；dw2*pi/1000； %dw为频率分辨率，将02p分为1000份Rp(min(db(1：wp/dw1)%检验通带波动Asround(max(db(ws/dw1：501)%检验最小阻带衰减作图局部省略。程序执行结果如下：N 45Rp0.0428As50程序运行结果如图24-3所示。图24-3 例24-3设计的数字低通滤波器特性由Rp、As数据和曲线可见，用哈明窗设计的结果能够满足设计指标要求。由N值可知，FIR数字低通滤波器的阶数一般都比较高，其公式不再用数学表达式列出。例24-4 用MATLAB信号处理箱提供的fir1子函数，设计一个FIR

15、数字低通滤波器，要求同例24-3：通带截止频率为wp0.3p，Rp0.05 dB；阻带截止频率为ws0.45p，As50 dB。解 查表24-1，选择哈明窗。程序如下：wp0.3*pi；ws0.45*pi； %输入设计指标deltawwswp；%计算过渡带的宽度N0ceil(6.6*pi/deltaw)；%按哈明窗计算滤波器长度N0NN0mod(N01，2)；%为实现FIR类型I偶对称滤波器，应确保N为奇数windowshamming(N)；%使用哈明窗，此句可省略TXXwc(wswp)/2/pi；%截止频率取归一化通阻带频率的平均值TXXbfir1(N1，wc，windows)%用fir1子

16、函数求系统函数系数，windows可省略TXXdb，mag，pha，grd，wfreqzm(b，1)； %求解频率特性n0：N1；dw2*pi/1000；%dw为频率分辨率，将02p分为1000份Rp(min(db(1：wp/dw1)%检验通带波动Asround(max(db(ws/dw1：501) %检验最小阻带衰减其中打横线的三条程序与例24-3不同，作图局部省略。程序运行结果如下：N 45Rp0.0428As50程序运行结果如图24-3所示。由Rp、As数据和曲线可见，用fir1子函数设计的结果与例24-3结果完全相同。4.用窗函数法设计FIR数字高通滤波器例24-5 选择适宜的窗函数设

17、计一个FIR数字高通滤波器，要求：通带截止频率为wp0.45p，Rp0.5 dB；阻带截止频率为ws0.3p，As20 dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解 查表24-1，选择三角形窗。程序如下：wp0.45*pi；ws0.3*pi； %输入设计指标deltawwpws；%计算过渡带的宽度N0ceil(6.1*pi/deltaw)；%按三角形窗计算滤波器长度N0NN0mod(N01，2)；%为实现FIR类型I偶对称滤波器，应确保N为奇数windows(triang(N)；%使用三角形窗，并将列向量变为行向量wc(wswp)/2；%截止频率取通阻带频率的平

18、均值hdideallp(pi，N)ideallp(wc，N)； %建立理想高通滤波器bhd.*windows；%求FIR系统函数系数db，mag，pha，grd，wfreqzm(b，1)； %求解频率特性n0：N1；dw2*pi/1000；%dw为频率分辨率，将02p分为1000份Rp(min(db(wp/dw1：501) %检验通带波动Asround(max(db(1：ws/dw1) %检验最小阻带衰减作图局部省略。程序运行结果如下：N 41Rp0.3625As25程序运行结果如图24-4所示。由Rp、As数据和曲线可见，用三角形窗设计的结果能够满足设计指标要求。图24-4 例24-5设计的

19、数字高通滤波器特性例24-6 用MATLAB信号处理箱提供的fir1子函数，设计一个FIR数字高通滤波器，要求：通带截止频率为fp450 Hz，Rp0.5 dB；阻带截止频率为fs300Hz，As20 dB；采样频率Fs2000 Hz。描绘滤波器的脉冲响应、窗函数及滤波器的幅频响应和相频响应曲线。解 查表24-1，选择三角形窗。程序如下：fs300；fp450；Fs2000； %输入设计指标wsfs/(Fs/2)*pi；wpfp/(Fs/2)*pi；%计算归一化角频率deltawwpws；%计算过渡带的宽度N0ceil(6.1*pi/deltaw)；%按三角形窗计算滤波器长度N0NN0mod(

20、N01，2)；%为实现FIR类型偶对称滤波器，应确保N为奇数windowstriang(N)；%使用三角形窗wc(wswp)/2/pi； %截止频率取归一化通阻带频率的平均值bfir1(N1，wc，high，windows)；%用fir1子函数求系统函数系数db，mag，pha，grd，wfreqzm(b，1)；%求解频率特性n0：N1；dw2*pi/1000； %dw为频率分辨率，将02p分为1000份Rp(min(db(wp/dw1：501)%检验通带波动Asround(max(db(1：ws/dw1)%检验最小阻带衰减作图局部省略。程序运行结果如下：N 41Rp0.3625As25图24

21、-5 例24-6设计的数字高通滤波器特性程序运行结果如图24-5所示。由Rp、As数据和曲线可见，用三角形窗设计的结果能够满足设计指标要求，且与例24-5结果相同，只是本例给出的是实际频率和采样频率指标，需先进行频率的归一化。作图程序也需进行相应的修改。如幅频响应和相频响应局部作图程序如下：subplot(2，2，3)，plot(w/2/pi*Fs，db)；axis(0，Fs/2，40，2)；title(幅频响应)；xlabel(f(Hz)；ylabel(H(ejomega)；set(gca，XTickMode，manual，XTick，0，fs，fp，Fs/2)；set(gca，YTickM

22、ode，manual，YTick，20，3，0)；gridsubplot(2，2，4)，plot(w/2/pi*Fs，pha)；axis(0，Fs/2，4，4)；title(相频响应)；xlabel(f(Hz)；ylabel(phi(omega)；set(gca，XTickMode，manual，XTick，0，fs，fp，Fs/2)；set(gca，YTickMode，manual，YTick，pi，0，pi)；grid5.用窗函数法设计FIR数字带通滤波器例24-7 选择适宜的窗函数设计一个FIR数字带通滤波器，要求：下阻带截止频率ws10.2p，As65 dB；通带低端截止频率wp10.

23、3p，Rp0.05 dB；通带高端截止频率wp20.7p，Rp0.05 dB；上阻带截止频率ws20.8p，As65 dB。描绘实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解 查表24-1，选择布莱克曼窗。程序如下：wp10.3*pi；wp20.7*pi； %输入设计指标ws10.2*pi；ws20.8*pi；wpwp1，wp2；wsws1，ws2；deltawwp1ws1；%计算过渡带的宽度N0ceil(11*pi/deltaw)；%按布莱克曼窗计算滤波器长度N0NN0mod(N01，2)；%为实现FIR类型偶对称滤波器，应确保N为奇数windows(blackman(

24、N)； %使用布莱克曼窗，并将列向量变为行向量wc1(ws1wp1)/2；wc2(ws2wp2)/2；%截止频率取通阻带频率的平均值hdideallp(wc2，N)ideallp(wc1，N)； %建立理想带通滤波器bhd.*windows；%求FIR系统函数系数db，mag，pha，grd，wfreqzm(b，1)；%求解频率特性n0：N1；dw2*pi/1000； %dw为频率分辨率，将02p分为1000份Rp(min(db(wp1/dw1：wp2/dw1) %检验通带波动ws01：ws1/dw1，ws2/dw1：501；%建立阻带频率样点数组Asround(max(db(ws0)%检验最

25、小阻带衰减subplot(2，2，1)，stem(n，b)；%作图axis(0，N，1.1*min(b)，1.1*max(b)；subplot(2，2，2)，stem(n，windows)；axis(0，N，0，1.1)；subplot(2，2，3)，plot(w/pi，db)；axis(0，1，150，10)；set(gca，XTickMode，manual，XTick，. 0，ws1/pi，wp1/pi，wp2/pi，ws2/pi，1)；set(gca，YTickMode，manual，YTick，. 100，65，20，3，0)；gridsubplot(2，2，4)，plot(w/pi，

26、pha)；axis(0，1，4，4)；set(gca，XTickMode，manual，XTick，. 0，ws1/pi，wp1/pi，wp2/pi，ws2/pi，1)；set(gca，YTickMode，manual，YTick，. pi，0，pi)；grid一条程序过长，需要转到下一行书写时，在中断行用“，.结尾，下一行继续书写程序，如上面程序中的set语句。程序运行结果如下：N 111Rp0.0033As73程序运行结果如图24-6所示。图24-6 例24-7设计的数字带通滤波器特性注意：带通滤波器的设计指标中，两边过渡带的宽度一般应该一致，deltaw可以由任意一边的带宽来确定。如果两

27、边过渡带的宽度不一致，那么应求两个过渡带中小的一个作为设计的依据(带阻滤波器的设计同)，此时deltaw一句改为deltawmin(wp1ws1)，(ws2wp2)；由Rp、As数据和曲线可见，用布莱克曼窗设计的结果不仅满足设计指标要求，且具有很小的通带波动(Rp0.0033 dB)，很高的阻带衰减值(最小阻带衰减As73 dB)，过渡带很窄。例24-8 用MATLAB信号处理箱提供的fir1子函数设计一个FIR数字带通滤波器，要求：下阻带截止频率fs1100 Hz，As65 dB；通带低端截止频率fp1150Hz，Rp0.05 dB；通带高端截止频率fp2350Hz，Rp0.05 dB；上阻

28、带截止频率fs2400Hz，As65 dB；采样频率Fs1000 Hz。描绘实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解 查表24-1，选择布莱克曼窗。程序如下：fp1150；fp2350； %输入设计指标fs1100；fs2400；Fs1000；ws1fs1/(Fs/2)*pi；ws2fs2/(Fs/2)*pi； %计算归一化角频率wp1fp1/(Fs/2)*pi；wp2fp2/(Fs/2)*pi；deltawwp1ws1；%计算过渡带的宽度N0ceil(11*pi/deltaw)；%按布莱克曼窗计算滤波器长度N0NN0mod(N01，2)%为实现FIR类型偶对称滤波

29、器，应确保N为奇数windowsblackman(N)；%使用布莱克曼窗wc1(ws1wp1)/2/pi；%截止频率取归一化通阻带频率的平均值wc2(ws2wp2)/2/pi；bfir1(N1，wc1，wc2，windows)；db，mag，pha，grd，wfreqzm(b，1)； %求解频率特性n0：N1；dw2*pi/1000； KG*2%dw为频率分辨率，将02p分为1000份Rp(min(db(wp1/dw1：wp2/dw1) %检验通带波动ws01：ws1/dw1，ws2/dw1：501；%建立阻带频率样点数组Asround(max(db(ws0)%检验最小阻带衰减作图局部省略。程

30、序运行结果如下：N 111Rp0.0033As73程序运行结果如图24-7所示。图24-7 例24-8用fir1设计的数字带通滤波器特性由Rp、As数据和曲线可见，用布莱克曼窗设计的结果完全能够满足设计指标要求。图24-6和图24-7的区别仅在于频率特性，一个使用了归一化的频率单位，一个使用了实际的频率单位。6.用窗函数法设计FIR数字带阻滤波器例24-9 选择适宜的窗函数设计一个FIR数字带阻滤波器，要求：下通带截止频率wp10.2p，Rp0.1 dB；阻带低端截止频率ws10.3p，As40 dB；阻带高端截止频率ws20.7p，As40 dB；上通带截止频率wp20.8p，Rp0.1 d

31、B。描绘实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解 查表24-1，选择汉宁窗。程序如下：wp10.2*pi；wp20.8*pi； %输入设计指标ws10.3*pi；ws20.7*pi；wpwp1，wp2；wsws1，ws2；deltawws1wp1；%计算过渡带的宽度N0ceil(6.2*pi/deltaw)；%按汉宁窗计算滤波器长度N0NN0mod(N01，2)%为实现FIR类型偶对称滤波器，应确保N为奇数windows(hanning(N)；%使用汉宁窗，并将列向量变为行向量wc1(ws1wp1)/2；wc2(ws2wp2)/2；%截止频率取通阻带频率的平均值hd

32、ideallp(wc1，N)ideallp(pi，N)ideallp(wc2，N)；%建立理想带阻bhd.*windows；%求FIR系统函数系数db，mag，pha，grd，wfreqzm(b，1)； %求解频率特性n0：N1；dw2*pi/1000；%dw为频率分辨率，将02p分为1000份wp01：wp1/dw1，wp2/dw1：501； %建立通带频率样点数组Rp(min(db(wp0) %检验通带波动Asround(max(db(ws1/dw1：ws2/dw1) %检验最小阻带衰减作图局部省略。程序运行结果如下：N 63Rp0.0888As44程序运行结果如图24-8所示。由Rp、A

33、s数据和曲线可见，用汉宁窗设计的结果完全能够满足设计指标要求。图24-8 例24-9的数字带阻滤波器特性例24-10 用凯塞窗设计一个长度为75的FIR数字带阻滤波器，要求：下通带截止频率wp10.2p，Rp0.1 dB；阻带低端截止频率ws10.3p，As60 dB；阻带高端截止频率ws20.7p，As60 dB；上通带截止频率wp20.8p，Rp0.1 dB。描绘实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。解 凯塞窗参数b0.112(As8.7)。用fir1子函数编写的程序如下：N75；As60； %输入设计指标wp10.2*pi；wp20.8*pi；ws10.3*pi；ws20.7*pi；beta0.1102*(As8.7)%计算b值windowskaiser(N，beta)；%使用凯塞窗wc1(ws1wp1)/2/p