7章++FIR滤波器设计

第7章 FIR滤波器设计

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

7.1 FIR滤波器原理概述及滤波函数

7.1.1 FIR滤波器原理及设计方法分类

根据第 6 章对数字滤波器的介绍，我们知道FIR滤波器的传递函数为：

X(z)n?0可得FIR滤波器的系统差分方程为：

y(n)?b(0)x(n)?b(1)x(n?1)???b(N?1)x(n?N?1)N?1 ??b(m)x(n?m)?b(n)?x(n)m?0 H(z)?Y(z)??h(n)zN?1?n (7-1)

因此，FIR滤波器又称为卷积滤波器。根据第 4 章中所描述的系统频率响应，FIR滤波器的频率响应表达式为：

He????b(n)ej?n?0N?1?jn? （7-2）

信号通过FIR滤波器不失真条件与（6-6）式所描述的相同，即滤波器在通带内具有恒

定的幅频特性和线性相位特性。理论上可以证明（这里从略）：当FIR滤波器的系数满足下列中心对称条件：

b(n)?b(N?1?n)或b(n)??b(N?1?n) （7-3）

时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个 N 阶的线性相位FIR滤波器，群延迟为常数，即滤波后的信号简单地延迟常数个时间步长。这一特性使通带频率内信号通过滤波器后仍保持原有波形形状而无相位失真。

本章主要介绍的FIR数字滤波器设计方法及 MATLAB 信号处理工具箱提供的FIR数字滤波器设计函数，见表7-1。由于篇幅所限，本章我们主要介绍窗函数法和最优化设计方法。

表7-1 FIR滤波器设计的主要方法

函数设计方法 说明 理想滤波器加窗处理 平方误差最小化逼近理工具函数 fir1(单频带) ， fir2(多频带) ， kaiserord firls ， remez,remezord 窗函数法 最优化设计 想幅频响应或Park-McClellan 算法产生等波纹滤波器 约束最小二乘逼近 在满足最大误差限制条件下使整个频带平方误差最小化 具有光滑、正弦过渡带的低通滤波器设计 fircls,fircls1 升余弦函数 Fircos 7.1.2 FIR数字滤波器滤波函数

相对于IIR 滤波器的滤波函数，FIR数字滤波器滤波函数除了dimpulse和dstep仅适用于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=filter(b,a,x)。

7.2 FIR滤波器的窗函数设计

7.2.1 窗函数的基本原理

FIR滤波器设计的主要任务是根据给定的性能指标确定滤波器的系数b，即系统单位脉冲序列h(n)，它是一个有限长序列。

FIR滤波器的理想频率响应，可写成复数形式的Fourier级数形式：

Hde????h?n?ej?dn?????j?n (7-4)

式中，hd(n)是对应的单位脉冲响应序列。这说明滤波器的频率响应和单位脉冲响应互为Fourier变换对。因此其单位脉冲响应可由下式求得，

hd?n??求得序列hd?n?后，通过z变换，可得到Hd?z?

H?e??e?d? (7-5) ?2???jjn?d1 Hd(z)?(7-6)

n????h?d(n)z?n

注意，这里hd?n?为无限长序列，因此Hd?z?是物理上不可实现的。如何变成物理上可实现呢？一个自然的想法是只取其中的某些项，即只截取hd?n?中的一部分，比如n=0,?,N-1，N为正整数。这种处理相当于将hd?n?,n=-∞~∞与函数w(n)相乘，w(n)具有下列形式：

?0,w(n)???1,n?0,n?N0?n?N

w(n)相当于一个矩形，我们称之为矩形窗。即我们可采用矩形窗函数w(n)将无限脉冲响应hd?n?截取一段h(n)来近似为hd?n?，这种截取在数学上表示为：

h(n)= hd?n?w(n) (7-7) 这里应该强调的是，加窗函数不是可有可无的，而是将设计变为物理可实现所必须的。 截取之后的滤波器传递函数变为：

H(z)??h(n)zn?0N?1?n (7-8)

式中，N为窗口宽度，H(z)是物理可实现系统。

为了获得线性相位，FIR滤波器h(n)必须满足中心对称条件（即7-3式），序列h(n)的延迟为???N?1?/2。

这种方法的基本原理是用一定宽度的矩形窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，从而得到FIR滤波器的脉冲响应，故称为FIR滤波器的窗函数设计法。

经过加矩形窗后所得的滤波器实际频率响应能否很好地逼近理想频率响应呢？图 7-1 示意给出了理想滤波器加矩形窗后的情况。理想低通滤波器的频率响应如图中左上角图，矩形窗的频率响应为左下角图。时间域内的乘积（7-7）式要求实际频率响应为这两个频率响应函数在频域内的卷积（卷积定理），即得到图形为图7-1（右图）。

图 7-1 FIR滤波器理想与实际频率响应

由图可看出，加矩形窗后使实际频率响应偏离理想频率响应，主要影响有三个方面：

（1）理想幅频特性陡直边缘处形成过渡带，过渡带宽取决于矩形窗函数频率响应的主瓣宽度。

（2）过渡带两侧形成肩峰和波纹，这是矩形窗函数频率响应的旁瓣引起的，旁瓣相对值越大，旁瓣越多，波纹越多。

（3）随窗函数宽度N的增大，矩形窗函数频率响应的主瓣宽度减小，但不改变旁瓣的相对值。

为了改善FIR滤波器性能，要求窗函数的主瓣宽度尽可能窄，以获得较窄的过渡带；旁瓣相对值尽可能小，数量尽可能少，以获得通带波纹小,阻带衰减大，在通带和阻带内均平稳的特点，这样可使滤波器实际频率响应更好地逼近理想频率响应。

这里我们明确两个概念：截断和频谱泄漏。信号是无限长的，而在进行信号处理时只能采取有限长信号，所以需要将信号“截断”。在信号处理中， “截断”被看成是用一个有限长的“窗口”看无限长的信号，或者从分析的角度是无限长的信号x(t)乘以有限长的窗函数w(t)。由傅立叶变换性质可知,时间域内的乘积对应于频率域的卷积,即

