《FIR滤波器设计.pdf》由会员分享，可在线阅读，更多相关《FIR滤波器设计.pdf（66页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第 7 章 FIR 滤波器设计 第六章我们介绍了无限冲激响应（IIR）滤波器的设计方法。其中最常用的由模拟滤波器转换为数字滤波器的方法为双线性变换法，因为这种方法无混叠效应，效果较好。但通过前面的例子我们看到，IIR数字滤波器相位特性不好（非线性，如图 6-11、图6-13、图6-15），也不易控制。然而在现代信号处理中，例如图像处理、数据传输、雷达接收以及一些要求较高的系统中对相位特性要求较为严格，这种滤波器就无能为力了。改善相位特性的方法是采用有限冲激响应滤波器。本章首先对 FIR 滤波器原理及其使用函数作基本介绍，然后重点介绍窗函数法设计 FIR 滤波器，并对最优滤波器设计函数进行介绍。
2、FIR 滤波器原理概述及滤波函数 FIR 滤波器原理及设计方法分类 根据第 6 章对数字滤波器的介绍，我们知道 FIR 滤波器的传递函数为：10)()()()(NnnznhzXzYzH (7-1)可得 FIR 滤波器的系统差分方程为：10)()()()()1()1()1()1()()0()(NmnxnbmnxmbNnxNbnxbnxbny 因此，FIR 滤波器又称为卷积滤波器。根据第 4 章中所描述的系统频率响应，FIR 滤波器的频率响应表达式为：10)(NnjnjenbeH （7-2）信号通过 FIR 滤波器不失真条件与（6-6）式所描述的相同，即滤波器在通带内具有恒定的幅频特性和线性相位特
3、性。理论上可以证明（这里从略）：当 FIR 滤波器的系数满足下列中心对称条件：)1()()1()(nNbnbnNbnb或 （7-3）时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位 FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个 N 阶的线性相位 FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。本章主要介绍的 FIR 数字滤波器设计方法及 MATLAB 信号处理工具箱提供的 FIR 数字滤波器设计函数，见表 7-1。由于篇幅所限，本章我们主要介绍窗函数法和最优
4、化设计方法。表 7-1 FIR 滤波器设计的主要方法 函数设计方法 说明 工具函数 窗函数法 理想滤波器加窗处理 fir1(单 频 带)，fir2(多频带)，kaiserord 最优化设计 平方误差最小化逼近理想幅频响应或Park-McClellan 算法产生等波纹滤波器 firls，remez,remezord 约束最小二乘逼近 在满足最大误差限制条件下使整个频带平方误差最小化 fircls,fircls1 升余弦函数 具有光滑、正弦过渡带的低通滤波器设计 Fircos FIR 数字滤波器滤波函数 相对于 IIR 滤波器的滤波函数，FIR 数字滤波器滤波函数除了 dimpulse 和 dst
5、ep 仅适用于 IIR 滤波器外，其他各种函数可直接应用于 FIR 滤波器，只是输入的分母多项式向量a=1。另外，MATLAB 还提供了一个函数 fftfilt,该函数利用效率高的基于 FFT 算法实现对数据的滤波，该函数只适用于 FIR 滤波器，调用形式为：y=fftfilt(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=filte
6、r(b,a,x)。FIR 滤波器的窗函数设计 窗函数的基本原理 FIR 滤波器设计的主要任务是根据给定的性能指标确定滤波器的系数 b，即系统单位脉冲序列 h(n)，它是一个有限长序列。FIR 滤波器的理想频率响应，可写成复数形式的 Fourier 级数形式：nnjdjdenheH (7-4)式中，hd(n)是对应的单位脉冲响应序列。这说明滤波器的频率响应和单位脉冲响应互为Fourier 变换对。因此其单位脉冲响应可由下式求得，deeHnhnjjdd21 (7-5)求得序列 nhd后，通过 z 变换，可得到 zHd nnddznhzH)()(7-6)注意，这里 nhd为无限长序列，因此 zHd是
7、物理上不可实现的。如何变成物理上可实现呢一个自然的想法是只取其中的某些项，即只截取 nhd中的一部分，比如 n=0,N-1，N 为正整数。这种处理相当于将 nhd,n=-与函数 w(n)相乘，w(n)具有下列形式：NnNnnnw0,1,0,0)(w(n)相当于一个矩形，我们称之为矩形窗。即我们可采用矩形窗函数 w(n)将无限脉冲响应 nhd截取一段 h(n)来近似为 nhd，这种截取在数学上表示为：h(n)=nhdw(n)(7-7)这里应该强调的是，加窗函数不是可有可无的，而是将设计变为物理可实现所必须的。截取之后的滤波器传递函数变为：10)()(NnnznhzH (7-8)式中，N 为窗口宽
8、度，H(z)是物理可实现系统。为了获得线性相位，FIR 滤波器 h(n)必须满足中心对称条件（即 7-3 式），序列 h(n)的延迟为2/1 N。这种方法的基本原理是用一定宽度的矩形窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，从而得到 FIR 滤波器的脉冲响应，故称为 FIR 滤波器的窗函数设计法。经过加矩形窗后所得的滤波器实际频率响应能否很好地逼近理想频率响应呢图 7-1 示意给出了理想滤波器加矩形窗后的情况。理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图。时间域内的乘积（7-7）式要求实际频率响应为这两个频率响应函数在频域内的卷积（卷积定理），即得到图形为图 7
9、-1（右图）。图 7-1 FIR 滤波器理想与实际频率响应 由图可看出，加矩形窗后使实际频率响应偏离理想频率响应，主要影响有三个方面：（1）理想幅频特性陡直边缘处形成过渡带，过渡带宽取决于矩形窗函数频率响应的主瓣宽度。（2）过渡带两侧形成肩峰和波纹，这是矩形窗函数频率响应的旁瓣引起的，旁瓣相对值越大，旁瓣越多，波纹越多。（3）随窗函数宽度 N 的增大，矩形窗函数频率响应的主瓣宽度减小，但不改变旁瓣的相对值。为了改善 FIR 滤波器性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣相对值尽可能小，数量尽可能少，以获得通带波纹小,阻带衰减大，在通带和阻带内均平稳的特点，这样可使滤波器实际
10、频率响应更好地逼近理想频率响应。这里我们明确两个概念：截断和频谱泄漏。信号是无限长的，而在进行信号处理时只能采取有限长信号，所以需要将信号“截断”。在信号处理中，“截断”被看成是用一个有限长的“窗口”看无限长的信号，或者从分析的角度是无限长的信号 x(t)乘以有限长的窗函数 w(t)。由傅立叶变换性质可知,时间域内的乘积对应于频率域的卷积,即 )()()()(fWfXtwtx （7-9）这里，x(t)是频宽有限信号，而 w(t)是频宽无限信号，表示互为 Fourier 变换对。截断后的信号也必须是频宽无限信号，这样就是有限频带的信号分散到无限频带中去，这样就产生了所谓频谱泄漏。从能量的角度来看
11、，频谱泄漏也是能量的泄漏，因为加窗后使原来信号集中的窄频带内的能量分散到无限的频带宽度范围内。频谱泄漏是不可避免的，但要尽量减小。上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，能否改变窗的形状而达到这个目的呢回答是肯定的。其实数字信号处理的前驱者们设计了不同于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足不同的要求。为此，用窗函数法设计 FIR 数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度 N 和窗函数 w(n)。下面我们介绍窗函数。MATLAB 信号处理中提供的窗函数（1）矩形窗：前面分析中所用的矩形窗可用下面函数
12、来实现 w=boxcar(N)，N 为窗的长度（以下函数与此同），w 为返回的窗函数序列。（2）汉宁窗：w=hanning(N)汉宁窗的表达式为：NkNkkw,.,1,12cos15.0)(7-10)（3）哈明窗：w=hamming(N)哈明窗的表达式为：1,.,1,0,12cos46.054.0)1(NkNkkw (7-11)（4）Bartlett 窗：w=bartlett(N)Bartlett 窗的表达式为：当 N 为奇数时，NkNNkNkNkkw21,1)1(22211,1)1(2)(7-12)当 N 为偶数时，NkNNkNNkNkkw12,1)(2221,1)1(2)(7-13)（5）
13、Blackman 窗：w=blackman(N)Blackman 窗的表达式为：NkNkNkkw,.,1,11408.0112cos5.042.0)(7-14)Blackman 窗比其他相同尺寸窗(哈明窗,汉宁窗)具有主瓣较宽和旁瓣泄漏较小的特点。（6）三角窗：w=triang(N)三角窗的表达式为：当 N 为奇数时，NkNNkNNkNkkw21,1)1(2211,12)(7-15)当 N 为偶数时，NkNNkNNkNkkw12,)1(2221,112)(7-16)三角窗和 Bartlett 窗十分类似。三角窗的两端值不为零，而 Bartlett 窗则为零，这一点可从例 7-1 中看出。（7）
14、Kaiser 窗：w=kaiser(n,beta)其中，beta 是 Kaiser 窗参数,影响窗旁瓣幅值的衰减率。Kaiser 窗表达式：0201211)(INkIkw (7-17)式中，I0.是修正过的零阶 Bessel 函数。Kaiser 窗用于滤波器设计时，若旁瓣幅值为d，则 21,02150,2107886.0215842.050,7.81102.04.0 （7-18）（8）Chebyshev 窗：w=chebwin(n,r)式中,r 是窗口的旁瓣幅值在主瓣以下的分贝数。Chebyshev 窗具有主瓣宽度最小，而旁瓣等高,高度可调整的特点。下面我们在 MATLAB 观看各种窗函数的形
15、状和频率域图象来验证上述所讲特点。【例 7-1】用 MATLAB 编程绘制各种窗函数的形状。窗函数的长度为 21。%Samp7_1 clf Nwin=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=bartlett(Nwin);%Bartlett 窗 stext=Bartelett 窗;end
16、 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 figure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin);%Blackman 窗 stext=Blackman 窗;case 2 w=triang(Nwin);%三角窗 stext=三角窗;case 3 w=kaiser(Nwin,4);
17、%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 程序运行结果见图 7-2。可以看到各种窗函数的形状。0510152000.20.40.60.81nw(n)矩 形 窗0
18、510152000.20.40.60.81nw(n)汉 宁 窗0510152000.20.40.60.81nw(n)哈 明 窗0510152000.20.40.60.81nw(n)Bartelett窗 0510152000.51nw(n)Blackman窗0510152000.20.40.60.81nw(n)三 角 窗0510152000.20.40.60.81nw(n)Kaiser窗(Beta=4)0510152000.20.40.60.81nw(n)Chebyshev窗(r=40)图 7-2 各种窗函数的时间域形状【例 7-2】用 MATLAB 编程，采用 512 个频率点绘制各种窗函数的
19、幅频特性。%Samp7_2 clf;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(Nwin);%Bartlett 窗 stext=Bartelett 窗;end y,f=freqz(w,1,Nf);%求解窗函数的幅频特性，窗函数相当于一个数字滤
20、波器 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;end figure(2)for ii=1:4 switch ii case 1 w=blackman(Nwin);%Blackman 窗 stext=Blackman 窗;case 2 w=triang(Nwin);%三角窗 stext=三角窗;case 3 w=kaiser(N
21、win,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*log10(mag/max(mag);%绘制幅频响应 xlabel(归一化频率);ylabel(振幅/dB);title(stext);grid on;end
22、程序运行结果见图 7-3。可以看到各种窗函数的幅频形状。对照该图可知这些窗函数具有上面所分析的窗函数的特征。00.51-70-60-50-40-30-20-100归 一 化 频 率振幅/dB矩 形 窗00.51-120-100-80-60-40-200归 一 化 频 率振幅/dB汉 宁 窗00.51-120-100-80-60-40-200归 一 化 频 率振幅/dB哈 明 窗00.51-100-80-60-40-200归 一 化 频 率振幅/dBBartelett 窗 00.51-140-120-100-80-60-40-200归 一 化 频 率振幅/dBBlackman窗00.51-200
23、-150-100-500归 一 化 频 率振幅/dB三 角 窗00.51-100-80-60-40-200归 一 化 频 率振幅/dBKaiser窗(Beta=4)00.51-100-80-60-40-200归 一 化 频 率振幅/dBChebyshev窗(r=40)图 7-3 各种窗函数的幅频形状 由图 7-3 可见，各种窗函数都具有明显的主瓣（Mainlobe）和旁瓣（Sidelobe）。主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰减为 13 dB）；Blackman 窗具有最大的旁瓣衰减，但也具有最宽的主瓣宽度。不同窗函数在这
24、两方面的特点是不同的，因此应根据具体的问题进行选择。通常来讲，哈明窗和汉宁窗的主瓣具有较小的旁瓣和较大的衰减速度，是较为常用的窗函数。表 7-2 总结了各种窗函数主瓣和旁瓣的特征（理论分析可参考其他的数字信号处理教材），大家可对照窗函数的幅频形状（图 7-3）认真理解体会。表 7-2 各种窗函数的特点 窗函数 主瓣宽 第一旁瓣相对主瓣衰减（分贝）矩形窗 N4-13 汉宁窗 N8-31 哈明窗 N8-41 Bartlett 窗 N8-25 Blackman 窗 N12-57 三角窗 N8-25 Kaiser 窗 可调整 可调整 Chebyshev 窗 可调整 可调整 主旁瓣频率宽度还与窗函数长度
25、 N 有关。增加窗函数长度 N 将减小窗函数的主瓣宽度，但不能减小旁瓣幅值衰减的相对值（分贝数），这个值是由窗函数决定的。这个特点可由下面的例子清楚地看出。【例 7-3】绘 制 矩 形 窗 函 数 的 幅 频 响 应，窗 长 度 分 别 为：(1)N=10;(2)N=20;(3)N=50;(4)N=100.%Samp7_3 clf;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);%用
26、不同的窗长度求得复数频率特性 mag=abs(y);%求得幅频特性 posplot=2,2,int2str(ii);%指定绘图位置 subplot(posplot);plot(f/pi,20*log10(mag/max(mag);%绘出幅频形状 xlabel(归一化频率);ylabel(振幅/dB);stext=N=int2str(Nwin);%给出标题，指出所用的数据个数 title(stext);grid on;end 00.51-70-60-50-40-30-20-100归 一 化 频 率振幅/dBN=1000.51-70-60-50-40-30-20-100归 一 化 频 率振幅/dB
27、N=2000.51-80-60-40-200归 一 化 频 率振幅/dBN=5000.51-80-60-40-200归 一 化 频 率振幅/dBN=100 图 7-4 数据长度不同的矩形窗的幅频形状 程序运行结果见图 7-4。显然，随着 N 的增大，主瓣和旁瓣都变窄，但第一旁瓣相对主瓣的幅值下降分贝数相同，第二旁瓣相对第一旁瓣幅值下降的分贝数也相同。这样，随着 N的增大，旁瓣也得到抑制，有力地减少了频谱泄漏，但不能完全消除。减少主瓣宽度和抑制旁瓣是一对矛盾，不可兼得，只能根据不同用途折衷处理。运用窗函数设计数字滤波器 用于信号分析中的窗函数可根据用户的不同要求选择。用于滤波器的窗函数，一般要求
28、窗函数的主瓣宽度窄，以获得较好的过渡带；旁瓣相对值尽可能少，增加通带的平稳度和增大阻带的衰减。基于窗函数的 FIR 数字滤波器设计的算法十分简单，其主要步骤为：(1)对滤波器理想频域幅值响应进行傅立叶逆变换获得理想滤波器的单位脉冲响应hd(n)。一般假定理想低通滤波器的截止频率为c，其幅频特性满足 ccjeH001 (7-19)则根据傅立叶逆变换，单位脉冲响应为：nndenhcnjdccsin21)(,n (7-20)其中，为信号延迟。（2）由性能指标（阻带衰减的分贝数）根据表 7-2 第 3 列的值确定满足阻带衰减的窗函数类型 w(n)。滤波器的阶数越高，滤波器的幅频特性越好，但数据处理的费
29、用也越高，因此像 IIR滤波器一样，FIR 滤波器也要确定满足性能指标的滤波器最小阶数。由前面的讨论（图 7-1）可知，滤波器的主瓣宽度相当于过渡带宽，因此，使过渡带宽近似于窗函数主瓣宽（表 7-2中的第二列）可求得满足性能指标的窗口长度 N，此时，信号延迟为(N-1)/2。（3）求实际滤波器的单位脉冲响应 h(n)：根据 h(n)=hd(n)*w(n)。（4）检验滤波器的性能。可设定一些信号采用 节指出的函数或节所给的函数进行滤波。下面采用实例说明如何根据上面步骤设计 FIR 滤波器。【例 7-4】用窗函数设计一个线性相位 FIR 低通滤波器，并满足性能指标：通带边界的归一化频率 wp=,阻
30、带边界的归一化频率 ws=，阻带衰减不小于 30dB，通带波纹不大于 3dB。假设一个信号，其中 f1=5Hz,f2=20Hz。信号的采样频率为 50Hz。试将原信号与通过滤波器的信号进行比较。由题意，阻带衰减不小于 30dB,根据表 7-2,选取汉宁窗，因为汉宁窗的第一旁瓣相对主瓣衰减为 31dB，满足滤波要求。在窗函数设计法中，要求设计的频率归一化到 0区间内，Nyquist 频率对应于，因此通带和阻带边界频率为和。程序如下%Samp7_4 wp=*pi;ws=*pi;%滤波器边界频率 wdelta=ws-wp;%过渡带宽 N=ceil(8*pi/wdelta)%根据过渡带宽等于表 7-2
31、 中汉宁窗函数主瓣宽求得滤波器所用窗函数的最小长度 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 的采样频率绘出该滤波器的幅频和相频响应 subplot(2,1,1)
32、,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=;t=0:dt:3;%采样间隔和检测信号的时间序列 x=sin(2*pi*f1*t)
33、+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。该例设计通带边界wp=,阻带边界频率ws=，对应于50Hz的采样频率通带边界频率为 fp=Fs/2*Fn
34、ormal=50/2*=,fs=50/2*=,其中 Fs 为采样频率，Fnormal 为归一化频率。由图 7-5 上图可以看到，在小于的频段上，几乎看不到下降，即满足通带波纹不大于 3dB 的要求。在大于的频段上，阻带衰减大于 30dB，满足设计要求。由图 7-5 下图可见，在通带范围内，相位频率为一条直线，表明该滤波器为线性相位。图 7-6给出了滤波器的输入信号和输出信号，输入信号包括 3Hz 和 20Hz 的信号，由图 7-5 可知，20Hz 的信号不能通过该滤波器，通过滤波器后只剩下 3Hz 的信号，输出结果也证明了这一点。但要注意，由于 FIR 滤波器所需的阶数较高，信号延迟(N-1)
35、/2 也较大，输出信号前面有一段直线就是延迟造成的。上述程序显示的 N 取 50 才能满足设计要求。这样相位延迟为(N-1)/2*1/Fs=，可以看到输出信号前面一段直线的距离大约为。验证了 FIR 滤波器相位延迟的理论。在输出信号的前部，有一些小信号，这是截断信号边界所致，后面的部分就没有了这种信号。若采用零相位的 filtfilt 函数（说明见第六章第三节）输出，则可最大限度地减小边界的影响。0510152025-150-100-50050频 率/Hz振幅/dB0510152025-3000-2000-10000频 率/Hz相位/o 图 7-5 例 7-4 所设计滤波器的幅频响应（上图）和
36、相频响应（下图）00.511.522.53-2-1012输 入 信 号00.511.522.53-2-1012时 间/s输 出 信 号 图 7-6 例 7-4 所设计滤波器的输入和输出信号 标准型 FIR 滤波器 节给出了运用理想脉冲响应与窗函数乘积的方法给出了滤波器传递函数的设计方法。其实 MATLAB 已将上述方法复合成一个函数，提供基于上述原理设计标准型 FIR 滤波器的工具函数。fir1 就是采用经典窗函数法设计线性相位 FIR 数字滤波器的函数，且具有标准低通，带通，高通,带阻等类型。函数调用格式为：b=fir1(n,wn,ftype,window)式中，n 为 FIR 滤波器的阶数
37、，对于高通,带阻滤波器，n 需取偶数；wn 为滤波器截止频率，范围为 01（归一化频率）。对于带通,带阻滤波器，wn=w1,w2(w1w2)；对于多带滤波器，如 wn=w1，w2，w3,w4,频率分段为：0ww1，w1ww2，w2ww3，w3ww4。ftype为滤波器的类型：缺省时为低通或带通滤波器；high为高通滤波器；stop为带阻滤波器，DC-1为第一频带为通带的多带滤波器；DC-0为第一频带为阻带的多带滤波器。window 为窗函数列向量，其长度为 n+1。缺省时，自动取哈明窗。MATLAB 提供的窗函数有 boxcar、hanning、hamming、bartlett、blackma
38、n、kaiser、chebwin，调用方式见上节。b 为 FIR 滤波器系数向量，长度为 n+1。FIR 滤波器的传递函数具有下列形式：nznbzbzbbzb)1()3()2()1()(21 (7-21)用函数 fir1 设计的 FIR 滤波器的群延迟为 n/2。考虑到 n 阶滤波器系数个数为 N，即n+1，这里的延迟与前面所讲的(N-1)/2 的延迟一致。注意这里的滤波器的最小阶数比窗函数的长度少 1。【例 7-5】用窗函数设计一个线性相位 FIR 低通滤波器，技术指标同上节例 7-4。%Samp7_5 wp=*pi;ws=*pi;%滤波器的边界频率 wdelta=ws-wp;%过渡带宽度
39、N=ceil(8*pi/wdelta);%求解滤波器的最小阶数，根据汉宁窗主瓣宽 Wn=+*pi/2;%截止频率取通带和阻带边界频率的中点 b=fir1(N,Wn/pi,hanning(N+1);%设计 FIR 滤波器,注意 fir1 要求输入归一化频率 H,f=freqz(b,1,512,50);%采用 50Hz 的采样频率求出频率响应 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
40、);ylabel(相位/o);grid on;程序运行与上例中的图 7-5 一致。【例 7-6】设计一个 48 阶 FIR 带通滤波器，通带边界的归一化频率为和。假设一个信号，其中含有 f1=1Hz，f2=10Hz,f3=20Hz，三种频率成分信号的采样频率为 50Hz。试将原信号与通过滤波器的信号进行比较。%Samp7_6 wp=;N=48;%通带边界频率（归一化频率）和滤波器阶数 Fs=50;b=fir1(N,wp);%设计 FIR 带通滤波器 figure(1)H,f=freqz(b,1,512,Fs);%以 50Hz 为采样频率求出滤波器频率响应 subplot(2,1,1),plot
41、(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;f1=1;f2=10;f3=20;%输入信号的三种频率成分 t=0:1/50:3;%时间序列 x=sin(2*pi*f1*t)+*cos(2*pi*f2*t)+*sin(2*pi*f3*t);%输入信号%y=filter(b,1,x);%可以采用过滤器进行滤波 y=fftfilt(b,x);%采用 fftfilt 对输入信号滤波 f
42、igure(2)subplot(2,1,1),plot(t,x),title(输入信号)%绘出输入信号波形 subplot(2,1,2),plot(t,y)%绘出输出信号波形 hold on;plot(N/2*ones(1,2),ylim,r)%绘制延迟到的时刻 title(输出信号),xlabel(时间/s)程序运行结果见图 7-7 和图 7-8。通带归一化频率和对应于采样频率为 50Hz 的通带范围为和(采用 6-19 式计算)。由图 7-7 上图可见满足这一设计要求。在这个频带范围内的相位满足线性相位，符合 FIR 滤波器的一般特点。图 7-8 为检测滤波器的输入信号和输出信号。输入信号
43、中含有 1Hz,10Hz 和 20Hz 的信号。根据图 7-7 上图可知，1Hz 和 20Hz 的频率在阻带范围内，不能通过该滤波器，只有 10Hz 的信号可以通过该滤波器，输出信号表明了这一点。滤波器的相位延迟根据 N/2*=得到，输出信号前面的直线部分大体为这个时间延迟,另外滤波后周期为 10Hz 的信号相位（红线开始部分），跟滤波前的相位（信号开始部分）也一致，说明通过该滤波器滤波后没有改变信号的相位。0510152025-100-50050频 率/Hz振幅/dB0510152025-2000-1500-1000-5000500频 率/Hz相位/o 图 7-7 例 7-6 设计滤波器的幅
44、频特性（上图）和相频特性（下图）00.511.522.53-2-1.5-1-0.500.511.5输 入 信 号00.511.522.53-0.4-0.200.20.40.6输 出 信 号时 间/s 图 7-8 例 7-6 滤波器的输入信号和输出信号 【例 7-7】FIR 低通滤波器阶数为 40，截止频率为 200Hz，采样频率为 Fs=1000Hz。试设计此滤波器并对信号 x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t)滤波，f1=50Hz,f2=250Hz，选取滤波器输出的第 81 个采样点到第 241 个采样点之间的信号并与对应的输入信号进行比较。由于采样频率为 10
45、00Hz，所以该滤波器的归一化频率的 1 对应于 Nyquist 频率 500Hz，因此归一化截止频率为 200/500(参看（6-19）式)。%Samp7_7 clf;N=1000;Fs=1000;%数据总数和采样频率 fc=200;n=0:N-1;t=n/Fs;%时间序列 f1=50;f2=250;x=sin(2*pi*f1*t)+sin(2*pi*f2*t);%输入信号 b=fir1(40,fc*2/Fs);%设计 40 阶的低通滤波器，归一化截止频率据 6-19 式 yfft=fftfilt(b,x,256);%对数据进行滤波 n1=81:241;t1=t(n1);%选择采样点间隔 x
46、1=x(n1);%与采样点对应的输入信号 subplot(2,1,1);plot(t1,x1);grid on;%绘制输入信号 title(输入信号);n2=n1-40/2;t2=t(n2);%输出信号，扣除了相位延迟 N/2 y2=yfft(n2);subplot(2,1,2);plot(t2,y2);%绘制输出信号 title(输出信号);grid on;xlabel(时间/s)程序输出结果见图 7-9。可见经过滤波器的滤波，完全滤去了 250Hz 的高频成分，只剩下 50Hz 的低频成分。0.080.10.120.140.160.180.20.220.24-2-1012输 入 信 号0.
47、060.080.10.120.140.160.180.20.22-1-0.500.51输 出 信 号时 间/s 图 7-9 例 7-7设计滤波器的输入信号和输出信号 【例 7-8】设计采样频率为 1000Hz,阻带频率从 100Hz200Hz 的 100 阶的带阻 FIR 滤波器。对信号 x(t)=sin(2*pi*f1*t)+sin(2*pi*f2*t)滤波，f1=50Hz，f2=150Hz，并与对应的输入信号进行比较。%Samp7_8 clf;N=300;Fs=1000;%数据总数和采样频率 Order=100;%滤波器阶数 n=0:N-1;t=n/Fs;%时间序列 wn=100 200/
48、(Fs/2);%边界频率转换为归一化频率，据 6-19 式 b=fir1(Order,wn,stop);%设计 100 阶的带阻滤波器 figure(1)H,f=freqz(b,1,512,Fs);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;f1=50;f2=150;%输入信号频率 x=sin(2*pi*f1*t)+sin(2*pi*f2
49、*t);%输入信号 y=fftfilt(b,x,256);%对数据进行滤波 figure(2)subplot(2,1,1);plot(t,x);grid on;%绘制输入信号 title(输入信号);subplot(2,1,2);plot(t,y);%绘制输出信号 hold on;plot(Order/2/Fs*ones(1,2),ylim,r)%绘制延迟到的时刻 title(输出信号);grid on;xlabel(时间/s)程序输出结果见图 7-10 和图 7-11。由图 7-10 上图可见，阻带范围为 100200Hz，完全符合设计要求。在通带的相位是线性的。由图 7-11 可见，滤波器
50、滤除了 150Hz（在阻带范围内）的信号，保留了 50Hz 的信号。相位延迟 100/2/Fs=，与图形一致。050100150200250300350400450500-100-50050频 率/Hz振幅/dB050100150200250300350400450500-8000-6000-4000-20000频 率/Hz相位/o 图 7-10 例 7-8 设计带阻滤波器的幅频（上图）和相频特性（下图）00.050.10.150.20.250.30.35-2-1012输 入 信 号00.050.10.150.20.250.30.35-1.5-1-0.500.511.5输 出 信 号时 间/s