毕业设计——基于matlab的FIR数字滤波器设计之欧阳德创编

毕业设计任务书

时间：2021.03.07 创作：欧阳德 设计题目：基于MATLAB的IIR数字滤波器设计

专业：通信工程 班级学号： 姓名： 指导教师：

设计期限：2012年3月 5日开始

2012年5月20

日结束

院、系：信息工程学院

2012

年3月7日 一、毕业设计的目的

1、通过毕业设计把自己在大学中所学的知识应用到实践当中。

2、深入了解利用Matlab设计FIR数字滤波器的基本方法。

3、在毕业设计的过程中基本掌握了Matlab编译程序的基本方法。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

4、提高自己的自学能力和动手能力。

5、锻炼自己通过网络及各种资料解决实际问题的能力。

二、主要设计内容 利用窗函数法、频率抽样法设计FIR滤波器，绘制出滤波器的特性图。利用所设计的滤波器对多个频带叠加的正弦信号进行处理，对比滤波前后的信号时域和频域图，验证滤波器的效果。最后找一段语音信号，并对此信号进行采样和加噪，绘制出采样后语音信号的时域波形和频谱图，然后用所设计的滤波器对加噪后的信号进行滤波，绘制出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化。 三、重点研究问题

基于Matlab的FIR数字滤波器的基本设计方法，能够根据性能指标要求灵活地进行滤波器的设计。

四、主要技术指标或主要设计参数

（1）滤波器类型

（2）滤波器阶数和采样频率 （3）通带和阻带截止频率 （4）通带和阻带衰减

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

五、设计成果要求

1、完成毕业设计书文档 2、完成程序的编译和调试 3、对程序主要语句做出注释

本科生毕业设计（论文）开题报告

2012年3月20日

学生姓名 题目名称 课题来源 学号 专业 通信工程 基于Matlab的FIR数字滤波器设计 导师提供 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

数字滤波技术是数字信号处理的一个重要组成部分,滤波器的设计是信号处理的核心问题之一。数字滤波器是一个离散的系统,它可以对输入的离散信号进行一系列运算处理，从输入的信号中获得所需要的信息。数字滤波器的系统函数通常表示为 jbzjMH(z)j01aizii1N 数字滤波器分为有限冲激响应数字滤波器，即FIR数字滤波器和无限冲激响应，即IIR数字滤波器。从公式的角度来看，FIR数字滤波器的始终为零；IIR数字滤波器aiai 至少有一个非零。 实现数字滤波器的方法一般有两种：一种方法是吧滤波器所要完成的运算编成程序并让计算机执行，也就是采用计算机软件来实现；另一种方法是设计专用的数字硬件、专用的数字信号处理器或采用通用的数字信号主要内容 处理器来实现。 本设计根据 FIR滤波器的设计原理,提出了Matlab环境下FIR滤波器的窗函数法、频率抽样法, Matlab环境为设计FIR滤波器提供了一个可靠而有效的工作平台。 Matlab软件以矩阵运算为基础，把计算、可视化及程序设计有机融合到交互式工作环境中，并且为数字滤波的研究和应用提供了一个直观、高效、便捷的利器。工程人员可以直观方便地进行科学研究与工程应用。Matlab是美国MathWorks 公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括 Matlab和Simlink两大部分。 在数字信号处理中 ,由于信号中经常混有各种复杂成分，所以很多信号分析都是基于滤波器而进行的， 因此数字滤波器占有极其重要的地位 。数字滤波器是具有一定传输选择特性的数字信号处理装置,其输入与输出均为数字信号,实质上是一个由有限精度算法实现的线性时不变离散欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

系统。它的基本工作原理是利用离散系统特性对系统输入信号进行加工和变换,改变输入序列的频谱或信号波形,让有用频率的信号分量通过,抑制无用的信号分量输出。数字滤波器和模拟滤波器有着相同的滤波概念,根据其频率响应特性可分为低通、高通、带通、带阻等类型。与模拟滤波器相比,数字滤波器除了具有数字信号处理固有优点外,还有滤波精度高、稳定性好、灵活性强等优点。FIR滤波器可以得到严格的线性相位,但它的传递函数的极点固定在原点,只能通过改变零点位置来改变性能,为了达到高的选择性,必须用较高的阶数,对于同样的滤波器设计指标,FIR滤波器要求的阶数可能比IIR滤波器高5～10倍。 在设计中，我将利用窗函数法、频率抽样法设计FIR滤波器，绘制出滤波器的特性图。利用所设计的滤波器对多个频带叠加的正弦信号进行处理，对比滤波前后的信号时域和频域图，验证滤波器的效果。最后找一段语音信号，并对找到的信号进行采样和加噪，绘制出采样后语音信号的时域波形和频谱图，然后用所设计的滤波器对加噪后的信号进行滤波，绘制出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化。 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

由于FIR数字滤波器具有严格的相位特性，对于信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率抽样法和优化设计方法。我们本文主要介绍窗函数法和频率抽样设计法。 其中利用窗函数法设计FIR滤波器的基本思路是：先给定频率响应函数，利用IDFT求出理想滤波器的单位响应在时域表达设计 h(n)逼近理想 hd(n)hd(n),从时域出发。我们要设计的是FIR滤波器，其h（n）必然hd(n)是有限长的，所以要用有限长的h（n）来逼近无限长的方法是将采 取 的 主 要 技 术 路 线 或 方 法 hd(n)。最有效的进行截断，或者说，是用一个有限长度的窗函数序列w，即利用h(n)hd(n)(n)（n）来截取hd(n)截取为有限长因果序列。按照线性相位滤波器的要求，线性相位FIR数字低通滤波器的单位抽样响应h(n)必须是偶对称的。矩形窗设计的FIR低通滤波器，最大相对肩峰值为8.95%，N增加钾时，2pi/N减小，故起伏振荡变密，最大肩峰则总是8.95%，这种现象称为吉布斯（Gibbs）现象。为了消除吉布斯效应，一般采用其他类型的窗函数，Matlab设计FIR滤波器有多种方法和对应的函数。窗函数设计法不仅在数字滤波器的设计中占有重要的地位，同时可以用于功率谱的估计，从根本上讲，使用窗函数的目的就是消除由无限序列的截短而引起的Gibbs现象所带来的影响。 利用频率抽样法的基本思路是：设所需滤波器的频率响应为Hd(ej)。现要求设计一个M阶的FIR滤波器h[k]，使得Hd(ej)Hd(ej)在M+1相个抽样点上，FIR滤波器的频率响应等，即 与所需的频率响应H(ej)Hd(ej)H(ejm)h[k]ejkm,m0,1,...,Mk0M Hd(ej)由设计要求给定，h[k]需要通过设计来确定。如果M+1个方程是线性无关的，则可以通过求解M+1阶的线性方程得出FIR滤波器的欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

h[k]。对m的一些特殊抽样法，上式方程的解可以直接由IDFT得到。由Hd(ej)于要求设计出的滤波器是实系数的线性相位FIR滤波器，所以抽样值还需要满足线性相位滤波器的约束条件。 的欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

1、毕业设计书文档 2、程序的编译和调试 3、程序主要语句注释 预 期 的 成 果 及 形 式 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

第3周 根据任务书，查阅相关资料 第4周 完成开题报告 第5周 复习数字信号处理中关于FIR数字滤波器的内容 第6周 复习matlab课程中的笔记，熟悉matlab的使用 第7周 开始利用matlab进行程序的编译 第8周 继续利用matlab进行程序的编译，并修改程序中所出现的错误 时间安排 第9周 完成外文翻译 第10周 继续利用matlab进行程序的编译和修改 第11周 继续利用matlab进行程序的编译和修改，并完善程序注释 第12周 完成毕业论文 第13周 制作PPT，准备答辩 第14周 答辩 签 名： 年 月 日 指导教师意见 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

备注 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

基于Matlab的FIR数字滤波器设计

摘要：在数字信号处理中 ,由于信号中经常混有各种复杂成分，所

以很多信号分析都是基于滤波器而进行的， 因此数字滤波器占有极其重要的地位。在数字控制系统中输入信号中所含的干扰对系统的性能会产生很大的影响，因此需要对输入信号进行处理，以提取有用信号。有限长冲激响应（FIR）滤波器在数字信号处理中发挥着重要作用，采用Matlab软件对FIR数字滤波器进行仿真设计，简化了设计中繁琐的计算。本文采用窗函数法，频率采样法通过调用Matlab函数设计FIR数字滤波器。绘制对应的幅频特性曲线。最后用基于Matlab函数设计的FIR数字滤波器进行语音滤波处理，通过滤波前后信号的频谱图和生成的声音文件的对比，分析不同滤波器的滤波效果。

关键词：FIR数字滤波器，仿真，窗函数法，频率抽样法，Matlab

