基于窗函数法FIR数字滤波器的设计（毕业设计）

??

重庆三峡学院毕业设计（论文）

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

专 业 电子信息工程 年 级 2 0 0 6 级 学生姓名 郝 海 学生学号 200615190154

指导教师 徐正坤 职称 助教

完成毕业设计（论文）时间2009年12月

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

郝海

重庆三峡学院应用技术学院电子信息工程（仪器仪表）专业06级 重庆万州 404100

摘要 简述了数字滤波器中的有限长单位冲激响应(FIR)滤波器的原理,对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 FIR数字滤波器的优化 ------------------------------------------------------------- 15 4.1均方误差最小化准则 ----------------------------------------------------------- 15 4.2切比雪夫最佳一致逼近定理 ----------------------------------------------------- 16 4.3利用切比雪夫逼近理论设计FIR数字滤波器 --------------------------------------- 16 4.4瑞米兹算法 ------------------------------------------------------------------- 17 5 MATLAB简介与数字滤波器的MATLAB实现 --------------------------------------------- 19 5.1 MATLAB简介 ------------------------------------------------------------------ 19 5.2 MATLAB程序 ------------------------------------------------------------------ 20 结论 ------------------------------------------------------------------------------ 26 谢辞 ------------------------------------------------------------------------------ 26 参考文献 -------------------------------------------------------------------------- 26 附录 ------------------------------------------------------------------------------ 27

2010届电子信息工程专业毕业设计（论文）

1 引言

数字信号处理(DSP,digital signal processing)是从20世纪60年代以来，随着信息学科和计算机的高速发展而迅速发展起来的一门新兴学科。数字信号处理是把信号用数字或符号表示的序列，通过计算机或通用(专用)信号处理设备，用数字的数值计算方法处理(例如滤波、变换、压缩、增强、估计、识别等)，以达到提取有用信息便于应用的目的。

数字滤波是数字信号处理的一部分。数字滤波器按照单位取样响应h(n)的时域特性可以分为无限脉冲响应(IIR)系统和有限脉冲响应(FIR)系统。FIR 数字滤波器的优点在于它可以做成具有严格线性相位,而同时可以具有任意的幅度特性;它的传递函数没有极点;这保证了设计出的FIR 数字滤波器一定是平稳的。

所谓数字滤波器设计，简单地说，就是要找到一组能满足特定滤波要求的系数向量a和b。而滤波器设计完成后还需要进一步考虑如何将其实现，即选择什么样的滤波器结构来完成滤波运算。FIR数字滤波器的设计方法很多，其中较为常用的是窗函数设计法、频率采样设计法和最优化设计法。本文讨论利用窗函数法、均方误差最小化法和等波纹切比雪夫逼近法(调用remez函数)来分别实现各种FIR滤波器的设计。

窗函数法设计的基本思想是把给定的频率响应通过IDTFT(Inverse Discrete Time Fourier Transform)，求得脉冲响应，然后利用加窗函数对它进行截断和平滑，以实现一个物理可实现且具有线性相位的 FIR 数字滤波器的设计目的。其核心是从给定的频率特性，通过加窗确定有限长单位取样响应h(n)。

均方误差最小化法是等波纹切比雪夫逼近法的基础，这一准则是使误差能量最小，但是由于吉布斯效应，窗谱的肩峰过大，造成所设计出的滤波器通带起伏不均匀且过大，而阻带衰减则过小，不能满足要求。

等波纹切比雪夫逼近法，一致逼近法的原理即为切比雪夫最佳一致逼近法则，也可称为等波纹逼近。其目的是在所需要的区间内，使误差函数E(x)较均匀一致，并且通过合理地选择多项式，使E(x)的最大值达到最小。通俗的讲，就是使最大误差最小化。应用切比雪夫理论，提出了一种FIR数字滤波器的计算机辅助设计方法。这种方法可获得很好的通带和阻带性能，并能准确地指定通带和阻带的边缘，是一种有效的设计方法。最后利用 MATLAB 提供的Remes函数实现Remes算法，设计滤波器逼近理想频率响应。

