基于MATLAB的FIR+数字低通滤波器设计 - 图文

毕业设计（论文）任务书

题目：基于MATLAB的FIR 数字低通滤波器设计系 名 信息工程系 专 业 电子信息工程 学 号 6009202371 学生姓名 马成 指导教师 李晓峰 职 称 讲师

2012年 12 月 15 日

一、原始依据（包括设计或论文的工作基础、研究条件、应用环境、工作目

的等。）

研究条件：在大学四年专业学习的基础上，阅读有关MATLAB软件使用方法以及数字滤波器设计等方面的书籍；掌握MATLAB编程语言，熟练利用计算机进行MATLAB仿真设计。

应用环境：数字滤波器在现实生活中与人们息息相关，广泛使用于各种声音、图像以及文字等处理系统中。将MATLAB强大的运算处理能力有机融入数字滤波器设计中可实现对于数字滤波器的快速设计以及各种处理变换。

工作目的：本课题的主要任务就是利用MATLAB软件中的数字信号处理工具箱实现FIR低通数字滤波器的设计。

二、参考文献

[1]徐明远，刘增力，《MATLAB仿真在信号处理中的应用》[M].西安：西安电子科技大学出版社，2007.11.

[2]陈桂明，张明照，《应用MATLAB语言处理信号与数字图像》[M].北京：科学

出版社，2000.5.

[3]飞思科技产品研发中心.MATLAB基础与提高[M].北京：电子工业出版社，2005.

[4]高西全，丁玉美.数字信号处理[M].西安：西安电子科技大学出版社，第三版，2008.

[5]李亚奇，张雅琦.线性相位FIR数字滤波器[J].电子测量技术，2005（6）：35-37.

[6]郭德才.基于Matlab的FIR低通滤波器的设计与仿真[J].通化师范学院学报，2009, 30（8）：38-41.

[7]赵刚.基于数字滤波器设计的讨论[J].南开大学学报（自然科学版），2003（3）：101-103.

[8]陈明军.改进窗函数在FIR数字滤波器设计中的应用[J].继电器，2007（13）：65-67.

三、设计（研究）内容和要求（包括设计或研究内容、主要指标与技术参数，并根据课题性质对学生提出具体要求。）

1、研究内容：研究目前利用MATLAB应用于数字信号处理方面设计数字滤波器的内容，熟练掌握MATLAB语言和数字滤波器设计方法，实现基于MATLAB的数字低通滤波器设计。

2、主要技术指标及设计具体要求：

本设计要求在熟练应用MATLAB软件的基础上，采用目前常用的窗函数法、频率采样法和等波纹最佳逼近法等FIR 数字滤波器的设计方法设计一个数字低通滤波器，并对设计结果进行比较分析，研究它们各自的优缺点及适用对象。

指导教师（签字）

年 月 日

审题小组组长（签字）年

月 日

天津大学仁爱学院本科生毕业设计（论文）开题报告

课题名称 系 名 学生姓名 基于MATLAB的FIR 数字低通滤波器设计 信息工程系 马成 专 业 指导教师 电子信息工程 李晓峰 一、课题来源及意义 数字滤波技术作为数字信号处理的基本分支之一，就是提取信号的有用分量，削弱无用分量的技术，被广泛应用于数据处理，图像处理、雷达、声纳信号处理、地址石油勘探等很多领域，越来越受到人们的关注。由于单位冲击响应的不同数字滤波器有两种类型：有限冲击响应（Finite Impulse Response，FIR）数字滤波器和无限冲击响应（Infinite Impulse Response，IIR）数字滤波器。两种类型滤波器相比而言，对于同样的滤波器设计指标，虽然FIR 滤波器成本较高，信号延迟较大并且FIR 滤波器没有现成的计算公式（必须要用计算机辅助设计软件（如MATLAB）来计算），但是FIR 滤波器可以采用FFT算法，运算速度较快；精度高，具有严格的线性相位等特点优于IIR 数字滤波器已被广泛应用。本课题设计的就是基于MATLAB的FIR 数字低通滤波器设计。 二、国内外发展现状 无论是在理论研究上还是在如语音、数字音频、图像处理、通讯、雷达、军事、航空航天、医疗等实际应用上都有着美好的技术前景和巨大的使用价值。 采用数字技术则避免很多类似的难题，如模拟电路元件对温度的敏感性等等 数字滤波器在其他方面也有很多突出的优点都是模拟技术所不能及的，所以采用数字滤波器对信号进行处理是目前的发展方向。 三、研究目标 主要研究基于MATLAB的FIR低通滤波器。具体是采用窗函数设计法、频率采样法、等波纹逼近法进行设计，并用MATLAB软件编写程序进行仿真

