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

第9章FIR数字滤波器设计教学目标通过本章的学习，理解FIR数字滤波器的设计机理；掌握FIR数字滤波器的基本设计方法，包括：窗函数法、频率采样法等；掌握使用MATLAB实现以上设计方法的基本函数和设计流程。9.1引言IIR数字滤波器虽然有许多优良的特性，但是也有一些缺点，例如，设计IIR滤波器必须注意稳定性问题，再者若想利用快速傅里叶变换技术进行快速卷积实现滤波器，则必须要求单位脉冲响应是有限长的。此外，IIR滤波器的优异幅度响应，一般是以相位的非线性为代价的，非线性相位或引起频率色散。而FIR滤波器具有严格的线性相位特性，并且FIR滤波器的单位脉冲响应()hn仅含有有限个非零值，是因果有限长序列，因此稳定性能够得以保证。本章首先介绍较易掌握的窗函数法，基于对指定的频率响应的傅里叶级数进行截断来设计FIR数字滤波器；而后，介绍一种等波纹线性相位FIR数字滤波器的设计，主要采用加权切比雪夫逼近方法；其次，介绍一种频率采样方法，从FIR数字滤波器系统函数之间存在频率抽样关系的角度设计滤波器。最后，应用MATLAB对以上算法加以实现，并针对生物医学信号的特点进行了滤波器的设计。9.2FIR数字滤波器设计的基本方法与IIR数字滤波器的设计不同，FIR滤波器的设计与模拟滤波器的设计没有任何联系。故前一章所讨论的各种变换法对FIR滤波器不适用，因为它们给出的是IIR函数。FIR滤波器的设计基于对给定幅度响应的直接逼近，同时还需要考虑线性相位响应。一个长度为N的FIR滤波器传输函数()Hz是1z的1N次多项式，其对应的频率响应函数为：10()()NjjnnHehne(9.2.1)由于任何长度为N的有限长度序列()xn的特性可由其离散时间傅里叶变换()jXe的2N个样本来描述，故要设计一个长度为N的FIR滤波器，即可利用冲激响应序列()hn或其频率响应()jHe的N个样本来进行。此外，同样也需要估计FIR数字滤波器的阶数。下面介绍两种常用的估计数字滤波器阶数的最小值方程。Kaiser方程：-20lg1314.6()/2psspN(9.2.2)其中，p为归一化通带边界角频率，s为归一化阻带边界角频率，p为峰值通带波纹，s为峰值阻带波纹。Bellanger方程：2lg(10)13()/2psspN(9.2.3)以上两种阶数估计方法同样可以应用于估算高通、带通和带阻FIR滤波器的阶数。下面介绍几种具体的用于设计FIR数字滤波器的方法及MATLAB实现过程。9.3窗函数法设计FIR数字滤波器窗函数法，也称为傅里叶级数法。一般是先给定所要求的理想滤波器频率响应()jdHe，要求设计一个FIR数字滤波器，使其频率响应10()()NjjnnHehne无限逼近()jdHe。由于该设计过程通常都是在时域进行的，因而先由()jdHe的傅里叶逆变换得到()dhn，即()()jwjwnddnHehne(9.3.1)1()()2jjnddhnHeed(9.3.2)一般来说，理想的选频滤波器的()jdHe是逐段恒定的，且在频带边界处有不连续点，因此序列()dhn是无限长的，这时不能用(9.3.1)式中的傅氏级数系数来设计滤波器。因为第一，滤波器的单位取样响应()dhn是无限长的，n从到无法求和；第二由于()dhn是从开始的，所以是非因果的，且不能用有限的延时来实现它。第9章FIR数字滤波器设计3为解决上述问题，可采取如下办法(1)用有限项和来逼近无限项和。由傅里叶级数理论可知，(9.3.1)式级数的有限部分和(1)/2(1)/2()NjwndnNhne是()jdHe在均方意义下的最优逼近。这时用到的单位脉冲响应显然是()dhn的一段，是有限长的。(2)将有限长的()dhn进行(N-1)/2的有限延时，从而由非因果的系统得到了因果的系统。幅频特性完全不被时延所影响。因为时域的延时，在频域相应于线性相移，往往可以忽略。所以要用有限长的()来逼近无限长的()，最有效的方法就是截断()，也就是用一个有限长度的窗函数序列()来截取()，即()()()dhnhnn(9.3.3)由此可见，窗函数序列的形状及长度的选择是该方法设计FIR滤波器的关键。加窗处理对理想矩形频率响应会产生如下影响：第一，使理想频率特性不连续点处边沿加宽，形成一个过渡带，其宽度等于窗的频率响应()RW的主瓣宽度4N。第二，在截止频率c的两边2cN的地方(即过渡带的两边)，()H出现最大的肩峰值，肩峰的两侧形成起伏振荡，其振荡幅度取决于旁瓣的相对幅度，而振荡的多少，则取决于旁瓣的多少。第三，增加截取长度N，则在主瓣附近的窗的频率响应为：sin()sin()sin22()sin()22RNNxWNx(9.3.4)其中2xN。可见，改变N，只能改变窗谱的主瓣宽度、改变坐标的比例以及改变()RW的绝对值大小，但是不能改变主瓣与旁瓣的相对比例，这个相对比例是由sinxx决定的，也就是说只由窗函数的形状来决定。表9.3.1归纳了几种常用窗函数的主要性能，供设计FIR滤波器时参考。4表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)获取所要求的频率响应函数()jdHe；(2)利用式(9.3.2)求()dhn；(3)由过渡带带宽及最小衰减的要求，选定窗()n的形状及N的大小，一般N要通过几次试探而最后确定；(4)根据式(9.3.3)，求得所设计的FIR滤波器的单位抽样响应；(5)求()[()]jHeDTFThn，并检验其是否满足设计要求。下面举两个具体的生物医学信号的例子，来具体演示使用窗函数法设计FIR滤波器的实现过程。【例9.3.1】用窗函数法设计线性相位数字低通滤波器，对脉搏波信号进行去噪，要求通带截止频率15Hz，阻带截止频率25Hz，通带最大衰减1pdB，阻带最小衰减为40dB，抽样频率为100Hz，要求滤波器阶数尽量低。解：(1)确定数字滤波器过渡带宽。由已知条件可知，数字滤波器的数字通带截止频率为：2/215/1000.3ppffs。数字滤波器的阻带截止频率为:2/225/1000.5sstffs。数字滤波器的过渡带宽:第9章FIR数字滤波器设计50.2csp。数字滤波器的阻带最小衰减为：40sdB。选择窗函数，计算窗函数长度N，满足阻带衰减表的窗函数包括汉宁窗、汉明窗、布莱克曼窗和凯塞窗，根据过渡带宽，可以计算出前三种窗对窗长度要求分别为31N，33N，33N，采样凯塞窗，此时需要窗长度为22.3N，很显然，只有用凯塞窗设计出的阶数最低，根据题意，选择凯塞窗。这里N取24。(2)加窗，确定系统单位脉冲响应()hnsin[()sin[0.4()()()()cdwnnhnnn，12N。sin[0.4()()()()()()dnhnhnwnwnn，12N。(3)MATLAB对以上算法的实现代码如下：clearall;closeall;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(211);plot(data,'r');title('原脉搏波形');subplot(212);plot(newpule,'g');title('去噪后波形');采用脉搏波信号验证该滤波器，其仿真结果如图9.3.1所示：6图9.3.1例9.3.1的仿真结果【例9.3.2】编程设计一个50阶多通带FIR滤波器，通带分别为[0,0.1]、[0.3,0.5]和[0.7,0.9]，要求采用汉明窗实现。解：其MATLAB实现代码如下：clc;clearall;closeall;wc=[0.10.30.50.70.9];%归一化截止频率N=51;hn=fir1(N-1,wc,'DC-1',hamming(N));%采用汉明窗设计低通滤波器omega=linspace(0,pi,512);%频率采样512个点mag=freqz(hn,1,omega);%计算频率响应magdb=20*log10(abs(mag));%计算对数幅度频率响应subplot(121);stem([0:N-1],hn,'.');gridon;xlabel('n');ylabel('h(n)');title('单位抽样响应');subplot(122);plot(omega/pi,magdb);gridon;00.511.522.580100120140原脉搏波形t(s)Pressure(mmHg)00.511.522.580100120140去噪后波形t(s)Pressure(mmHg)第9章FIR数字滤波器设计7xlabel('\omega/\pi');ylabel('db');title('幅度频率响应');set(gca,'xtick',[00.10.30.50.70.91]);set(gca,'xticklabel',[00.10.30.50.70.91]);axis([01-8050]);set(gca,'ytick',[-80-44-60]);set(gca,'yticklabel',[-80-44-60]);其仿真结果如图9.3.2所示：图9.3.2例9.3.2的仿真结果9.4频率采样法设计FIR滤波器窗函数法是从时域出发，用窗函数截取理想的()dhn获得()hn，以有限长序列()hn去逼近理想的()dhn这样得到的频率响应()jwHe以频域方差最小的方式逼近于理想的()jwdHe。频率抽样法是从频域出发，以频域抽样点内插的方式来逼近理想()jwdHe。因为有限长序列()hn又可以用离散傅里叶变换()Hk来唯一确定，()Hk与所要求的FIR数字滤波器系统函数()dhn之间存在着频率抽样关系。即()dHz在z平面单位圆上按角度等分的抽样值等于()dHk的各相应值，以此()dHk作为实际FIR数字滤波器频率响应的抽样值()Hk，或者说()Hk正是所要求的频率响应的N个等间隔的抽样值。我们知道，一个有限长序列，同样可以用N个频域的抽样值来唯一确定。在前面我们讨论过用抽样值表达z函数的内插公式为：1101()()1NNkkNzHkHzNWz(9.4.1)式中，()Hk是频率抽样值2()()()kNjkNzWHkHzHe(9.4.2)01020304050-0.100.10.20.30.40.50.6nh(n)单位抽样响应00.10.30.50.70.91-80-70-60-50-40-30-20-10010/dB幅度频率响应8这就为设计FIR滤波器开创了另外一条通道，即直接从频域出发，对理想频率特性()jdHe抽样22/()()()jkjNdkNddHeHeHk(9.4.3)以此()dHk作为实际FIR数字滤波器的频率响应的抽样值()Hk,即令()()dHkHk，0,1,,1kN(9.4.4)知道()Hk后，由傅里叶逆变换定义，可以用这N个抽样值()Hk来唯一确定有限长序列()hn，即2110011()()()NNjnknkNNkkhnHkWHkeNN，0,1,2,