MATLAB软件是由美国Math works公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB新的版本集中了日常数学处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。在MATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。 MATLAB提供了一个人机交互的数学系统环境，该系统的基本数据结构是矩阵，在生成矩陈对象时，不要求作明确的维数说明。与利用c语言或FORTRAN语言作数值计算的程序设计相比，利用MATLAB可以节省大量的编程时间。在工程技术界，MATLAB被用来解决一些实际课题和数学模型问题。典型的应用包括数值计算、算法预设计与验证，以及一些特殊的短阵计算应用，如自动控制理论、统计、数字信号处理(时间序列分拆)等。

1

2010届电子信息工程专业毕业设计（论文）

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变换为：

N?1H(z)??h(n)zm?0?m

在有限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)偶对称

2

2010届电子信息工程专业毕业设计（论文）

k?1H(?)??h(n)cos[(n?0N?12?n)?]

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

(2) h(n)奇对称

k?1H(?)??n?0h(n)sin[()??N?12?n)?]

?(?)??(N?12?2相位函数仍是线性,但在零频率(ω=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为奇数、偶对称时的相位函数、冲激响应、幅度函数波形图

3

2010届电子信息工程专业毕业设计（论文）

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

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

四种线性相位FIR滤波器的特性可以总结如下： 第一种情况，偶对称、奇数点，四种滤波器都可设计；

第二种情况，偶对称、偶数点，可设计低、带通滤波器，不能设计高通和带阻； 第三种情况，奇对称、奇数点，只能设计带通滤波器，其它滤波器都不能设计； 第四种情况，奇对称、偶数点，可设计高、带通滤波器，不能设计低通和带阻。

2.3 FIR数字滤波器的设计原理

一个截止频率为?c(rad/s)的理想数字低通滤波器，其传递函数的表达式是：

?j???e,???c?j? (式2.3.1) Hd(e)????0,?c????由式2.3.1可以看出，这个滤波器在物理上是不可实现的，因为冲激响应具有无限性和因果性。为了产生有限长度的冲激响应函数，我们取样响应为h(n)，长度为N，其系数函数为H(z)：

N?1H(z)?

?h(n)zn?0?n (式2.3.2)

用h(n)表示截取hd(n)后冲激响应，即h(n)??(n)hd(n)，式子中?(n)为窗函数，长度为N。当τ=(N-1)/2时，截取的一段h(n)对(N-1)/2对称，可保证所设计的滤波器具有线性相位。

一般来说，FIR数字滤波器输出y(n)的Z变换形式Y(z)与输入x(n)的Z变换形式之间的关系如下：

Y(z)?H(z)X(z)?(h(0)?h(1)z?1???h(n)z?n)X(z) (式2.3.3)

从上面的Z变换和结构图可以很容易得出FIR滤波器的差分方程表示形式。对式2.3.3进行反Z变换，可得：

4

2010届电子信息工程专业毕业设计（论文）

y(n)?h(1)x(n)?h(2)x(n?1)??h(n)x(1) (式2.3.4)



图2-5卷积型滤波器

式(2.3.4)为FIR数字滤波器的时域表示方法，其中x(n)是在时间n的滤波器的输入抽样值。根据式(2.3.4)即可对滤波器进行设计。从上面的公式我们可以看出，在对滤波器实际设计时，整个过程的运算量很大。设计完成后对已设计的滤波器的频率响应进行校核，运算量也很大。并且在数字滤波器设计的过程中，要根据设计要求和滤波效果不断地调整，以达到设计的最优化。在这种情况下，要进行大量复杂的运算，单纯靠公式计算和编制简单的程序很难在短时间内完成。而利用MATLAB工具进行计算机辅助设计，则可以快速有效地设计数字滤波器，大大的减少了计算量。

2.4 数字滤波器的性能指标

我们在进行滤波器设计时，需要确定其性能指标。一般来说，滤波器的性能要求往往以频率响应的幅度特性的允许误差来表征。以低通滤波器特性为例，频率响应有通带、过渡带及阻带三个范围。

在通带内： 1- AP≤H(ej?) ≤1 ?c≤?c 在阻带中： H(ej?) ≤Ast ?st≤?≤?c

其中?c为通带截止频率, ?st为阻带截止频率，Ap为通带误差, Ast为阻带误差。

图2-6 低通滤波器的幅度特性

与模拟滤波器类似，数字滤波器按频率特性划分为低通、高通、带通、带阻、全通等类型，由于数字滤波器的频率响应是周期性的，周期为2π。

由于频率响应的周期性，频率变量以数字频率?来表示，所以数字滤波器设计中必须给出抽样频率。图2-7为各种数字滤波器理想幅度，可以看出：

1、 一个高通滤波器相当于一个全通滤波器减去一个低通滤波器。 2、 一个带通滤波器相当于两个低通滤波器相减。

3、 一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。

5

2010届电子信息工程专业毕业设计（论文）

这里的相加相减都是相当于并联结构。

图2-7中所示的各种数字滤波器理想频率响应只表示了正频率部分，这样的理想频率响应是不可能实现的，原因是频带之间幅度响应是突变的，因而其单位抽样响应是非因果的。因此要给出实际逼近容限。数字滤波器的系统函数H(z)，它在z平面单位圆上的值为滤波器频率响应

H(ej?)，表征数字滤波器频率响应特征的三个参量是幅度平方响应、相位响应和群延时响应。

低通

Hd(ej?)

Hd(e高通

j?)

Hd(e带通

j?)

Hd(e带阻

j?

)

Hd(e全通

j?)

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

3窗函数设计法

3.1窗函数设计原理分析

设数字滤波器的传输函数为H(e数。

N?1j?),hd(n)是与其对应的单位脉冲响应, H(z) 为系统函

H(ej?)?1?h(n)en?0?j?n (式3.1.1)

hd(n)?2???H??d(ej?)ej?nd? (式3.1.2)

6

2010届电子信息工程专业毕业设计（论文）

N?1H(z)??h(n)zn?0?n (式3.1.3)

一般说来, hd(n)是无限长的,需要求对Hd(ej?)的一个逼近。采用窗函数设计法时,可通过对理想滤波器的单位采样响应加窗设计滤波器

h(n)??(n)hd(n) (式3.1.4)

其中, ?(n)是一个长度有限的窗,在区间0 ≤ n ≤ N外值为0 ,且关于中间点对称

?(n)??(N?1?n) (式3.1.5)

频率响应根据(式3.1.5) ,由卷积定理得出

H(ej?)?12?

Hd(ej?)??(ej?) (式3.1.6)

理想的频率响应被窗函数的离散时间傅立叶变换?(ej?)“平滑”了。

采用窗函数设计法设计出来的滤波器的频率响应对理想响应Hd(ej?)的逼近程度,由两个因素决定:①?(ej?)主瓣的宽度;②?(ej?)旁瓣的幅度大小。

理想的情况是?(ej?)主瓣的宽度窄,旁瓣的幅度小。但对于一个长度固定的窗函数来说,这些不能独立地达到最小。窗函数的一些通用性质为:

1、窗函数的长度N增加,主瓣的宽度减小,使得过渡带变小。

关系为:NB = C其中:B是过渡带的宽度;C是取决于窗函数的一个参数。如矩形窗为4π。调整N可以有效地控制过渡带的宽度,但N的改变不改变主瓣和旁瓣的相对比例。随着N值增加,过渡带变窄,波动频率也随着增加,虽然总的幅度有所减少,但截止频率附近的肩峰并不减少,而只是随着N值的增加,肩峰被抑制在愈来愈小的范围内,使肩峰宽度变窄。

2、窗函数的旁瓣的幅度大小取决于窗函数的选择。选择恰当的窗函数使主瓣包含更多的能量,相应旁瓣的幅度就减小。旁瓣幅度的减小,可以减少通带和阻带的波动,使通带尽可能趋近水平,阻带尽可能达到最大衰减。但通常此时过渡带会变宽。

3、取不同的窗函数对幅度特性的整形效果比单纯的增加窗口长度要强得多。 3.2设计方法

这种方法也叫傅里叶级数法。一般是先给出所要求的理想的滤波器的频率响应Hd(eN?1j?)，要

求设计一个FIR滤波器频率响应H(ej?j?)??h(n)en?0?j?n来逼近Hd(ej?)。设计是在时域进行的，

因而先由Hd(e)的傅里叶反变换导出hd(n)，即

hd(n)?12?????Hd(ej?)ej?nd? (式3.2.1)

由于Hd(ej?)是矩形频率响应特性，故hd(n)一定是无限长序列，且是非因果的，而FIR滤波器

的h(n)必然是有限长的，所以要用有限长的h(n)来逼近无限长的hd(n)，最有效的方法是截断

hd(n)或者说用一个有限长度的窗口函数序列?(n)来截取hd(n)，即

7

2010届电子信息工程专业毕业设计（论文）

h(n)??(n)hd(n) (式3.2.2)

因而窗函数序列的形状及长度的选择就是关键。

我们以一个截止频率为?c的线性相位的理想矩形幅度特性的低通滤波器为例来讨论。设低通特性的群延时为?，即

Hd(ej??e?j??,??c????c (式3.2.3) )???0,?c????,???????c这表明，在通带?≤?c范围内，Hd(ej?)的幅度是均匀的，其值为1，相位是???。利用(1)式可得

12??cc hd(n)????e?j??ej?nd???csin??c(n??)???c(n??) (式3.2.4)

hd(n)是中心点在?的偶对称无限长非因果序列，要得到有限长的h(n)，一种最简单的方法就是

取矩形窗RN(n)，即

?(n)?RN(n)

但是按照线形相位滤波器的约束，h(n)必须是偶对称的，对称中心应为长度的一半(N-1)/2，因而必须?=(N-1)/2，所以有

??hd(n),0?n?N?1h(n)?h(n)?(n)???d??0,n为其他 (式3.2.5) ?N-1????2?将(式3.2.4)代入(式3.25)，可得

?N?1??sin?(n?)????c2???ch(n)???,0?n?N?1 (式3.2.6) N?1?c(n?)?2??0,n为其他值此时，一定满足h(n)?h(N?1?n)这一线性相位的条件。

下面求h(n)的傅里叶变换，也就是找出待求FIR滤波器的频率特性，以便能看出加窗处理后究竟对频率响应有何影响。

按照复卷积公式，在时域是相乘、频域上是周期性卷积关系，即

H(ej?)?12?????Hd(ej?)ej(???)d? (式3.2.7)

8

2010届电子信息工程专业毕业设计（论文）

因而H(ej?)逼近Hd(ej?)的好坏，完全取决于窗函数的频率特性W(ej?)。 窗函数?(n)的频率特性W(ej?)为

N?1W(ej?)???(n)en?0?j?n (式3.2.8)

对矩形窗RN(n)，则有

sin(sin(?N2N2N?1WR(ej?)??en?0?j?n?e?j?N?12) (式3.2.9) )也可表示成幅度函数与相位函数

?j(N?12)?WN(ej?)?WR(?)e (式3.2.10)

其中

sin(WR(?)?sin(WR(ej??N2N2)) (式3.2.11)

)就是频域抽样内插函数，其幅度函数WR(?)在???2?/N之内为一个主瓣，两侧形成

许多衰减振荡的旁瓣，如果将理想频率响应也写成

j??j(N?12)?Hd(e)?Hd(?)e (式3.2.12)

则其幅度函数为

??1,???cHd(?)?? (式3.2.13)

??0,?c????3.3窗函数介绍

实际应用的窗函数，可分为以下主要类型：

1、幂窗--采用时间变量某种幂次的函数，如矩形、三角形、梯形或其它时间(t)的高次幂； 2、三角函数窗--应用三角函数，即正弦或余弦函数等组合成复合函数，例如汉宁窗、海明窗等； 3、指数窗--采用指数时间函数，如e?st形式，例如高斯窗等。

下面介绍几种常用窗函数的性质和特点。

9

2010届电子信息工程专业毕业设计（论文）

（1）矩形窗

矩形窗属于时间变量的零次幂窗，函数形式为:

?1?,t?T ?(t)??T (式 3.3.1)

?0,t?T?相应的窗谱为：

2sin?TW(?)??T (式 3.3.2)

矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。

图3-1矩形窗的时域及频域波形

（2）三角窗亦称费杰(Fejer)窗，是幂窗的一次方形式，其函数形式是：

?1t?(1?),t?T?(t)??TT?0,t?T? (式 3.3.3)

三角窗与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣，如图3-2所示。

10

2010届电子信息工程专业毕业设计（论文）

图3-2三角窗的时域及频域波形

（3）汉宁(Hanning)窗

汉宁窗又称升余弦窗，其时域表达式为：

相应的窗谱为：

?t?111(?cos),t?T?T?(t)??T22?0,t?T? (式3.3.4)

W(?)?sin?T

?T?1?sin(?T??)sin(?T??)???2??T???T???? (式3.3.5)

由此式可以看出，汉宁窗可以看作是3个矩形时间窗的频谱之和，或者说是 3个 sin(t)型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 π/T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗。但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降。 （4）海明(Hamming)窗

海明窗也是余弦窗的一种，又称改进的升余弦窗，其时间函数表达式为：

其窗谱为：

?t?1(0.54?0.4cos),t?T?TT?(t)???0,t?T? (式3.3.6)

W(?)?1.08sin?T

?T?sin(?T??)sin(?T??)??0.46???T?????T??? (式3.3.7)

11

2010届电子信息工程专业毕业设计（论文）

海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更小。分析表明，海明窗的第一旁瓣衰减为-42dB。海明窗的频谱也是由3个矩形窗的频谱合成，但其旁瓣衰减速度为20dB／(10oct)，这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。

（5）高斯窗

高斯窗是一种指数窗。其时域函数为：

?1?at2,t?T?e?(t)??T?0,t?T? (式3.3.8)

式中a为常数，决定了函数曲线衰减的快慢。a值如果选取适当，可以使截断点(T为有限值)处的函数值比较小，则截断造成的影响就比较小。高斯窗谱无负的旁瓣，第一旁瓣衰减达一55 dB。高斯窗的主瓣较宽，故而频率分辨率低。高斯窗函数常被用来截断一些非周期信号，如指数衰减信号等。不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。图3-3是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。

图3-3几种常用的窗函数的时域和频域波形

对于窗函数，还有一些要求：

1）3dB带宽B，它是主瓣归一化的幅度下降到-3dB时的带宽。当数据长度为N时，矩形窗主瓣两个过零点之间的宽度为4π/N。

12

2010届电子信息工程专业毕业设计（论文）

M?[Hg(?)?k?0]Ckcos(?)-cos(?k)?kM?cos(?)-cos(?k?0?k (式4.4.3)

k)把Hg(?)代入E(e)?W(e)|Hg(e)?Hd(e)|式,可求得误差函数E(?)。

如果在子集F上，对所有的频率?，都有|E(?)|?|?|，这说明，E(?)是纹波的极值，初始猜测的?0，?1，?,?M?1恰是交错点组。这时，设计工作即可以结束。当然，对第一次猜测的位置，不会恰好如此。一般，在某些频率处，总有E(?)??，这说明，需要交换上次猜测的交错点组中的某些点，得到一组新的交错点组。

第二步：对上次确定的交错点组?0，?1，?,?M?1中的每一个点，都在其附近检查时候在某一个频率处有E(?)??，如若有，再在该点附近找出局部极值点，用这一局部极值点代替原来的点。待这M+2个点都检查过后，便得到一组新的交错点组?0，?1，?,?M?1，再利用上式是求出?，Hg(?)，E(?)，这样就完成一次迭代，也就是完成了一次交错点组的交换。通过交换算法，使得这一次的交错点组中的每一个?i都是由上一次的交错点组所产生的E(?)，局部极值频率点，因此，用这次的交错点组求出的?将增大。

第三步：利用和第二步相同的方法，把在各频率处使E(?)??的点作为新的局部极值点，从而又得到一组新的交错点组。重复上述步骤，因为新的交错点组的选择都是作为每一次求出的

E(?)的局部极值点，因此，在迭代中，每次的|?|都是递增的。?最后收敛到自己的上限。也

j?j?j?j?即Hg(?)最佳的一致逼近Hd(?)的解。因此，若再迭代一次，新的误差曲线E(?)的峰值将不会大于|?|，这时迭代可以结束。由最后的交错点组可按

M?[ Hg(?)?k?0]Ckcos(?)-cos(?k)?kM?cos(?)-cos(?k?0?k (式4.4.4)

k)得到Hg(?)，将Hg(?)再附上

MHg(ej?)??a(n)cos(n?) (式4.4.5)

n?0线性相位后作逆变换，便可以得到单位抽样响应h(n)。最后求出的?是阻带的峰值偏差?2，而k?2便是通带的峰值偏差。交错点数?0，?1，?,?M?1是限制在通带和阻带内的交错点组插值产生的。对通带和阻带内的逼近误差E(?)不需要事先指定，而是用切比雪夫最佳一致逼近理论保证了逼近的最大偏差为最小。

通过上述的方法和步骤，最后我们得到了FIR数字滤波器的优化设计的计算机实现结果，如下：

remez(remezord): fedge=[800 1000];

mval=[1 0];dev=[0.0559 0.01]; fs=4000;

18

2010届电子信息工程专业毕业设计（论文）

[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); [h,w]=freqz(b,1,256);

plot(w*2000/pi,20*log10(abs(h))); grid;

xlabel('频率/Hz') ylabel('幅度/dB')

图4-5 瑞米兹设计实现的FIR滤波器

预先知道极值点的最大数目很重要，因为有些设计方法只能设计极点数目有最大可能值的最多波纹滤波器。瑞米兹算法则可以设计任何最优线性相位FIR滤波器。它是一种最为实用的最优化算法。

5 MATLAB简介与数字滤波器的MATLAB实现

5.1 MATLAB简介

MATLAB软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB是英文MATrix LABoratory(短阵实验室)的缩写。它的第1版(DoS版本1.0)发行于1984年，经过十几年不断改进，现今已推出它的WINDWS XP版本(2009b版)。新的版本集中了日常数学处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。在MATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。MATLAB提供了一个人机交互的数学系统环境，该系统的基本数据结构是矩阵，在生成矩陈对象时，不要求作明确的维数说明。与利用C语言或FORTRAN语言作数值计算的程序设计相比，利用MATLAB可以节省大量的编程时间。在工程技术界，MATLAB被用来解决一些实际课题和数学模型问题。典

19

2010届电子信息工程专业毕业设计（论文）

型的应用包括数值计算、算法预设计与验证，以及一些特殊的短阵计算应用，如自动控制理论、统计、数字信号处理(时间序列分拆)等。进入了20世纪90年代，MATLAB已经成为国际控制界公认的标准计算软件。MATLAB因其提供计算、图形、文字处理的统一环境而深受欢迎。

MTALAB系统由五个主要部分组成，下面分别加以介绍。

(1)MATALB语言是高层次的矩阵／数组语言。具有条件控制、函数调用、数据结构、输入输出、面向对象等程序语言特性。利用它既可以进行小规模编程，完成算法设计和算法实验的基本任务，也可以进行大规模编程，开发复杂的应用程序。

(2)MATLAB工作环境这是对MATLAB提供给用户使用的管理功能的总称。包括管理工作空间中的变量据输入输出的方式和方法，以及开发、调试、管理M文件的各种工具。

(3)MATLAB图形系统的基础，包括完成2D和3D数据图示、图像处理、动画生成、图形显示等功能的高层MATLAB命令，也包括用户对图形图像等对象进行特性控制的低层MATLAB命令，以及开发GUI应用程序的各种工具。

(4)MATLAB数学函数库这是对MATLAB使用的各种数学算法的总称。包括各种初等函数的算法，也包括矩阵运算、矩阵分析等高层次数学算法。

(5)MATLAB应用程序接口(API) 这是MATLAB为用户提供的一个函数库，使得用户能够在MATLAB环境中使用C程序或Fortran程序，包括从MATLAB中调用于程序(动态链接)，读写MAT文件的功能。 可以看出MATLAB是一个功能十分强大的系统，是集数值计算、图形管理、程序开发为一体的环境。

除此之外，MATLAB还具有根强的功能扩展能力，与它的主系统一起，可以配备各种各样的工具箱，以完成一些特定的任务。用户可以根据自己的工作任务，开发自己的工具箱。 5.2 MATLAB程序

我们设计的低通FIR数字滤波器的理想性能指标是：通带截止频率Wp=0.3π,阻带截止频率Ws=0.5π， 阻带衰减At不小于40dB， 通带衰减不大于3dB, fs=4000。 boxcar:

wp=0.3*pi; ws=0.5*pi; wdelta=ws-wp; N=ceil(8*pi/wdelta); Wn=(0.3+0.5)*pi/2;

b=firl(N,Wn/pi,hanning(N)); freqz(b,1,512)

20

2010届电子信息工程专业毕业设计（论文）

图5-1 矩形窗的幅频特性曲线

这种窗的主瓣比较集中，应用很普遍，但是旁瓣比较高，并且有负旁瓣。截取长度N增加，会减少过渡带宽4π/N，不会改变肩峰的相对值。 hanning:

wp=0.3*pi; ws=0.5*pi; wdelta=ws-wp; N=ceil(8*pi/wdelta); Wn=(0.3+0.5)*pi/2;

b=firl(N,Wn/pi,hanning(N+1)); freqz(b,1,512)

21

2010届电子信息工程专业毕业设计（论文）

图5-2 汉宁窗的幅频特性曲线

此窗旁瓣互相抵消，能量更集中在了主瓣。但是，它的代价是主瓣宽度比矩形窗的主瓣宽度增加了一倍。 kaiser:

wn=kaiser(28,3.44); nn=[0:1:27]; alfa=(28-1)/2 ;

hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa)); h=hd.*wn';

