FIR滤波器的设计.doc

______________________________________________________________________________________________________________
基于窗函数法FIR数字滤波器的设计

陈敏敏
宁夏大学电子信息工程 专业09级 12009243656
摘要 简述了数字滤波器中的有限长单位冲激响应(FIR)滤波器的原理,对FIR滤波器的窗函数设计方法进行了研究。窗函数法在FIR 数字滤波器的设计中有着广泛的应用。通过MATLAB 的仿真实现。传统的数字滤波器设计方法繁琐且结果不直观,本文利用MATLAB具有强大的科学计算和图形显示这一优点,与窗函数法设计理论相结合共同设计FIR数字滤波器,不但使设计结果更加直观,而且提高了滤波器的设计精度,从而更好地达到预期效果。

关键字 FIR数字滤波器 窗函数 MATLAB仿真

目 录

1 引言 1
2 FIR数字滤波器的介绍 2
2.1 FIR数字滤波器的特点 2
2.2线性相位FIR数字滤波器的特点 2
2.2.1　单位冲激响应h(n)的特点 2
2.2.2　线性相位的条件 2
2.2.3　线性相位特点和幅度函数的特点 2
2.3 FIR数字滤波器的设计原理 4
2.4 数字滤波器的性能指标 5
3窗函数设计法 6
3.1窗函数设计原理分析 6
3.2设计方法 7
3.3窗函数介绍 9
3.4窗函数法设计步骤 13
3.5设计实例 13
3.6窗函数法计算中的主要问题 14
4 MATLAB简介与数字滤波器的MATLAB实现 19
4.1 MATLAB简介 19
4.2 MATLAB程序 20
结论 26
谢辞 26
参考文献 26
附录 27

精品资料
1 引言

数字信号处理(DSP,digital signal processing)是从20世纪60年代以来，随着信息学科和计算机的高速发展而迅速发展起来的一门新兴学科。数字信号处理是把信号用数字或符号表示的序列，通过计算机或通用(专用)信号处理设备，用数字的数值计算方法处理(例如滤波、变换、压缩、增强、估计、识别等)，以达到提取有用信息便于应用的目的。
数字滤波是数字信号处理的一部分。数字滤波器按照单位取样响应h(n)的时域特性可以分为无限脉冲响应(IIR)系统和有限脉冲响应(FIR)系统。FIR 数字滤波器的优点在于它可以做成具有严格线性相位,而同时可以具有任意的幅度特性;它的传递函数没有极点;这保证了设计出的FIR 数字滤波器一定是平稳的。
所谓数字滤波器设计，简单地说，就是要找到一组能满足特定滤波要求的系数向量a和b。而滤波器设计完成后还需要进一步考虑如何将其实现，即选择什么样的滤波器结构来完成滤波运算。FIR数字滤波器的设计方法很多，其中较为常用的是窗函数设计法、频率采样设计法和最优化设计法。本文讨论利用窗函数法来实现各种FIR滤波器的设计。
窗函数法设计的基本思想是把给定的频率响应通过IDTFT(Inverse Discrete Time Fourier Transform)，求得脉冲响应，然后利用加窗函数对它进行截断和平滑，以实现一个物理可实现且具有线性相位的 FIR 数字滤波器的设计目的。其核心是从给定的频率特性，通过加窗确定有限长单位取样响应h(n)。
MATLAB软件是由美国Math works公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB新的版本集中了日常数学处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。在MATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。 MATLAB提供了一个人机交互的数学系统环境，该系统的基本数据结构是矩阵，在生成矩陈对象时，不要求作明确的维数说明。与利用c语言或FORTRAN语言作数值计算的程序设计相比，利用MATLAB可以节省大量的编程时间。在工程技术界，MATLAB被用来解决一些实际课题和数学模型问题。典型的应用包括数值计算、算法预设计与验证，以及一些特殊的短阵计算应用，如自动控制理论、统计、数字信号处理(时间序列分拆)等。

2 FIR数字滤波器的介绍