Design of the MATLAB-basedFIR digital

filter

Abstract:In digital signal processing, because the signal is often mixed

with a variety of complex composition, so a lot of signal analysis are based on the filter, digital filter occupies an extremely important position.In digital control system, interference, which is mixed in the input signal, has a great effect on performance of the system. Therefore, processing of input signal has to be done to get useful signal. Finite impulse response (FIR) filter plays an important role in the processing of digital signal. Designing the FIR filter by Matlab can simplify the complicated computation in simulation and improve the performance. By using the methods of window function, frequency sampling ,the design of FIR digital filter has been processed in Matlab. In the view of the

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

designed program of Matlab and the figure of the amplitude-frequency characterization. At last, by using the FIR digital filters which have been designed to process the sound signal based on the Matlab function, the filtering effect of different digital filters is analyzed by comparing the signal’s spectrum viewers and the sound files which have been generated. The experimental results show that the FIR filters designed in this paper are effective.

Key words: FIR digital filter, simulation, windowing method, frequency sampling

method,

Matlab

欧阳德创编 2021.03.07

目录

欧阳德创编 2021.03.07

第一章 绪论

Matlab是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括Matlab和Simulink两大部分。

1.1 Matlab简介

1.1.1 MATLAB的发展

MATLAB是英文MATrix LABoratory（矩阵实验室）的缩写。早期的MATLAB是用FORTRAN语言编写的，尽管功能十分简单，但作为免费软件，还是吸引了大批使用者。经过几年的校际流传，在John Little。Cleve Moler和Steve Banger合作，于1984年成立MathWorks公司，并正式推出MATLAB第一版版。从这时起，MATLAB的核心采用C语言编写，功能越来越强大，除原有的数值计算功能外，还新增了图形处理功能。

MathWorks公司于1992年推出了具有划时代意义的4.0版；1994年推出了4.2版扩充了4.0版的功能，尤其在图形界面设计方面提供了新方法；1997年春5.0版问世，5.0版支持了更多的数据结构，使其成为一种更方便、更完善的编程语言；1999年初推出的MATLAB5.3版在很多方面又进一步改进了MATLAB语言的功能，随之推出的全新版本的最优化工具箱和Simulink3.0达到了很高水平；2000年10月，MATLAB6.0版问世，在操作页面上有了很大改观，为用户的使用提供了很大方便，在计算机性能方面，速度变的更快，性能也更好，在图形界面设计上更趋合理，与C语言接口及转换的兼容性更强，与之配套的Simulink4.0版的新功能也特别引人注

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

目；2001年6月推出的MATLAB6.1版及Simulink4.1版，功能已经十分强大；2002年6月推出的MATLAB6.5版及Simulink5.0版，在计算方法、图形功能、用户界面设计、编程手段和工具等方面都有了重大改进；2004年，MathWorks公司推出了最新的MATLAB7.0版，其中集成了最新的MATLAB7编译器、Simumlink6.0仿真软件以及很多工具箱。这一版本增加了很多新的功能和特性，内容相当丰富。

Matlab主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，代表了当今国际科学计算软件的先进水平。

1.1.2 Matlab的主要功能

1. 数值计算和符号计算功能 2. 绘图功能 3. 语言体系

4. MATLAB的工具箱

1.2 Matlab的优势及特点

1.2.1 MATLAB的优势

(1) 工作平台编程环境十分友好 （2）编程语言简单易用

（3）数据的计算处理能力十分强大 （4）图像处理能力强大

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

（5）模块集合工具箱应用广泛 （6）程序的接口和发布平台很实用 （7）可以开发用户界面

1.2.2 Matlab 语言的特点

MATLAB语言被称为计算机语言，其利用丰富的函数资源，使程序员从繁琐的程序代码中出来，其最突出的特点就是简洁。MATLAB用更直观的、符合人们思维习惯的代码，代替了C和FORTRAN语言的冗长代码，给用户带来最直观、最简洁的程序开发环境，下面简单介绍一下MATLAB的主要特点。

语言简洁紧凑，使用方便，库函数十分丰富。MATLAB程序书写的形式自由，利用丰富的库函数避开了繁琐的子程序编程任务，由于库函数都是由本领域的专家编写，所以不必担心函数的可靠性。

高效方便的矩阵和数组运算，MATLAB语言不需要定义数组的维数，并给出了矩阵函数、特殊矩阵函数、特殊矩阵专门的库函数，使得在求解信号处理、建模、系统识别、优化和控制等领域的问题时，显得大为简洁、方便、高效，这是其他高级语言所不能的。

MATLAB既具有结构化的控制语句，又具有面向对象编程的特性。

MATLAB语法不严格，程序设计自由度大，通过建立M后缀名文件的形式，与用户已经编好的FORTRAN、C语言成语混合编程，方便地调用有关的FORTRAN、C语言的子程序。可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上面运行。

MATLAB的图形功能强大。在C和FORTRAN语言里，绘图都很不容易，但在MATLAB里，数据的可视化非常简单。此外，

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

MATLAB还具有较强的编辑图形界面的能力。

MATLAB拥有功能强大的工具箱，主要用来扩充其符号计算功能、图示建模仿真功能、文字处理功能以及与硬件实施交互功能。

源程序的开放性强。除内部函数以外，所有MATLAB的核心文件和工具箱文件都是可读可改变的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。

MATLAB软件自1984年推向市场以来，历经十几年的发展和竞争，现已成为国际公认的最优秀的科技应用软件。它功能强大、界面友好、语言自然、开放性强，很快成为应用学科计算机辅助分析、设计、仿真、教学乃至科技文字吹不可缺少的基础软件。

第二章 数字滤波器

2.1 数字滤波器简介

数字滤波器是一个离散的系统。它可以对输入的离散信号进行一系列运算处理，从输入的信号中获得所需要的信息。数字滤波器的系统函数通常表示为

bzjMjH(z)j01aizii1N (1-1)

数字滤波器分为有限冲激响应数字滤波器，即FIR数字滤波器和无限冲激响应，即IIR数字滤波器。从公式的角度来看，FIR数字滤波器的

aia 始终为零；IIR数字滤波器i至少有一个非零。

实现数字滤波器的方法一般有两种：一种方法是吧滤波器所要完成的运算编成程序并让计算机执行，也就是采用计算机软件来实

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

现；另一种方法是设计专用的数字硬件、专用的数字信号处理器或采用通用的数字信号处理器来实现。实现一个数字滤波器一般需要三个基本的运算单元：加法器、单位延时器和常数乘法器。这些基本的单元可以有两种表示方法—方框图法和信号流程图法。

设计一个数字滤波器的一般步骤为： （1）按所给要求确定滤波器的性能

（2）用一个因果稳定的离散线性时不变的系统函数逼近此性能的要求

（3）利用算法来实现这个系统函数 （4）利用计算机仿真或硬件来实现

2.2 IIR数字滤波器

无限长单位冲激响应滤波器，即IIR数字滤波器具有下面几个特点：

（1） 系统的单位冲激响应h(n)为无限长的; （2） 系统函数H(z)在有限z平面上有极点存在；

（3） 结构上存在着输出到输入的反馈，也就是结构上是递归

型的。

IIR滤波器的设计就是在给定的技术指标下去确定滤波器的阶数N和系数 {

ai，

bi}。在已满足给定的技术指标下，应选用阶数尽可能低的滤

波器，因为滤波器的阶数越低，在实现时成本就越低。

在设计IIR滤波器时，最常用的方法是利用模拟滤波器来设计数字滤波器。其原因为：

（1） 模拟滤波器的设计技术相对成熟，可以广泛利用； （2） 模拟滤波器有大量的参考程序和表格； （3） 它的解可以为闭合形式的。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

2.3 FIR数字滤波器

有限长单位冲激响应滤波器，即FIR数字滤波器具有下面几个特点：

（1） 系统的单位冲激响应h(n)在有限个n值处不为零; （2） 系统函数H(z)在|z|>0处收敛，在|z|<0处只有零，即

有限z平面上只有零点，儿全部极点都在z=0处（因果系统）；

（3） 没有输出到输入的反馈，但有些结构中（例如频率抽样

结构）也包含有反馈的递归部分。，

FIR滤波器是指在有限范围内系统的单位脉冲响应h[k]仅有非零值的滤波器。M阶FIR滤波器的系统函数H(z)为

H(z)h[k]zkk0kM (1-2)

其中H(z)是z的M阶多项式，在有限的z平面内H(z)有M个零点，在z平面原点z=0有M个极点.

j FIR滤波器的频率响应H(e)为

H(e)h[k]ejkjk0M (1-3)

它的另外一种表示方法为

H(ej)H(ej)ej()其中

(1-4)