四、研究内容 1、研究FIR滤波器的定义、分类、应用以及设计方法。 2、了解FIR滤波器窗函数的设计方法及原理。 3、了解FIR滤波器频率采样的设计方法及原理。 4、了解FIR滤波器等波纹逼近的设计方法及原理。 5、确定滤波器的技术指标，并进行MATLAB仿真 五、研究方法与手段 本课题主要研究的是基于Matlab的FIR低通滤波器设计,因此在研究过程中主要是采用理论分析结合软件仿真来实现。总的研究方法是在课题研究之前.掌握一定的理论基础,了解FIR滤波器是如何进行信号各个波段的滤波的,在理论中寻找适合本课题的设计方案,然后了解Matlab软件应用,使用Matlab软件设计数字低通滤波器并进行仿真优化,实现数字低通滤波器的滤波。具体措施如下： 1. 首先阅读广泛资料,对FIR滤波器的定义、设计方法做一下介绍,分析FIR低通滤波器的研究意义及作用。 2．Matlab作为工作平台和开发工具,熟悉Matlab软件程序设计的基本原理和实验环境,实现仿真。 3．针对FIR低通滤波器的滤波基本原理进行研究,查阅相关资料书籍。并利用窗函数的设计方法进行滤波设计,并运用Matlab进行模拟仿真。 4．查漏补缺六、进度安排 1、2012.12.15-2013.03.03 查找国内外相关资料，完成开题报告 2、2013.03.07-2013.04.07 重点学习基于Matlab的FIR低通滤波器的设计方案.实现数字低通滤波器的设计与仿真。 3、2013.04.08-2013.05.08 对设计方案的软件进行仿真，确定最终的方案 4、2013.05.08-2013.06.01 完成毕业论文，准备答辩. 八、主要参考文献 [1]徐明远，刘增力，《MATLAB仿真在信号处理中的应用》[M].西安：西安电子科技大学出版社，2007.11. [2]陈桂明，张明照，《应用MATLAB语言处理信号与数字图像》[M].北京：科学出版社，2000.5. [3]飞思科技产品研发中心.MATLAB基础与提高[M].北京：电子工业出版社，2005. 进行论文的改进。

[4]高西全，丁玉美.数字信号处理[M].西安：西安电子科技大学出版社，第三版，2008. [5]李亚奇，张雅琦.线性相位FIR数字滤波器[J].电子测量技术，2005（6）：35-37.[6]郭德才.基于Matlab的FIR低通滤波器的设计与仿真[J].通化师范学院学报，2009, 30（8）：38-41. [7]赵刚.基于数字滤波器设计的讨论[J].南开大学学报（自然科学版），2003（3）：101-103. [8]陈明军.改进窗函数在FIR数字滤波器设计中的应用[J].继电器，2007（13）：65-67. [9]丁玉美, 高西全.数字信号处理 [M]. 第三版. 西安: 西安电子科技大学出版社, 2008, 6. [10]闫胜利. FIR滤波器及设计原理[J]. 长春工程学院学报(自然科学版), 2003, 6, 4(1): 21-24. [11]姚齐国. 基于MATLAB的数字滤波器的设计[J]. 江西理工大学学报, 2006, 2, 27(1): 50-52. [12]杨守卫. FIR数字滤波器应用分析探讨[J]. 河北省工程咨询院学报, 2011, 7, 29(15): 47-49. [13]朱敏. MATLAB数字信号处理工具箱的开发和应用[J]. 信息与电脑, 2010, 2, 26(8): 154-155. [14]姚海燕. FIR数字滤波器设计窗函数法与频率抽样法比较[J].安阳工学院学报,2007,6, 12(6): 51-53. [15]刘波. MATLAB信号处理[M]. 北京: 电子工业出版社, 2006, 7. 选题是否合适： 是□ 否□ 课题能否实现： 能□ 不能□ 指导教师（签字） 年 月 日 选题是否合适： 是□ 否□ 课题能否实现： 能□ 不能□ 审题小组组长（签字） 年 月 日