2.1 FIR数字滤波器的特点
数字信号处理主要是研究用数字或符号的序列来表示信号波形，并用数字的方式去处理这些序列，把它们改变成在某种意义上更为希望的形式，以便估计信号的特征参量，或削弱信号中的多余分量和增强信号中的有用分量。
有限长单位冲激响应(FIR)数字滤波器可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。此外，FIR滤波器的单位抽样响应是有限长的，因而滤波器一定是稳定的。再有，只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而总能用因果系统来实现。最后，FIR滤波器由于单位冲激响应是有限长的，可以用快速傅立叶变换(FFT)算法来实现过滤信号，从而可大大提高运算效率。但是，要取得很好的衰减特性，FIR滤波器H(z)的阶次比IIR滤波器的要高。
2.2线性相位FIR数字滤波器的特点
2.2.1　单位冲激响应h(n)的特点
FIR滤波器的单位冲激响应h(n)是有限长(0≤n≤N-1),其Z变换为：

在有限Z平面有(N-1)个零点,而它的(N-1)个极点均位于原点z=0处。
2.2.2　线性相位的条件
如果FIR滤波器的单位抽样响应h(n)为实数而且满足以下任一条件:
偶对称：h(n)=h(N-1-n)
奇对称：h(n)=-h(N-1-n)
其对称中心在n=(N-1)/2处,则滤波器具有准确的线性相位。
2.2.3　线性相位特点和幅度函数的特点
(1) h(n)偶对称

幅度函数H(ω)包括正负值,相位函数是严格线性相位,说明滤波器有(N-1)/2个抽样的延时,它等于单位抽样响应h(n)长度的一半。图2-1中，线性相位无90°附加相移,幅度函数在π处存在零点,且对ω=π呈奇对称,因此不适合作高通滤波器。图2-2所示：线性相位无90°附加相移,幅度函数对在ω=0、π、2π呈偶对称,因此适合作低通、高通滤波器。
(2) h(n)奇对称

相位函数仍是线性,但在零频率(ω=0)处有π/2的截距。不仅有(N-1)个抽样的延时,还产生一个π/2的相移。
图2-3中,线性相位有90°附加相移,幅度函数在0、2π处为零点,且对ω=0、2π呈奇对称,对ω=π呈偶对称。
图2-4中,线性相位有90°附加相移,幅度函数在0、π、2π处为零,且对ω=0、π、2π呈奇对称。图2-3、图2-4所示的滤波器均适合在微分器和90°移相器中应用。

图2-1 长度N为偶数、偶对称时的相位函数、冲激响应、幅度函数波形图

图2-2 长度N为奇数、偶对称时的相位函数、冲激响应、幅度函数波形图

图2-3 长度N为偶数、奇对称时的相位函数、冲激响应、幅度函数波形图

图2-4 长度N为奇数、奇对称时的相位函数、冲激响应、幅度函数波形图

四种线性相位FIR滤波器的特性可以总结如下：
第一种情况，偶对称、奇数点，四种滤波器都可设计；
第二种情况，偶对称、偶数点，可设计低、带通滤波器，不能设计高通和带阻；
第三种情况，奇对称、奇数点，只能设计带通滤波器，其它滤波器都不能设计；
第四种情况，奇对称、偶数点，可设计高、带通滤波器，不能设计低通和带阻。
2.3 FIR数字滤波器的设计原理
一个截止频率为(rad/s)的理想数字低通滤波器，其传递函数的表达式是：　　
(式2.3.1)
由式2.3.1可以看出，这个滤波器在物理上是不可实现的，因为冲激响应具有无限性和因果性。为了产生有限长度的冲激响应函数，我们取样响应为，长度为N，其系数函数为：
　　 (式2.3.2)
用表示截取后冲激响应，即，式子中为窗函数，长度为N。当τ=(N-1)/2时，截取的一段对(N-1)/2对称，可保证所设计的滤波器具有线性相位。
一般来说，FIR数字滤波器输出的Z变换形式Y(z)与输入的Z变换形式之间的关系如下：
　　 (式2.3.3)
从上面的Z变换和结构图可以很容易得出FIR滤波器的差分方程表示形式。对式2.3.3进行反Z变换，可得：
　　 (式2.3.4)


