FIR数字滤波器设计

今天给大侠带来FIR数字滤波器设计，由于篇幅较长，分三篇。今天带来第三篇，FIR数字滤波器设计，包括窗函数法设计FIR滤波器、频率采样法设计FIR滤波器以及基于firls函数和remez函数的最优化方法设计FIR滤波器。话不多说，上货。

这里也超链接了上、中两篇，方便参考学习。

FIR数字滤波器设计（上）

FIR数字滤波器设计（中）

数字滤波器的输入输出均为数字信号，信号通过数字滤波器后，可以改变频率成分的相对比例或滤除某些频率成分。数字滤波器可以分为IIR数字滤波器和FIR数字滤波器。

本篇只介绍FIR数字滤波器的设计，可以根据所给定的频率特性直接设计FIR数字滤波器。FIR数字滤波器在保证幅度特性满足要求的同时，能够做到严格的线性特性。

本篇采用了窗函数法、频率采样法以及基于firls函数和remez函数的最优化方法设计FIR滤波器。对FIR滤波器进行了详细的理论分析，并且对应于每种方法都给出了设计实例。通过编写MATLAB语言程序，运行程序，得到幅频和相频特性图。

对于窗函数和firls函数设计的滤波器，还通过建立Simulink系统模块进行仿真，观察滤波器滤波情况。

FIR数字滤波器设计

FIR数字滤波器的设计方法有窗函数法、频率采样法和基于firls函数和remez函数的最优化方法。MATLAB语言中的数字信号处理工具箱，提供了一些滤波器的函数，使FIR滤波器的运算更加方便和快捷。

在MATLAB中提供的滤波函数有fir1()，此函数以经典的方法实现加窗线性相位FIR数字滤波器设计，可以设计出低通、高通、带通和带阻滤波器；fir2函数设计的FIR滤波器，其滤波的频率特性由矢量f和m决定，f和m分别为滤波器的期望幅频响应的频率相量和幅值相量。

Firls()和remez()的基本格式用于设计I型和II型线性相位FIR滤波器，I型和II型的区别是偶函数还是奇函数。freqz()用于求数字滤波器的频率响应。并且提供了各种窗函数的函数，比如，hamming()是海明窗函数，hanning()是汉宁窗函数，kaiser()是凯泽窗函数，使在设计的过程中，不用自己重新设计窗函数。

Simulink是MATLAB众多工具包中的一员，对于建模，Simulink提供了一个图形化的用户界面（GUI）。Simulink包括一个复杂的由接收器、信号源、线性和非线性组件以及连接件组成的模块库。定义完一个模型后，就可以通过Simulink的菜单或者在MATLAB的命令窗口输入命令对它进行仿真。使用Scopes或者其他的显示模块，可以在运行仿真时观察到仿真的结果。另外，还可以在仿真时改变参数，并且立即就可看到变化。

一、窗函数法设计FIR滤波器

设我们所要设计的FIR滤波器的传输函数是 ， 是与其对应的单位脉冲响应，因此

（3-1）

 （3-2)

如果我们能够在已知的情况下，求出，经过Z变换可得到滤波器的系统函数。通常情况下理想数字滤波器的单位脉冲相应是无限长的，且是非因果序列。获得有限脉冲响应滤波器的一种可能方法是对截取一段来近似代替，可是这样会改变原来的滤波器指标，出现吉布斯效应误差。

窗函数法就是用被称为窗函数的有限加权序列w(n)来修正式(1)的傅里叶基数以求得要求的有限脉冲响应序列，即：

（3-3）

w(n)是有限长序列，当n<0或n>N-1时，w(n)=0。

这种方法的重点在于选择某种合适的窗函数。要求窗函数主瓣宽度尽可能窄，以获得最小的过渡带；旁瓣相对值尽可能小，以使得通带波纹小，并且阻带衰减大。

下面介绍几种常用的窗函数：

1、矩形窗(Rectangle Window)

（3-4）

其频率函数为：

（3-5）

2、三角形窗(Bartlett Window)

（3-6）

其频率函数为：

（3-7)

3、汉宁(Hanning)窗，又称升余弦窗

（3-8）

利用傅里叶变换得到频率函数为：

（3-9）

当  时，，所以窗函数的幅度函数为：

 (3-10)

4、汉明(Hamming)窗，又称改进的升余弦窗

（3-11）

其幅度函数为：

（3-12）

5、布莱克曼(Blankman)窗，又称二阶升余弦窗