[h1,w1]=freqz(h,1);

plot(w1/pi,20*log10(abs(h1))); axis([0,1,-80,101]); grid;

xlabel('归一化频率/π') ylabel('幅度/dB')

22

2010届电子信息工程专业毕业设计（论文）

图5-3 凯塞窗的频率特性曲线

这是一种适应性较强的窗。它拥有一个可自由选择的参数β，可以同时调整主瓣宽度与旁瓣电平。因此，改变β就能对主瓣宽度与旁瓣衰减进行选择。

23

2010届电子信息工程专业毕业设计（论文）

输入滤波器性能要求 带通 微分器 希尔伯特变换器 WATE加权函数 集F的密集格点 EFF要求的幅度 用公式表示逼近问题 瑞米兹交换算法 计算单位冲激响应h(n) 打印出最佳误差及h(n)

图5-4 最佳线性相位FIR滤波器设计算法框图

24

2010届电子信息工程专业毕业设计（论文）

(r+1)个极值频率的初始猜测值

N 最优逼近 Y 极值点位置是否变化？ N 极值点是否多于（r+1）个？ Y 舍去ω=0，π中较小的一个局部极值，保留（r+1）个最大极值 计算E（ω）并求出| E（ω）|>δ处的局部极大值 通过r点内插得到p(ω) 计算极值频率集上的最优