H(ej)和()分别为系统的幅度响应和相位响应。

若系统的相位响应()满足下面的条件

() (1-5)

即系统的群延迟是一个与没有关系的常数，称为系统H(z)具有

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

严格线性相位。由于严格线性相位条件在数学层面上处理起来较为困难，因此在FIR滤波器设计中一般使用广义线性相位。

jH(e)可以表示为 如果一个离散系统的频率响应

H(ej)A()ej()(1-6)

其中和是与无关联的常数，A()是可正可负的实函数，则称系统是广义线性相位的。

如果M阶FIR滤波器的单位脉冲响应h[k]是实数，则可以证明系统是线性相位的充要条件为

h[k]h[Mk](1-7)

当h[k]满足h[k]=h[M-k],称h[k]偶对称。当h[k]满足h[k]=-h[M-k],称h[k]奇对称。按阶数h[k]又可分为M奇数和M偶数，所以线性相位的FIR滤波器可以有四种类型。

四种线性相位FIR滤波器的性质如表1-1所示

表1-1 四种线性相位FIR滤波器的特性

类型 阶数M h[k]的对称性 A()关于I 偶数 偶对称 偶对称 偶对称 II 奇数 偶对称 偶对称 奇对称 III 偶数 奇对称 奇对称 奇对称 IV 奇数 奇对称 奇对称 偶对称 0的对性 A()关于的对性 A()的周期 2 0 任意 任意 4 0 任意 0 2 0.5 0 0 4 0.5 0 任意  A(0) A() 欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

可适用的滤波器类型 LP,HP,BP,SP 微分器，变换LP,BP 器，Hilbert 微分器，变换器，Hilbert，HP 2.4 IIR与FIR数字滤波器的比较

（1） IIR滤波器h（n）无限长，FIR滤波器h（n）有限长。 (2)在技术指标相同的条件下，IIR滤波器的输出对输入有反馈，所以可以用比FIR少的阶数来满足要求，存储单元少，运算次数也少，经济实惠。

（3）FIR滤波器的相位是严格线性的，而IIR滤波器做不到这一点，IIR滤波器的选择性越好，其相位的非线性越严重。

（4）FIR滤波器主要采用非递归结构， 有限精度的运算误差很小。而IIR滤波器在运算中会产生寄生振荡。

（5）FIR滤波器可以使用快速傅里叶变换算法，而IIR滤波器不能这样。

（6）IIR滤波器可以利用模拟滤波器的公式、数据和表格，计算量小。FIR滤波器设计时往往要借助计算机。

（7）IIR滤波器极点位于z平面任意位置，而FIR滤波器极点固定在原点。

（8）IIR滤波器用于设计规范化的选频滤波器，FIR滤波器可设计各种幅度特性和相频特性的滤波器。

第三章 FIR数字滤波器的设计

在数字信号处理中 ,由于信号中经常混有各种复杂成分，所以很多信号分析都是基于滤波器而进行的， 因此数字滤波器占有极其重要的地位 。数字滤波器是具有一定传输选择特性的数字信号处理

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

装置,其输入与输出均为数字信号,实质上是一个由有限精度算法实现的线性时不变离散系统。它的基本工作原理是利用离散系统特性对系统输入信号进行加工和变换,改变输入序列的频谱或信号波形,让有用频率的信号分量通过,抑制无用的信号分量输出。数字滤波器和模拟滤波器有着相同的滤波概念,根据其频率响应特性可分为低通、高通、带通、带阻等类型。与模拟滤波器相比,数字滤波器除了具有数字信号处理固有优点外,还有滤波精度高、稳定性好、灵活性强等优点。FIR滤波器可以得到严格的线性相位,但它的传递函数的极点固定在原点,只能通过改变零点位置来改变性能,为了达到高的选择性,必须用较高的阶数,对于同样的滤波器设计指标,FIR滤波器要求的阶数可能比IIR滤波器高5～10倍。

由于FIR数字滤波器具有严格的相位特性，对于信号处理和数据传输是很重要的。目前FIR滤波器的设计方法主要有三种：窗函数法、频率抽样法和优化设计方法。

我们本章主要介绍窗函数法、频率抽样法。

3.1 窗函数法设计FIR滤波器

窗函数设计法又称为傅里叶级数法。这种方法首先给出

Hd(ej)，

Hd(ej)表示要逼近的理想滤波器的频率响应，则由

IDTFT可得出滤波器的单位脉冲响应为

1hd[k]2Hd(ej)ejkd (3-1)

hd[k]由于是理想滤波器，故是无限长序列。但是我们所要设计

的FIR滤波器，其h[k]是有限长的。为了能用FIR滤波器近似理想滤波器，需将理想滤波器的无线长单位脉冲响应行截断。当截断后的单位脉冲响应

hd[k]hd[k]分别从左右进

不是因果系统的时候，可将

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

其右移从而获得因果的FIR滤波器。 另一种设计方案是将线性相位因子e的频率响应中，然后利用IDTFT计算出

j(0.5M)加入到理想滤波器

hd[k]hd[k]后，取在0≦k≦M

范围的值为FIR滤波器单位脉冲响应。 理想滤波器的频率响应

Hd(ej)Hd(ej)和设计出的滤波器的频率响应

的积分平方误差定义为

1222也可以表示为

Hd(ej)H(ej)d(3-2)

22khd[k]h[k]2M2(3-3)

2k1hd[k]hd[k]h[k]k0kM1hd[k]2

上式中的第一项和第三项与所设计出的滤波器参数是没有关系的，为了使上式中的第二项达到最小，可选择

h[k]hd[k],0kM佳滤波器。

(3-4)

所以用上面的方法得出的滤波器是在积分平方误差最小意义下的最

Gibbs现象就是理想滤波器的单位脉冲响应

hd[k]截断获得的

FIR滤波器的幅度函数A()在通带和阻带都呈现出振荡现象。随着滤波器阶数的增加，幅度函数在通带和阻带振荡的波纹数量也随之增加，波纹的宽度随之减小，然而通带和阻带最大波纹的幅度与滤波器的阶数M无关。窗函数的主瓣宽度决定了度，窗函数长度N增大，过渡带减小。

下面介绍一些常用的窗函数，用N=M+1表示窗函数的长度。

Hd(ej)过渡带的宽

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

（1） 矩形窗

10kM[k]0otherwise (3-5)

矩形窗的主瓣宽度为宽度近似为

1.8N。

4N。用矩形窗设计的FIR滤波器过渡带

（2） Hanning窗

0.50.5cos(2kM)(0kM)[k]0(otherwise)(3-6)

Hanning窗的主瓣宽度为

8N。由Hanning窗的定义可知，

Hanning窗在其两个端点的值为零，这就使得在实际的应用中不能利用两个端点的数据。我们可将N+2点的Hanning窗除去两个端点来定义长度为N的Hanning窗。修改后的长度为N的Hanning窗定义为

