FIR滤波器设计

《FIR滤波器设计》由会员分享，可在线阅读，更多相关《FIR滤波器设计（51页珍藏版）》请在装配图网上搜索。

1、word第7章 FIR滤波器设计第六章我们介绍了无限冲激响应IIR滤波器的设计方法。其中最常用的由模拟滤波器转换为数字滤波器的方法为双线性变换法，因为这种方法无混叠效应，效果较好。但通过前面的例子我们看到，IIR数字滤波器相位特性不好非线性，如图 6-11、图6-13、图6-15 ，也不易控制。然而在现代信号处理中，例如图像处理、数据传输、雷达接收以与一些要求较高的系统中对相位特性要求较为严格，这种滤波器就无能为力了。改善相位特性的方法是采用有限冲激响应滤波器。本章首先对FIR滤波器原理与其使用函数作根本介绍，然后重点介绍窗函数法设计FIR滤波器，并对最优滤波器设计函数进展介绍。7.1 FIR

2、滤波器原理概述与滤波函数7.1.1 FIR滤波器原理与设计方法分类根据第 6 章对数字滤波器的介绍，我们知道FIR滤波器的传递函数为： (7-1)可得FIR滤波器的系统差分方程为：因此，FIR滤波器又称为卷积滤波器。根据第 4 章中所描述的系统频率响应，FIR滤波器的频率响应表达式为： 7-2信号通过FIR滤波器不失真条件与6-6式所描述的一样，即滤波器在通带内具有恒定的幅频特性和线性相位特性。理论上可以证明这里从略：当FIR滤波器的系数满足如下中心对称条件： 7-3时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不

3、变的。对于一个 N 阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。本章主要介绍的FIR数字滤波器设计方法与 MATLAB 信号处理工具箱提供的了FIR数字滤波器设计函数，见表7-1。由于篇幅所限，本章我们主要介绍窗函数法和最优化设计方法。表7-1 FIR滤波器设计的主要方法函数设计方法说明工具函数窗函数法理想滤波器加窗处理fir1(单频带) ， fir2(多频带) ， kaiserord最优化设计平方误差最小化逼近理想幅频响应或Park-McClellan 算法产生等波纹滤波器firls ，

4、remez,remezord约束最小二乘逼近在满足最大误差限制条件下使整个频带平方误差最小化fircls,fircls1升余弦函数具有光滑、正弦过渡带的低通滤波器设计Fircos7.1.2 FIR数字滤波器滤波函数相对于IIR 滤波器的滤波函数，FIR数字滤波器滤波函数除了dimpulse和dstep仅适用于IIR滤波器外，其他各种函数可直接应用于FIR滤波器，只是输入的分母多项式向量a=1。另外，MATLAB还提供了一个函数fftfilt,该函数利用效率高的基于FFT算法实现对数据的滤波，该函数只适用于FIR滤波器，调用形式为：y=fftfilt(b,x,n)式中，b为FIR滤波器的系数向量

5、；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、带宽度X围内。频谱泄漏是不可防止的，但要尽量减小。上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，能否改变窗的形状而达到这个目的呢？回答是肯定的。其实数字信号处理的前驱者们设计了不同于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足不同的要求。为此，用窗函数法设计FIR数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度N和窗函数w(n)。下面我们介绍窗函数。7.2.2 MATLAB信号处理中提供的窗函数1矩形窗：前面分析中所用的矩形窗可用下面函数来实现w=boxcar (N)，N 为窗的长度以下函数与此同，w为返回的窗函数序列

11、。2汉宁窗：w=hanning(N)汉宁窗的表达式为： (7-10)3哈明窗：w=hamming(N)哈明Hamming窗的表达式为： (7-11)4Bartlett窗：w=bartlett(N)Bartlett 窗的表达式为：当 N 为奇数时， (7-12)当 N 为偶数时， (7-13)5 Blackman 窗：w= blackman(N)Blackman 窗的表达式为： (7-14)Blackman 窗比其他一样尺寸窗 (哈明Hamming窗,汉宁Hanning窗) 具有主瓣较宽和旁瓣泄漏较小的特点。6三角窗：w=triang(N)三角窗的表达式为：当 N 为奇数时， (7-15)当 N