图5-5 瑞米兹算法的步骤

25

2010届电子信息工程专业毕业设计（论文）

结论

数字信号处理系统具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等突出优点，这些都是模拟信号处理系统所无法比拟的。

与传统的通过调整窗函数来进行优化设计的思想相比，使最后所设计出的滤波器的频率响应同理想频率响应之间的逼近误差最小。由于其中的逼近误差函数可灵活定义，算法能够适用于具有不同频域技术指标要求的滤波器的设计。

比较上述设计的结果可以得出以下结论：

(1).窗函数法是从时域出发，通过一定的窗函数截取有限长的单位脉冲响应来逼近理想单位脉冲响应；频率抽样法则是从频域出发对理想脉冲响应的逼近。

(2).采用切比雪夫最佳一致逼近法时,滤波器既有好的衰减特性又有好的边缘频率．滤波器的阶次也可以预先估计。

(3).通过MATLAB软件提供的Remez 函数实现Remez算法，设计滤波器逼近理想频率响应，从而达到FIR数字滤波器的最优化设计。应用Matlab对DSP进行系统级的设计极大地改进了传统的设计方法。Matlab系统级的设计环境,有助于在设计早期发现错误和应对系统复杂性不断增加的挑战,方便了复杂DSP应用系统的设计。所设计的数字滤波器和传统设计的滤波器比较,速度提高了很多,设计过程得到了很大的简化,方便了开发设计。

