第九章 FIR数字滤波器设计20150907

1、第9章 FIR数字滤波器设计第9章 FIR数字滤波器设计教学目标通过本章的学习，理解FIR数字滤波器的设计机理；掌握FIR数字滤波器的基本设计方法，包括：窗函数法、频率采样法等；掌握使用MATLAB实现以上设计方法的基本函数和设计流程。9.1 引言IIR数字滤波器虽然有许多优良的特性，但是也有一些缺点，例如，设计IIR滤波器必须注意稳定性问题，再者若想利用快速傅里叶变换技术进行快速卷积实现滤波器，则必须要求单位脉冲响应是有限长的。此外，IIR滤波器的优异幅度响应，一般是以相位的非线性为代价的，非线性相位或引起频率色散。而FIR滤波器具有严格的线性相位特性，并且FIR滤波器的单位脉冲响应仅含有有

2、限个非零值，是因果有限长序列，因此稳定性能够得以保证。本章首先介绍较易掌握的窗函数法，基于对指定的频率响应的傅里叶级数进行截断来设计FIR数字滤波器；而后，介绍一种等波纹线性相位FIR数字滤波器的设计，主要采用加权切比雪夫逼近方法；其次，介绍一种频率采样方法，从FIR数字滤波器系统函数之间存在频率抽样关系的角度设计滤波器。最后，应用MATLAB对以上算法加以实现，并针对生物医学信号的特点进行了滤波器的设计。9.2 FIR数字滤波器设计的基本方法与IIR数字滤波器的设计不同，FIR滤波器的设计与模拟滤波器的设计没有任何联系。故前一章所讨论的各种变换法对FIR滤波器不适用，因为它们给出的是IIR函

3、数。FIR滤波器的设计基于对给定幅度响应的直接逼近，同时还需要考虑线性相位响应。一个长度为的FIR滤波器传输函数是的次多项式，其对应的频率响应函数为： (9.2.1)由于任何长度为的有限长度序列的特性可由其离散时间傅里叶变换的个样本来描述，故要设计一个长度为的FIR滤波器，即可利用冲激响应序列或其频率响应的个样本来进行。此外，同样也需要估计FIR数字滤波器的阶数。下面介绍两种常用的估计数字滤波器阶数的最小值方程。Kaiser方程： (9.2.2)其中，为归一化通带边界角频率，为归一化阻带边界角频率，为峰值通带波纹，为峰值阻带波纹。Bellanger方程： (9.2.3)以上两种阶数估计方法同样

4、可以应用于估算高通、带通和带阻FIR滤波器的阶数。下面介绍几种具体的用于设计FIR数字滤波器的方法及MATLAB实现过程。9.3 窗函数法设计FIR数字滤波器窗函数法，也称为傅里叶级数法。一般是先给定所要求的理想滤波器频率响应，要求设计一个FIR数字滤波器，使其频率响应无限逼近。由于该设计过程通常都是在时域进行的，因而先由的傅里叶逆变换得到，即 (9.3.1) (9.3.2)一般来说，理想的选频滤波器的是逐段恒定的，且在频带边界处有不连续点，因此序列是无限长的，这时不能用(9.3.1)式中的傅氏级数系数来设计滤波器。因为第一，滤波器的单位取样响应是无限长的，n从到无法求和；第二由于是从开始的，

5、所以是非因果的，且不能用有限的延时来实现它。 为解决上述问题，可采取如下办法(1) 用有限项和来逼近无限项和。由傅里叶级数理论可知，(9.3.1)式级数的有限部分和是在均方意义下的最优逼近。这时用到的单位脉冲响应显然是的一段，是有限长的。(2) 将有限长的进行(N-1)/2的有限延时，从而由非因果的系统得到了因果的系统。幅频特性完全不被时延所影响。因为时域的延时，在频域相应于线性相移，往往可以忽略。所以要用有限长的hn来逼近无限长的hdn，最有效的方法就是截断hdn，也就是用一个有限长度的窗函数序列n来截取hdn，即 (9.3.3)由此可见，窗函数序列的形状及长度的选择是该方法设计FIR滤波器

