熊欣
（武汉理工大学信息工程学院，湖北武汉 430070）
摘要：目前常用的FIR 滤波器的设计方法主要有三种：窗函数法、频率采样法和最优化设计法。本文针对三种设计方法基于MATLAB 进行FIR各种形式滤波器的设计与仿真，并比较了三种方法的特点。结果表明，在同样的设计指标下，利用等波纹切比雪夫逼近法则的设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。 关键词：FIR滤波器；MATLAB; 窗函数法；频率采样法；等波纹切比雪夫逼近法
1 引言
数字滤波器按照单位取样响应h(n)的时域特性可以分为无限脉冲响应（IIR）系统和有限脉冲响应（FIR）系统。FIR 数字滤波器的优点在于它可以做成具有严格线性
相位，而同时可以具有任意的幅度特性；它的传递函数没有极点；这保证了设计出的FIR
数字滤波器一定是平稳的。
所谓数字滤波器设计，简单地说，就是要到一组能满足特定滤波要求的系数向量a和b。而滤波器设计完成后还需要进一步考虑如何将其实现，即选择什么样的滤波器结构来完成滤波运算。FIR数字滤波器的设计方法很多，其中较为常用的是窗函数设计法、频率采样设计法和最优化设计法。本文讨论利用窗函数法、频率采样法和等波纹切比雪夫逼近法（调用remez函数）来分别实现各种FIR滤波器的设计。
窗函数法设计的基本思想是把给定的频率响应通过IDTFT（Inverse Discrete Time Fourier Transform）[1]，求得脉冲响应，然后利用加窗函数对它进行截断和平滑，以实现一个物理可实现且具有线性相位的 FIR 数字滤波器的设计目的。其核心是从给定的频率特性，通过加窗确定有限长单位取样响应h(n)；频率采样法设计的基本思想是把给出的理想频率响应进行取样，通过 IDFT 从频谱样点直接求得有限脉冲响应；等波纹切比雪夫逼近法，利用 MATLAB 提供的 remez 函数实现Remez算法，设计滤波器逼近理想频率响应。
2 FIR滤波器设计与仿真
FIR 滤波器设计的任务是选择有限长度的h(n)，使传输函数H(e jw)满足一定的幅度特性和线性相位要求。由于FIR 滤波器很容易实现严格的线性相位，所以FIR 数字滤波器设计的核心思想是求出有限的脉冲响应来逼近给定的频率响应[2]。
目前在MATLAB信号处理工具箱中，
一共有10种FIR数字滤波器设计函数，对
应5种设计方法，如表1：
表1 FIR滤波器设计函数
函数设计方法
fir1,fir2,kaiserord 窗函数法
firls,remez,remezord等波纹最小平方误差
设计
fircls,fircls1 最小二乘约束设计
cremez 任意频响“复滤波器”
设计
firrcos 升余弦设计
2.1 窗函数法金刚石碎片
设计FIR数字滤波器的最简单的方法是窗
函数法，通常也称之为傅立叶级数法[3]。
FIR数字滤波器的设计首先给出要求的理
想滤波器的频率响应H d(e jw)，设计一个FIR
数字滤波器频率响应H(e jw)，去逼近理想的
滤波响应H d(e jw)。然而，窗函数法设计FIR
数字滤波器是在时域进行的，因而必须由理
想的频率响应H d(e jw)推导出对应的单位
取样响应h d(n)，再设计一个FIR数字滤波
器的单位取样响应h(n)去逼近h d(n)。设计过
程如下：
H d (e )
加窗的作用是通过把理想滤波器的无限长脉冲响应h d (n)乘以窗函数w(n)来产生一个被截断的脉冲响应，即h(n)= h d (n)w(n)，并且对频率响应进行平滑。MATLAB 工具箱提供的窗函数有：矩形窗(Rectangular window)、三角窗(Triangular window)、布拉克曼窗(Blackman window)、汉宁窗(Hanning window)、海明窗(Hamming window)、凯塞窗(Kaiser window)、巴特里特窗(Bartlett window)、切比雪夫窗(Chebyshev window)，性能对比如表2。
表
2 几种窗函数的性能比较
窗函数主要用来减少序列因截断而产生的Gibbs 效应。但当这个窗函数为矩形时，得到的FIR 滤波器幅频响应会有明显的Gibbs 效应，并且任意增加窗函数的长度（即FIR 滤波器的抽头数）Gibbs 效应
也不能得到改善。为了克服这种现象，窗函数应该使设计的滤波器：（1）频率特性的主瓣宽度应尽量窄，且尽可能将能量集中在主瓣内；(2)窗函数频率特性的旁瓣ω趋于π 的过程中，其能量迅速减小为零。
下面用矩形窗和海明窗分别设计FIR DF 并分析两者的区别，给出信号的采样频率是1000HZ ，数字滤波器的截止频率是200HZ ，阶数为60。程序如下： passrad=0.4*pi; w1=boxcar(61); w2=hamming(61) n=1:1:61;
hd=sin(passrad*(n-31))./(pi*(n-31));
hd(31)=passrad/pi; h1=hd.*rot90(w1); h2=hd.*rot90(w2); [mag1,rad]=freqz(h1); [mag2,rad]=freqz(h2); subplot(2,2,1);
plot(rad,20*log10(abs(mag1))); grid on;
title('designed by Rectangular window'); subplot(2,2,2);
plot(rad,20*log10(abs(mag2))); grid on;
title('designed by Hamming window'); [h1,w1]=freqz(h1,1,100,2); subplot(2,2,3);
plot(w1,unwrap(angle(h1))); grid on;
[h2,w2]=freqz(h2,1,100,2); subplot(2,2,4);
窗函数
旁瓣峰值幅度（dB ）
过渡带宽度（π/N ） 阻带最小衰减
（dB ）
矩形窗 13 4 21 三角窗 25 8 25 汉宁窗 31 8 44 海明窗 41 8 53 凯塞窗（β=7.856） 57 10 80 布拉克曼窗
57 12 74
plot(w2,unwrap(angle(h2))); grid on;
图1 窗函数设计的FIR 低通滤波器频率响应
可以看出，采用特殊的窗函数如
Hamming 窗，可以减小这种Gibbs 效应，但同时也会使滤波器的过度带变宽。
窗函数法还可以用MATLAB 信号处理工具箱中的FIR 滤波器设计函数fir1、fir2和kaiserord 来实现。函数fir1用来设计线性相位的低通、高通、带通、带阻FIR 滤波器（默认为低通），它使用标准的窗函数法进行设计（默认窗函数为Hamming ），滤波器的阶数由参数n 指定（fir1在设计高通和带阻滤波器时得到的滤波器的阶数总是偶阶），截止频率由参数Wn （归一化截止频率，对应于频响曲线上的-6dB 点）定义；函数fir2用来设计多通带任意响应FIR 滤波器，该滤波
器的幅频特性由向量对f 和m 确定，f 为归一化频率向量，m 为对应频率点上的幅度。当设计的滤波器在频率为π的幅度响应不
是0时，滤波器的阶数n 为偶数；当函数kaiserord 中各参数得到后可以利用返回值n 和beta 设计凯塞窗函数，然后利用返回值Wn 和ftype 传输给fir1进行滤波器的设计。三种函数比较完整的语法形式如下： b= fir1(n, Wn , ‘ftype’, window)； b= fir2(n, f, m, npt, lap, window)；
[n, Wn , beta, ftype] = kaiserord (f,a,dev, fs)； 下面分别用这三种函数设计：
I 用fir1设计一个27阶FIR 低通滤波器，
其中通带截止频率为0.2π，用汉宁窗函数。 n=27;
Wn=0.2;
window=hanning(28);
b=fir1(n,Wn,window);
freqz(b);
图1 fir1设计的FIR 低通滤波器频率响应金属槽筒
II 用fir2设计一个60阶的FIR 滤波器，要求滤波器0到π/4的幅度响应为0 ，π/4到π/2的幅度响应为1/4，π/2到3π/4的幅度响应为0，3π/4到1的幅度响应为1。 n=60;
f=[0 0.25 0.25 0.50 0.50 0.75 0.75 1];
m=[0 0 1/4 1/4 0 0 1 1];
%对幅频响应插值时插值点的个数
npt=1024;
%插值时不连续点转变成连续时的点数 lap=50; %衰减为30dB 的切比雪夫窗函数
window=chebwin(61,30);
b=fir2(n,f,m,npt,lap,window);
freqz(b);
图2 fir2设计的FIR 多通带滤波器频率响应III 用kaiserord 函数设计一个带通FIR 滤波器，通带范围是1500HZ 到2500HZ ，通带的波纹最大为0.03，阻带范围是0HZ 到1000HZ 和3500HZ 到5000HZ ，阻带的波纹最大为0.01 ，信号的采样频率为10KHZ. f=[1000 1500 2500 3500]; a=[0 1 0]; fs=10000;
dev=[0.01 0.03 0.01];
[n,Wn,beta,ftype]=kaiserord(f,a,dev, fs);
hh=fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale') freqz(hh);
查看Workspace 得到滤波器的阶数为45，截止频率为0.25π和0.6π，凯塞窗函数的β值为3.3953。
图3 kaiserord 设计的FIR 带通滤波器频率响应 2.2 频率采样法频率采样法是从频域出发，根据频域采样定理，对给定的理想滤波器的频率响应H d (e jw )加以等间隔的抽样[4]，得到H d (k)：
H d (k)= H d (e jw
)∣w=(2π)k/N k=0,1,…,N-1 (1) 再利用H d (k)可求得FIR 滤波器的系统函数)
及频率响应H(e jw
):
H(Z
(2) (3) 其中，)(w φ是一个内插函数：
(4) 从以上公式可以看出，在每个采样频率
点w k =2πk/N 处，滤波器的实际频率响应是
严格地和理想频率响应数值相等，即：
(5)
防老剂rd而在各采样点间的频率响应则是其的
加权内插函数延伸叠加的结果。但对于一个
无限长的序列，用频率采样法必然有一定的逼近误差，误差的大小取决于理想频响曲线的形状, 理想频响特性变换越平缓, 则内插函数值越接近理想值，误差越小。为了提高逼
地磁指数预报近的质量，可以通过在频率相应的过渡带内
插入比较连续的采样点，扩展过渡带使其比
较连续，从而使得通带和阻带之间变换比较
缓慢，以达到减少逼近误差的目的。
选取w ∈[0,2π]内N 个采样点的约束条
件为：
(6) 下面用频率采样法设计一个FIR DF 通带截止频率为0.2π，通带波纹最大为0.04，阻带截止频率为0.3π，阻带波纹最大为0.02，滤波器的阶数通过remezord 函数估算，程序如下： f=[0.2 0.3]; a=[1 0];
dev=[0.04 0.02];
%给出滤波器的参数
[n f0 a0 w]=remezord(f,a,dev); N=n;
alpha=(N-1)/2; k=0:N-1;
wp=0.2*pi; ws=0.3*pi;
%计算理想低通滤波器的截止频率
wc=(wp+ws)/2;
m=fix(wc*N/(2*pi)+1);
%在两边过渡带取值为0.5的采样点
黄曲霉毒素测定T = 0.5; Hrs=[ones(1,m),T,zeros(1,N-2*m-1),T,ones(1,
m-1)]; k1 = 0:floor(alpha); k2 = floor(alpha+1):N-1;
phai=[-alpha*(2*pi)/N*k1,alpha*(2*pi)/N*(N
-k2)];
H =Hrs.*exp(j*phai); %计算单位冲激响应 h =ifft(H,N);
[h1,w1] = freqz(h,1,256,1);
hr = abs( h1); h1 = 20* log10(hr);
%画出FIR DF 的单位取样响应 figure(1);车载广告
k=0:N-1; stem(k,h,'k.')
axis([0,N-1,1.1*min(real(h)),1.1*max(real(h))
]);
xlabel('n'); ylabel('h(n)');
grid on;
%画出FIR DF 的低通衰减幅频特性
figure(2);
plot(w1,h1);
xlabel('Normalized Frequency(×πrad/sample)'); ylabel('Magnitude(dB)');
grid on; 由Workspace 查出滤波器参数为n=27, f0=[0 0.2 0.31], a0=[1 1 0 0 ],w=[1 2 ]。在理想滤波器的设计中，若不增加过渡点，阻带和通带之间的衰减约为21dB 。如果在通带和阻带之间增加一个采样点，阻带的最小衰减可以达到65dB （衰减程度由选择的采样点决定）。本程序在过渡带取值为0.5的采样点，阻带衰减约为40dB .
图4 FIR 的单位取样响应
图5 FIR的低通衰减幅频特性
2.3 等波纹切比雪夫逼近法
尽管窗函数法与频率采样法在FIR数字滤波器的设计中有着广泛的应用, 但两者都不是最优化的设计[5]。通常线性相位滤波器在不同的频带内逼近的最大容许误差要求
不同。等波纹切比雪夫逼近准则就是通过对通带和阻带使用不同的加权函数，实现在不同频段(通常指的是通带和阻带) 的加权误差最大值相同，从而实现其最大误差在满足性能指标的条件下达到最小值，即使得
H d(e jw) 和H(e jw)之间的最大绝对误差最小。
等波纹切比雪夫逼近是采用加权逼近误差E（e jw），它可以表示为：
E（e jw）=W（e jw）（H d(e jw)-H(e jw)）（7）其中，W（e jw）为逼近误差加权函数，在误差要求高的频段上，可以取较大的加权值，否则，应当取较小的加权值。
尽管按照FIR数字滤波器单位取样响应h(n)的对称性和N的奇、偶性，FIR数字滤波器可以分为4种类型，但滤波器的频率响应可以写成统一的形式：
H(e jw)=e-j(N-1)w/2e j(π/2)k H(ω) (8)
其中，k∈{0 ,1} , H (ω)为幅度函数，且是一个纯实数，表达式也可以写成统一的形式：
H(e jw)=Q（ω）P（ω）(9)
其中，Q（ω）为ω的固定函数，P（ω）为M个余弦函数的线性组合。若令：W(ω)= W（e jw）Q（ω），H d(ω)= H d(e jw)/ Q（ω），因此，由式（9）、（10）将E（e jw）改写成：E（e jw）=W(ω)[ H d(ω)- P（ω）] (10) 故等波纹切比雪夫逼近法设计FIR数字滤波器的步骤是：
①给出所需的频率响应H d(e jw)，加权函数W （e jw）和滤波器的单位取样响应h ( n)的长度N。
②由①中给定的参数来形成所需的W(ω)、
H d(ω)和P（ω）的表达式。
③根据Remez 算法，求解逼近问题。
④利用傅立叶逆变换计算出单位取样响应h ( n) 。
Remez 算法是由Parks 和McClellan 等人在1972年推导出来的。它是将FIR数字滤波器中的五个参数( N ,σ1 , σ2 ,ωp ,ωs) 中的N ,ωp , ωs 和σ1/σ2 固定，而视σ1 (或σ2) 为变量的一种迭代方法。在MATLAB工具箱中可以直接调用remez函数（采用Remez 算法），来进行FIR数字滤波器的设计。其具体算法有几种，常见的一种算法格式为：
b=remez (n, f, a, w, ‘ftype’);
下面用remez函数设计一个27阶的FIR 低通滤波器，其通带截止频率为0.2π，通带波纹最大为0.04，阻带截止频率为0.3π，阻带波纹最大为0.02，程序如下：
n=27;
f=[0 0.2 0.3 1];
a=[1 1 0 0];
w=[0.04 0.02];
b=remez(n,f,a,w);
freqz(b);
图6 remez设计的FIR低通滤波器频率响应
3 结论
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议