谢辞

在这里，对所有参与此次设计并给予帮助的老师和同学表示感谢。

首先，要感谢学院给予我们一次这样好的机会。让我们可以通过设计动手动脑，既巩固了过去四年中学习过的相关专业知识，又学习到了以前课本上未曾接触到的知识或技术。同时，也要感谢领导们的悉心照顾，为我们参加设计的同学提供良好的实验环境和上机环境，给我们带来非常大的便利。

其次，毕业设计的顺利完成，离不开指导老师、同学及朋友的关心和帮助。在整个的毕业设计中，老师、同学及朋友提供有利于设计的宝贵建议和意见，在他们的帮助下，毕业设计得以不断的完善，最终帮助我完整的写完了整个论文。

在这里，尤其要对这次设计的指导老师徐老师，表示诚挚的感谢。徐老师指引我的设计的方向和架构，并对本论文初稿进行逐字批阅，指正出其中误谬之处，使我有了思考的方向。他的循循善诱的教导和不拘一格的思路给予我无尽的启迪；他的严谨细致、一丝不苟的作风，将一直是我工作、学习中的榜样。在此，谨向徐老师表示崇高的敬意和衷心的感谢！

参考文献

[1]胡广书. 数字信号处理—理论、算法与实现[M]. 北京:清华大学出版社. 1997 [2]程佩清. 数字信号处理(第二版) [M]. 北京:清华大学出版社. 2001 [3]张志涌. 精通MATLAB6.5[M]. 北京航空航天出版社. 2003