6、的关键。加窗处理对理想矩形频率响应会产生如下影响：第一，使理想频率特性不连续点处边沿加宽，形成一个过渡带，其宽度等于窗的频率响应的主瓣宽度。第二，在截止频率的两边的地方(即过渡带的两边)，出现最大的肩峰值，肩峰的两侧形成起伏振荡，其振荡幅度取决于旁瓣的相对幅度，而振荡的多少，则取决于旁瓣的多少。第三，增加截取长度N，则在主瓣附近的窗的频率响应为： (9.3.4)其中。可见，改变N，只能改变窗谱的主瓣宽度、改变坐标的比例以及改变的绝对值大小，但是不能改变主瓣与旁瓣的相对比例，这个相对比例是由决定的，也就是说只由窗函数的形状来决定。表9.3.1归纳了几种常用窗函数的主要性能，供设计FIR滤波器时参

7、考。表9.3.1 几种常用窗函数的主要性能参数 窗函数 主瓣宽度 精确过渡带旁瓣峰值衰减 /dB阻带最小衰减/dB 矩形窗4/N1.8/N-13-21 三角窗8/N6.1/N-25-25 汉宁窗8/N6.2/N-31-44 海明窗8/N6.6/N-41-53 布拉克曼窗 12/N11/N-57-74 贝塞尔窗 自由选择10/N-57-80窗函数法设计数字滤波器的基本步骤如下：(1) 获取所要求的频率响应函数；(2) 利用式(9.3.2)求；(3) 由过渡带带宽及最小衰减的要求，选定窗的形状及N的大小，一般N要通过几次试探而最后确定；(4) 根据式(9.3.3)，求得所设计的FIR滤波器的单位抽

8、样响应；(5) 求，并检验其是否满足设计要求。下面举两个具体的生物医学信号的例子，来具体演示使用窗函数法设计FIR滤波器的实现过程。【例9.3.1】用窗函数法设计线性相位数字低通滤波器，对脉搏波信号进行去噪，要求通带截止频率15Hz，阻带截止频率25Hz，通带最大衰减，阻带最小衰减为，抽样频率为100Hz，要求滤波器阶数尽量低。解：(1)确定数字滤波器过渡带宽。由已知条件可知，数字滤波器的数字通带截止频率为：。数字滤波器的阻带截止频率为:。数字滤波器的过渡带宽:。数字滤波器的阻带最小衰减为：。选择窗函数，计算窗函数长度，满足阻带衰减表的窗函数包括汉宁窗、汉明窗、布莱克曼窗和凯塞窗，根据过渡带宽