毕业设计（论文）说明书

题目：

MATLAB的FIR数字低通滤波器设计系 名 信息工程系 专 业 电子信息工程 学 号 6009202371 学生姓名 马成 指导教师 李晓峰

2013 年 6 月 9 日

基于

摘 要

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

关键词：FIR数字滤波器；窗函数法;频率抽样法;

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 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；windowing method；frequency;sampling;method；

目 录

第一章

绪论 .............................................................................................................. 1

1.1 课题来源及意义 ............................................................................................. 1 1.2 国内外发展现状 ............................................................................................. 1 1.3 研究目标 ............................................................................................................. 1 1.4 研究内容 ............................................................................................................. 1

第二章 数字滤波器线性相位条件 ........................................................... 2

2.1 FIR数字滤波器概念 .................................................................................... 2 2.2 FIR数字滤波器的线性线性相位定义 ............................................... 3 2.3 FIR数字滤波器线性相位时域约束条件 .......................................... 3

第三章 MATLAB简介 ............................................................................................ 4

3.1 MATLAB基本功能 ............................................................................................ 4 3.2 MATLAB的优势及特点.................................................................................. 4

第四章 FIR数字滤波器的设计 .................................................................. 6

4.1 窗函数法设计FIR数字滤波器 .............................................................. 6 4.2 利用频率采样法设计FIR数字滤波器 .............................................. 9 4.3 利用等波纹最佳逼近法设计FIR数字滤波器 ............................ 11 5.1 窗函数法的MATLAB实现 ......................................................................... 13 5.2 频率抽样法的MATLAB实现 .................................................................... 17 5.3 利用滤波器处理加有噪声的音频波形 ............................................ 20

第五章 利用MATLAB实现FIR滤波器设计..................................... 13

第六章 总结与展望 .......................................................................................... 24 附录 ..................................................................................................................................... 26 外文资料 中文译文 致谢

参考文献 .......................................................................................................................... 25

1

天津大学仁爱学院2013届本科生毕业生设计（论文）

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

4.1 窗函数法设计FIR数字滤波器

要设计出的滤波器的理想频率响应函数为

1 hd(n)?2?，则对应的单位脉冲响应为

??H??d(ej?)ej?nd? (4-1)

窗函数设计法的基本原理是用有限长单位脉冲响应序列h(n)逼近hd(n)。由于hd(n)往往是无限长序列，且是非因果的，所以用窗函数ω(n)将hd(n)截断，并进行加权处理，得到

N?1n?0 （4-2）

h(n)就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数 H(e)??h(n)e?j?n (4-3)

j?为（4-3）式中，N为所选窗函数ω(n)的长度。

用窗函数法设计的滤波器性能取决于窗函数ω(n)的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。 出

这样选定窗函数类型和长度N后，求出单位脉冲响应

，并求

。是否满足要求，要进行验算。一般在h(n)尾部加零使长度满足2的整

数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。不满足要求，则要重新选择窗函数类型和长度N，再次验算，直至满足要求。

如果要求线性相位特性，则h(n)还必须满足：h(n)=+/-h(N-1-n)，根据式中的正负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类[8]。 下面介绍下典型的窗函数：

常见的窗函数有矩形窗（Rectangle Window）、三角形窗（Bartlerr Window）、汉宁（Hanning）窗——升余弦窗、哈明（Hamming）窗——改进的升余弦窗、布莱克曼（Blackman）窗、凯塞—贝塞尔窗（Kaiser-Basel Window） 1矩形窗（Rectangle Window）

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

6

天津大学仁爱学院2013届本科生毕业生设计（论文）

N?1?1n?? ?R(n)??2 （4-5）

??0otherwise其频谱的幅度函数为：

WRk(?)sin(?N/2)sin(?/2) （4-6）

2三角形窗（Bartlett Window）

三角窗亦称费杰（Fejer）窗，是幂窗的一次方形式。与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。 三角形窗的窗函数为：

