1、FIR滤波器设计第7章 FIR滤波器设计之迟辟智美创作第六章我们介绍了无限冲激响应（IIR）滤波器的设计方法.其中最经常使用的由模拟滤波器转换为数字滤波器的方法为双线性变换法，因为这种方法无混叠效应，效果较好.但通过前面的例子我们看到，IIR数字滤波器相位特性欠好（非线性，如图 6-11、图6-13、图6-15 ），也不容易控制.然而在现代信号处置中，例如图像处置、数据传输、雷达接收以及一些要求较高的系统中对相位特性要求较为严格，这种滤波器就无能为力了.改善相位特性的方法是采纳有限冲激响应滤波器.本章首先对FIR滤波器原理及其使用函数作基本介绍，然后重点介绍窗函数法设计FIR滤波器，并对最优滤
2、波器设计函数进行介绍.7.1 FIR滤波器原理概述及滤波函数7.1.1 FIR滤波器原理及设计方法分类根据第 6 章对数字滤波器的介绍，我们知道FIR滤波器的传递函数为： (7-1)可得FIR滤波器的系统差分方程为：因此，FIR滤波器又称为卷积滤波器.根据第 4 章中所描述的系统频率响应，FIR滤波器的频率响应表达式为： （7-2）信号通过FIR滤波器不失真条件与（6-6）式所描述的相同，即滤波器在通带内具有恒定的幅频特性和线性相位特性.理论上可以证明（这里从略）：当FIR滤波器的系数满足下列中心对称条件： （7-3）时，滤波器设计在迫近平直幅频特性的同时，还能获得严格的线性相位特性.线性相位
3、FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的.对一个 N 阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长.这一特性使通带频率内信号通过滤波器后仍坚持原有波形形状而无相位失真.本章主要介绍的FIR数字滤波器设计方法及 MATLAB 信号处置工具箱提供的FIR数字滤波器设计函数，见表7-1.由于篇幅所限，本章我们主要介绍窗函数法和最优化设计方法.表7-1 FIR滤波器设计的主要方法函数设计方法说明工具函数窗函数法理想滤波器加窗处置fir1(单频带) ， fir2(多频带) ， kaiserord最优化设计平方误差最小化迫近理想幅频响应或Park-McCl
4、ellan 算法发生等波纹滤波器firls ， remez,remezord约束最小二乘迫近在满足最年夜误差限制条件下使整个频带平方误差最小化fircls,fircls1升余弦函数具有光滑、正弦过渡带的低通滤波器设计Fircos7.1.2 FIR数字滤波器滤波函数相对IIR 滤波器的滤波函数，FIR数字滤波器滤波函数除dimpulse和dstep仅适用于IIR滤波器外，其他各种函数可直接应用于FIR滤波器，只是输入的分母多项式向量a=1.另外，MATLAB还提供了一个函数fftfilt,该函数利用效率高的基于FFT算法实现对数据的滤波，该函数只适用于FIR滤波器，调用形式为：y=fftfilt
5、(b,x,n)式中，b为FIR滤波器的系数向量；x为输入数据；n为FFT长度，缺省时，函数选用最佳的FFT长度，y为滤波器的输出.该函数执行下面的把持：n=length(x);y=ifft(fft(x).*fft(b,n)./fft(a,n);应注意，y=fftfilt(b,x)等价于y=filter(b,a,x).7.2 FIR滤波器的窗函数设计7.2.1 窗函数的基来源根基理FIR滤波器设计的主要任务是根据给定的性能指标确定滤波器的系数b，即系统单元脉冲序列h(n)，它是一个有限长序列.FIR滤波器的理想频率响应，可写成复数形式的Fourier级数形式： (7-4)式中，hd(n)是对应的
6、单元脉冲响应序列.这说明滤波器的频率响应和单元脉冲响应互为Fourier变换对.因此其单元脉冲响应可由下式求得， (7-5)求得序列后，通过z变换，可获得 (7-6)注意，这里为无限长序列，因此是物理上不成实现的.如何酿成物理上可实现呢？一个自然的想法是只取其中的某些项，即只截取中的一部份，比如n=0,N-1，N为正整数.这种处置相当于将,n=-与函数w(n)相乘，w(n)具有下列形式：w(n)相当于一个矩形，我们称之为矩形窗.即我们可采纳矩形窗函数w(n)将无限脉冲响应截取一段h(n)来近似为，这种截取在数学上暗示为： h(n)= w(n) (7-7)这里应该强调的是，加窗函数不是可有可无的
7、，而是将设计酿成物理可实现所必需的.截取之后的滤波器传递函数酿成： (7-8)式中，N为窗口宽度，H(z)是物理可实现系统.为了获得线性相位，FIR滤波器h(n)必需满足中心对称条件（即7-3式），序列h(n)的延迟为.这种方法的基来源根基理是用一定宽度的矩形窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，从而获得FIR滤波器的脉冲响应，故称为FIR滤波器的窗函数设计法.经过加矩形窗后所得的滤波器实际频率响应能否很好地迫近理想频率响应呢？图 7-1 示意给出了理想滤波器加矩形窗后的情况.理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图.时间域内的乘积（7-7）式要求实际频
8、率响应为这两个频率响应函数在频域内的卷积（卷积定理），即获得图形为图7-1（右图）.图 7-1 FIR滤波器理想与实际频率响应由图可看出，加矩形窗后使实际频率响应偏离理想频率响应，主要影响有三个方面：（1）理想幅频特性陡直边缘处形成过渡带，过渡带宽取决于矩形窗函数频率响应的主瓣宽度.（2）过渡带两侧形成肩峰和波纹，这是矩形窗函数频率响应的旁瓣引起的，旁瓣相对值越年夜，旁瓣越多，波纹越多.（3）随窗函数宽度N的增年夜，矩形窗函数频率响应的主瓣宽度减小，但不改变旁瓣的相对值.为了改善FIR滤波器性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣相对值尽可能小，数量尽可能少，以获得通带波纹
9、小,阻带衰减年夜，在通带和阻带内均平稳的特点，这样可使滤波器实际频率响应更好地迫近理想频率响应.这里我们明确两个概念：截断和频谱泄漏.信号是无限长的，而在进行信号处置时只能采用有限长信号，所以需要将信号“截断”.在信号处置中，“截断”被看成是用一个有限长的“窗口”看无限长的信号，或者从分析的角度是无限长的信号x(t)乘以有限长的窗函数w(t).由傅立叶变换性质可知,时间域内的乘积对应于频率域的卷积,即 （7-9）这里，x(t)是频宽有限信号，而w(t)是频宽无限信号，暗示互为Fourier变换对.截断后的信号也必需是频宽无限信号，这样就是有限频带的信号分散到无限频带中去，这样就发生了所谓频谱泄
10、漏.从能量的角度来看，频谱泄漏也是能量的泄漏，因为加窗后使原来信号集中的窄频带内的能量分散到无限的频带宽度范围内.频谱泄漏是不成防止的，但要尽量减小.上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，能否改变窗的形状而到达这个目的呢？回答是肯定的.其实数字信号处置的前驱者们设计了分歧于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足分歧的要求.为此，用窗函数法设计FIR数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度N和窗函数w(n).下面我们介绍窗函数.7.2.2 MATLAB信号处置中提供的窗函数（1）矩形窗：前面分析中所
11、用的矩形窗可用下面函数来实现w=boxcar (N)，N 为窗的长度（以下函数与此同），w为返回的窗函数序列.（2）汉宁窗：w=hanning(N)汉宁窗的表达式为： (7-10)（3）哈明窗：w=hamming(N)哈明窗的表达式为： (7-11)（4）Bartlett窗：w=bartlett(N)Bartlett 窗的表达式为：当 N 为奇数时， (7-12)当 N 为偶数时， (7-13)（5） Blackman 窗：w= blackman(N)Blackman 窗的表达式为： (7-14)Blackman 窗比其他相同尺寸窗 (哈明窗,汉宁窗) 具有主瓣较宽和旁瓣泄漏较小的特点.（6）
12、三角窗：w=triang(N)三角窗的表达式为：当 N 为奇数时， (7-15)当 N 为偶数时， (7-16)三角窗和Bartlett窗十分类似.三角窗的两端值不为零，而Bartlett窗则为零，这一点可从例7-1中看出.（7）Kaiser窗：w=kaiser(n,beta)其中，beta是Kaiser窗参数,影响窗旁瓣幅值的衰减率.Kaiser窗表达式： (7-17)式中， I0.是修正过的零阶 Bessel 函数.Kaiser窗用于滤波器设计时，若旁瓣幅值为，则 （ 7-18 ）（8） Chebyshev窗：w=chebwin(n,r)式中, r 是窗口的旁瓣幅值在主瓣以下的分贝数.Ch
13、ebyshev窗具有主瓣宽度最小，而旁瓣等高,高度可调整的特点.下面我们在MATLAB观看各种窗函数的形状和频率域图象来验证上述所讲特点.【例7-1】 用MATLAB编程绘制各种窗函数的形状.窗函数的长度为21.%Samp7_1clfNwin=21;n=0:Nwin-1;%数据总数和序列序号figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin);%矩形窗 stext=矩形窗; case 2 w=hanning(Nwin);%汉宁窗 stext=汉宁窗; case 3 w=hamming(Nwin);%哈明窗 stext=哈明窗; case 4 w
14、=bartlett(Nwin);%Bartlett窗 stext=Bartelett窗; end posplot=2,2,int2str(ii);%指定绘制窗函数的图形位置subplot(posplot); stem(n,w);%绘出窗函数 hold on plot (n ,w,r);%绘制包络线 xlabel(n); ylabel(w(n); title(stext); hold off; grid onendfigure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin);%Blackman 窗 stext=Blackman窗; case 2 w=
15、triang(Nwin);%三角窗 stext=三角窗; case 3 w=kaiser(Nwin,4);%Kaiser窗 stext=Kaiser窗(Beta=4); case 4 w=chebwin(Nwin,40);%Chebyshev 窗 stext=Chebyshev窗(r=40); end posplot=2,2,int2str(ii); subplot(posplot); stem(n,w);%绘出窗函数 hold on plot (n,w,r);%绘出包络线 xlabel(n);ylabel(w(n);title(stext); hold off;grid on;end法式运行
16、结果见图 7-2 .可以看到各种窗函数的形状.图 7-2 各种窗函数的时间域形状【例 7-2】 用 MATLAB 编程，采纳512个频率点绘制各种窗函数的幅频特性.%Samp7_2clf;Nf=512;%窗函数复数频率特性的数据点数Nwin=20;%窗函数数据长度figure(1)for ii=1:4 switch ii case 1 w=boxcar(Nwin);%矩形窗 stext=矩形窗; case 2 w=hanning(Nwin);%汉宁窗 stext=汉宁窗; case 3 w=hamming (Nwin);%哈明窗 stext=哈明窗; case 4 w=bartlett(Nwi
17、n);%Bartlett窗 stext=Bartelett窗; end y,f=freqz(w,1,Nf);%求解窗函数的幅频特性，窗函数相当于一个数字滤波器 mag=abs(y);%求得窗函数幅频特性 posplot=2,2,int2str(ii); subplot(posplot); plot(f/pi,20* log10(mag/max(mag);%绘制窗函数的幅频特性xlabel(归一化频率);ylabel(振幅/dB); title(stext);grid on;endfigure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin);%Bla
18、ckman 窗 stext=Blackman窗; case 2 w=triang(Nwin);%三角窗 stext=三角窗; case 3 w=kaiser(Nwin,4);%Kaiser窗 stext=Kaiser窗(Beta=4); case 4 w=chebwin(Nwin,40);%Chebyshev 窗 stext=Chebyshev窗(r=40); end y,f=freqz(w,1,Nf);%以 Nf点数求解窗函数的幅频响应 mag=abs(y);%求得窗函数幅频响应 posplot=2,2,int2str(ii); subplot(posplot);plot(f/pi,20*
19、log10(mag/max(mag); %绘制幅频响应 xlabel(归一化频率);ylabel(振幅/dB); title(stext);grid on;end法式运行结果见图 7-3 .可以看到各种窗函数的幅频形状.对比该图可知这些窗函数具有上面所分析的窗函数的特征.图 7-3 各种窗函数的幅频形状由图 7-3 可见，各种窗函数都具有明显的主瓣（Mainlobe）和旁瓣（Sidelobe）.主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场所.矩形窗具有最窄的主瓣，但也有最年夜的旁瓣峰值（第一旁瓣衰减为 13 dB）；Blackman窗具有最年夜的旁瓣衰减，但也具有最宽的主瓣宽度.分歧窗函数
20、在这两方面的特点是分歧的，因此应根据具体的问题进行选择.通常来讲，哈明窗和汉宁窗的主瓣具有较小的旁瓣和较年夜的衰减速度，是较为经常使用的窗函数.表7-2总结了各种窗函数主瓣和旁瓣的特征（理论分析可参考其他的数字信号处置教材），年夜家可对比窗函数的幅频形状（图7-3）认真理解体会.表7-2 各种窗函数的特点窗函数主瓣宽第一旁瓣相对主瓣衰减（分贝）矩形窗-13汉宁窗-31哈明窗-41Bartlett窗-25Blackman 窗-57三角窗-25Kaiser窗可调整可调整Chebyshev 窗可调整可调整主旁瓣频率宽度还与窗函数长度N有关.增加窗函数长度N将减小窗函数的主瓣宽度，但不能减小旁瓣幅值衰
21、减的相对值（分贝数），这个值是由窗函数决定的.这个特点可由下面的例子清楚地看出.【例7-3】绘制矩形窗函数的幅频响应，窗长度分别为：(1)N=10;(2)N=20;(3)N=50;(4)N=100.%Samp7_3clf;Nf=512;for ii=1:4 switch ii case 1 Nwin=10; case 2 Nwin=20; case 3 Nwin=50; case 4 Nwin=100; end w=boxcar(Nwin);%矩形窗 y,f=freqz(w,1,Nf);%用分歧的窗长度求得复数频率特性 mag=abs(y);%求得幅频特性 posplot=2,2,int2st
22、r(ii);%指定绘图位置 subplot(posplot); plot (f/pi,20*log10(mag/max(mag);%绘出幅频形状 xlabel(归一化频率);ylabel(振幅/dB); stext=N=int2str(Nwin);%给出题目，指出所用的数据个数 title(stext);grid on;end图 7-4 数据长度分歧的矩形窗的幅频形状法式运行结果见图7-4.显然，随着N的增年夜，主瓣和旁瓣都变窄，但第一旁瓣相对主瓣的幅值下降分贝数相同，第二旁瓣相对第一旁瓣幅值下降的分贝数也相同.这样，随着N的增年夜，旁瓣也获得抑制，有力地减少了频谱泄漏，但不能完全消除.减少主
23、瓣宽度和抑制旁瓣是一对矛盾，不成兼得，只能根据分歧用途折衷处置.7.2.3 运用窗函数设计数字滤波器用于信号分析中的窗函数可根据用户的分歧要求选择.用于滤波器的窗函数，一般要求窗函数的主瓣宽度窄，以获得较好的过渡带；旁瓣相对值尽可能少，增加通带的平稳度和增年夜阻带的衰减.基于窗函数的FIR数字滤波器设计的算法十分简单，其主要步伐为：(1)对滤波器理想频域幅值响应进行傅立叶逆变换获得理想滤波器的单元脉冲响应hd(n).一般假定理想低通滤波器的截止频率为，其幅频特性满足 (7-19)则根据傅立叶逆变换，单元脉冲响应为： (7-20)其中，为信号延迟.（2）由性能指标（阻带衰减的分贝数）根据表7-2
24、第3列的值确定满足阻带衰减的窗函数类型w(n).滤波器的阶数越高，滤波器的幅频特性越好，但数据处置的费用也越高，因此像IIR滤波器一样，FIR滤波器也要确定满足性能指标的滤波器最小阶数.由前面的讨论（图7-1）可知，滤波器的主瓣宽度相当于过渡带宽，因此，使过渡带宽近似于窗函数主瓣宽（表7-2中的第二列）可求得满足性能指标的窗口长度N，此时，信号延迟为(N-1)/2.（3）求实际滤波器的单元脉冲响应h(n)：根据h(n)=hd(n)*w(n).（4）检验滤波器的性能.可设定一些信号采纳 7.1.2 节指出的函数进行滤波.下面采纳实例说明如何根据上面步伐设计FIR滤波器.【例 7-4】 用窗函数设
25、计一个线性相位FIR低通滤波器，并满足性能指标：通带鸿沟的归一化频率wp=0.5,阻带鸿沟的归一化频率ws=0.66，阻带衰减不小于30dB，通带波纹不年夜于3dB.假设一个信号，其中f1=5Hz,f2=20Hz.信号的采样频率为50Hz.试将原信号与通过滤波器的信号进行比力.由题意，阻带衰减不小于30dB,根据表7-2,选取汉宁窗，因为汉宁窗的第一旁瓣相对主瓣衰减为31dB，满足滤波要求.在窗函数设计法中，要求设计的频率归一化到0区间内，Nyquist频率对应于.法式如下%Samp7_4wp=0.5*pi;ws=0.66*pi;%滤波器鸿沟频率wdelta=ws-wp; %过渡带宽N=cei
26、l(8*pi/wdelta) %根据过渡带宽即是表 7-2中汉宁窗函数主瓣宽求得滤波器所用窗函数的最小长度Nw=N;wc=(wp+ws)/2;%截止频率在通带和阻带鸿沟频率的中点n=0:N-1;alpha=(N-1)/2;%求滤波器的相位延迟m=n-alpha+eps; %eps为MATLAB系统的精度hd=sin(wc*m)./(pi*m);%由（7-20）式求理想滤波器脉冲响应win=hanning(Nw);%采纳汉宁窗h=hd.*win;%在时间域乘积对应于频率域的卷积b=h; figure(1)H,f=freqz(b,1,512,50);%采纳 50 Hz 的采样频率绘出该滤波器的幅频
27、和相频响应subplot(2,1,1),plot(f,20*log10(abs(H)xlabel(频率/Hz);ylabel(振幅/dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;%impz(b,1);%可采纳此函数给出滤波器的脉冲响应%zplane(b,1);%可采纳此语句给出滤波器的零极点图%grpdelay(b,1);%可采纳此函数给出滤波器的群延迟f1=3;f2=20;%检测输入信号含有两种频率成份dt=0.02; t=0:dt:3; %采样间隔和检测信
28、号的时间序列x=sin(2*pi*f1* t)+cos(2* pi*f2* t);%检测信号%y=filter(b,1,x);%可采纳此函数给出滤波器的输出y=fftfilt(b,x);%给出滤波器的输出figure(2)subplot(2,1,1), plot(t,x),title(输入信号)%绘输入信号subplot(2,1,2),plot(t,y)%绘输出信号hold on; plot(1 1*(N-1)/2*dt,ylim, r) %绘出延迟到的时刻xlabel(时间/s),title(输出信号)法式运行结果见图7-5和图7-6.该例设计通带鸿沟Hz, 其中Fs为采样频率，Fnorma
29、l为归一化频率.由图7-5上图可以看到，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不年夜于3dB的要求.在年夜于16.5Hz的频段上，阻带衰减年夜于30dB，满足设计要求.由图7-5下图可见，在通带范围内，相位频率为一条直线，标明该滤波器为线性相位.图7-6给出了滤波器的输入信号和输出信号，输入信号包括3Hz和20Hz的信号，由图7-5可知，20Hz的信号不能通过该滤波器，通过滤波器后只剩下3Hz的信号，输出结果也证明了这一点.但要注意，由于FIR滤波器所需的阶数较高，信号延迟(N-1)/2也较年夜，输出信号前面有一段直线就是延迟造成的.上述法式显示的N取50才华满足设计要求.
30、这样相位延迟为s.验证了FIR滤波器相位延迟的理论.在输出信号的前部，有一些小信号，这是截断信号鸿沟所致，后面的部份就没有了这种信号.若采纳零相位的filtfilt函数（说明见第六章第三节）输出，则可最年夜限度地减小鸿沟的影响.图 7-5 例7-4所设计滤波器的幅频响应（上图）和相频响应（下图）图7-6 例7-4所设计滤波器的输入和输出信号7.2.4标准型FIR滤波器7.2.3节给出了运用理想脉冲响应与窗函数乘积的方法给出了滤波器传递函数的设计方法.其实MATLAB已将上述方法复合成一个函数，提供基于上述原理设计标准型FIR滤波器的工具函数.fir1就是采纳经典窗函数法设计线性相位FIR数字滤波器的函数，且具有标准低通，带通，高通,带阻等类型.函数调用格式为：b=fir1(n,wn,ftype,window)式中，n为FIR滤波器的阶数，对高通,带阻滤波器，n需取偶数；wn为滤波器截止频率，范围为01（归一化频率）.对带通,带阻滤波器，wn=w1,w2