图2-5卷积型滤波器
式(2.3.4)为FIR数字滤波器的时域表示方法，其中是在时间n的滤波器的输入抽样值。根据式(2.3.4)即可对滤波器进行设计。从上面的公式我们可以看出，在对滤波器实际设计时，整个过程的运算量很大。设计完成后对已设计的滤波器的频率响应进行校核，运算量也很大。并且在数字滤波器设计的过程中，要根据设计要求和滤波效果不断地调整，以达到设计的最优化。在这种情况下，要进行大量复杂的运算，单纯靠公式计算和编制简单的程序很难在短时间内完成。而利用MATLAB工具进行计算机辅助设计，则可以快速有效地设计数字滤波器，大大的减少了计算量。
2.4 数字滤波器的性能指标
我们在进行滤波器设计时，需要确定其性能指标。一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波器特性为例，频率响应有通带、过渡带及阻带三个范围。
在通带内：  1- AP≤ ≤1 ≤
在阻带中：       ≤ ≤≤
其中为通带截止频率, 为阻带截止频率，Ap为通带误差, 为阻带误差。

图2-6 低通滤波器的幅度特性
与模拟滤波器类似，数字滤波器按频率特性划分为低通、高通、带通、带阻、全通等类型，由于数字滤波器的频率响应是周期性的，周期为2π。
由于频率响应的周期性，频率变量以数字频率来表示，所以数字滤波器设计中必须给出抽样频率。图2-7为各种数字滤波器理想幅度，可以看出：
1、 一个高通滤波器相当于一个全通滤波器减去一个低通滤波器。
2、 一个带通滤波器相当于两个低通滤波器相减。
3、 一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。
这里的相加相减都是相当于并联结构。
图2-7中所示的各种数字滤波器理想频率响应只表示了正频率部分，这样的理想频率响应是不可能实现的，原因是频带之间幅度响应是突变的，因而其单位抽样响应是非因果的。因此要给出实际逼近容限。数字滤波器的系统函数，它在z平面单位圆上的值为滤波器频率响应，表征数字滤波器频率响应特征的三个参量是幅度平方响应、相位响应和群延时响应。
低通
高通
带通
带阻
全通

图2-7 各种理想数字滤波器的幅度频率响应.

3窗函数设计法