9、，可以计算出前三种窗对窗长度要求分别为，采样凯塞窗，此时需要窗长度为，很显然，只有用凯塞窗设计出的阶数最低，根据题意，选择凯塞窗。这里取24。(2) 加窗，确定系统单位脉冲响应 ， 。 ， 。(3) MATLAB对以上算法的实现代码如下：clear all; close all; clc; fs=100;wp=0.3; %归一化通带截止频率ws=0.5; %归一化阻带截止频率wc=0.4;N=24;data=load(pulse.txt);hn=fir1(N-1,wc,low,kaiser(N); %采用凯塞窗设计低通滤波器newpule=filtfilt(hn,1,data);subplot

10、(211);plot(data,r);title(原脉搏波形);subplot(212);plot(newpule,g);title(去噪后波形);采用脉搏波信号验证该滤波器，其仿真结果如图9.3.1所示：图9.3.1 例9.3.1的仿真结果【例9.3.2】编程设计一个50阶多通带FIR滤波器，通带分别为、和，要求采用汉明窗实现。解：其MATLAB实现代码如下：clc; clear all; close all;wc=0.1 0.3 0.5 0.7 0.9; %归一化截止频率 N=51;hn=fir1(N-1,wc,DC-1,hamming(N); %采用汉明窗设计低通滤波器omega=lin

11、space(0,pi,512); %频率采样512个点mag=freqz(hn,1,omega); %计算频率响应magdb=20*log10(abs(mag); %计算对数幅度频率响应subplot(121);stem(0:N-1,hn,.);grid on;xlabel(n);ylabel(h(n);title(单位抽样响应);subplot(122);plot(omega/pi,magdb);grid on;xlabel(omega/pi);ylabel(db);title(幅度频率响应);set(gca,xtick,0 0.1 0.3 0.5 0.7 0.9 1);set(gca,xt

12、icklabel,0 0.1 0.3 0.5 0.7 0.9 1);axis(0 1 -80 50);set(gca,ytick,-80 -44 -60);set(gca,yticklabel,-80 -44 -60);其仿真结果如图9.3.2所示：图9.3.2 例9.3.2的仿真结果9.4 频率采样法设计FIR滤波器窗函数法是从时域出发，用窗函数截取理想的获得，以有限长序列去逼近理想的这样得到的频率响应以频域方差最小的方式逼近于理想的。频率抽样法是从频域出发，以频域抽样点内插的方式来逼近理想。因为有限长序列又可以用离散傅里叶变换来唯一确定，与所要求的FIR数字滤波器系统函数之间存在着频率抽样

13、关系。即在z平面单位圆上按角度等分的抽样值等于的各相应值，以此作为实际FIR数字滤波器频率响应的抽样值，或者说正是所要求的频率响应的N个等间隔的抽样值。我们知道，一个有限长序列，同样可以用N个频域的抽样值来唯一确定。在前面我们讨论过用抽样值表达z函数的内插公式为： (9.4.1) 式中，是频率抽样值 (9.4.2) 这就为设计FIR滤波器开创了另外一条通道，即直接从频域出发，对理想频率特性抽样 (9.4.3) 以此作为实际FIR数字滤波器的频率响应的抽样值,即令， (9.4.4) 知道后，由傅里叶逆变换定义，可以用这N个抽样值来唯一确定有限长序列，即， (9.4.5)式中，为带设计的滤波器的单

14、位冲激响应。再将进行离散时间傅里叶变换，便得到设计的滤波器的频率响应为 (9.4.6)因此，可以得出： (9.4.7)即 (9.4.8)式中 (9.4.9)称为内插函数。由此可见，设计指标所要求的频率响应与实际设计的频率响应在采样点上完全相等，逼近误差等于零。但在采样点之间的频响是由各取样点的内插函数的延伸叠加形成，其逼近误差是有限的。误差的大小取决于频率响应曲线的圆滑程度和取样点的密度。频率响应曲线越圆滑，误差越小；取样点数越多，即取样频率越高，误差越小。以上就是频率采样法设计滤波器的基本原理。当设计线性相位FIR滤波器时，其采样值的幅度和相位一定要满足表9.4.1线性相位滤波器的约束条件。

15、表9.4.1线性相位滤波器的约束条件偶对称单位脉冲响应 h(n)=h(N-1-n) 情 况 相位响应 N为奇数 情 况 N为偶数奇对称单位脉冲响应 h(n)= - h(N-1-n) 情 况 相位响应 N为奇数 情 况 N为偶数(1) 对于偶对称，N为奇数时，其的表达式为 (9.4.10)其中，幅度函数关于为偶对称，即 (9.4.11)对在之间等间隔采样N点，得到，将其幅度(纯标量)与相角的形式表示 (9.4.12)根据式(9.4.10)和式(9.4.11)可知 ， (9.4.13)即满足偶对称。(2) 对于偶对称，N为偶数时，则其的表达式仍为式(9.4.10)。但是其幅度函数关于是奇对称的，关

16、于为偶对称，即 (9.4.14)所以，这时的也应满足奇对称要求。 (9.4.15)(3) 对于奇对称，N为奇数时，其的表达式为 (9.4.16)其中，幅度函数关于奇对称，即 (9.4.17)由式(9.4.12)和式(9.4.13)可知 , (9.4.18)即满足奇对称要求。(4) 对于奇对称，N为偶数，则其的表达式仍为 (9.4.19) (9.4.20)但是，其幅度函数关于是偶对称的，关于为奇对称，即 (9.4.21)所以，这时的也应该满足对称要求：， (9.4.22)频率采样法设计滤波器的频率特性如图9.4.1所示。用该滤波器特性逼近目标滤波器，通带和阻带出现波动，过渡带加宽而使实际的与理想

17、的相比产生误差。图9.4.1频率采样的响应从时域分析，若希望设计逼近理想低通滤波器，对应的单位脉冲响应为： (9.4.23)根据采样域采样定理，在频域之间等间隔采样N点并作离散傅里叶逆变换，得到的应是以N为周期，周期性延拓乘以，即 (9.4.24)如果有间断点，那么响应单位脉冲响应应是无限长的。由于时域混叠，引起所设计的和有偏差。要消除偏差，希望在频域的采样点数N加大。N越大，设计出的滤波器越逼近带设计的滤波器。从频域分析，有采样定理知，频率域等间隔采样,经过离散傅里叶逆变换得到，其z变换和的关系为 (9.4.25)将代入到上式，得到 (9.4.26)式中 (9.4.27)上式表明，在采样点。

18、在采样点处与相等，逼近误差为0.在两相邻采样点间，由有限项的之和形成，特性愈平滑，误差愈小，特性曲线间断点处误差最大。表现形式为间断点用倾斜线取代，且间断点附近形成振荡特性，使阻带衰减减小，往往导致不能满足技术指标要求。最直观的想法是增加采样点数，即加大N值，这时过渡带就等于采样间隔，即 (9.4.28)所以加大N，可使过渡带变窄，但增加要适当，否则会增加滤波器的体积与成本。因为是理想矩形，无论怎样增多频率采样的点数，在通、阻带交界处，幅值总是从1突变到0，必然会引起较大的起伏振荡。因此得出结论，增加N并不会改善滤波器的阻带衰减特性。改进在通带边缘由于采样点的徒然变化而引起的起伏振荡特性，考虑

19、在不连续点的边缘加上一些过渡的采样点，减小频带边缘的突变，也就能减小起伏振荡，增大阻带最小衰减。这些采样点上的取值不同，效果也就不同。因为每一个频率采样值都产生一个与常数成正比且在频率上位移的频率响应，而FIR数字滤波器的频率响应就是各与相应的内插函数相乘后的线性组合。设计适当的过渡带采样值，就可使它的有用频带(通带、阻带)的波纹得以减小，设计出性能较好的滤波器。一般过渡带取一、二、三点采样值即可得到满意结果，如图9.4.2所示。低通滤波器设计中不加过渡点时，阻带最小衰减为-20dB，加一个过渡点(采用最优化设计值)，阻带最小衰减可提高-40-54dB，加两个过渡点(最优设计)可达-60-75

20、dB，加三个过渡点(最优设计)则可达-80-95dB。图9.4.2理想低通滤波器增加过渡点增加过度点可使阻带衰减明显提高，代价是过渡带变宽。增加一个过渡点，过渡带宽为，增加两个过渡点，过度带宽为，式(9.4.28)修正为： (9.4.29)其中，m为增加的过渡点数。若过渡带不满足要求，可通过加大N来调整。现归纳频率采样法设计FIR数字滤波器的具体步骤如下：(1) 根据阻带最小衰减选择过渡带采样点个数m；(2) 确定过渡带宽度，估计频域采样点数N (即滤波器长度)。如果增加m个过度采样点，则过渡带宽度近似变成。当N确定后，m越大，过渡带越宽。如果给定过渡带宽度，则要求，滤波器长度N必须满足如下估

21、计公式： (9.4.30)(3) 构造一个希望逼近的频率响应函数 (9.4.31)设计标准常数特性的FIR数字滤波器时，一般构造幅频响应函数为相应的理想频率响应特性。(4) 按照式(9.3.12)和式(9.3.13)进行频率采样，加入过渡采样点。过渡带采样值可以设置为经验值，或用累试法确定，也可以采用优化算法计算。(5) 对进行N点IDFT，得到第一类线性相位FIR数字滤波器的单位脉冲响应。 检验设计结果。如果阻带最小衰减值未达到指标要求，则要改变过渡带采样值，直到满足指标要求为止。如果滤波器边界频率未达到指标要求，则需要微调的边界频率。具体流程如图9.4.3所示。图9.4.3频率采样法的设计

22、步骤下面举两个具体的例子，来演示频率采样法设计FIR数字滤波器的过程，及其MATLAB实现方法。【例9.4.1】给定一个理想低通FIR滤波器的频率特性现希望用函数抽样法设计该滤波器，要求具有线性相位。假定滤波器的长度为31。 解： 利用频率抽样法设计数字滤波器时，需要将理想的连续频率响应转换为理想的离散频率响应。 因为N为奇数，所以， 或令然后根据。即可求出所要设计的滤波器，且上述幅值方法保证了所设计的滤波器具有线性相位。其响应的MATLAB程序如下：%用频率采样法设计FIR数字滤波器fl=pi/6;N=31;factor=(N+1)/2;m=factor*fl/pi;m=floor(m);j

23、=sqrt(-1);for k=0:(N-1) if(k=N-m) a=-j*(N-1)*k*pi/N; H(k+1)=exp(a); else H(k+1)=0; endend;for i=0(N-1) sum=0; for k=0(N-1) a=j*2*pi*i*k/N; sum=sum+H(k+1)*exp(a); end h(i+1)=sum/N;endh=real(h);h1=freqz(h,1,N);t=0:0.5/N:0.5-0.5/N;figure(color,w);stem(n,h,.);xlabel(n);ylabel(h(n);grid on;figure(color,w

24、);plot(t,abs(h1);xlabel(w/2pi);ylabel(Amplitude);grid on;设计出滤波器的频谱抽样序列图和幅频响应图分别如图9.4.4和图9.4.5所示。 图9.4.4低通滤波器的频谱抽样序列图 图9.4.5低通滤波器的幅频响应图【例9.4.2】给定以理性带阻滤波器的频率特性抽样点数N=33，现希望用频率抽样法来完成带阻滤波器设计。解：其响应的MATLAB程序如下： %用频率采样法设计FIR数字滤波器fl=pi/4;fh=3*pi/4;N=99;factor=(N+1)/2;m=factor*fl/pi;m=floor(m);m2=factor*fh/pi

25、;m2=ceil(m2);j=sqrt(-1);n=0:N-1;for k=0:(N-1)if(k=N-m)|( (k=m2)&(k=(N-m2) a=-j*(N-1)*k*pi/N; H(k+1)=exp(a); else H(k+1)=0; endend;for i=0(N-1) sum=0; for k=0(N-1) a=j*2*pi*i*k/N; sum=sum+H(k+1)*exp(a); end h(i+1)=sum/N;endh=real(h);h1=freqz(h,1,N);t=0:0.5/N:0.5-0.5/N;figure(color,w);stem(n,h,.);xlab

26、el(n);ylabel(h(n);grid on;figure(color,w);plot(t,abs(h1);xlabel(w/2pi);ylabel(Amplitude);grid on;设计出滤波器的频谱抽样序列图和幅频响应图分别如图9.4.6和图9.4.7所示。 图9.4.6带阻滤波器的频谱抽样序列图 图9.4.7带阻滤波器的幅频响应图本 章 小 结本章主要讲述了FIR数字滤波器的基本设计方法和MATLAB实现过程，主要介绍了如下内容：1. FIR数字滤波器的基本设计原理及过程； 2. 使用窗函数法设计FIR数字滤波器的方法；3. 使用频率采样法设计FIR数字滤波器的方法。其中，以生物医学信号处理过程中的数字滤波器设计为主要研究对象进行举例，并附以MATLAB源程序，供读者实践使用。习 题9.1 窗函数的主要类型有哪些？有哪些频域指标？9.2 什么是Gibbs现象？为什么选用矩形窗会出现Gibbs现象？9.3 频率抽样法设计F