2(k1))(0kM)0.50.5cos(M[k]0(otherwise)(3-7)

在Matlab信号处理工具箱中所采用的就是这种修改后的定义方式。

（3） Hamming窗

对升余弦加以改进，可以得到旁瓣更小的效果，窗形式为

0.0.46cos(2k/M)(0kM)[k]0(otherwise)(3-8)

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

与汉宁窗相比，主瓣宽度相同为瓣峰值小于主瓣峰值的1%。 （4） Blackman窗

8N，但旁瓣幅度更小，旁

为了更进一步抑制旁瓣，可再加上余弦的二次谐波分量，得到Blackman窗

0.420.5cos(2kM)0.08cos(4kM)(0kM)w[k]0(otherwise)(3-9)

Blackman窗的主瓣宽度为（5） Kaiser窗

此种窗是一种应用广泛的可调节窗，它可以通过改变窗函数的形状来控制窗函数旁瓣的大小，从而在设计中可用滤波器的衰减指标来确定窗函数的形状。长度为N的Kaiser窗定义为

12N。

[k]I0(1(12kI0()M)2),0kM(3-10)

其中是一个可调参数，可以通过改变的值来调整窗函数的形状，从而达到不同的阻带衰减要求。上式中的正贝塞尔函数。可用幂级数表示为

I0(x)是零阶第一类修

(x/2)n2I0(x)1[]n!n1(3-11)

对于任意的一个实变量x，函数的值都是正的。在实际计算中，上式的求和一般取20项就能达到所需精度。随着参数的增加，Kaiser窗在两端的衰减是逐渐加大的。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

3.2 频率抽样法设计线性相位FIR滤波器

频率抽样法是从频域出发，在频域直接设计，把给定的理想频率响应

Hd(ej)加以等间隔抽样，并以此作为实际FIR滤波器的频率

Hd(ej)响应。设所需滤波器的频率响应为。

Hd(ej)现要求设计一个M阶的FIR滤波器h[k]，使得个抽样点上，FIR滤波器的频率响应

H(ej)Hd(ej)在M+1

与所需的频率响应

相等，即

jmHd(e)H(eHd(ej)j)h[k]ejkm,m0,1,...,Mk0M(3-12)

由设计的要求给定，h[k]需要通过设计来确定。如果M+1个

方程是线性无关的，则可以通过求解M+1阶的线性方程来得出FIR滤波器的h[k]。对

m的一些特殊抽样方法，上述方程的解可以直

接由IDFT得到。由于要求设计出的滤波器是实系数的线性相位FIR滤波器，所以件。

I型和II型线性相位滤波器的0，III型和IV型线性相位滤波

Hd(ej)的抽样值还需要满足线性相位滤波器的约束条

2Hd(ej)器的。为了使设计出的滤波器具有线性相位，在M+1个抽样点上的值

Hd[m]j应为

Hd[m]Hd(e13)

2mM1)ejejMmM1Ad(2m(M1)) (3-

ejejMmM1Ad[m]

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

下面分别讨论四种线性相位滤波器在抽样点

Hd[m]上的值：

I型（M为偶数，h[k]偶对称）线性相位FIR滤波器在M+1个抽样点值为

MjmM1(0mM2)Hd[m]Ad[m]e*Hd[M1m](M21mM)(3-14)

上式表明I型线性相位FIR滤波器值可由

Hd[m]Hd[m]Hd[m]在M21mM的

H[m]在1mM2的值确定。在d的值确定后，对

做M+1点的IDFT即可得到I型线性相位滤波器的h[k]。

II型（M为奇数，h[k]偶对称）线性相位FIR滤波器在M+1个抽样点值为

MjmM1(0m(M1)2)Ad[m]eHd[m]0(m(M1)2)H*[M1m]((M3)21mM)d(3-

15)

上式表明II型线性相位FIR滤波器

Hd[m]在

(M3)21mM的值可由Hd[m]在1m(M1)2的值确定。

III型（M为偶数，h[k]奇对称）线性相位FIR滤波器在M+1个抽样点值为

0(m0)MjmM1Hd[m]jAd[m]e(1mM2)*Hd[M1m](M21mM)(3-16)

上式表明III型滤波器线性相位FIR滤波器

Hd[m]在M21mM欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

的值可由

Hd[m]在1mM2的值确定。

IV型（M为奇数，h[k]奇对称）线性相位FIR滤波器在M+1个抽样点值为

0(m0)MjmM1Hd[m]jAd[m]e(1m(M1)2)*Hd[M1m]((M3)2mM)(3-17)

上式表明IV型线性相位FIR滤波器可由

Hd[m]Hd[m]在(M3)2mM的值

在1m(M1)2的值确定。

进行频率抽样，就是在z平面单位圆上的N个等间隔

对

Hd(ej)点上抽样出频率响应值。在单位圆上可以有两种抽样方式，第一种是第一个抽样点在w=0处，第二种是第一个抽样点在w=pi/M处，每种方式可分为M为偶数与M为奇数两种。

为了提高逼近质量，使逼近误差更小，也就是减小在通带边缘由于抽样点的徒然变化而引起的起伏变化（这种起伏振荡使阻带内最小衰减变小，例如从衰减30dB变小为衰减20dB）。和窗口法的平滑截断一样，这里是使理想频率响应的不连续点的边缘加上一些过渡的抽样点（在这些点上抽样的最佳值由计算机算出），从而增加过渡带，减小频带边缘的突变，也就是减小了起伏振荡，增大了阻带最小衰减。这些抽样点上的取值不同，效果也就不同。如果精心设计过渡带的抽样值，就有可能使它的游泳频带的博文减小，从而设计出较好的滤波器。一般过渡带取一、二、三点抽样值即可得到满意结果。

在理想低通滤波器的设计中，若不增加过渡点，阻带和通带之间的衰减约为-21dB，如果在通带和阻带之间增加一个采样点，阻带的最小衰减可以提高到-65dB，如果增加两个采样点，阻带的最小衰减可以提高到-75dB，如果增加3个采样点，阻带的最小衰减可以提

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

高到-85dB至-95dB。

频率抽样法的优点是可以在频域直接设计，并且适合于最优化设计；缺点是抽样频率只能等于2pi/M的整数倍或等于2pi/M的整数倍上加上pi/M，因而不能确保截止频率Wc的自由取值。要想实现自由选择频率，则必须增加抽样点数M，但这种计算量加大。

第四章 利用Matlab实现FIR滤波器设计

在利用Matlab设计FIR滤波器时，分别采用窗函数法、频率抽样法和优化设计方法去设计所需的滤波器。在设计的过程中，用设计的滤波器对加有噪声的语音信号或不同频率叠加的正弦输入信号进行滤波，对比输入前后的图像，以此验证滤波器的性能。在程序绘制的图像中，有滤波器的特性图、输入信号的时域频域图和输出信号的时域频域图。

4.1 窗函数法的Matlab实现

在窗函数法的Matlab实现中，程序中经常使用的函数有fir1和kaiserord。

程序中fir1函数的用法：b=fir1(n,Wn,’ftype’,window) ①n为滤波器的阶数

②Wn为滤波器的截止频率，它是一个0到1的数。如果Wn是一个含有两个数的向量，则函数返回一个带通滤波器

③ftype为滤波器的类型，ftype=’high’时，设计的是高通滤波器；ftype=’stop’时，设计的是带阻滤波器；没有此参数时，设计的是低通滤波器

④window为指定的窗函数，矩形窗为boxcar(n)，汉宁窗为hanning(n)，海明窗为hamming(n)，布莱克曼窗为blackman(n)，凯

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

撒窗为kaiser(n,beta)，没有此参数时，默认为hamming窗 程序中kaiserord函数的用法：

[n,Wn,beta,ftype]=kaiserord(f,a,dev,Fs)

①f是一个向量，为设计滤波器过渡带的起始点和结束点 ②a是一个向量，指定频率段的幅度值

②dev是一个向量，长度和a相同，为各个通带和阻带内容许的幅度最大误差

④n为能够满足要求的滤波器的最小阶数 ⑤Wn为滤波器的截止频率

⑥ftype为根据待设计滤波器的要求得到的滤波器的类型 高通滤波器是容许高频信号通过、但减弱（或减少）频率低于截止频率信号通过的滤波器。对于不同滤波器而言,每个频率的信号的减弱程度不同。它有时被称为低频剪切滤波器；在音频应用中也使用低音消除滤波器或者噪声滤波器。低通滤波器与高通滤波器特性恰恰相反。 (1)

利用窗函数法设计低通滤波器

设计要求：①使用hamming窗，采样频率2000Hz

②通带截频0.1，阻带截频0.17

③通带衰减小于等于0.1dB，阻带衰减大于等于50dB

程序参见附录二中的1-（1）利用窗函数法设计低通滤波器

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

滤波器的增益响应200-20-40增益/分贝-60-80-100-120-1400100200300400500600频率/赫兹7008009001000图4-1 窗函数法设计低通滤波器的增益响应

从参考程序及图4-1可以得到所设计出滤波器的参数如下： ①滤波器的采样频率为2000Hz，滤波器的阶数为266

②滤波器的通带截频0.1，阻带截频0.17，过渡带宽0.07 ③通带衰减为0.019dB，阻带衰减为53dB

对比设计要求与所设计出滤波器的参数可知，其各项参数均满足设计指标，所设计出的滤波器即为设计所要求的滤波器。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波前时域图21幅度0-1-200.050.10.150.250.3时间/秒信号滤波前频域图0.20.350.40.450.5300200幅度10000100200300400500600频率/赫兹7008009001000

图4-2 信号滤波前的时域图和频域图

信号滤波后时域图21幅度0-1-20.20.250.30.350.4时间/秒信号滤波后频域图0.450.58060幅度402000100200300400500600频率/赫兹7008009001000图4-3 信号滤波后的时域图和频域图

从图4-2和图4-3的图像中可以看到：输入信号是由两个不同频

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

率的正弦信号叠加而成，信号频域图中位于滤波器通带内的频率分量保留了下来，位于滤波器阻带内的频率分量被滤除，滤波器的效果符合设计要求。

（2） 利用窗函数法设计带通滤波器

设计要求：①使用Kaiser窗，采样频率8000Hz

②通带截频0.325与0.5525，阻带截频0.25与

0.6025

③阻带衰减大于等于40dB，通带和阻带波纹0.01

程序参见附录二中的1-（2）利用窗函数法设计带通滤波器

滤波器的增益响应200-20增益/分贝-40-60-80-100-12005001000150020002500频率/赫兹300035004000

图4-4 窗函数法设计带通滤波器的增益响应

从参考程序及图4-4可以得到所设计出滤波器的参数如下： ①滤波器的采样频率为8000Hz，滤波器的阶数为90 ②滤波器的通带截频0.325与0.5525，阻带截频0.25与0.6025，过渡带宽0.075与0.05

③阻带衰减为40dB，通带和阻带的波纹均为0.01

对比设计要求与所设计出滤波器的参数可知，其各项参数均满足设计指标，所设计出的滤波器即为设计所要求的滤波器。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波前时域图42幅度0-2-400.0050.010.0150.020.025时间/秒信号滤波前频域图150100幅度50005001000150020002500频率/赫兹300035004000

图4-5 信号滤波前的时域图和频域图

信号滤波后时域图42幅度0-2-40.0050.010.015时间/秒信号滤波后频域图0.020.0258060幅度4020005001000150020002500频率/赫兹300035004000

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

图4-6 信号滤波后的时域图和频域图

从图4-5和图4-6的图像中可以看到：输入信号是由四个不同频率的正弦信号叠加而成，信号频域图中位于滤波器通带内的频率分量保留了下来，位于滤波器阻带内的频率分量被滤除，滤波器的效果符合设计要求。

(3) 利用窗函数法设计多通带滤波器 设计要求：①使用Kaiser窗，采样频率200Hz

②通带截频0.2、0.4、0.7、0.8 阻带截频0.1、0.5、0.6、0.9 ③阻带衰减大于等于30dB，通带和阻带波纹0.01

程序参见附录二中的1-（3）利用窗函数法设计多通带滤波器

滤波器的增益响应100-10-20-30增益/分贝-40-50-60-70-80-900102030405060频率/赫兹708090100

图4-7 窗函数法设计多通带滤波器的增益响应

从参考程序及图4-7可以得到所设计出滤波器的参数如下： ①滤波器的采样频率为200Hz，滤波器的阶数为46

②滤波器的通带截频0.2、0.4、0.7、0.8，阻带截频0.1、0.5、0.6、0.9，过渡带宽均为0.1 ③阻带衰减为38dB，通带和阻带的波纹均为0.01

对比设计要求与所设计出滤波器的参数可知，其各项参数均满

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

足设计指标，所设计出的滤波器即为设计所要求的滤波器。

信号滤波前时域图42幅度0-2-400.10.20.30.50.6时间/秒信号滤波前频域图0.40.70.80.91150100幅度5000102030405060频率/赫兹708090100

图4-8 信号滤波前的时域图和频域图

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波后时域图42幅度0-2-400.10.20.30.50.6时间/秒信号滤波后频域图0.40.70.80.91100幅度5000102030405060频率/赫兹708090100

图4-9 信号滤波后的时域图和频域图

从图4-8和图4-9的图像中可以看到：输入信号是由六个不同频率的正弦信号叠加而成，信号频域图中位于滤波器通带内的频率分量保留了下来，位于滤波器阻带内的频率分量被滤除，滤波器的效果符合设计要求。

4.2 频率抽样法的Matlab实现

(1)

利用频率抽样法设计低通滤波器

设计要求：①通带截频0.5，阻带截频0.6

②阻带衰减大于等于15dB

程序参见附录二中的2-（1）利用频率抽样法设计低通滤波器

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

滤波器的增益响应0-5-10-15增益/分贝-20-25-30-35-40-45-5000.10.20.30.40.50.6归一化角频率0.70.80.91

图4-10 频率抽样法设计低通滤波器的增益响应

从参考程序及图4-7可以得到所设计出滤波器的参数如下： ①滤波器的阶数为63

②滤波器的通带截频0.5，阻带截频0.6，过渡带宽为0.1 ③阻带衰减为17dB

对比设计要求与所设计出滤波器的参数可知，其各项参数均满足设计指标，所设计出的滤波器即为设计所要求的滤波器。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波前时域图21幅度0-1-200.050.10.150.20.25时间/秒信号滤波前频域图300200幅度10000100200300400500600频率/赫兹7008009001000

图4-11 信号滤波前的时域图和频域图

信号滤波后时域图21幅度0-1-20.20.2050.210.2150.2250.230.235时间/秒信号滤波后频域图0.220.240.2450.25300200幅度10000100200300400500600频率/赫兹7008009001000

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

图4-12 信号滤波后的时域图和频域图

从图4-11和图4-12的图像中可以看到：输入信号是由三个不同频率的正弦信号叠加而成，信号频域图中位于滤波器通带内的频率分量保留了下来，位于滤波器阻带内的频率分量被滤除，滤波器的效果符合设计要求。

（2）利用频率抽样法设计高通滤波器

设计要求：①通带截频0.5，阻带截频0.6

②阻带衰减大于等于15dB

程序参见附录二中的2-（2）利用频率抽样法设计高通滤波器

滤波器的增益响应0-5-10-15-20增益/分贝-25-30-35-40-45-5000.10.20.30.40.50.6归一化频率0.70.80.91

图4-13 频率抽样法设计高通滤波器的增益响应

从参考程序及图4-7可以得到所设计出滤波器的参数如下： ①滤波器的阶数为32

②滤波器的通带截频0.6，阻带截频0.5，过渡带宽为0.1 ③阻带衰减为18dB

对比设计要求与所设计出滤波器的参数可知，其各项参数均满足设计指标，所设计出的滤波器即为设计所要求的滤波器。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波前时域图42幅度0-2-400.050.10.150.20.25时间/秒信号滤波前频域图300200幅度10000100200300400500600频率/赫兹7008009001000

图4-14 信号滤波前的时域图和频域图

信号滤波后时域图21幅度0-1-20.20.2050.210.2150.2250.230.235时间/秒信号滤波后频域图0.220.240.2450.25300200幅度10000100200300400500600频率/赫兹7008009001000

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

图4-15 信号滤波后的时域图和频域图

从图4-14和图4-15的图像中可以看到：输入信号是由三个不同频率的正弦信号叠加而成，信号频域图中位于滤波器通带内的频率分量保留了下来，位于滤波器阻带内的频率分量被滤除，滤波器的效果符合设计要求。

4.3 利用滤波器处理加有噪声的音频波形

（1） 利用窗函数法设计的低通滤波器处理加有噪声的音频波形 程序参见附录二3-（1）利用窗函数法设计的低通滤波器处理加噪声的音频波形

加噪前音频语音波形的时域图10.5幅度(Y)0-0.5-10.050.10.20.25时间(t)加噪前音频波形的频域图0.150.30.35200150幅度(FY)1005000100200300400500600频率(f)7008009001000

图4-16 原始音频的时域与频域

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

加噪声后音频波形的时域图10.5幅度(Y1)0-0.5-10.050.10.20.25时间(t)加噪声后音频波形的频域图0.150.30.35200幅度(FY1)150100500400500600700800频率(f)900100011001200

图4-17 加噪声后音频的时域与频域波形

滤波器的增益响应200-20-40增益/分贝-60-80-100-120-1400200400600800频率/赫兹100012001400

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

图4-18 滤波器的增益响应

滤波后音频波形的时域图10.5幅度(sf)0-0.5-10.050.10.20.25时间(t)滤波后音频波形的频域图0.150.30.35200150100500400500600700800频率(f)900100011001200幅度(Fsf)

图4-19 滤波后音频的时域与频域波形

从参考程序及以上的四个图像中可以得到如下结论：

①从原始信号波形的频域图可以看到其频率分量主要在500到900Hz之间，噪声的频率分量主要集中在950Hz，利用通带截频为800Hz的低通滤波器可以滤除噪声。对比图4-16和图4-19滤波前后的波形和频谱，可以看到波形得到了重现

②滤波器的采样频率为22050Hz，滤波器的阶数为266

③滤波器的通带截频0.8，阻带截频0.82，过渡带宽0.02 ④通带衰减为0.019dB，阻带衰减约为53dB

（2） 利用频率抽样法设计的高通滤波器处理加有噪声的音频波形 程序参见附录二3-（2）

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

加噪前信号波形的时域图10.5幅度(Y)0-0.5-100.050.10.20.25时间(t)加噪前信号波形的频域图0.150.30.350.4200150幅度(FY)1005000100200300400500600频率(f)7008009001000

图4-20 加噪前信号的时域与频域波形

加噪声后信号波形的时域图10.5幅度(Y1)0-0.5-100.050.10.20.250.3时间(t)加噪声后信号波形的频域图0.150.350.4200幅度(FY1)1501005000100200300400500600频率(f)7008009001000

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

图4-21 加噪后信号的时域与频域波形

滤波器的增益响应200-20增益/分贝-40-60-80-10000.10.20.30.40.50.6归一化频率0.70.80.91

图4-22 滤波器的增益响应

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

信号滤波后时域图10.5幅度0-0.5-100.050.10.150.20.25时间/秒信号滤波后频域图0.30.350.4200150幅度1005000100200300400500频率/赫兹600700800900

图4-23 信号滤波后的时域图和频域图

从参考程序及以上的四个图像中可以得到如下结论：

①从原始信号波形的频域图可以看到其频率分量主要在500到900Hz之间，噪声的频率分量主要集中在250Hz，利用通带截频为300Hz的低通滤波器可以滤除噪声。对比图4-20和图4-23滤波前后的波形和频谱，可以看到波形得到了重现

②滤波器的采样频率为22050Hz，滤波器的阶数为266 ③滤波器的通带截频0.4，阻带截频0.3，过渡带宽0.1

结 论

论文正文主要简单介绍了Matlab、数字滤波器及利用matlab实现FIR滤波器的多种技术设计。

Matlab语言简洁紧凑，使用方便，库函数十分丰富。MATLAB

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

程序书写的形式自由，利用丰富的库函数避开了繁琐的子程序编程任务。

在数字信号处理中 ,由于信号中经常混有各种复杂成分，所以很多信号分析都是基于滤波器而进行的， 因此数字滤波器占有极其重要的地位 。数字滤波器分为有限冲激响应数字滤波器，即FIR数字滤波器和无限冲激响应，即IIR数字滤波器。我们主要介绍了FIR数字滤波器。

目前FIR滤波器的设计方法主要有三种：窗函数法、频率抽样法和优化设计方法。我们主要介绍前两种方法。

涉及FIR滤波器的多种技术设计。各种方法都有其优点和缺点,需根据不同的滤波器类型选择不同的方法。窗函数法在设计标准滤波器，例如低通、高通、带通，是很有用的。另一方面, 频率抽样法的优点是可以在频域直接设计，并且适合于最优化设计；缺点是抽样频率只能等于2pi/M的整数倍或等于2pi/M的整数倍上加上pi/M，因而不能确保截止频率Wc的自由取值。要想实现自由选择频率，则必须增加抽样点数M，但这种计算量加大。

本文实现了基于Matlab的数字低通和数字高通滤波器。 高通滤波器是容许高频信号通过、但减弱（或减少）频率低于截止频率信号通过的滤波器。对于不同滤波器而言,每个频率的信号的减弱程度不同。它有时被称为低频剪切滤波器；在音频应用中也使用低音消除滤波器或者噪声滤波器。低通滤波器与高通滤波器特性恰恰相反。

此次设计当中有很多问题困扰我，通过查阅资料、同学和邵霞老师的帮助逐步解决了问题，在此艰难的过程中让我懂得了很多。

万事开头难，不要畏惧，做好了开头也就成功了一半；其发现问题要立即解决问题；自己钻研不出来的，要敢于问问题；做事认真仔细，不要怕麻烦，否则只会更麻烦。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

致 谢

首先，在邵霞老师开的第一次会议上，邵老师把她自己的一些构想给我说了一下，这些构想使我对完成FIR滤波器的设计有了很大的信心，在此特别感谢邵老师对我的悉心指导。在与邵老师接触的过程中，她给我们的感觉就是我们就像朋友间的交流，毫无距离感。另外，邵老师思维缜密、知识渊博、生活上平易近人，使我从邵老师身上学到了很多专业的知识和做人的道理，这都是我以后做人做事的榜样。

其次，感谢通信教研室老师和领导，他们对我们在设计的过程中出现的问题同样给与了一些的建议，我也从他们身上学到了许多东西，这对我以后的生活与工作都有巨大的作用。

再则，我要感谢和我一组的同学，由于我们一组的几个同学用到的知识都有重叠，所以我们也在互帮互助同成长，在此过程中我深深体会到了团队的重要性，虽然我们不是一个团队，一个和谐团结的团队是一只非常可怕的队伍，众志成城，其利断金。

感谢和班里和系里的一些同学，他们也为我提供了无私的帮助，而且很重要，我们相互交流，相互学习，这个过程中我们一起解决了很多问题，同时也学到了很多新的知识。

最后我要感谢我的父母家人，他们对我的影响是最大的。 谢谢，谢谢你们！

参考文献

【1】陈后金.数字信号处理.高等教育出版社，2006

【2】罗军辉、罗勇江等. MATLAB7.0在数字信号处理中的应用

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

[M].北京:机械工业出版社，2005

【3】邹鲲，袁俊泉，龚享铱.Matlab 6.x信号处理.清华大学出版社.2002

【4】张明照，刘政波，刘斌.应用Matlab实现信号分析和处理.科学出版社.2006

【5】刘波，文忠.曾涯.Matlab信号处理.电子工业出版社.2006 【6】陈永春.MATLAB语言高级编程[M].北京：清华大学出版社，2003

【7】程佩青.数字信号处理教程.清华大学出版社.2007 【8】黄顺吉.数字信号处理及其应用.国防工业出版社.1982 【9】苏金明，张莲花，刘波.Matlab工具箱应用.电子工业出版社.2004

【10】邹理和.数字滤波器[M].北京：国防工业出版社，1979

附 录

附录一 外文原文及翻译

外文原文

FIR Filter Design Techniques

Abstract：

This report deals with some of the techniques used to design FIR filters. In the beginning, the windowing method and the frequency sampling methods are discussed in detail with their merits and demerits. Different optimization techniques involved in FIR filter design are also covered, including Rabiner’s method for FIR filter design. These optimization techniques reduce the error caused by frequency sampling technique at the non-sampled frequency points. A brief discussion of

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

some techniques used by filter design packages like Matlab are also included.

Introduction

FIR filters are filters having a transfer function of a polynomial in z and is an all-zero filter in the sense that the zeroes in the z-plane determine the frequency response magnitude characteristic.The z transform of a N-point FIR filter is given by

(1)

FIR filters are particularly useful for applications where exact linear phase response is required. The FIR filter is generally implemented in a non-recursive way which guarantees a stable filter. FIR filter design essentially consists of two parts (i) approximation problem (ii) realization problem

The approximation stage takes the specification and gives a transfer function through four steps. They are as follows:

(i) A desired or ideal response is chosen, usually in the frequency domain. (ii) An allowed class of filters is chosen (e.g.the length N for a FIR filters). (iii) A measure of the quality of approximation is chosen.

(iv) A method or algorithm is selected to find the best filter transfer function.

The realization part deals with choosing the structure to implement the transfer function which may be in the form of circuit diagram or in the form of a program.

There are essentially three well-known methods for FIR filter design namely:

(1) The window method

(2) The frequency sampling technique (3) Optimal filter design methods

The Window Method

In this method, [Park87], [Rab75], [Proakis00] from the desired

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

frequency response specification Hd(w), corresponding unit sample response hd(n) is determined using the following relation

(2)

(3)

In general, unit sample response hd(n) obtained from the above relation is infinite in duration, so it must be truncated at some point say n= M-1 to yield an FIR filter of length M (i.e. 0 to M-1). This truncation of hd(n) to length M-1 is same as multiplying hd(n) by the rectangular window defined as

w(n) = 1 0≦n≦M-1 (4)

0 otherwise

Thus the unit sample response of the FIR filter becomes

h(n) = hd(n) w(n) (5)

= hd(n) 0≦n≦M-1 = 0 otherwise

Now, the multiplication of the window function w(n) with hd(n) is equivalent to convolution of Hd(w) with W(w), where W(w) is the frequency domain representation of the window function

(6)

Thus the convolution of Hd(w) with W(w) yields the frequency response of the truncated FIR filter

(7)

The frequency response can also be obtained using the following relation

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

(8)

But direct truncation of hd(n) to M terms to obtain h(n) leads to the Gibbs phenomenon effect which manifests itself as a fixed percentage overshoot and ripple before and after an approximated discontinuity in the frequency response due to the non-uniform convergence of the fourier series at a discontinuity.Thus the frequency response obtained by using (8) contains ripples in the frequency domain. In order to reduce the ripples, instead of multiplying hd(n) with a rectangular window w(n), hd(n) is multiplied with a window function that contains a taper and decays toward zero gradually, instead of abruptly as it occurs in a rectangular window. As multiplication of sequences hd(n) and w(n) in time domain is equivalent to convolution of Hd(w) and W(w) in the frequency domain, it has the effect of smoothing Hd(w).

The several effects of windowing the Fourier coefficients of the filter on the result of the frequency response of the filter are as follows: (i) A major effect is that discontinuities in H(w) become transition bands between values on either side of the discontinuity.

(ii) The width of the transition bands depends on the width of the main lobe of the frequency response of the window function, w(n) i.e. W(w). (iii) Since the filter frequency response is obtained via a convolution relation , it is clear that the resulting filters are never optimal in any sense. (iv) As M (the length of the window function) increases, the mainlobe width of W(w) is reduced which reduces the width of the transition band, but this also introduces more ripple in the frequency response.

(v) The window function eliminates the ringing effects at the bandedge and does result in lower sidelobes at the expense of an increase in the width of the transition band of the filter.

Some of the windows [Park87] commonly used are as follows: 1. Bartlett triangular window:

W(n)=2(n+1)/N+1 n=0,1,2…….,(N-1)/2 (9)

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

=2-2(n+1)/N+1 n=(N-1)/2,……,N-1 = 0 otherwise 2 Generalized cosine windows

(Rectangular, Hanning, Hamming and Blackman)

W(n)=a-bcos(2p(n+1)/(N+1))+ccos(4p(n+1)/(N+1)) n=0,1….N-1 (10)

= 0 otherwise 3.Kaiser window with parameter β:

(11)

The general cosine window has four special forms that are commonly used. These are determined by the parameters a,b,c

TABLE I

Value of coefficients for a,b and c from [Park87]

Window Retangular Hanning a 1 0.5 b 0 0.5 c 0 0 0 Hamming 0. 0.46 Blackman 0.42 0.5 0.08 The Bartlett window reduces the overshoot in the designed filter but spreads the transition region considerably.The Hanning,Hamming and Blackman windows use progressively more complicated cosine functions to provide a smooth truncation of the ideal impulse response and a frequency response that looks better. The best window results probably come from using the Kaiser window, which has a parameter . that allows adjustment of the compromise between the overshoot reduction and transition region width spreading.

The major advantages of using window method is their relative simplicity as compared to other methods and ease of use. The fact that well defined equations are often available for calculating the window

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

coefficients has made this method successful.

There are following problems in filter design using window method: (i) This method is applicable only if Hd(w) is absolutely integrable i.e only if (2) can be evaluated. When Hd(w) is complicated or cannot easily be put into a closed form mathematical expression, evaluation of hd(n) becomes difficult.

(ii) The use of windows offers very little design flexibility e.g. in low pass filter design, the passband edge frequency generally cannot be specified exactly since the window smears the discontinuity in frequency. Thus the ideal LPF with cut-off frequency fc, is smeared by the window to give a frequency response with passband response with passband cutoff frequency f1 and stopband cut-off frequency f2.

(iii) Window method is basically useful for design of prototype filters like lowpass,highpass,bandpass etc. This makes its use in speech and image processing applications very limited. The Frequency Sampling Technique

In this method, [Park87], [Rab75], [Proakis00] the desired frequency response is provided as in the previous method. Now the given frequency response is sampled at a set of equally spaced frequencies to obtain N samples. Thus , sampling the continuous frequency response Hd(w) at N points essentially gives us the N-point DFT of Hd(2pnk/N). Thus by using the IDFT formula, the filter co-efficients can be calculated using the following formula

(12)

Now using the above N-point filter response, the continuous frequency response is calculated as an interpolation of the sampled frequency response. The approximation error would then be exactly zero at the sampling frequencies and would be finite in frequencies between them. The smoother the frequency response being approximated, the smaller will be the error of interpolation between the sample points.

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

One way to reduce the error is to increase the number of frequency samples [Rab75]. The other way to improve the quality of approximation is to make a number of frequency samples specified as unconstrained variables. The values of these unconstrained variables are generally optimized by computer to minimize some simple function of the approximation error e.g. one might choose as unconstrained variables the frequency samples that lie in a transition band between two frequency bands in which the frequency response is specified e.g. in the band between the passband and the stopband of a low pass filter.

There are two different set of frequencies that can be used for taking the samples. One set of frequency samples are at fk = k/N where k = 0,1,….N-1. The other set of uniformly spaced frequency samples can be taken at fk =(k+1/2)/N for k = 0,1,….N-1.

The second set gives us the additional flexibility to specify the desired frequency response at a second possible set of frequencies. Thus a given band edge frequency may be closer to type-II frequency sampling point that to type-I in which case a type-II design would be used in optimization procedure.

In a paper by Rabiner and Gold [Rabi70], Rabiner has mentioned a technique based on the idea of frequency sampling to design FIR filters. The steps involved in this method suggested by Rabiner are as follows: (i) The desired magnitude response is provided along with the number of samples,N . Given N, the designer determines how fine an interpolation will be used.

(ii) It was found by Rabiner that for designs they investigated, where N varied from 15 to 256, 16N samples of H(w) lead to reliable computations, so 16 to 1 interpolation was used.

(iii) Given N values of Hk , the unit sample response of filter to be designed, h(n) is calculated using the inverse FFT algorithm.

(iv) In order to obtain values of the interpolated frequency response two procedures were suggested by Rabiner. They are

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

(a) h(n) is rotated by N/2 samples(N even) or (N-1)/2 samples for N odd to remove the sharp edges of impulse response, and then 15N zero-valued samples are symmetrically placed around the impulse response. (b) h(n) is split around the N/2nd sample, and 15N zero-valued samples are placed between the two pieces of the impulse response.

Merits of frequency sampling technique

(i) Unlike the window method, this technique can be used for any given magnitude response.

(ii) This method is useful for the design of non-prototype filters where the desired magnitude response can take any irregular shape.

There are some disadvantages with this method i.e the frequency response obtained by interpolation is equal to the desired frequency response only at the sampled points. At the other points, there will be a finite error present.

Optimal Filter Design Methods

Many methods are present under this category. The basic idea in each method is to design the filter coefficients again and again until a particular error is minimized. The various methods are as follows: (i) (ii) (iii)

Least squared error frequency domain design

Nonlinear equation solution for maximal ripple FIR filters Polynomial interpolation solution for maximal ripple FIR filters

Least squared error frequency domain design

As seen in the previous method of frequency sampling technique there is no constraint on the response between the sample points, and poor results may be obtained.

The frequency sampling technique is more of an interpolation method rather than an approximation method. This method [Rab75], [Parks87] controls the response between the sample points by considering a number of sample points larger than the order of the filter.The purpose of most filters is to separate desired signals from undesired signals or noise. As the energy of the signal is related to the square of the signal, a squared error approximation criterion is appropriate to optimize the design of the FIR filters.

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

The frequency response of the FIR filter is given by (8) for a N-point FIR filter. An error function is defined as follows

(13)

where wk=(2*pi*k)/L and Hd(wk) are L samples of the desired response, which is the error measure as a sum of the squared differences between the actual and desired frequency response over a set of L frequency samples. The method consists of the following steps:

(i) First ‘L’ samples from the continuous frequency response are taken, where L>N(length of the impulse response of filter to be designed) (ii) Then using the following formula

(14)

the L-point filter impulse response is calculated.

(iii) Then the obtained filter impulse response is symmetrically truncated to desired length N.

(iv) Then the frequency response is calculated using the following relation

(15)

(v) The magnitude of the frequency response at these frequency points for wk =(2k)/L will not be equal to the desired ones , but the overall least square error will be reduced effectively this will reduce the ripple in the filter response.

To further reduce the ripple and overshoot near the band edges, a transition region will be defined with a linear transfer function. Then the L frequency samples are taken at wk =(2*pi*k)/L using which the first N samples of the filter are calculated using the above method.Using this method, reduces the ripple in the interpolated frequency response. Nonlinear Equation solution for maximal ripple FIR filters

The real part of the frequency response of the designed FIR filter can

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

be written as a(n)cos(wn) [Rab75] where limits of summation and a(n) vary according to the type of the filter. The number of frequencies at which H(w) could attain an extremum is strictly a function of the type of the linear phase filter i.e. whether length N of filter is odd or even or filter is symmetric or anti-symmetric. At each extremum, the value of H(w) is predetermined by a combination of the weighting function W(w), the desired frequency response, and a quantity that represents the peak error of approximation distributing the frequencies at which H(w) attains an extremal value among the different frequency bands over which a desired response was being approximated. Since these filters have the maximum number of ripples, they are called maximal ripple filters. This method is as follows:

At each of the Ne unknown external frequencies, E(w) attains the maximum value of either and E(w) or equivalently H(w) has zero derivative. Thus two Ne equations of the form

H(wi)= ±δ/W(w) + D(wi) (16) d/dw {H(wi)} at w=wi=0 (17)

are obtained.

These equations represent a set of 2Ne nonlinear equations in two Ne unknowns, Ne impulse response coefficients and Ne frequencies at which H(w) obtains the extremal value. The set of two Ne equations may be solved iteratively using nonlinear optimiation procedure.

An important thing to note is that here the peak error (δ) is a fixed quantity and is not minimized by the optimization scheme. Thus the shape of H(w) is postulated apriori and only the frequencies at which H(w) attains the extremal values are unknown.

The disadvantage of this method is that the design procedure has no way of specifying band edges for the different frequency bands of the filter. Thus the optimization algorithm is free to select exactly where the bands will lie.

Polynomial Interpolation Solution for Maximal Ripple FIR filters This algorithm [Rab75] is basically an iterative technique for producing

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

a polynomial H(w) that has extrema of desired values. The algorithm begins by making an initial estimate of the frequencies at which the extrema in H(w) will occur and then uses the well-known Lagrange interpolation formula to obtain a polynomial that alternatively goes through the maximum allowable ripple values at these frequencies. It has been experimentally found that the initial guess of extremal frequencies does not affect the ultimate convergence of the algorithm but instead affects the number of iterations required to achieve the desired result. Let us consider the case of design of a low pass filter using the above method.

Fig. 1. Iterative solution for a maximum ripple lowpass filter from [Rab75]

The Fig. 1 shows the response of a lowpass filter with N = 11.The number of extremal frequencies i.e. the frequencies where ripples occur are 6 in this case. They are divided into 3 passband extrema and 3 stopband

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

extrema. The filled dots indicate the initial guess as to the extremal frequencies of H(w). The solid line is the initial Lagrange polynomial obtained by choosing polynomial coefficients so that the values of the polynomial at the guessed set of frequencies are identical to the assigned extreme values. But this polynomial has extrema that exceeds the specified maxima values. The next stage of the algorithm is to locate the frequencies at which the extrema of the first Lagrange interpolation occur. These frequencies are now used as the new frequencies for which the extrema of the filter response occur. This second set of frequencies are indicated by open dots in Fig. 1. Now similarly the new set of frequencies are taken as those frequencies where the maximum exceeds the specified maxima. Thus the method is completely iterative in nature. Conclusions

The report has described the various techniques involved in the design of FIR filters. Every method has its own advantages and disadvantages and is selected depending on the type of filter to be designed.The window method is basically used for the design of prototype filters like the low-pass, high-pass, band-pass etc. They are not very suitable for designing of filters with any given frequency response. On the other hand, the frequency sampling technique is suitable for designing of filters with a given magnitude response. The ideal frequency response of the filter is approximated by placing appropriate frequency samples in the z -plane and then calculating the filter co-efficients using the IFFT algorithm. The disadvantage of the frequency sampling technique was that the frequency response gave errors at the points where it was not sampled. In order to reduce these erros the different optimization technique for FIR filter design were presented wherein the remaining frequency samples are chosen to satisfy an optimization criterion. An appendix consisting of the filter design methods used by the software package Matlab is also presented.

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

外文翻译

FIR滤波器设计技术

这份报告列举了一些设计FIR滤波器所使用的技术。首先讨论了窗函数法和频率抽样法的优点和缺点。FIR数字滤波器也包含了许多优化设计的方法，这些优化技术减少了在频率采样时非采样频率点的误差频率。对于用于设计数字滤波器的技术，例如matlab，进行了简明扼要的探讨。

介绍

FIR滤波器的系统函数是一个z的多项式，因FIR滤波器的频率响应是频率的实函数，也称其为零相位滤波器。N阶FIR滤波器的系统函数表示为

1 （1）

FIR滤波器是十分重要的，可应用于精确线性相位相应。FIR滤波器的实现方式保证了它是一个稳定的滤波器。

FIR滤波器的设计可分为两部分：近似问题、实现问题 解决近似问题，要通过四个步骤找出传递函数：

在频域内找出期望的或最理想的反应 选择滤波器的阶数（FIR滤波器的长度N） 选择近似结果中较好的

选择一种算法寻找最优的滤波器传递函数

选择部分结构处理实现传递函数的形式可能是线路图或程序。 本质上来说，有三种著名的FIR滤波器设计方法： 窗函数法、频率抽样法、滤波器的优化设计 窗函数法：

在该方法中，[Park87],[Rab75], [Proakis00]从理想的频率响

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

应Hd（w）出发，其对应的单位脉冲相应关系如下：

（2）

（3）

一般来说，单位脉冲相应hd（n）的持续时间是无限的，所以在某种程度上说，它必须截断。n=M-1约束着FIR滤波器的长度M。以M-1截断的hd（n）乘以窗函数就得到了滤波器的单位脉冲响应。

矩形窗口的定义为 w(n) （4）

= 0 其它

FIR滤波器的单位脉冲相应为 h(n) （5）

= hd(n) 0≦n≦M-1 = 0 其它

现在,多元化的窗函数w(n)与hd(n)相当于Hd(w)与W（w）的卷积，其中，W（w）是窗函数的频域表示

=

hd(n)

=

1

0

≦

n

≦

M-1

w(n)

(6)

因此Hd(w)与W(w)的卷积为FIR数字滤波器的截断后的频率响应

(7)

频率响应也可以利用以下的关系式

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

(8)

由于非均匀收敛的傅里叶级数的不连续性，其自身的波纹前后有一种近似于不连续的频率响应，因此直接截断的hd(n)来获得h(n)将导致吉布斯现象。与此同时，利用（8）得到的频率响应在频域内有波纹的振荡。为了减少波纹，hd(n)不是乘以一个矩形窗口w(n)，而是乘以一个含有圆锥和逐渐衰减到零的窗口。作为主体的序列的hd(n)和w(n)在时域内的卷积相当于其在频域内的乘积，其效果是平滑的。

滤波器的窗函数的傅里叶系数对滤波结果的频率响应的影响如下：

（i）一个主要的结果就是过渡带的不连续的两边出现中断 （ii）过渡带的宽度取决于窗函数的频率响应的主瓣宽度

（iii）滤波器的频率响应是通过卷积关系得到的，可以肯定的是，由产生的滤波器绝不是最佳的

（iv）随着M的增加，其主瓣宽度降低从而降低了过渡带的宽度，但是这也过滤掉了更多的脉冲频率响应。

（v）窗函数消除边缘响应引起的效果，并以较低的旁瓣代价增加过渡带的宽度

常用的窗函数如下[Park87]: 1.Bartlett 三角窗：

W(n)= 2(n+1)/N+1 n=0,1,2…….,(N-1)/2 (9)

= 2-2(n+1)/N+1 n=(N-1)/2,……,N-1

= 0 其它 2.广义余弦窗

(Rectangular, Hanning, Hamming and Blackman)

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

W(n)=a-bcos(2p(n+1)/(N+1))+ccos(4p(n+1)/(N+1))

n=0,1….N-1 (10)

=0 其它 3.Kaiser 窗

（11）

一般的余弦窗有四种常用的形式。以下就是确定的参数a,b,c

表一

从[Park87]得到的系数 Window Retangular Hanning Hamming a 1 b 0 c 0 0 0 0.5 0.5 0. 0.46 Blackman 0.42 0.5 0.08 Bartlett 窗函数的设计减少了信息的误差，但其过渡带较宽。Hanning,Hamming 和Blackman 窗的使用会更好，它们可以用于复杂的余弦函数，并可以提供理想的光滑截断的脉冲响应和频率响应。研究的结果表明，最佳的窗函数可能是有一个参数的Kaiser 窗，它可以实现衰减和过渡带宽度的妥协。

窗函数的主要优点是它们比起其它方法更加的简单，且易于使用。事实上，计算窗函数的明确的方程系数就可以成功的使用该方法。

在使用窗函数来设计滤波器时，会遇到三个问题：

1.该方法只适用于Hd(w)是绝对可积的情况，即只有（2）式可以评估。当Hd(w)是复杂的或不能轻易被评价的闭合形式，写出Hd(n)的数学表达式就变得困难了。

欧阳德创编 2021.03.07

欧阳德创编 2021.03.07

2.使用窗函数的灵活性比较差，例如，在低通滤波器的设计中，通频带的边缘频率一般不能用窗口完全掠过不连续区域。因此理想低通滤波器的截止频率，是通带截止频率f1和阻带截止频率f2相关的一个频率响应。

3.窗函数法在设计标准滤波器，例如低通、高通、带通，是很有用的。但这也使其在语音、图像处理的程序上的应用是十分有限的。

频率抽样技术:

在该方法中，[Rad75]、[Park87]、[Proakis00]是按照前面的方法提供理想的频率响应。现在，是在给定的频率响应中取一系列等间隔的频率，用以得到N的抽样。因此，采样频率的响应Hd(w)在其本质上是给了我们Hd(2pnk/N)。因此，利用该滤波器可以计算出下面的公式：

(12)

现在使用上述的N阶滤波器响应、连续性频率响应作为计算差值采样频率响应的方法。其近似误差就会完全接近于零点采样频率的误差，并将它们之间的频率有限化。越平滑的频率响应将会越近似，存在于样本点间的差值误差会很小。

时间：2021.03.07 创作：欧阳德 欧阳德创编 2021.03.07