3.1窗函数设计原理分析
设数字滤波器的传输函数为,是与其对应的单位脉冲响应, 为系统函数。
(式3.1.1)
(式3.1.2)
(式3.1.3)
一般说来, 是无限长的,需要求对的一个逼近。采用窗函数设计法时,可通过对理想滤波器的单位采样响应加窗设计滤波器
　　　 (式3.1.4)
其中, 是一个长度有限的窗,在区间0 ≤ n ≤ N外值为0 ,且关于中间点对称
　　　 (式3.1.5)
频率响应根据(式3.1.5) ,由卷积定理得出
　　　 (式3.1.6)
理想的频率响应被窗函数的离散时间傅立叶变换“平滑”了。
　采用窗函数设计法设计出来的滤波器的频率响应对理想响应的逼近程度,由两个因素决定:①主瓣的宽度;②旁瓣的幅度大小。
理想的情况是主瓣的宽度窄,旁瓣的幅度小。但对于一个长度固定的窗函数来说,这些不能独立地达到最小。窗函数的一些通用性质为:
1、窗函数的长度N增加,主瓣的宽度减小,使得过渡带变小。
关系为:NB = C其中:B是过渡带的宽度;C是取决于窗函数的一个参数。如矩形窗为4π。调整N可以有效地控制过渡带的宽度,但N的改变不改变主瓣和旁瓣的相对比例。随着N值增加,过渡带变窄,波动频率也随着增加,虽然总的幅度有所减少,但截止频率附近的肩峰并不减少,而只是随着N值的增加,肩峰被抑制在愈来愈小的范围内,使肩峰宽度变窄。
2、窗函数的旁瓣的幅度大小取决于窗函数的选择。选择恰当的窗函数使主瓣包含更多的能量,相应旁瓣的幅度就减小。旁瓣幅度的减小,可以减少通带和阻带的波动,使通带尽可能趋近水平,阻带尽可能达到最大衰减。但通常此时过渡带会变宽。
3、取不同的窗函数对幅度特性的整形效果比单纯的增加窗口长度要强得多。
3.2设计方法
这种方法也叫傅里叶级数法。一般是先给出所要求的理想的滤波器的频率响应，要求设计一个FIR滤波器频率响应来逼近。设计是在时域进行的，因而先由的傅里叶反变换导出，即
(式3.2.1)
由于是矩形频率响应特性，故一定是无限长序列，且是非因果的，而FIR滤波器的必然是有限长的，所以要用有限长的来逼近无限长的，最有效的方法是截断或者说用一个有限长度的窗口函数序列来截取，即
(式3.2.2)
因而窗函数序列的形状及长度的选择就是关键。
我们以一个截止频率为的线性相位的理想矩形幅度特性的低通滤波器为例来讨论。设低通特性的群延时为，即
(式3.2.3)
这表明，在通带≤范围内，的幅度是均匀的，其值为1，相位是。利用(1)式可得
(式3.2.4)
是中心点在的偶对称无限长非因果序列，要得到有限长的，一种最简单的方法就是取矩形窗，即

但是按照线形相位滤波器的约束，必须是偶对称的，对称中心应为长度的一半(N-1)/2，因而必须=(N-1)/2，所以有
(式3.2.5)
将(式3.2.4)代入(式3.25)，可得
(式3.2.6)
此时，一定满足这一线性相位的条件。
下面求的傅里叶变换，也就是找出待求FIR滤波器的频率特性，以便能看出加窗处理后究竟对频率响应有何影响。
按照复卷积公式，在时域是相乘、频域上是周期性卷积关系，即
(式3.2.7)
因而逼近的好坏，完全取决于窗函数的频率特性。
窗函数的频率特性为
(式3.2.8)
对矩形窗，则有
(式3.2.9)
也可表示成幅度函数与相位函数
(式3.2.10)
其中
(式3.2.11)
就是频域抽样内插函数，其幅度函数在之内为一个主瓣，两侧形成许多衰减振荡的旁瓣，如果将理想频率响应也写成
(式3.2.12)
则其幅度函数为
(式3.2.13)
3.3窗函数介绍
实际应用的窗函数，可分为以下主要类型：
1、幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间(t)的高次幂；
2、三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁窗、海明窗等；
3、指数窗--采用指数时间函数，如形式，例如高斯窗等。
    下面介绍几种常用窗函数的性质和特点。
（1）矩形窗
矩形窗属于时间变量的零次幂窗，函数形式为:
         (式 3.3.1)
相应的窗谱为：
(式 3.3.2)
矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。

图3-1矩形窗的时域及频域波形
（2）三角窗亦称费杰(Fejer)窗，是幂窗的一次方形式，其函数形式是：
         (式 3.3.3)
三角窗与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣，如图3-2所示。
图3-2三角窗的时域及频域波形
（3）汉宁(Hanning)窗
汉宁窗又称升余弦窗，其时域表达式为：
         (式3.3.4)
相应的窗谱为：
         (式3.3.5)
由此式可以看出，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是 3个 sin(t)型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗。但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降。
（4）海明(Hamming)窗
海明窗也是余弦窗的一种，又称改进的升余弦窗，其时间函数表达式为：
         (式3.3.6)
其窗谱为：
         (式3.3.7)
海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更小。分析表明，海明窗的第一旁瓣衰减为-42dB。海明窗的频谱也是由3个矩形窗的频谱合成，但其旁瓣衰减速度为20dB／(10oct)，这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。
（5）高斯窗
    高斯窗是一种指数窗。其时域函数为：
         (式3.3.8)