（7-9）

这里，x(t)是频宽有限信号，而w(t)是频宽无限信号，表示互为Fourier变换对。截断后的信号也必须是频宽无限信号，这样就是有限频带的信号分散到无限频带中去，这样就产生了所谓频谱泄漏。从能量的角度来看，频谱泄漏也是能量的泄漏，因为加窗后使原来信号集中的窄频带内的能量分散到无限的频带宽度范围内。频谱泄漏是不可避免的，但要尽量减小。

上边只考虑了矩形窗，如果我们使窗的主瓣宽度尽可能地窄，旁瓣尽可能地小，可以获得性能更好的滤波器，能否改变窗的形状而达到这个目的呢？回答是肯定的。其实数字信号处理的前驱者们设计了不同于矩形窗的很多窗函数，这些窗函数在主瓣和旁瓣特性方面各有特点，可满足不同的要求。为此，用窗函数法设计FIR数字滤波器时，要根据给定的滤波器性能指标选择窗口宽度N和窗函数w(n)。下面我们介绍窗函数。

7.2.2 MATLAB信号处理中提供的窗函数

（1）矩形窗：前面分析中所用的矩形窗可用下面函数来实现w=boxcar (N)，N 为窗的长度（以下函数与此同），w为返回的窗函数序列。

（2）汉宁窗：w=hanning(N)

汉宁窗的表达式为：

（3）哈明窗：w=hamming(N) 哈明窗的表达式为：

(7-10)

（4）Bartlett窗：w=bartlett(N) Bartlett 窗的表达式为：

当 N 为奇数时，

(7-11)

(7-12)

当 N 为偶数时，

(7-13)

（5） Blackman 窗：w= blackman(N) Blackman 窗的表达式为：

(7-14)

Blackman 窗比其他相同尺寸窗 (哈明窗,汉宁窗) 具有主瓣较宽和旁瓣泄漏较小的特点。

（6）三角窗：w=triang(N) 三角窗的表达式为： 当 N 为奇数时，

(7-15)

当 N 为偶数时，

(7-16)

三角窗和Bartlett窗十分类似。三角窗的两端值不为零，而Bartlett窗则为零，这一点可从例7-1中看出。

（7）Kaiser窗：w=kaiser(n,beta)

其中，beta是Kaiser窗参数,影响窗旁瓣幅值的衰减率。

Kaiser窗表达式：

式中， I0[.]是修正过的零阶 Bessel 函数。 (7-17)

Kaiser窗用于滤波器设计时，若旁瓣幅值为，则

（ 7-18 ）

（8） Chebyshev窗：w=chebwin(n,r)

式中, r 是窗口的旁瓣幅值在主瓣以下的分贝数。

Chebyshev窗具有主瓣宽度最小，而旁瓣等高,高度可调整的特点。

下面我们在MATLAB观看各种窗函数的形状和频率域图象来验证上述所讲特点。 【例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); ortlett 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 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); %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-2】 用 MATLAB 编程，采用512个频率点绘制各种窗函数的幅频特性。

%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); %求解窗函数的幅频特性，窗函数相当于一个数字滤波器 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(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('归一化频率');ylabel('振幅/dB'); title(stext);grid on; end

程序运行结果见图 7-3 。可以看到各种窗函数的幅频形状。对照该图可知这些窗函数具有上面所分析的窗函数的特征。

图 7-3 各种窗函数的幅频形状

由图 7-3 可见，各种窗函数都具有明显的主瓣（Mainlobe）和旁瓣（Sidelobe）。主瓣频宽和旁瓣的幅值衰减特性决定了窗函数的应用场合。矩形窗具有最窄的主瓣，但也有最大的旁瓣峰值（第一旁瓣衰减为 13 dB）；Blackman窗具有最大的旁瓣衰减，但也具有最宽的主瓣宽度。不同窗函数在这两方面的特点是不同的，因此应根据具体的问题进行选择。通常来讲，哈明窗和汉宁窗的主瓣具有较小的旁瓣和较大的衰减速度，是较为常用的窗函数。

表7-2总结了各种窗函数主瓣和旁瓣的特征（理论分析可参考其他的数字信号处理教材），大家可对照窗函数的幅频形状（图7-3）认真理解体会。

表7-2 各种窗函数的特点

窗函数 主瓣宽 第一旁瓣相对主瓣衰减（分贝） 矩形窗 -13 汉宁窗 -31 哈明窗 -41 Bartlett窗 -25 Blackman 窗 -57 三角窗 -25 Kaiser窗 Chebyshev 窗 可调整 可调整 可调整 可调整

主旁瓣频率宽度还与窗函数长度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); %用不同的窗长度求得复数频率特性 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

图 7-4 数据长度不同的矩形窗的幅频形状

程序运行结果见图7-4。显然，随着N的增大，主瓣和旁瓣都变窄，但第一旁瓣相对主瓣的幅值下降分贝数相同，第二旁瓣相对第一旁瓣幅值下降的分贝数也相同。这样，随着N的增大，旁瓣也得到抑制，有力地减少了频谱泄漏，但不能完全消除。减少主瓣宽度和抑制旁瓣是一对矛盾，不可兼得，只能根据不同用途折衷处理。

7.2.3 运用窗函数设计数字滤波器

用于信号分析中的窗函数可根据用户的不同要求选择。用于滤波器的窗函数，一般要求窗函数的主瓣宽度窄，以获得较好的过渡带；旁瓣相对值尽可能少，增加通带的平稳度和增大阻带的衰减。

基于窗函数的FIR数字滤波器设计的算法十分简单，其主要步骤为：

(1)对滤波器理想频域幅值响应进行傅立叶逆变换获得理想滤波器的单位脉冲响应hd(n)。一般假定理想低通滤波器的截止频率为，其幅频特性满足