?2n??N?1 ?B(n)???2?2n?N?1?0?n?1(N?1)21(N?1)?n?N?122 (4-7)

其频谱的幅度函数为：

2 WBg(?)?N?sin(?N/4)? （4-8） ???sin(?/2)?3汉宁（Hanning）窗------升余弦窗

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

汉宁窗的窗函数为：

2n???0.5?0.5cos ?H(n)??N?1??00?n?Notherwise (4-9)

4哈明（Hamming）窗————改进的升余弦窗

这种改进的升余弦窗，能量更加集中在主瓣中，主瓣的能量约占99.96%，瓣峰值幅度为40dB，但其主瓣宽度和汉宁窗的相同，仍为8π/N.可见哈明窗是一种高效窗函数，所以Matlab窗函数设计函数的默认窗函数就是哈明窗。 哈明窗的窗函数为：

7

天津大学仁爱学院2013届本科生毕业生设计（论文）

2n??0.54?0.64cos??H(n)??N?1??00?n?Notherwise (4-10)

其幅度函数为：

5布莱克曼（Blackman）窗

(4-11)

该窗函数位移不同，幅度函数也不同，会使旁瓣进一步抵消，主瓣宽度为 12

π/N。

布莱克曼窗的窗函数为：

2?n4?n?(0.42?0.5cos?0.8cos)? ?B(n)??N?1N?1??0N?12 (4-12) otherwisen?其频谱的幅度函数为：

WB(?)?0.42WR(?)?0.25[WR(??2?2?)?WR(??)]N?1N?1 (4-13) 4?4??0.04[WR(??)?WR(??)]N?1N?1

6凯撒——贝赛尔窗（Kaiser-Basel Window）

凯塞窗是一种最优窗函数，不同于前面五种窗函数，凯塞窗是一种参数可调的窗函数，其函数形式如下：

?K(n)?Io(?)Io(?)0?n?N?1 (4-14)

其中

2n ???1?(?1)2N?11?Io(?)?1??[?]2 (4-15)

2k?1k?一般取15-25项可以满足精度要求。参数可以控制窗的形状。一般越

大，主瓣越宽，而旁瓣幅度会随之减小，典型的是介绍6种窗函数的基本参数

数据在4到9之间。表4-1

8

天津大学仁爱学院2013届本科生毕业生设计（论文）

表4-1 6种窗函数的基本参数

窗函数类型 矩形窗 三角窗 汉宁窗 哈明窗 布莱克曼窗 凯赛窗

分瓣峰值（dB）

-13 -25 -31 -41 -57 -57 过度带宽度(P/N) 阻带最小衰减（dB）

4 -21 8 -25 8 -44 8 -53 12 -74 7.442 -80

4.2 利用频率采样法设计FIR数字滤波器

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

Hd(ej?)加以等间隔抽样，并以此作为实际FIR滤波器的频率响应。设所需要的

滤波器的频率响应为Hd(ej?)。

现要求设计一个M阶的FIR滤波器h[k]，使得Hd(ej?)在M+1个抽样点上，FIR滤波器的频率响应Hd(ej?)与所需的频率响应Hd(ej?)相等，即 Hd(e)?H(ej?j?m)??h[k]e?jk?,m?0,1,....,M (4-16)

k?0MHd(ej?)由设计的要求给定，h[k]需要通过设计来确定。如果M+1个方程是线性

无关的，则可以通过求解M+1阶线性方程来得出FIR滤波器的h[k]。对?m的一些特殊抽样方法，上述方程的解可以直接由IDFT得到。由于要求设计出的滤波器是实系数的线性相位FIR滤波器，所以Hd(ej?)的抽样值还需要满足线性相位滤波器的约束条件。

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

??。为了使设计出的滤波器具有线性相位，Hd(ej?)在M+1个抽样点上的值Hd[m]应为

Hd[m]?Hd(ej2?M?1π2)?eeAd[m]j??jM?mM?1Ad(

?ej?e?j2?m)M?1M?mM?1 (4-17)

下面分别讨论四种线性相位滤波器在抽样点Hd[m]上的值：

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

9

天津大学仁爱学院2013届本科生毕业生设计（论文）