12、 为偶数时， (7-16)三角窗和Bartlett窗十分类似。三角窗的两端值不为零，而Bartlett窗如此为零，这一点可从例7-1中看出。7Kaiser窗：w=kaiser(n,beta)其中，beta是Kaiser窗参数,影响窗旁瓣幅值的衰减率。Kaiser窗表达式： (7-17)式中， I0.是修正过的零阶 Bessel 函数。Kaiser窗用于滤波器设计时，假如旁瓣幅值为，如此 7-18 8 Chebyshev窗：w=chebwin(n,r)式中, r 是窗口的旁瓣幅值在主瓣以下的分贝数。Chebyshev窗具有主瓣宽度最小，而旁瓣等高,高度可调整的特点。下面我们在MATLAB观看各种

13、窗函数的形状和频率域图象来验证上述所讲特点。【例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=bartlett(Nwin);%Bbartlett窗 stext=Bartelett窗; e

14、nd 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=triang(Nwin);%三角窗 stext=三角窗; case 3 w=kaiser(

15、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程序运行结果见图 7-2 。可以看到各种窗函数的形状。图 7-2 各种窗函数的时间域形状【例 7-

16、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(Nwin);%Bartlett窗 stext=Bartelett窗; end y,f=freqz

17、(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);%Blackman 窗 stext=Blackman窗; case 2 w=triang(Nwi

18、n);%三角窗 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* log10(mag/max(mag); %绘制幅频响应 xlabel(归一化频率);yl

19、abel(振幅/dB); title(stext);grid on;end程序运行结果见图 7-3 。可以看到各种窗函数的幅频形状。对照该图可知这些窗函数具有上面所分析的窗函数的特征。图 7-3 各种窗函数的幅频形状由图 7-3 可见，各种窗函数都具有明显的主瓣Mainlobe和旁瓣Sidelobe。主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值第一旁瓣衰减为 13 dB；Blackman窗具有最大的旁瓣衰减，但也具有最宽的主瓣宽度。不同窗函数在这两方面的特点是不同的，因此应根据具体的问题进展选择。通常来讲，哈明Hamming窗和汉宁Hannin

20、g窗的主瓣，具有较小的旁瓣和较大的衰减速度，是较为常用的窗函数。表7-2总结了各种窗函数主瓣和旁瓣的特征理论分析可参考其他的数字信号处理教材，大家可对照窗函数的幅频形状图7-3认真理解体会。表7-2 各种窗函数的特点窗函数主瓣宽第一旁瓣相对主瓣衰减分贝矩形窗-13汉宁Hanning窗-31哈明Hamming窗-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,int2str(ii);%指定绘图位置 subplot(posplot); plo

22、t (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的增大，旁瓣也得到抑制，有力地减少了频谱泄漏，但不能完全消除。减少主瓣宽度和抑制旁瓣是一对矛盾，不可兼得，只能根据不同用途折衷处理。7.2.3

23、 运用窗函数设计数字滤波器用于信号分析中的窗函数可根据用户的不同要求选择。用于滤波器的窗函数，一般要求窗函数的主瓣宽度窄，以获得较好的过渡带；旁瓣相对值尽可能少，增加通带的平稳度和增大阻带的衰减。基于窗函数的FIR数字滤波器设计的算法十分简单，其主要步骤为：(1)对滤波器理想频域幅值响应进展傅立叶逆变换获得理想滤波器的单位脉冲响应hd(n)。一般假定理想低通滤波器的截止频率为，其幅频特性满足 (7-19)如此根据傅立叶逆变换，单位脉冲响应为： (7-20)其中，为信号延迟。2由性能指标阻带衰减的分贝数根据表7-2第3列的值确定满足阻带衰减的窗函数类型w(n)。滤波器的阶数越高，滤波器的幅频特性

24、越好，但数据处理的费用也越高，因此像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】 用窗函数设计一个线性相位FIR低通滤波器，并满足性能指标：通带边界的归一化频率wp=0.5,阻带边界的归一化

25、频率ws=0.66，阻带衰减不小于30dB，通带波纹不大于3dB。假设一个信号，其中f1=5Hz,f2=20Hz。信号的采样频率为50Hz。试将原信号与通过滤波器的信号进展比拟。由题意，阻带衰减不小于30dB,根据表7-2,选取汉宁hanning窗，因为汉宁Hanning窗的第一旁瓣相对主瓣衰减为31dB，满足滤波要求。在窗函数设计法中，要求设计的频率归一化到0区间内，Nyquist频率对应于。程序如下%Samp7_4wp=0.5*pi;ws=0.66*pi;%滤波器边界频率wdelta=ws-wp; %过渡带宽N=ceil(8*pi/wdelta) %根据过渡带宽等于表 7-2中汉宁Hann

26、ing窗函数主瓣宽求得滤波器所用窗函数的最小长度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),plot(f,20*l

27、og10(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; %采样间隔和检测信号的时间序列x=sin(2*pi*f1* t)+cos(2*

28、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=0.5,阻带边界频率ws=0.66，对应于50Hz的采样频率通带边界频率为fp=Fs/2*Fno

29、rmal=50/2*0.5=12.5Hz,Hz, 其中Fs为采样频率，Fnormal为归一化频率。由图7-5上图可以看到，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不大于3dB的要求。在大于16.5Hz的频段上，阻带衰减大于30dB，满足设计要求。由图7-5如下图可见，在通带X围内，相位频率为一条直线，明确该滤波器为线性相位。图7-6给出了滤波器的输入信号和输出信号，输入信号包括3Hz和20Hz的信号，由图7-5可知，20Hz的信号不能通过该滤波器，通过滤波器后只剩下3Hz的信号，输出结果也证明了这一点。但要注意，由于FIR滤波器所需的阶数较高，信号延迟(N-1)/2也较大，输

30、出信号前面有一段直线就是延迟造成的。上述程序显示的N取50才能满足设计要求。这样相位延迟为s。验证了FIR滤波器相位延迟的理论。在输出信号的前部，有一些小信号，这是截断信号边界所致，后面的局部就没有了这种信号。假如采用零相位的filtfilt函数说明见第六章第三节输出，如此可最大限度地减小边界的影响。图 7-5 例7-4所设计滤波器的幅频响应上图和相频响应如下图图7-6 例7-4所设计滤波器的输入和输出信号7.2.4标准型FIR滤波器节给出了运用理想脉冲响应与窗函数乘积的方法给出了滤波器传递函数的设计方法。其实MATLAB已将上述方法复合成一个函数，提供基于上述原理设计标准型FIR滤波器的工具

31、函数。fir1就是采用经典窗函数法设计线性相位FIR数字滤波器的函数，且具有标准低通，带通，高通,带阻等类型。函数调用格式为：b=fir1(n,wn,ftype,window)式中，n为FIR滤波器的阶数，对于高通,带阻滤波器，n需取偶数；wn为滤波器截止频率，X围为01归一化频率。对于带通,带阻滤波器，wn=w1,w2(w1w2)；对于多带滤波器，如wn=w1，w2，w3,w4,频率分段为：0ww1，w1ww2，w2ww3，w3ww4。 ftype为滤波器的类型：缺省时为低通或带通滤波器；high为高通滤波器；stop为带阻滤波器，DC-1为第一频带为通带的多带滤波器；DC-0为第一频带为阻

32、带的多带滤波器。window为窗函数列向量，其长度为n+1。缺省时，自动取Hamming哈明窗。MATLAB提供的窗函数有boxcar、hanning、hamming、bartlett、blackman、kaiser、chebwin，调用方式见上节。b为FIR滤波器系数向量，长度为n+1。FIR滤波器的传递函数具有如下形式： (7-21)用函数fir1设计的FIR滤波器的群延迟为n/2。考虑到n阶滤波器系数个数为N，即n+1，这里的延迟与前面所讲的(N-1)/2的延迟一致。注意这里的滤波器的最小阶数比窗函数的长度少1。【例7-5】 用窗函数设计一个线性相位FIR低通滤波器，技术指标同上节例7-

33、4。%Samp7_5wp=0.5*pi;ws=0.66*pi;%滤波器的边界频率wdelta=ws-wp;%过渡带宽度N=ceil(8* pi/wdelta);%求解滤波器的最小阶数，根据汉宁Hanning 窗主瓣宽Wn=(0.5+0.66)*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(振幅/

34、dB);grid on;subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)xlabel(频率/Hz);ylabel(相位/o);grid on;程序运行与上例中的图7-5一致。【例 7-6】 设计一个48阶FIR带通滤波器，通带边界的归一化频率为0.35和0.65。假设一个信号，其中含有f1=1Hz，f2=10Hz,f3=20Hz，三种频率成分信号的采样频率为50Hz。试将原信号与通过滤波器的信号进展比拟。%Samp7_6wp=0.35 0.65;N=48;%通带边界频率归一化频率和滤波器阶数Fs=50;b=fir1(N,wp);%设计FIR带通滤波器fi

35、gure(1)H,f=freqz(b,1,512,Fs);%以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);ylabel(相位/o);grid on;f1=1;f2=10;f3=20;%输入信号的三种频率成分t=0:1/50:3;%时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t)

36、;%输入信号%y=filter(b,1,x);%可以采用过滤器进展滤波y=fftfilt(b,x);%采用fftfilt对输入信号滤波figure(2)subplot(2,1,1), plot(t,x),title(输入信号)%绘出输入信号波形subplot(2,1,2),plot(t,y)%绘出输出信号波形hold on;plot(N/2*0.02*ones(1,2),ylim, r) %绘制延迟到的时刻title(输出信号),xlabel(时间/s)程序运行结果见图7-7和图7-8。通带归一化频率对应于采样频率为50Hz的通带XHz(采用6-19式计算)。由图7-7上图可见满足这一设计要求

37、。在这个频带X围内的相位满足线性相位，符合FIR滤波器的一般特点。图7-8为检测滤波器的输入信号和输出信号。输入信号中含有1Hz,10Hz和20Hz的信号。根据图7-7上图可知，1Hz和20Hz的频率在阻带X围内，不能通过该滤波器，只有10Hz的信号可以通过该滤波器，输出信号明确了这一点。滤波器的相位延迟根据N/2*0.02s=0.48s得到，输出信号前面的直线局部大体为这个时间延迟,另外滤波后周期为10Hz的信号相位红线开始局部，跟滤波前的相位信号开始局部也一致，说明通过该滤波器滤波后没有改变信号的相位。图 7-7 例 7-6 设计滤波器的幅频特性上图和相频特性如下图图 7-8 例 7-6

38、滤波器的输入信号和输出信号【例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个采样点之间的信号并与对应的输入信号进展比拟。由于采样频率为1000Hz，所以该滤波器的归一化频率的1对应于Nyquist频率500Hz，因此归一化截止频率为200/500(参看6-19式)。%Samp7_7clf;N=1000;Fs=1000;%数据总数和采样频率fc=200;n=0:N-1;t=n/Fs;%

39、时间序列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);%选择采样点间隔x1=x(n1);%与采样点对应的输入信号subplot(2,1,1);plot(t1,x1); grid on;%绘制输入信号title(输入信号);n2=n1-40/2;t2=t(n2); %输出信号，扣除了相位延迟N/2y2=yfft(n2);subplot(2,1,2);plot

40、(t2,y2);%绘制输出信号title(输出信号);grid on; xlabel(时间/s)程序输出结果见图7-9。可见经过滤波器的滤波，完全滤去了250Hz的高频成分，只剩下50Hz的低频成分。图 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_8clf;N=300;Fs=1000;%数据总数和采样频率Order=100; %滤波器阶数n

41、=0:N-1;t=n/Fs;%时间序列wn=100 200/(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=

42、sin(2*pi*f1*t)+sin(2*pi*f2*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上图可见，阻带X围为100200Hz，完全符合设计要求。在通带的

43、相位是线性的。由图7-11可见，滤波器滤除了150Hz在阻带X围内的信号，保存了50Hz的信号。相位延迟100/2/Fs=0.05s，与图形一致。图 7-10 例7-8设计带阻滤波器的幅频上图和相频特性如下图图 7-11例7-8设计滤波器的输入和输出信号该小节只给出了FIR低通，带通和带阻滤波器的例子，请大家在课下自己设计高通，第一频带为通频带和第一频带为阻频带的多带滤波器的例子，以加深对该函数的理解。7.2.5 多频带FIR滤波器除了设计标准型FIR滤波器外，MATLAB信号处理工具箱还提供另一种基于窗函数滤波器设计的工具函数fir2，用于设计具有任意形状频率响应的FIR滤波器，其调用格式为

44、：b=fir2n，f,m,npt,window式中，n为滤波器的阶数；f和m分别为滤波器期望幅频响应的频率向量和幅值向量，取值X围为01(归一化频率)。m，f具有一样的长度,window为窗函数，得到列向量，长度必须为n+1。缺省时自动取hamming哈明窗；npt为对频率响应进展内插的点数，缺省时为512。b为FIR滤波器系数向量，长度为n+1，滤波器具有与7-21式一样的形式。【例 7-9】 用窗函数设计一个多频带的FIR滤波器，滤波器阶数分别为10和100，滤波器的特性同上一章例6-12，即f=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0,m=0 0

45、 1 1 0 0 1 1 1 0 0,比拟理想和实际滤波器的幅频响应。假设一个信号，其中f1=12Hz,f2=36Hz。信号的采样频率为100Hz。试将原信号与通过滤波器的信号进展比拟。%Samp7_9clff=0:0.1:1;%归一化频率点数m=0 0 1 1 0 0 1 1 1 0 0;%幅频特性值Order=10;% 滤波器的阶数b=fir2(Order,f,m,hamming(Order+1);%设计滤波器h,w=freqz(b,1,128);%计算滤波器的频率响应subplot(2,1,1)plot(f,m,w/pi,abs(h),r:)%绘制理想幅频响应和设计的滤波器幅频响应leg

46、end(理想特性, 实际设计) %给出图例title(Order=10);xlabel(归一化频率);ylabel(振幅);Order=100;b=fir2(Order,f,m,hamming(Order+1);%设计阶数为100的滤波器h,w=freqz(b,1,128);%计算滤波器的频率响应subplot(2,1,2),plot(f,m,w/pi,abs(h),r:);%绘制理想幅频响应和设计的幅频响应ylim(0 1)legend(理想特性, 实际设计) %给出图例title(Order=100);xlabel(归一化频率);ylabel(振幅);f1=12;f2=36;% 输入信号的

47、两种频率成分t=0:1/100:2;% 时间序列x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);%输入信号y=fftfilt(b,x);%对输入信号进展滤波figure(2)subplot(2,1,1), plot(t,x),title(输入信号)%绘制输入信号subplot(2,1,2),plot(t,y)%绘制输出信号hold on;plot(Order/2/100*ones(1,2),ylim, r) %绘制延迟到的时刻title(输出信号),xlabel(时间/s)程序输出结果见图7-12和图7-13。由该例可知，只有取100阶时，实际滤波器的幅频响应才逼近理想

48、滤波器的幅频响应。与第6章例6-12的输出比拟可知，一样性能的FIR滤波器阶数比IIR滤波器要高得多。另外，该例中，信号含有两种频率成分，15Hz。同理，归一化频率中的40Hz。可见12Hz和36Hz的波均在通带的X围内，因此均可通过。该例的输出结果与例6-12的输出结果比拟可知，FIR滤波器的相位延迟比IIR滤波器的相位延迟大得多。图 7-12 例7-9所设计滤波器的幅频特性上图：阶数为10；如下图：阶数为100图 7-13 例7-9所设计滤波器的输入和输出信号最优FIR滤波器设计MATLAB信号处理工具箱提供了比基于窗函数法FIR滤波器设计工具函数fir1和fir2更为通用的函数firls

49、和remez。它们采用不同的优化方法设计最优的标准多频带FIR数字滤波器。函数remez实现Park-McClellan算法，这种算法利用Remez交换算法和Chebyshev近似理论来设计滤波器，使实际频率响应拟合期望频率响应达到最优。从实际和理想频响之间最大误差最小化的观点来看，函数remez设计的滤波器是最优的，因此，又称之为最优滤波器。在频率域内，滤波器呈现等波纹特点，因此又称之为等波纹滤波器。Park-McClellan滤波器设计方法是FIR滤波器设计中最流行的，应用最广的设计方法。函数firls和remez的调用格式语法规如此一样，只是优化算法不同。7.3.1 根本形式最优滤波器函

50、数根本调用格式为：b=firls(n,f,a)b=remez(n,f,a)式中，n为滤波器阶数；f为滤波器期望频率特性归一化频率向量，X围为01，为递增向量，允许定义重复频率点；a为滤波器期望频率特性的幅值向量，向量a和f必须为同长度,且为偶数；b为返回的滤波器系数，长度为n+1，且具有偶对称的关系，即b(k)=b(n+2-k)。假如滤波器的阶数为奇数，如此在Nyquist频率处对应于归一化频率1，幅频响应必须为0。滤波器的阶数为偶数如此无此限制。【例7-10】分别用firls和函数remez设计一个20阶FIR低通滤波器，通带边界频率为0.4,幅值为1，阻带边界频率为0.5，幅值为零。输入一

51、个采样频率为50Hz,频率为5Hz和15Hz的合成振动，比拟运用两种设计方法的输出的差异。%Samp7_10clf;n=20; %滤波器阶数f=0 0.4 0.5 1; %频率向量a=1 1 0 0; %振幅向量b=firls(n,f,a); %采用firls设计滤波器h,w1=freqz(b); %计算其频率响应bb=remez(n,f,a); %采用remez设计滤波器hh,w2=freqz(bb); %计算滤波器的频率响应figure(1)plot(w1/pi,abs(h),w2/pi,abs(hh),r:); %绘制滤波器的幅频响应xlabel(归一化频率);ylabel(振幅);le

52、gend(firls,remez); %给出图例grid on;figure(2)fs=50;t=0:1/fs:2; %采样频率和时间序列f1=5;f2=15; %输入信号频率x1=sin(2*pi*f1*t)+8.*cos(2*pi*f2*t);subplot(2,1,1),plot(t,x1)title(原始信号)y1=filter(b,1,x1);y2=filter(bb,1,x1);subplot(2,1,2),plot(t,y1,t,y2,r:)legend(firls,remez); %给出图例title(输出信号)xlabel(时间/s)图 7-14 例 7-10 所设计滤波器的

53、幅频响应图 7-15 例 7-10 所设计滤波器的输入和输出信号的比拟程序运行结果见图7-14和图7-15。比拟两种设计方法所设计的滤波器幅频响应可见，用firls设计的滤波器在整个频率X围内包括通带和阻带均具有较好的响应，但理想频率响应和实际频率响应的误差在带区内不均匀，且在边界频率处误差较大；而函数remez设计的滤波器在通带内内更接近于理想频响。由图7-15可见，输入5Hz和15Hz的两个频率的振动信号，5Hz的信号对应的采样频率为50Hz的归一化频率为0.25/(50/2)，15Hz的信号对应的归一化频率为0.615/(50/2)，对应于该滤波器的幅频响应，归一化频率为0.2的振动能通

54、过滤波器，而归一化为0.6的振动不能通过该滤波器，输出结果也看到这一点，只有5Hz的振动通过了滤波器。但通过两种类型的滤波器滤波，输出结果略有不同，这是由于两种滤波器的幅频特性的略微不同造成的。运用firls函数设计幅频响应在归一化为0.2的频率处略小于remez函数的设计，因此运用firls函数设计滤波器的输出的归一化频率为0.2的振幅也略小于remez函数设计的输出。运用firls函数设计幅频响应在归一化为0.6的频率处比remez函数设计的振幅响应衰减更大，因此运用firls函数设计滤波器的输出更为平滑，即高频成分更小，而remez函数设计的滤波器的输出具有较多的高频成分，使得输出结果不

55、对称。这完全可以根据滤波器的幅频响应分析出来。函数firls和remez可用于设计低通，高通,带通和带阻等一般类型的滤波器，这可由函数中给定的理想幅频响应的频率向量f和幅值向量确定。如设计一个带通滤波器，幅频响应定义为:f=0 0.3 0.4 0.7 0.8 1， a=0 0 1 1 0 0;如此，该理想滤波器幅频响应定义为：阻带频率00.8。设计一个高通滤波器，如理想幅频响应向量对的给定形式为： f=0 0.7 0.8 1,a=0 0 1 1,如此，该理想滤波器的幅频响应定义为：阻带频率00.8。设计一个带阻滤波器，如理想幅频响应向量对的给定形式为：f=0 0.3 0.4 0.5 0.6 1

56、,a=1 1 0 0 1 1，如此，该0.5，通带频率：00.3,0.6此外，函数firls和remez还可以设计多带滤波器，下面给出一个例子。【例7-11】用函数firls和remez设计一个50阶多通带滤波器，滤波器理想频率响应对为：f=0 0.1 0.15 0.25 0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1,a=1 1 0 0 1 1 0 0 1 1 0 0 1 1;将设计的滤波器的幅频响应和理想滤波器幅频响应进展比拟，绘制remez函数设计滤波器的脉冲响应。%Samp7_11clf;n=50;% 滤波器的阶数f=0 0.1 0.15 0.25

57、0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1;% 频率向量a=1 1 0 0 1 1 0 0 1 1 0 0 1 1;%振幅向量b=firls(n, f,a);% 采用 firls 设计滤波器h,w1=freqz(b);% 计算滤波器的频率响应bb=remez(n, f,a);%采用 remez 设计滤波器hh,w2=freqz(bb);%计算滤波器的频率响应figure(1)plot (w1/pi,abs(h),b.,w2/pi, abs(hh),r:,f,a,m);%绘滤波器幅频响应xlabel(归一化频率);ylabel(振幅);legend(fir

58、ls,remez,理想特性);%给出图例figure(2)impz(bb,1),title (脉冲响应)%给出滤波器的脉冲响应xlabel(样本数);ylabel(幅度)程序的输出结果为图7-16和图7-17。可见，firls所设计的滤波器的通带和阻带具有较小波纹,但在整个频带内不一致。而remez函数设计的滤波器具有较大的通带和阻带波纹，但在整个频带内较为一致。图7-17给出了remez函数设计滤波器的脉冲响应，且具有偶对称的关系，即b(k)=b(n+2-k)。大家还可以观看firls函数设计滤波器的脉冲响应是否具有偶对称的关系。图 7-16例 7-11所设计滤波器的幅频响应与理想滤波器幅频响应的比拟图 7-17 例7-11所设计滤波器的脉冲响应函数firls和remez还能设计具有任意线性过渡带连接阻带和通带，使过渡带具有更广阔平滑的过渡区间。如理想幅频响应按如下频响对给出：f