（3-13)

其幅度函数为：

(3-14）

6、凯泽(Kaiser)窗

（3-15）

其中：β是一个可自由选择的参数，是第一类修正零阶贝塞尔函数。

上述窗函数的基本参数如下表：

窗函数法设计滤波器的步骤：

1）根据技术要求确定等待求值滤波器的单位取样响应。

2）根据对过渡带和阻带衰减的要求，选择窗函数的形式，并估计窗口长度N。

3）计算滤波器的单位取样响应：

（3-16)

上式中， 是前面所选择好的窗函数。

4）检验技术指标是否满足要求。根据下式计算：

(3-17)

如果不满足要求，根据具体情况重复步骤(2)(3)(4)步，直到满足要求为止。

本篇以一个FIR滤波器的设计为例说明如何使用MATLAB设计数字滤波器。

设计实例：用窗函数法设计线性相位FIR低通数字滤波器，要求通带截止频率Wp=0.4*，阻带截止频率Ws=0.5*, 通带衰减不大于3db，阻带衰减不小于40db。

1、程序如下：

 

Wp=0.4*pi;Ws=0.5*pi;Wdel=Ws-Wp;N=ceil(8*pi/Wdel);Wn=(0.4+0.5)*pi/2;window=hanning(N+1);b=fir1(N,Wn/pi,window);freqz(b,1,512)

 

程序执行后所得幅频和相频如下图所示：

2、用Simulink仿真对窗函数设计的低通滤波器进行滤波：

FIR低通滤波器的Simulink仿真模块图如下所示：

输入信号：s(t)=sin(30πt)+sin(56πt)+sin(70πt)

数字滤波器参数设置：已知滤波器的阶数n=80。

首先在Filter Type中选择lowpass；

在Design Method选项中选择FIR Window，接着在Window选项中选取hann；指定Filter Order项中的Specify order为80；

采样频率Fs=100Hz,截止频率Fc=22.5HZ。

设置完以后点击窗口下方的Design Filter，在窗口上方就会看到所设计滤波器的幅频响应和相频响应曲线。

设计界面如下图所示：

信号源参数设置：三个信号源的频率分别是15Hz、28Hz、35Hz ：

运行仿真，得波形如下图所示：

15HZ信号波形如下：

28HZ信号波形如下：           

35HZ信号波形如下：

滤波前的叠加波形如下：         

滤波后波形如下：   

由上图可知，15Hz、28Hz、35Hz叠加后通过截止频率为22.5HZ的FIR低通滤波器，输出与15Hz信号波形相似，其它的高频被阻隔。

二、频率采样法设计FIR滤波器

1、对理想滤波器的系统函数进行频率采样以得到系统的理想频响的等间隔采样值H(k)。

H(k)实际上是所要求的滤波器的单位采样响应(h(n))的离散傅里叶变换（DFT），如下关系式：

 (3-18)

  (3-19)

为了减小H（k）的通带边缘由于抽样点的变化而引起的起伏振荡，可以增加过渡点，加宽过渡带以减小通带的起伏。

每一个抽样值产生一个与 成正比，并位移的频率响应，而H（k）与内插函数的线性组合就是FIR滤波器的频率响应，增加一点过渡可以使阻带衰减提高到-44~54dB，二点过渡衰减-65~75，三点过渡衰减-85~95dB。

如果不能使过渡带太宽，同时要求增大阻带衰减，可以增加取样点数N，但这样会增加计算量、延时和误差。

频率取样型FIR滤波器设计步骤：

(1)给定理想滤波器频率响应。

(2)根据过渡带宽和阻带衰减确定过渡点数和h(n)的长度N。

 (3-20)

(3)由IFFT计算IDFT得到：             

 (3-21)

设计实例：

频率采样法设计一个带通滤波器，

满足：低阻带边缘：w1s=0.2*；低通带边缘：w1p=0.35*；高通带边缘：w2p=0.65*；高阻带边缘：w2s=0.8*。设计过渡带中的频率样本值为t1和t2，取t1=0.109021，t2=0.59417456。

设计程序如下：

 

M=40;al=(M-1)/2;l=0:M-1;t1=0.109021;t2=0.59417456;Hrs=[zeros(1,5),t1,t2,ones(1,7),t2,t1,zeros(1,9),t1,t2,ones(1,7),t2,t1,zeros(1,4)];k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1;angh=[-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2)];H=Hrs.*exp(j*angh);h=real(ifft(H,M));freqz(h,1,512,1000)

 