?M?j?M?1?(0?m?M/2) (4-18) Hd[m]??Ad[m]e??Hd[M?1?m](M/2?1?M)上式表明I型线性相位FIR滤波器Hd[m]在M/2?1?m?M的值可由

Hd[m]在1?m?M/2的值确定。在Hd[m]的值确定以后，对Hd[m]做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) (4-19)

?H[M?1?m](M?3)/2?1?m?M)?d?上式表明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?1 Hd[m]??jAd[m]e(1?m?M/2) (4-20)

?H[M/2?1?m?M)d??上式表明III型滤波器线性相位FIR滤波器Hd[m]在M/2?1?m?M的值可由

Hd[m]在1?m?M/2的值来确定。

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

?0(m?0)??M?jm?M?1(1?m?(M?1)/2) (4-21) Hd[m]??jAd[m]e?H[M?1?m](M?3)/2?m?M)d??上式表明IV型线性相位FIR滤波器Hd[m]在的值可(M?3)/2?m?M由Hd[m]在1?m?(M?1)/2的值确定。

对Hd(ej?)进行频率抽样，就是在z平面单位圆上的N个等间隔点上抽样出频率响应值。在单位圆上可以有两种抽样方式，第一种是第一个抽样点在w=0处，第二种是第一个抽样点在w=?/M处，每种方式可分为M为偶数与M为奇数

10

天津大学仁爱学院2013届本科生毕业生设计（论文）

π(2)N=33，Wc=rad，用4种窗函数设计线性相位低通滤波器。

4绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，

比较4种窗函数对滤波特性的影响。

图5-3 矩形窗设计线性相位低通滤波器

图5-4 汉宁窗设计线性相位低通滤波器

16

天津大学仁爱学院2013届本科生毕业生设计（论文）

图5-5 哈明窗设计线性相位低通滤波器

图5-6 布莱克曼设计线性相位低通滤波器

从图中可以看出用各种窗函数设计的FIR滤波器的阻带最小衰减及过渡带均符合要求，而且在通带内均为严格线性相位。

5.2 频率抽样法的MATLAB实现

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

17

天津大学仁爱学院2013届本科生毕业生设计（论文）

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

②阻带衰减大于等于15dB

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

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

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

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

信号滤波前和滤波后的时域图和频率图如图5-8、5-9所示

18

天津大学仁爱学院2013届本科生毕业生设计（论文）

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

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

19

天津大学仁爱学院2013届本科生毕业生设计（论文）

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

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

调用信号产生函数xtg产生具有加性噪声的信号xt，并自动显示xt及其频谱，如图5-10所示；

图5-10 信号加噪声的波形和频谱

（1）设计一个低通滤波器，要求从高频噪声中提取xt中的单频调幅信号，要求信号幅频失真小于0.1dB，将噪声频谱衰减60dB。先观察xt的频谱，确定滤波器指标参数。

根据图5-10（b）和实验要求，可选择滤波器指标参数：通带截止频率fp=120Hz，阻带截至频率fs=150Hz，换算成数字频率，通带截止频率

?p?2?fp??0.24?，通带最大衰为0.1dB，阻带截至频率?s?2?fs??0.3?，阻带最

小衰为60dB。

调用信号产生函数xtg产生具有加性噪声的信号xt程序框图如图5-11所示：

20

天津大学仁爱学院2013届本科生毕业生设计（论文）

Fs=1000，T=1/Fs xt=xtg 产生信号xt, 并显示xt及其频谱 用窗函数法或等波纹最佳逼近法 设计FIR滤波器hn 对信号xt滤波：yt=fftfilt(hn,xt) 1、计算并绘图显示滤波器损耗函数 2、绘图显示滤波器输出信号yt End

图5-11 实验程序框图

提取单频信号后加噪声的波形和频谱如图5-12所示：

图5-12 提取单频信号后加噪声的波形和频谱

21

天津大学仁爱学院2013届本科生毕业生设计（论文）

通带截止频率fp=120Hz，阻带截至频率fs=150Hz。代入采样频率Fs=1000Hz，换算成数字频率，通带截止频率?p?2?fpT?0.24?，通带最大衰为0.1dB，阻带截至频率?s?2?fsT?0.3?，阻带最小衰为60dB。所以选取blackman窗函数。与信号产生函数xtg相同，采样频率Fs=1000Hz。