式中a为常数，决定了函数曲线衰减的快慢。a值如果选取适当，可以使截断点(T为有限值)处的函数值比较小，则截断造成的影响就比较小。高斯窗谱无负的旁瓣，第一旁瓣衰减达一55 dB。高斯窗的主瓣较宽，故而频率分辨率低。高斯窗函数常被用来截断一些非周期信号，如指数衰减信号等。不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。图3-3是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。

图3-3几种常用的窗函数的时域和频域波形
对于窗函数，还有一些要求：
1）3dB带宽B，它是主瓣归一化的幅度下降到-3dB时的带宽。当数据长度为N时，矩形窗主瓣两个过零点之间的宽度为4π/N。
2）最大边瓣峰值A(dB)。
3）边瓣谱峰渐进衰减速度D(dB/cot)。
所以，理想的窗函数应当具有最小的B和A，和最大的D。
3.4窗函数法设计步骤
1、首先是给定所要求的频率响应函数；
2、其次，求单位冲激响应；
3、再次，有过渡带宽及阻带最小衰减的要求，查表选定窗函数及N的大小，一般N的大小要通过几次试探而后确定；
4、求得所设计的FIR滤波器的单位冲激响应；
，n=0，1，…,N-1；
5、 求，检验是否满足设计要求，如不满足，则需要重新设计。
3.5设计实例
线性相位FIR低通滤波器的设计(用窗函数法)。指标要求:通带截止频率:0.2π,阻带起始频率:0.4π,阻带最小衰减:-50dB。

(1)设为理想线性相位滤波器

由所需低通滤波器的过渡带求出理想低通滤波器的截止数字频率ω=0.3π,得出:

(2)由阻带衰减确定窗函数,由过渡带宽确定N值。阻带最小衰减50dB,比对6种窗函数基本参数选定窗函数为海明窗。
所要求的过渡带宽:
Δω=0.4π-0.2π=0.2π
N=6.6π/0.2π=33,τ=(N-1)/2=16
(3)由海明窗函数确定FIR滤波器的h(n)。

得出：

(4)仿真检验各项指标,得出结论:满足设计要求。
取N=33,偶对称,
得:过渡带宽Δω:0.3476563π,
第一通带波纹:0.020837dB,
第一阻带最小衰减:60.9159dB。
我们设计的低通FIR数字滤波器的理想性能指标是：通带截止频率Wp=0.3π,阻带截止频率Ws=0.5π， 阻带衰减At不小于40dB， 通带衰减不大于3dB, fs=4000。
Matlab 程序
wp=0.3*pi;
ws=0.5*pi;
wdelta=ws-wp;
N=ceil(8*pi/wdelta);
Wn=(0.3+0.5)*pi/2;
b=fir1(N,Wn/pi,hanning(N+1));
freqz(b,1,512)
N=ceil(8/0.15);
n=0:N-1;
window=hanning(N);
[h1,w]=freqz(window,1);
figure(1);
stem(window);
axis([0 60 0 1.2]);
grid;
xlabel('n');
title('Hanning窗函数');
figure(2);
plot(w/pi,20*log(abs(h1)/abs(h1(1))));
axis([0 1 -350 0]);
grid;
xlabel('w/pi');
ylabel('幅度(db)');
title('Hanning窗函数的频谱');
hn=fir1(N-1,wc,hanning(N));
[h2,w]=freqz(hn,1,512);
figure(3);
stem(n,hn);
axis([0 60 -0.25 0.25]);
grid;
xlabel('n');
ylabel('h(n)');
title('Hanning窗函数的单位脉冲响应');
figure(4);
plot(w/pi,20*log(abs(h2)/abs(h2(1))));
grid;
xlabel('w/pi');
ylabel('幅度（db）');
>>
运行结果：

Welcome To
Download !!!

欢迎您的下载，资料仅供参考！