实验得幅频相频特性如下图所示：

三、最优化方法设计FIR数字滤波器

MATLAB信号处理工具箱提供了通用的工具函数firls和remez，采用不同的优化方法可以设计最优滤波器和多频带滤波器。

firls是fir1和fir2函数的扩展，其基本设计准则是利用最小二乘法使期望的频率响应和实际的频率响应之间的整体误差最小。

Remez函数则利用Parks-McClellan算法设计线性相位FIR滤波器。

该算法能使期望频率响应和实际频率响应之间的最大误差最小。Firls函数设计的滤波器在频率响应上出现等波纹，因此也叫等波纹滤波器。

函数firls、remez调用的语法规则相同。

1） firls。

其调用格式为：

b=firls(n,f,a)

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

b=firls(n,f,a,’ftype’)

b=firls(n,f,a,w,’ftype’)

该函数用于设计n阶FIR滤波器，其幅频特性由f和a向量确定，f是频率点向量，其范围为0~1，是递增向量；a是指定频率点的幅度响应，与f必须同长度；b为函数返回的滤波器系数，长度为n+1,且具有偶对称关系：b(k)=-b(n+2-k),  k=1,2,…,n+1；

b=firls(n,f,a,w)则使用权系数w给误差加权。

b=firls(n,f,a,w,’ftype’)调用形式中，参数’ftype’用于指定所设计的滤波器类型，ftype=Hilbert，为奇对称的线性相位滤波器。

返回的滤波器系数满足b(k)=-b(n+2-k),  k=1,2,…,n+1；ftype=differentiator，则采用特殊加权技术，生成奇对称的线性相位滤波器，使低频段误差大大小于高频段误差。

2） 函数remez的基本调用格式为：b=remez(n,f,a)。

其中各项意义同函数firls所述。

设计实例：

分别使用函数firls和remez函数设计20阶带通滤波器，阻带频率为ws1=0.3*，ws2=0.7*，通带频率为wp1=0.4*,wp2=0.6*，并画出幅频特性曲线。

1、程序如下所示：

 

f=[0 0.3 0.4 0.6 0.7 1];m=[0 0 1 1 0 0];n=20;b=firls(n,f,m); %firls函数设计FIR数字滤波器[h,w]=freqz(b); bb=remez(n,f,m); %remez函数设计FIR数字滤波器[hh,w]=freqz(bb);axes('position',[0.2 0.2 0.5 0.5]); plot(w/pi,abs(h),'b:',w/pi,abs(hh),'r-');xlabel('frequency');ylabel('magnitude');grid on;

 

实验得幅频相频特性如下图所示：

2、用Simulink仿真对用firls函数设计的带通滤波器进行滤波：

仿真模块如下：

输入信号：s(t)=sin(20πt)+sin(50πt)+sin(90πt)

数字滤波器参数设置：已知滤波器的阶数n=20。

首先在Filter Type中选择bandpass；

在Design Method选项中选择FIR least-squares；

指定Filter Order项中的Specify order为20；

采样频率Fs=100Hz,Fstop1=15HZ，Fpass1=20HZ，Fpass2=30HZ，Fstop2=35HZ。

设置完以后点击窗口下方的Design Filter，在窗口上方就会看到所设计滤波器的幅频响应曲线。

设计界面如下图所示：

信号源参数设置：三个信号源的频率分别是10Hz、25Hz、45Hz ：

运行仿真，得波形如下图所示：

（1）10HZ信号波形如下所示：

（2） 25HZ信号波形如下所示：        

（3）45HZ信号波形如下所示：

（4）滤波前波形如下所示：

（5）滤波后波形如下所示：

由以上波形可知，10Hz、25Hz、45Hz叠加后通过通频带为20Hz~30Hz的FIR带通滤波器，输出与25Hz信号相似，低频和高频部分被阻隔。

第三篇就到这里，本篇在数字滤波器的设计过程中，采用的设计方法是基于MATLAB的数字滤波器的设计。设计过程中充分发挥了MATLAB语言编程效率高的特点以及Simulink仿真能力强的特点，完成各种计算与图形绘制工作，使得数字滤波器的设计变得简便快捷。通过FIR数字滤波器的设计实例，说明如何利用MATLAB来完成数字滤波器的设计，设计结果的各项性能指标均达到指定要求，Simulink仿真验证了设计方法的正确性。本篇所述设计过程简便易行，具有很好的实际运用意义。

审核编辑 ：李倩