（2）①根据滤波器指标选择合适的窗函数，计算窗函数的长度N，调用MATLAB函数fir1设计一个FIR低通滤波器。并编写程序，调用MATLAB快速卷积函数fftfilt实现对xt的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性图和时域波形图。

②滤波器指标不变，但改用等波纹最佳逼近法，调用MATLAB函数remezord和remez设计FIR数字滤波器。并比较两种设计方法设计的滤波器阶数。

图5-13 窗函数法滤除噪声后的信号波形

图5-14 等波纹逼近法滤除噪声后的信号波形

22

天津大学仁爱学院2013届本科生毕业生设计（论文）

用窗函数法设计滤波器，滤波器长度 Nb=184。滤波器损耗函数和滤波器输出yw(nT)分别如图5-13(a)和（b）所示。

用等波纹最佳逼近法设计滤波器，滤波器长度 Ne=83。滤波器损耗函数和滤波器输出ye(nT)分别如图5-14(c)和（d）所示。

两种方法设计的滤波器都能有效地从噪声中提取信号，但等波纹最佳逼近法设计的滤波器阶数低得多，当然滤波实现的运算量以及时延也小得多，从图5-13(b)和5-14（d）可以直观地看出时延差别。

23

天津大学仁爱学院2013届本科生毕业生设计（论文）

第六章 总结与展望

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

Matlab语言简洁紧凑，使用方便，库函数十分丰富。MATLAB程序书写的形式自由，利用丰富的库函数避开了繁琐的子程序编程任务。

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

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

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

本设计实现了基于MATLAB的数字低通滤波器。

低通滤波器是容许低频信号通过、但减弱（或减少）频率低于截止频率信号通过的滤波器。对于不同滤波器而言,每个频率的信号的减弱程度不同。它有时被称为低频剪切滤波器；在音频应用中也使用低音消除滤波器或者噪声滤波器。

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

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

24

天津大学仁爱学院2013届本科生毕业生设计（论文）

参考文献

[1]徐明远，刘增力，《MATLAB仿真在信号处理中的应用》[M].西安：西安电子科技大学出版社，2007.11.

[2]陈桂明，张明照，《应用MATLAB语言处理信号与数字图像》[M].北京：科学出版社，2000.5. [3]飞思科技产品研发中心.MATLAB基础与提高[M].北京：电子工业出版社，2005. [4]高西全，丁玉美.数字信号处理[M].西安：西安电子科技大学出版社，第三版，2008. [5]李亚奇，张雅琦.线性相位FIR数字滤波器[J].电子测量技术，2005（6）：35-37. [6]郭德才.基于Matlab的FIR低通滤波器的设计与仿真[J].通化师范学院学报，2009, 30（8）：38-41.

[7]赵刚.基于数字滤波器设计的讨论[J].南开大学学报（自然科学版），2003（3）：101-103. [8]陈明军.改进窗函数在FIR数字滤波器设计中的应用[J].继电器，2007（13）：65-67. [9]丁玉美, 高西全.数字信号处理 [M]. 第三版. 西安: 西安电子科技大学出版社, 2008, 6.

[10]闫胜利. FIR滤波器及设计原理[J]. 长春工程学院学报(自然科学版), 2003, 6, 4(1):

21-24.

[11]姚齐国. 基于MATLAB的数字滤波器的设计[J]. 江西理工大学学报, 2006, 2, 27(1):

50-52.

[12]杨守卫. FIR数字滤波器应用分析探讨[J]. 河北省工程咨询院学报, 2011, 7, 29(15):

47-49.

[13]朱敏. MATLAB数字信号处理工具箱的开发和应用[J]. 信息与电脑, 2010, 2, 26(8):

154-155.

[14]姚海燕. FIR数字滤波器设计窗函数法与频率抽样法比较[J].安阳工学院学报,2007,6,

12(6): 51-53.

[15]刘波. MATLAB信号处理[M]. 北京: 电子工业出版社, 2006, 7.

25

天津大学仁爱学院2013届本科生毕业生设计（论文）

附录