26

2010届电子信息工程专业毕业设计（论文）

[4]陆心如. 数字信号处理原理及应用[M]. 北京:电子工业出版社. 1996 [5]黄大卫. 数字滤波器[M]. 北京:中国铁道出版社. 1991

[6]张雄伟. DSP 芯片的原理与开发应用[M] . 北京:电子工业出版社. 1997 [7]刘树堂,黄建国. 西安:西安交通大学出版社. 2001

[8]丁美玉,高西全. 数字信号处理[M].西安:西安电子科技大学出版社. 2001 [9]常迥. 数字信号处理. 北京：地质出版社，1980

[10]楼顺天,李博菡. 基于MATLAB 的系统分析与设计--信号处理[M]. 西安:西安电子科技大学

出版社. 2000

附录

在本次设计的最后，对研究中的比较重要和繁难的问题进行一下归纳和总结。主要解释的是在FIR数字滤波器设计的计算机实现过程中，所应用到的MATLAB软件的知识点和相关的低通滤波器的MATLAB程序。具体内容如下：

因为，低通FIR数字滤波器的理想性能指标要求为：通带截止频率Wp=0.3π,阻带截止频率Ws=0.5π，阻带衰减At不小于40dB，通带衰减不大于3dB, fs=4000。

所以，我们需要通过几种窗函数计算机实现的图象来进行对比。基于对窗函数的一些特性和公式的了解，我们选择凯塞窗、汉宁窗、矩形窗。在经过以上三种窗函数的对比分析之后，我们利用瑞米兹算法对滤波器进行了优化设计。 MATLAB程序： 凯塞窗：

