FIR滤波器设计 下载本文

内容发布更新时间 : 2024/10/22 23:48:24星期一 下面是文章的全部内容请认真阅读。

第7章 FIR滤波器设计

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

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

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

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

Y(z)N?1??h(n)z?n (7-1) H(z)?X(z)n?0可得FIR滤波器的系统差分方程为:

y(n)?b(0)x(n)?b(1)x(n?1)???b(N?1)x(n?N?1)

??b(m)x(n?m)?b(n)?x(n)m?0N?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滤波器设计的主要方法

— 192 —

函数设计方法 说明 理想滤波器加窗处理 平方误差最小化逼近理想幅频响应或Park-McClellan 算法产生等波纹滤波器 在满足最大误差限制条件下使整个频带平方误差最小化 具有光滑、正弦过渡带的低通滤波器设计 工具函数 fir1(单频带) , fir2(多频带) , kaiserord firls , remez,remezord 窗函数法 最优化设计 约束最小二乘逼近 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??12?????Hdej?ej?nd? (7-5)

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

— 193 —

Hd(z)?n????h?d(n)z?n (7-6)

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

?0,n?0,n?N w(n)??1,0?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(右图)。

— 194 —

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

由图可看出,加矩形窗后使实际频率响应偏离理想频率响应,主要影响有三个方面: (1)理想幅频特性陡直边缘处形成过渡带,过渡带宽取决于矩形窗函数频率响应的主瓣宽度。

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

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

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

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

x(t)?w(t)?X(f)?W(f) (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)

汉宁窗的表达式为:

w(k)?0.5?1?cos?2?????k????,k?1,...,N (7-10) N?1??(3)哈明窗:w=hamming(N) 哈明窗的表达式为:

w(k?1)?0.54?0.46cos?2?(4)Bartlett窗:w=bartlett(N)

— 195 —

??k??,k?0,1,...,N?1 (7-11) N?1?Bartlett 窗的表达式为:

当 N 为奇数时,

N?1?2(k?1),1?k??N?12 w(k)?? (7-12)

2(k?1)N?1?2?,?k?NN?12?当 N 为偶数时,

2(k?1)N?,1?k??N?12 w(k)?? (7-13)

2(N?k)N?2?,?1?k?NN?12?(5) Blackman 窗:w= blackman(N)

Blackman 窗的表达式为: w(k)?0.42?0.5cos?2???k?1?k?1????0.08?4??,k?1,...,N (7-14) N?1?N?1??Blackman 窗比其他相同尺寸窗 (哈明窗,汉宁窗) 具有主瓣较宽和旁瓣泄漏较小的特

点。

(6)三角窗:w=triang(N) 三角窗的表达式为: 当 N 为奇数时,

2kN?1?,1?k??N?12 w(k)?? (7-15)

2(N?k?1)N?1?,?k?N2?N?1当 N 为偶数时,

2k?1N?,1?k??N?12 w(k)?? (7-16)

2(N?k?1)N?2?,?1?k?NN2?三角窗和Bartlett窗十分类似。三角窗的两端值不为零,而Bartlett窗则为零,这一

点可从例7-1中看出。

(7)Kaiser窗:w=kaiser(n,beta)

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

Kaiser窗表达式:

2?2k???I0??1??1???N?1?????? w(k)? (7-17)

I0???式中, I0[.]是修正过的零阶 Bessel 函数。

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

— 196 —