1)用升余弦窗设计一线性相位低通FIR数字滤波器程序 N=15; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=hanning(N); hn=hdn.*(wn');

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(1);subplot(2,2,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('N=15时设计汉宁窗h(n)'); subplot(2,2,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(2,2,3); plot(w,angle(H));

xlabel('w');ylabel('φ(w)');title('h(n)的相位谱');axis([0,3,-4,4]); subplot(2,2,4); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,-110,0]); N=33; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=hanning(N); hn=hdn.*(wn');

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(2);subplot(2,2,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('N=33时设计汉宁窗h(n)');

26

天津大学仁爱学院2013届本科生毕业生设计（论文）

subplot(2,2,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(2,2,3); plot(w,angle(H));

xlabel('w');ylabel('φ(w)');title('h(n)的相位谱');axis([0,3,-4,4]); subplot(2,2,4); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,-110,0]);

?4种窗函数设计线性相位低通滤波器程序 (2)N=33，Wc= rad，用4N=33; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=boxcar(N); hn=hdn.*(wn');

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(1);subplot(1,3,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('N=33时设计矩形窗h(n)'); subplot(1,3,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(1,3,3); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,-110,0]); N=33; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=hanning(N); hn=hdn.*(wn');

27

天津大学仁爱学院2013届本科生毕业生设计（论文）

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(2);subplot(1,3,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('汉宁窗函数设计h(n)'); subplot(1,3,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(1,3,3); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,-110,0]); N=33; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=hamming(N); hn=hdn.*(wn');

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(3);subplot(1,3,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('哈明窗函数设计h(n)'); subplot(1,3,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(1,3,3); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,60,0]); N=33; wc=pi/4; a=(N-1)/2; n=0:(N-1); m=n-a+eps;

hdn=sin(wc*m)./(pi*m); wn=blackman(N);

28

天津大学仁爱学院2013届本科生毕业生设计（论文）

hn=hdn.*(wn');

[H,w]=freqz(hn,[1],1024,'whole'); dbH=20*log10((abs(H)+eps)/max(abs(H))); figure(4);subplot(1,3,1); stem(n,hn,'.');

xlabel('n');ylabel('h(n)');title('布莱克曼函数设计h(n)'); subplot(1,3,2); plot(w,abs(H));

xlabel('w');ylabel('H(jw)');title('h(n)的幅度谱');axis([0,3,0,1.5]); subplot(1,3,3); plot(w/pi,dbH);

xlabel('w/pi');ylabel('dB');title('损耗特性');axis([0,1,-100,0]); 5.2 利用频率抽样法设计低通滤波器程序

M=63; Wp=0.5*pi;%所需频率采样点个数及通带截止频率

m=0:(M+1)/2; Wm=2*pi*m./(M+1);%通频带上的采样点及阻带截止频率 mtr=floor(Wp*(M+1)/(2*pi))+2;%向负方向入floor(3.5)=3;floor(-3.2)=-4 Ad=[Wm<=Wp];Ad(mtr)=0.38;

Hd=Ad.*exp(-j*0.5*M*Wm);%构造频域采样向量H(k) Hd=[Hd conj(fliplr(Hd(2:(M+1)/2)))];

%fliplr函数实现矩阵的左右翻转conj是求复数的共轭 h=real(ifft(Hd));%h（n）=IDFT[H(k)]

w=linspace(0,pi,1000);%用于产生0,pi之间的1000点行矢量 H=freqz(h,[1],w);%滤波器的幅频特性图 figure(1)

plot(w/pi,20*log10(abs(H)));%参数分别是归一化频率与幅值

xlabel('归一化角频率');ylabel('增益/分贝');title('滤波器的增益响应'); axis([0 1 -50 0.5]);

f1=100;f2=300;f3=700;fs=2000;%待滤波正弦信号频率及采样频率 figure(2) subplot(211)

t=0:1/fs:0.25;%定义时间范围和步长

29

天津大学仁爱学院2013届本科生毕业生设计（论文）

s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%滤波前信号 plot(t,s);%滤波前的信号图像

xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图'); subplot(212)

Fs=fft(s,512); AFs=abs(Fs);%将信号变换到频域及信号频域图的幅值 f=(0:255)*fs/512;%频率采样

plot(f,AFs(1:256));%滤波前的信号频域图

xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图'); figure(3)

sf=filter(h,1,s);%使用filter函数对信号进行滤波 subplot(211)

plot(t,sf)%滤波后的信号图像

xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图'); axis([0.2 0.25 -2 2]);%限定图像坐标范围 subplot(212)

Fsf=fft(sf,512); AFsf=abs(Fsf);%滤波后的信号频域图及信号频域图的幅值 f=(0:255)*fs/512;%频率采样

plot(f,AFsf(1:256))%滤波后的信号频域图

xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图'); 5.3（1）程序 function xt=xtg(N)

%实验五信号x(t)产生,并显示信号的幅频特性曲线

%xt=xtg(N) 产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率Fs=1000Hz

%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz. Fs=1000;T=1/Fs;Tp=N*T; t=0:T:(N-1)*T;

fc=Fs/10;f0=fc/10; %载波频率fc=Fs/10，单频调制信号频率为f0=Fc/10; mt=cos(2*pi*f0*t); %产生单频正弦波调制信号mt，频率为f0 ct=cos(2*pi*fc*t); %产生载波正弦波信号ct，频率为fc xt=mt.*ct; %相乘产生单频调制信号xt nt=2*rand(1,N)-1; %产生随机噪声nt

30

天津大学仁爱学院2013届本科生毕业生设计（论文）

%=======设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声=======

fp=150; fs=200;Rp=0.1;As=70; fb=[fp,fs];m=[0,1];

% 滤波器指标

% 计算remezord函数所需参数f,m,dev

dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];

[n,fo,mo,W]=remezord(fb,m,dev,Fs); % 确定remez函数所需参数

hn=remez(n,fo,mo,W); % 调用remez函数进行设计,用于滤除噪声nt中的低频成分

yt=filter(hn,1,10*nt); %滤除随机噪声中低频成分，生成高通噪声yt %================================================================ xt=xt+yt; %噪声加信号 fst=fft(xt,N);k=0:N-1;f=k/Tp;

subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)'); axis([0,Tp/5,min(xt),max(xt)]);title('(a) 信号加噪声波形')

subplot(3,1,2);plot(f,abs(fst)/max(abs(fst)));grid;title('(b) 信号加噪声的频谱')

axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度') 5.2（2）程序

%《数字信号处理(第三版）学习指导》第10章实验5程序exp5.m % FIR数字滤波器设计及软件实现 clear all;close all;

%==调用xtg产生信号xt, xt长度N=1000,并显示xt及其频谱,========= N=1000;xt=xtg(N);

fp=120; fs=150;Rp=0.2;As=60;Fs=1000; % 输入给定指标 % (1) 用窗函数法设计滤波器

wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于pi归一化） B=2*pi*(fs-fp)/Fs; %过渡带宽度指标 Nb=ceil(11*pi/B); %blackman窗的长度N hn=fir1(Nb-1,wc,blackman(Nb)); Hw=abs(fft(hn,1024));

% 求设计的滤波器频率特性

ywt=fftfilt(hn,xt,N); %调用函数fftfilt对xt滤波

%以下为用窗函数法设计法的绘图部分（滤波器损耗函数，滤波器输出信号波形） %省略

% (2) 用等波纹最佳逼近法设计滤波器 fb=[fp,fs];m=[1,0];

% 确定remezord函数所需参数f,m,dev

31

天津大学仁爱学院2013届本科生毕业生设计（论文）

dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; [Ne,fo,mo,W]=remezord(fb,m,dev,Fs); hn=remez(Ne,fo,mo,W); Hw=abs(fft(hn,1024));

% 确定remez函数所需参数

% 调用remez函数进行设计 % 求设计的滤波器频率特性

yet=fftfilt(hn,xt,N); % 调用函数fftfilt对xt滤波

%以下为用等波纹设计法的绘图部分（滤波器损耗函数，滤波器输出信号yw(nT)波形） %省略

32

外文资料

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 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

1

(3) Optimal filter design methods The Window Method

In this method, [Park87], [Rab75], [Proakis00] from the desired 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

2

(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)

=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)

3

= 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.54 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 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.

4

(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.

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

5

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

(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.

6

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.

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

7

response.

Nonlinear Equation solution for maximal ripple FIR filters

The real part of the frequency response of the designed FIR filter can 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 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

8

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