wn=kaiser(28,3.44); nn=[0:1:27]; alfa=(28-1)/2 ;

hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa)); h=hd.*wn';

[h1,w1]=freqz(h,1);

plot(w1/pi,20*log10(abs(h1))); axis([0,1,-80,101]); grid;

xlabel('归一化频率/π')

ylabel('幅度/dB') 汉宁窗：

hanning: wp=0.3*pi; ws=0.5*pi; wdelta=ws-wp; N=ceil(8*pi/wdelta);

27

2010届电子信息工程专业毕业设计（论文）

Wn=(0.3+0.5)*pi/2;

b=firl(N,Wn/pi,hanning(N+1)); freqz(b,1,512)

矩形窗：

wp=0.3*pi; ws=0.5*pi; wdelta=ws-wp; N=ceil(8*pi/wdelta); Wn=(0.3+0.5)*pi/2;

b=firl(N,Wn/pi,hanning(N)); freqz(b,1,512)

瑞米兹(瑞米兹沃特)：

fedge=[800 1000];

mval=[1 0];dev=[0.0559 0.01]; fs=4000;

[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); [h,w]=freqz(b,1,256);

plot(w*2000/pi,20*log10(abs(h))); grid;

xlabel('频率/Hz') ylabel('幅度/dB')

28

2010届电子信息工程专业毕业设计（论文）

The Design of FIR Digital Filters Based on Window function Method

Hao Hai

Electrical & Information Engineering, Chongqing Three Gorges University, Wanzhou 404100

Abstract Description of the digital filter unit of the finite length impulse response (FIR) filter

principle, the FIR filter window function design method is studied. Window function method in the design of FIR digital filter has a wide range of applications. Such as corrugated paper, a Chebyshev approximation criteria for the design method of FIR digital filter through the MATLAB simulation to achieve. Conventional digital filter design method is cumbersome and results are not intuitive, this paper, MATLAB has a powerful scientific computing and graphics display the advantages, with the window function method combines design theory to design FIR digital filter design results not only more intuitive, but also improve the accuracy of the filter design to better achieve the desired results.

Key words FIR digital filter；window function；equivalent ripple Chebyshev approximation；

MATLAB simulation

29

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