则根据傅立叶逆变换，单位脉冲响应为： (7-19)

(7-20)

其中，为信号延迟。 （2）由性能指标（阻带衰减的分贝数）根据表7-2第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 节指出的函数或6.3.2节所给的函数进行滤波。

下面采用实例说明如何根据上面步骤设计FIR滤波器。

【例 7-4】 用窗函数设计一个线性相位FIR低通滤波器，并满足性能指标：通带边界的归一化频率wp=0.5,阻带边界的归一化频率ws=0.66，阻带衰减不小于30dB，通带波纹不大于3dB。假设一个信号，其中f1=5Hz,f2=20Hz。信号的采样频率为50Hz。试将原信号与通过滤波器的信号进行比较。

由题意，阻带衰减不小于30dB,根据表7-2,选取汉宁窗，因为汉宁窗的第一旁瓣相对主瓣衰减为31dB，满足滤波要求。在窗函数设计法中，要求设计的频率归一化到

0~区间内，Nyquist频率对应于，因此通带和阻带边界频率为0.5和0.66。程序如下

%Samp7_4 wp=0.5*pi;ws=0.66*pi; %滤波器边界频率 wdelta=ws-wp; %过渡带宽

N=ceil(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 的采样频率绘出该滤波器的幅频和相频响应 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; %采样间隔和检测信号的时间序列

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。该例设计通带边界wp=0.5,阻带边界频率ws=0.66，对应于50Hz的采样频率通带边界频率为fp=Fs/2*Fnormal=50/2*0.5=12.5Hz, fs=50/2*0.66=16.5Hz, 其中Fs为采样频率，Fnormal为归一化频率。由图7-5上图可以看到，在小于12.5Hz的频段上，几乎看不到下降，即满足通带波纹不大于3dB的要求。在大于16.5Hz的频段上，阻带衰减大于30dB，满足设计要求。由图7-5下图可见，在通带范围内，相位频率为一条直线，表明该滤波器为线性相位。图7-6给出了滤波器的输入信号和输出信号，输入信号包括3Hz和20Hz的信号，由图7-5可知，20Hz的信号不能通过该滤波器，通过滤波器后只剩下3Hz的信号，输出结果也证明了这一点。但要注意，由于FIR滤波器所需的阶数较高，信号延迟(N-1)/2也较大，输出信号前面有一段直线就是延迟造成的。上述程序显示的N取50才能满足设计要求。这样相位延迟为(N-1)/2*1/Fs=0.49s，可以看到输出信号前面一段直线的距离大约为0.49s。验证了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为滤波器截止频率，范围为0~1（归一化频率）。对于带通,带阻滤波器，wn=[w1,w2](w1

‘ ftype'为滤波器的类型：缺省时为低通或带通滤波器；'high'为高通滤波器；‘stop'为带阻滤波器，'DC-1'为第一频带为通带的多带滤波器；'DC-0'为第一频带为阻带的多带滤波器。

window为窗函数列向量，其长度为n+1。缺省时，自动取哈明窗。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-4。

%Samp7_5

wp=0.5*pi;ws=0.66*pi; %滤波器的边界频率 wdelta=ws-wp; %过渡带宽度

N=ceil(8* pi/wdelta);%求解滤波器的最小阶数，根据汉宁窗主瓣宽 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('振幅/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_6

wp=[0.35 0.65];N=48; %通带边界频率（归一化频率）和滤波器阶数 Fs=50;

b=fir1(N,wp); %设计FIR带通滤波器 figure(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);%输入信号 %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。通带归一化频率0.35和0.65对应于采样频率为50Hz的通带范围为8.75Hz和16.25Hz(采用6-19式计算)。由图7-7上图可见满足这一设计要求。在这个频带范围内的相位满足线性相位，符合FIR滤波器的一般特点。图7-8为检测滤波器的输入信号和输出信号。输入信号中含有1Hz,10Hz和20Hz的信号。根据图7-7上图可知，1Hz和20Hz的频率在阻带范围内，不能通过该滤波器，只有10Hz的信号可以通过该滤波器，输出信号表明了这一点。滤波器的相位延迟根据N/2*0.02s=0.48s得到，输出信号前面的直线部分大体为这个时间延迟,另外滤波后周期为10Hz的信号相位（红线开始部分），跟滤波前的相位（信号开始部分）也一致，说明通过该滤波器滤波后没有改变信号的相位。

图 7-7 例 7-6 设计滤波器的幅频特性（上图）和相频特性（下图）

图 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个采样点之间的信号并与对应的输入信号进行比较。

由于采样频率为1000Hz，所以该滤波器的归一化频率的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); %选择采样点间隔 x1=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的低频成分。

图 7-9 例7-7设计滤波器的输入信号和输出信号

【例7-8】设计采样频率为1000Hz,阻带频率从100Hz~200Hz的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]/(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*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上图可见，阻带范围为100~200Hz，完全符合设计要求。在通带的相位是线性的。由图7-11可见，滤波器滤除了150Hz（在阻带范围内）的信号，保留了50Hz的信号。相位延迟100/2/Fs=0.05s，与图形一致。

图 7-10 例7-8设计带阻滤波器的幅频（上图）和相频特性（下图）

图 7-11例7-8设计滤波器的输入和输出信号

该小节只给出了FIR低通，带通和带阻滤波器的例子，请大家在课下自己设计高通，第一频带为通频带和第一频带为阻频带的多带滤波器的例子，以加深对该函数的理解。

7.2.5 多频带FIR滤波器

除了设计标准型FIR滤波器外，MATLAB信号处理工具箱还提供另一种基于窗函数滤波器设计的工具函数fir2，用于设计具有任意形状频率响应的FIR滤波器，其调用格式为：

b=fir2（n，f,m[[,npt],window]）

式中，n为滤波器的阶数；f和m分别为滤波器期望幅频响应的频率向量和幅值向量，取值范围为0~1(归一化频率)。m，f具有相同的长度,window为窗函数，得到列向量，长度必须为n+1。缺省时自动取哈明窗；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 1 1 0 0 1 1 1 0 0],比较理想和实际滤波器的幅频响应。假设一个信号，其中f1=12Hz,f2=36Hz。信号的采样频率为100Hz。试将原信号与通过滤波器的信号进行比较。

%Samp7_9 clf

f=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:') %绘制理想幅频响应和设计的滤波器幅频响应 legend('理想特性', '实际设计') %给出图例

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; % 输入信号的两种频率成分

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阶时，实际滤波器的幅频响应才逼近理想滤波器的幅频响应。与第6章例6-12的输出比较可知，相同性能的FIR滤波器阶数比IIR滤波器要高得多。另外，该例中，信号含有两种频率成分，12Hz和36Hz，由于信号的采样频率为100Hz，因此Nyquist频率为50Hz，则归一化频率的0.2对应于10Hz，0.3对应于15Hz，因此，0.2~0.3的通带对应于10~15Hz。同理，归一化频率中的0.6~0.8的通带对应于30~40Hz。可见12Hz和36Hz的波均在通带的范围内，因此均可通过。该例的输出结果与例6-12的输出结果比较可知，FIR滤波器的相位延迟比IIR滤波器的相位延迟大得多。

图 7-12 例7-9所设计滤波器的幅频特性 上图：阶数为10；下图：阶数为100

图 7-13 例7-9所设计滤波器的输入和输出信号

7.3 最优FIR滤波器设计

MATLAB信号处理工具箱提供了比基于窗函数法FIR滤波器设计工具函数fir1和fir2更为通用的函数firls和remez。它们采用不同的优化方法设计最优的标准多频带FIR数字滤波器。

函数remez实现Park-McClellan算法，这种算法利用Remez交换算法和Chebyshev近似理论来设计滤波器，使实际频率响应拟合期望频率响应达到最优。从实际和理想频响之间最大误差最小化的观点来看，函数remez设计的滤波器是最优的，因此，又称之为最优滤波器。在频率域内，滤波器呈现等波纹特点，因此又称之为等波纹滤波器。Park-McClellan滤波器设计方法是FIR滤波器设计中最流行的，应用最广的设计方法。

函数firls和remez的调用格式语法规则相同，只是优化算法不同。

7.3.1 基本形式最优滤波器

函数基本调用格式为：

b=firls(n,f,a)

b=remez(n,f,a)

式中，n为滤波器阶数；f为滤波器期望频率特性归一化频率向量，范围为0~1，为递增向量，允许定义重复频率点；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，幅值为零。输入一个采样频率为50Hz,频率为5Hz和15Hz的合成振动，比较运用两种设计方法的输出的差异。

%Samp7_10

clf;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('振幅');

legend('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 所设计滤波器的幅频响应

图 7-15 例 7-10 所设计滤波器的输入和输出信号的比较

程序运行结果见图7-14和图7-15。比较两种设计方法所设计的滤波器幅频响应可见，用firls设计的滤波器在整个频率范围内（包括通带和阻带）均具有较好的响应，但理想频率响应和实际频率响应的误差在带区内不均匀，且在边界频率处误差较大；而函数remez设计的滤波器在通带内具有等波纹特性，在边界频率0.4和0.5处及过渡带内更接近于理想频响。由图7-15可见，输入5Hz和15Hz的两个频率的振动信号，5Hz的信号对应的采样频率为50Hz的归一化频率为0.2（5/(50/2)），15Hz的信号对应的归一化频率为0.6（15/(50/2)），对应于该滤波器的幅频响应，归一化频率为0.2的振动能通过滤波器，而归一化为0.6的振动不能通过该滤波器，输出结果也看到这一点，只有5Hz的振动通过了滤波器。但通过两种类型的滤波器滤波，输出结果略有不同，这是由于两种滤波器的幅频特性的略微不同造成的。运用firls函数设计幅频响应在归一化为0.2的频率处略小于remez函数的设计，因此运用firls函数设计滤波器的输出的归一化频率为0.2的振幅也略小于remez函数设计的输出。运用firls函数设计幅频响应在归一化为0.6的频率处比remez函数设计的振幅响应衰减更大，因此运用firls函数设计滤波器的输出更为平滑，即高频成分更小，而remez函数设计的滤波器的输出具有较多的高频成分，使得输出结果不对称。这完全可以根据滤波器的幅频响应分析出来。

函数firls和remez可用于设计低通，高通,带通和带阻等一般类型的滤波器，这可由函数中给定的理想幅频响应的频率向量f和幅值向量确定。

如设计一个带通滤波器，幅频响应定义为:

f=[0 0.3 0.4 0.7 0.8 1]， a=[0 0 1 1 0 0];

则，该理想滤波器幅频响应定义为：阻带频率0~0.3,0.8~1，通带频率0.4~0.7,过渡带:0.3~0.4,0.7~0.8。

设计一个高通滤波器，如理想幅频响应向量对的给定形式为： f=[0 0.7 0.8 1],a=[0 0 1 1],则，该理想滤波器的幅频响应定义为：阻带频率0~0.7，通带频率：0.8~1.0，过渡带 0.7~0.8。

设计一个带阻滤波器，如理想幅频响应向量对的给定形式为：f=[0 0.3 0.4 0.5 0.6 1],a=[1 1 0 0 1 1]，则，该阻带频率：0.4~0.5，通带频率：0~0.3,0.6~1.0,过渡带：0.3~0.4,0.5~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_11

clf;n=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]; %振幅向量

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('firls','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=[0 0.4 0.42 0.48 0.5 1], a=[1 1 0.8 0.2 0 0]; 这里，过渡带0.4~0.5给出多个响应值设计出的具有线性过渡带FIR滤波器的频率响应。大家可在课下练习。

7.3.2 加权最优滤波器

函数firls和remez还可以增加输入参数设置权向量w。在不同频段设置不同权值，使不同频段的误差值最小化得到不同程度的重视。具有权向量输入的函数firls和remez实现滤波器每个频率段加权处理。其调用格式为：

b=firls(n,f,a,w) b=remez(n,f,a,w)

式中，w为权向量，为f和a向量长度的一半，一个频带必须对应一个权值。

【例7-12】设计一个30阶的低通等波纹滤波器，通带边界频率0.4，阻带边界频率0.5，均为归一化频率，阻带波纹约为通带波纹的1/10。

%Samp7_12

clf;n=30; %滤波器的阶数 f=[0 0.4 0.5 1]; a=[1 1 0 0];

weit=[1 10];

b=firls(n,f,a,weit); %采用firls设计滤波器

[h,w1]=freqz(b); %计算滤波器的频率响应

bb=remez(n,f,a,weit); %采用remez设计滤波器 [hh,w2]=freqz(bb); % 计算滤波器的频率响应 figure(1)

plot(w1/pi,abs(h),'b.',w2/pi, abs(hh),'r:',f,a,'m');%绘滤波器幅频响应 xlabel('归一化频率');ylabel('振幅');

legend('firls','remez','理想特性'); %给出图例

程序运行结果见图7-18。该滤波器对阻带和通带波纹进行了控制，得到了符合要求的滤波器。由于remez要求等波纹的特点，因此，其在通带内的振动振幅较大，在阻带内的振动振幅确实较小，约为通带内振动振幅的1/10，符合设计要求。如果滤波器的阶数足够高，则可以获得更为理想的幅频响应，请大家自己修改程序进行试验。另外，大家可以自己设计一个输入信号，并观看其输出，试验这两种滤波器输出结果是否有差异。

图 7-18 例7-12所设计滤波器的幅频响应

7.3.3 反对称FIR滤波器——赫尔伯特变换器

利用函数firls和remez可设计滤波器系数为奇对称的滤波器，即：

（ 7-22 ）

这类滤波器要求滤波器的零频响应为0，若滤波器阶数为偶数，则还要求Nyquist频率（归一化频率为1）处的频率响应为零。该滤波器除了能对信号保持线性相位滤波外，还能实现信号的赫尔伯特(Hilbert)变换，故又称赫尔伯特变换器。所谓赫尔伯特变换就是使信号通过赫尔伯特变换后负频率作+90o的相移，正频率作-90o的相移，而振幅不发生改变。在地震学研究中，地震射线通过焦散面后会发生90o相移的畸变，通常采用赫尔伯特变换进行校正来正确识别地震震相。

函数调用格式为

b=firls(n,f,a,'h') b=firls(n,f, a,w,'h')

或

b=remez(n, f,a,'h') b=remez(n,f ,a ,w,'h')

式中，'h'为选择项，表示设计的滤波器是奇对称线性相位滤波器。

【例 7-13】利用函数remez和firls设计一个高通反对称线性相位滤波器，并绘制其频率特征图。设计一个采样频率为50Hz，频率为5Hz的振动作为输入信号，试验对应的数据点是否满足相位相差90o的特点。

%Samp7_13

clf;n=21; %滤波器的阶数为奇数

f=[0 0.1 0.2 1];a=[0 0 1 1]; %理想滤波器的幅频特性

b=firls(n,f,a,'h'); %采用 firls 设计奇对称系数滤波器

[h,w]=freqz(b,1,512); %计算滤波器的脉冲响应

bb=remez(n, f,a,'h'); %采用 remez 设计奇对称系数滤波器 [hh,w]=freqz(b,1,512); %计算滤波器的脉冲响应 figure(1)

plot (w/pi,abs(h),'b.',w/pi, abs(hh),'r--',f,a,'m'); %绘滤波器幅频响应 xlabel('归一化频率');ylabel('振幅');

legend('firls','remez','理想特性'); %给出图例 t=0:1/50:3;x=sin(2*pi*5*t); %时间序列和输入信号 figure(2)

subplot(2,1,1),plot(t(1:100),x(1:100)) %绘制输入信号的前100个样本 title('输入信号')

y1=filter(b,1,x); %运用firls设计的滤波器进行滤波 y2=filter(bb,1,x); %运用remez设计的滤波器进行滤波 subplot(2,1,2),

plot(t([1:100]+20/2),y1([1:100]+20/2),t([1:100]+20/2),y2([1:100]+20/2),'r:')

%绘制与输入信号对应的输出信号，考虑了延迟效应 legend('firls', 'remez') xlabel('时间/s') title('输出信号')

程序运行结果见图7-19和图7-20。该滤波器对阻带和通带波纹进行了控制，得到了符合要求的滤波器。由于两种滤波器的幅频响应非常一致，因此采用两种滤波器的输出信号的

o

振幅也很接近。但是跟输入信号比较可知，输出信号比输入信号前移了90的相位，这正是

o

赫尔伯特变换的结果。如果将输入信号改为-5Hz，可以看到输出的信号后移90的相位，请大家修改程序进行验证。

图 7-19 例7-13所设计滤波器的幅频响应

图 7-20 例7-13所设计滤波器的输入信号（上图）和输出信号（下图）

7.3.4 微分FIR滤波器

对信号时间域的微分等价于信号的傅立叶变换和一个虚数单位斜坡函数的乘积。也就是

说，信号的微分相当于让该信号通过一个频率响应为的滤波器。函数firls和remez可用于设计这种具有微分作用的FIR滤波器，调用格式为： b=firls(n,f,a,'d') b=remez(n,f,a,'d')

式中，'d'是选择项，表示设计的滤波器是具有微分器的作用。

【例7-14】利用remez和firls设计一个FIR微分器，频率在0~0.9范围内，绘制其频率特性图。

%Samp7_14

clf;n=30; %滤波器阶数

f=[0 0.9];a=[0 0.9]; %理想滤波器响应

b=firls(n,f,a,'d'); %采用firls设计滤波器 [h,w]=freqz(b); %计算滤波器频率响应

bb=remez(n,f,a,'d'); % 采用 remez 设计滤波器 [hh,w]=freqz(bb); %计算滤波器脉冲响应 figure(1)

plot(w/pi,abs(h),'b',w/pi, abs(hh),'r--',f,a,'m:'); %绘幅频响应 xlabel('归一化频率');ylabel('振幅');

legend('firls', 'remez','理想特性',4); %绘制图例，4表示图例位置在右下角 grid on ;

f1=5; %输入信号频率 t=0:1/1000:1; %时间序列 x=sin(2* pi*f1* t); %输入正弦信号 y=fftfilt(bb,x); %滤波 figure(2)

subplot(2,1,1), plot(t,x),title ('输入信号')% 绘输入信号 subplot(2,1,2), plot(t,y), %绘输出信号 ylim([-0.01 0.01]),xlabel('时间/s'),title('输出信号') 程序的输出为图7-21和图7-22。可见滤波器的幅频响应明显地表现为微分器的传递函数，由图7-22可见，程序输入一个正弦波后，输出确实为余弦波，表现为对原来信号进行了微分。

图7-21 例7-14所设计滤波器的频率响应

图7-22 例7-14的输入信号和输出信号

7.4 有限冲激响应数字滤波器的应用举例

7.4.1 淹没在噪声信号中地震波的提取

应用长春台2003年08月08日00点07分记录到的发生在吉林省内的2.9级地震。由于震级很小，在原始波形图上无法辨别，需要设计带通滤波器，把低频段和高频段干扰滤掉。由于地震波形数据的地震波频率成分较低。本例选择0.8Hz-5Hz的频段成分通过，即可满足要求,因此我们设计通带边界频率为0.8Hz--5Hz,阻带衰减为30dB，过渡带宽为0.5Hz。由于数字化采样率为50Hz，其Nyquist频率为25Hz，因此对应的归一化频率为0.8/25~5/25，因此通带边界频率为0.032*pi~0.2*pi,过渡带宽为0.025*pi。由于阻带衰减为30dB，因此我们

选用Hanning窗即可满足要求，根据表1中Hanning窗的主瓣宽度的定义可求出时的过渡带宽为0.025*pi的滤波器，最小阶数应为320。我们设计该滤波器并对信号进行滤波，源程序如下：

%Appl7_1 clear all;

wp=[0.032 0.2];N=320; %通带边界频率（归一化频率）和滤波器阶数 dt=0.02; %中国数字地震台网的采样间隔为0.02s,采样频率为50Hz b=fir1(N,wp,hanning(N+1)); %设计FIR带通滤波器 figure(1)

[H,f]=freqz(b,1,512,1/dt); %求出频率特性 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; load ChangChun.txt %加载地震波形记录

y=filtfilt(b,1,ChangChun); %产生零相位输出，采用filtfilt对输入信号滤波 t=[0:length(ChangChun)-1]*dt; figure(2)

subplot (2,1,1), plot (t,ChangChun),title ('输入信号')%绘出输入信号波形 subplot (2,1,2),plot (t,y) %绘出输出信号波形 title('输出信号'),xlabel('时间/s')

程序的输出为图7-23和图7-24，图7-23为滤波器的频率特性，图7-24为滤波器输入信号和输出信号的比较，可以看到，在滤波前，地震信号淹没在长周期的扰动之中，根本分不清地震信号的到时，滤波后可清楚地分辩地震波各个震相的到时，起到了预想的效果。

图7-23 设计的FIR滤波器的幅频特性（上图）和相频特性（下图）

图7-24 滤波前（上图）及滤波后（下图）的波形

第二个应用例子采用喀什地震台日常监测中记录到的一个地震信号的记录图，发震时刻2003年07月24日10时10分，震中距喀什地震台121km。为了更清楚地了解信号中的频率成分，我们采用下面的程序进行分析：

%Appl7_2

%用快速Fourier变换（FFT）对数字信号进行频谱分析的程序 load grbx3.txt; %读取数据序列

Xt=grbx3; %把数据赋值给变量 Fs=50; %设定采样率 单位(Hz) dt=1/Fs; %求采样间隔 单位(s) N=length(Xt); %得到序列的长度

Xf=fft(Xt); %对信号进行快速Fourier变换(FFT)

subplot(2,1,1),plot([0:N-1]/Fs,Xt); %绘制原始值序列 xlabel('时间/s'),title('时间域'); ylabel('振幅'); grid on;

subplot(2,1,2),plot([0:N-1]/(N*dt),abs(Xf)*2/N);%绘制信号的振幅谱 xlabel('频率/Hz'),title('频率域'); ylabel('振幅');

xlim([0 Fs/2]); %频率轴只画出Nyquist频率之前的部分 grid on;

运行上面的程序得到图7-25。可以看出原始地震信号中夹杂着多种干扰，完全掩盖了地震的信息。通过快速Fourier变换（FFT）分析可知，波形数据中有多组高频震动干扰。

图7-25 喀什地震台记录的2003年7月24日地震的原始信号及振幅谱

根据地震波干扰的频谱范围和特点，我们选择FIR带阻滤波器，地震信号采样频率Fs=50Hz，带阻开始频率为Fcs1 =3Hz，带阻结束频率Fcs2 =23Hz，采用400阶的FIR滤波器。由以下程序给出滤波器幅频相频特性和滤波前后数据的比较：

%Appl7_3

%用FIR数字滤波器实现对数字信号的滤波 load grbx3.txt %读取数据序列

Xt=grbx3; %把数据赋值给变量 Fs=50; %设定采样率 dt=1/Fs; %计算采样间隔 n=1:length(Xt);

Nn=length(Xt); %序列长度 t=n/Fs; %时间序列

Fcs1=3; Fcs2=23; %设置的通带和阻带边界频率，

Ws1=Fcs1/(Fs/2); Ws2=Fcs2/(Fs/2); %转换为标准频率 Wn=[Ws1 Ws2]; %通带、阻带频率(为标准频率) N=400; %滤波器阶数为400

pa=(N-1)/2/Fs; %计算相位延迟

b=fir1(N,Wn,'stop'); % 设计FIR带阻滤波器 figure(1);

[H,f]=freqz(b,1,Nn,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; figure(2);

subplot(2,1,1),plot(t,Xt); %绘制原始波形序列 xlabel('时间/s');ylabel('振幅');title('滤波前信号') grid on;

Yt=filter(b,1,Xt); %对信号进行滤波 t=t-pa; %计算相位延迟 subplot(2,1,2),plot(t,Yt); %绘制滤波后波形序列 ylim([-600 200])

xlabel('时间/s');ylabel('振幅');title('滤波后信号'); xlim([0 max(t)+1]);grid on;

运行上面的程序得到图7-26和图7-27。信号通过滤波器后滤除了高频振动干扰，滤波之前地震波形完全被高频干扰所掩盖，通过滤波后可以很容易分辨出地震中的各种震相。并且滤波后信号的相位无失真，滤波器达到了预期的要求。大家可以采用FFT分析滤波后的振幅谱，其中阻带内的频率成分完全被滤除。

图7-26 设计的滤波器的频率特性

图7-27 喀什地震台原始信号与滤波后信号的比较

第三个例子采用福建03311951.evt里 DOS台记录的地震图。截取的是NS向，截取时间从0—300秒。从原始图上可看出高频成分太多，把地震给淹没了，所以设计低通滤波器。因此我们设计通带边界频率为1Hz，阻带起始频率为2Hz，通带边界归一化频率为1/25，阻带边界归一化频率为2/25。采用通用的Hamming窗。程序如下：

%Appl7_4 figure(1)

load fjdos0351.txt; %读取数据序列 Xt= fjdos0351; %把数据赋值给变量

Fs=50; dt=0.02; %设定采样率(单位Hz)和采样间隔(单位s) N=length(Xt); t=[0:N-1]*dt; %得到序列的长度和时间序列 subplot(2,1,1),plot(t,Xt); %绘制原始值序列 title('滤波前'); ylabel('振幅');

Fcp=1; Fcs=2; %设置的通带边界频率，

Wp=Fcp/(Fs/2); Ws=Fcs/(Fs/2); %转换为标准频率 Wn=(Wp+Ws)/2; %给出设计时用的边界频率 wdelta=Ws-Wp; %过渡带宽

N=ceil(8*pi/wdelta); %按哈明窗求最小阶数 pa=(N-1)/2/Fs; %计算相位延迟 b=fir1(N,Wn); % 设计FIR带阻滤波器

Yt=filtfilt(b,1,Xt); %对信号进行滤波 t=t-pa; %计算相位延迟 subplot(2,1,2),plot(t,Yt); %绘制滤波后波形序列 xlim([0 300])

xlabel('时间/s');ylabel('振幅');title('滤波后信号'); figure(2)

[H,f]=freqz(b,1,128,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;

运行结果为图7-28和图7-29。由图7-28可见滤波器滤除了高频与一部分低频，使得地震波信息显现出来。滤波器的频率特性为图7-29，完全符合设计要求。

图7-28 滤波前和滤波后数据的比较

图7-29 滤波器的幅频特性

7.4.2 爆破信号的滤除

地震台观测的数据中除了地震震源激发的地震波外，还包含了如地方爆破等高频干扰数据。呼和浩特记录的2003年7月25日一远震UD分向。在P波中叠加了一个爆破记录。爆破频率相对地震波的频率明显较高，为了去掉该爆破，需要设计低通滤波器，把高频爆破滤掉。由于远震波形数据的地震波频率成分较低。我们选择低于2Hz的频率成分即可满足要求,因此我们设计通带边界频率为2Hz,阻带边界频率为2.5Hz,阻带衰减为40dB，可选择哈明窗，归一化通带边界频率为2/25，根据阻带宽度为0.5/25*pi以及哈明窗的主瓣宽度，可求出滤波器的最小阶数为400，设计的程序如下：

%Appl7_5 clear all;

wp=2/25;N=400; %通带边界频率（归一化频率）和滤波器阶数 dt=0.02; %中国数字地震台网的采样间隔为0.02s,采样频率为50Hz b=fir1(N,wp,hamming(N+1)); %设计FIR低通滤波器 figure(1)

[H,f]=freqz(b,1,512,1/dt); %求出频率特性 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; load cz02.txt %加载数据 x=cz02'; %输入信号

y=filtfilt(b,1,x); %产生零相位输出，采用filtfilt对输入信号滤波 t=[0:length(x)-1]*dt; figure(2)

subplot (2,1,1), plot (t,x),title ('输入信号')%绘出输入信号波形 subplot (2,1,2),plot (t,y) %绘出输出信号波形 title('输出信号'),xlabel('时间/s')

输出结果为图7-30和图7-31。可见滤波器的设计完全符合设计要求，滤波器完全滤除了地震信号前面的爆破信号。

图7-30 设计的滤除爆破信号的滤波器的频率特性

图7-31 原始信号和通过滤波器的信号比较

7.4.3 低频干扰的滤除

我们选择云南地震记录文件yn9705080600.ev2的bss台站记录，运用firls和remez函数对该地震记录进行处理。采用地震波截取软件得到地震数据的ASCII文件bssud.txt。经分析知该地震含有高频干扰和低频干扰，我们设计带通滤波器保留1.6-2.5Hz的频率成分，程序如下：

%Appl7_6

load bssud.txt %加载数据 x=bssud; %输入信号 n=320; %滤波器阶数

dt=0.02; %中国数字地震台网的采样间隔为0.02s,采样频率为50Hz t=[0:length(x)-1]*dt; %时间序列

f=[0 0.5/25 0.6/25 1.5/25 1.6/25 1]; %频率向量 a=[0 0 1 1 0 0]; %振幅向量

b=firls(n,f,a); %采用firls设计滤波器 [h,w]=freqz(b); %计算其频率响应 bb=remez(n,f,a); %采用remez设计滤波器 [hh,w]=freqz(bb); figure(1)

subplot (2,1,1),plot (t,x),title ('输入信号') %绘出输入信号波形

subplot (2,1,2),plot(w/pi,abs(h),w/pi,abs(hh),'r:'); %绘制滤波器的幅频响应 xlabel('归一化频率');ylabel('振幅');

legend('firls','remez'); %给出图例 %滤波

y1=filtfilt(b,1,x); %产生零相位输出，采用filtfilt对输入信号滤波 y2=filtfilt(bb,1,x); %产生零相位输出，采用filtfilt对输入信号滤波 t=[0:length(x)-1]*dt; figure(2)

subplot (2,1,1),plot (t,y1,'b') %绘出输出信号波形 title('输出信号'),xlabel('时间/s')

legend('firls'); %给出图例 subplot (2,1,2),plot (t,y2,'r') %绘出输出信号波形 title('输出信号'),xlabel('时间/s')

legend('remez'); %给出图例

运行结果为图7-32和图7-33。可见通过滤波器滤波后使得地震波更清楚地显现出来，滤波器的幅频特性也符合设计要求。

图7-32 云南1997年5月8日地震波原始数据及设计滤波器的幅频特性

输出信号1000firls5000-500-1000020406080100时间/s输出信号1201401601802001000remez5000-500-1000020406080100时间/s120140160180200 图7-33 运用firls和remez滤波后的地震图

7.5 无限冲激响应数字滤波器和有限冲激响应数字滤波器的比较

本章和第6章讨论了许多种设计无限冲激响应和有限冲激响应数字滤波器的方法。这样自然会产生一个问题：究竟哪类滤波器最好？是无限冲激响应还是有限冲激响应？为什么要给出这么多的设计方法？对这些问题的回答是：我们之所以讨论有限冲激响应和无限冲激响应两类滤波器，讨论它们的各种设计方法，是因为没有一种滤波器，也没有一种设计方法在所有的情况下都是最佳的。

选择有限冲激响应滤波器还是选择无限冲激响应滤波器取决于人们对每一类滤波器优缺点的权衡比较结果。例如，无限冲激响应滤波器的优点是可以利用现成的设计公式。就是说一旦规定了用哪种已知类型的滤波器（例如比特沃思、切比雪夫或椭圆滤波器），则直接把设计指标代入设计方程组，就可以得到待求数字滤波器的所有系数（或极点和零点）。若只要求设计少数几阶滤波器，或计算设备能力有限，则无限冲激响应滤波器是一种可选方案。

有限冲激响应滤波器没有现成的设计公式。经常采用迭代算法来满足预定技术指标，因

此设计这些滤波器需要功能较强、容量较大的计算机。而设计无限冲激响应滤波器只要用模拟滤波器设计参数表就可以了。然而设计过程上的简单性是牺牲了滤波器响应的灵活性换来的。此外这些设计中一般忽略了滤波器的相位响应。因此，我们虽然可以用比较简单的设计方法得到幅度特性极好的椭圆低通滤波器，但它的相位响应具有严重的非线性（尤其是在频带边缘上）。

与此相反，有限冲激响应有精确的线性相位特性。同时窗函数法和其他大多数算法都能够逼近更加任意的频率响应。因此，看起来设计有限冲激响应滤波器似乎比设计无限冲激响应滤波器灵活自如得多，因为设计有限冲激响应滤波器有最佳性定理，该定理在许多实际情况下是有意义的。

最后，实现数字滤波器时还要考虑经济性问题，经常用硬件的复杂程度和计算机速度来衡量它。这两个因素或多或少同滤波器满足指标所需的阶数有关系。如果我们撇开相位问题不予考虑，一般来说，无限冲激响应滤波器能有效地满足给定的频率响应指标。然而在许多情况下，为得到线性相位特性，采用有限冲激响应滤波器花费高昂的代价也是值得的。

第七章 习题

1. 用Hanning窗函数设计一个带阻FIR滤波器，性能指标为阻带边界频率0.4?~0.6?,

通带边界频率为0.3?和0.7?。绘制其频率特性图、脉冲响应图、零极点图，并设计一个序列进行滤波，检验此滤波器的滤波效果。

2. 用Hamming窗函数设计一个带通FIR滤波器，性能指标为通带边界频率0.4?~0.6?,

低通边界频率为0.3? 和0.7?。绘制其频率特性图、脉冲响应图、零极点图，并设计一个序列进行滤波，检验此滤波器的滤波效果。

3. 用窗函数法设计一个多频带FIR滤波器，理想频率响应为：频率[0 0.1 0.2 0.3 0.4 0.5

0.6 0.7 0.8 0.9 1] ?,幅值[0 0 1 1 1 0 0 1 1 0 0]，绘制其理想幅频图并与实际滤波器幅频图相比较。设计一个例子说明此滤波器满足要求。绘制其群延迟图。 4. 用MATLAB函数firls和remez设计一个带通滤波器，理想频率响应对为：频率[0 0.3 0.4

0.7 0.8 1] ?,幅值[0 0 1 1 0 0]，绘制其理想幅频图并与实际滤波器幅频图相比较。设计一个例子说明此滤波器满足要求，绘制其群延迟图。

5. 已知周期信号x(t)=0.75+3.4cos2?ft?2.7cos4?ft?1.5sin3.5?ft?2.5sin7?f,其

中，f=25/16Hz，若截断时间长度为信号最大周期的0.9和1.1倍，试绘制和比较采用下面窗函数提取x(t)的频谱：（1）矩形窗；（2）汉宁窗；（3）哈明窗；（4）巴特利特窗；（5）Blackmann窗；（6）三角窗；（7）Kaiser窗；（8）Chebyshev窗。

本文来源：https://www.bwwdw.com/article/q57p.html