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

1、重庆三峡学院毕业设计(论文 题 目 基于窗函数法 FIR 数字滤波器的设计专 业 电子信息工程年 级 2 0 0 6 级学生姓名 郝 海学生学号 200615190154指导教师 徐正坤 职称 助教完成毕业设计(论文时间 2009年 12月基于窗函数法 FIR 数字滤波器的设计郝海重庆三峡学院应用技术学院电子信息工程(仪器仪表专业 06级 重庆万州 404100 摘要 简述了数字滤波器中的有限长单位冲激响应 (FIR滤波器的原理 , 对 FIR 滤波器的窗函数设 计方法进行了研究。 窗函数法在 FIR 数字滤波器的设计中有着广泛的应用。 介绍了一种基于等波 纹切比雪夫逼近准则的 FIR 数字滤

2、波器的设计方法 , 通过 MATLAB 的仿真实现。 传统的数字滤波器 设计方法繁琐且结果不直观 , 本文利用 MATLAB 具有强大的科学计算和图形显示这一优点 , 与窗函 数法设计理论相结合共同设计 FIR 数字滤波器 , 不但使设计结果更加直观 , 而且提高了滤波器的设 计精度 , 从而更好地达到预期效果。关键字 FIR 数字滤波器 窗函数 等波纹切比雪夫逼近 MATLAB仿真目 录1 引言 - 12 FIR数字滤波器的介绍 - 22.1FIR 数字滤波器的特点 - 2 2.2线性相位 FIR 数字滤波器的特点 - 2 2.2.1 单位冲激响应 h(n的特点 - 2 2.2.2 线性相位

3、的条件 - 2 2.2.3 线性相位特点和幅度函数的特点 - 2 2.3FIR 数字滤波器的设计原理 - 4 2.4数字滤波器的性能指标 - 53窗函数设计法 - 6 3.1窗函数设计原理分析 - 6 3.2设计方法 - 7 3.3窗函数介绍 - 9 3.4窗函数法设计步骤 - 13 3.5设计实例 - 133.6窗函数法计算中的主要问题 - 144 FIR数字滤波器的优化 - 154.1均方误差最小化准则 - 15 4.2切比雪夫最佳一致逼近定理 - 16 4.3利用切比雪夫逼近理论设计 FIR 数字滤波器 - 16 4.4瑞米兹算法 - 175 MATLAB简介与数字滤波器的 MATLAB

4、 实现 - 195.1MATLAB 简介 - 19 5.2MATLAB 程序 - 20结论 - 26 谢辞 - 26 参考文献 - 26 附录 - 271 引言数字信号处理 (DSP,digital signal processing 是从 20世纪 60年代以来,随着信息学科和 计算机的高速发展而迅速发展起来的一门新兴学科。数字信号处理是把信号用数字或符号表示的 序列, 通过计算机或通用 (专用 信号处理设备, 用数字的数值计算方法处理 (例如滤波、 变换、 压 缩、增强、估计、识别等 ,以达到提取有用信息便于应用的目的。数字滤波是数字信号处理的一部分。 数字滤波器按照单位取样响应 h(n的

5、时域特性可以分为 无限脉冲响应 (IIR系统和有限脉冲响应 (FIR系统。 FIR 数字滤波器的优点在于它可以做成具有 严格线性相位 , 而同时可以具有任意的幅度特性 ; 它的传递函数没有极点 ; 这保证了设计出的 FIR 数字滤波器一定是平稳的。所谓数字滤波器设计,简单地说,就是要找到一组能满足特定滤波要求的系数向量 a 和 b 。 而滤波器设计完成后还需要进一步考虑如何将其实现,即选择什么样的滤波器结构来完成滤波运 算。 FIR 数字滤波器的设计方法很多,其中较为常用的是窗函数设计法、频率采样设计法和最优 化设计法。 本文讨论利用窗函数法、 均方误差最小化法和等波纹切比雪夫逼近法 (调用

6、remez 函数 来分别实现各种 FIR 滤波器的设计。窗函数法设计的基本思想是把给定的频率响应通过 IDTFT(Inverse Discrete Time Fourier Transform ,求得脉冲响应,然后利用加窗函数对它进行截断和平滑,以实现一个物理可实现且 具有线性相位的 FIR 数字滤波器的设计目的。 其核心是从给定的频率特性, 通过加窗确定有限长 单位取样响应 h(n。均方误差最小化法是等波纹切比雪夫逼近法的基础,这一准则是使误差能量最小,但是由于 吉布斯效应,窗谱的肩峰过大,造成所设计出的滤波器通带起伏不均匀且过大,而阻带衰减则过 小,不能满足要求。等波纹切比雪夫逼近法,一致

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

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

9、处理 (时间序列分拆 等。2 FIR数字滤波器的介绍2.1 FIR数字滤波器的特点数字信号处理主要是研究用数字或符号的序列来表示信号波形,并用数字的方式去处理这些 序列,把它们改变成在某种意义上更为希望的形式,以便估计信号的特征参量,或削弱信号中的 多余分量和增强信号中的有用分量。有限长单位冲激响应 (FIR数字滤波器可以做成具有严格的线性相位, 同时又可以具有任意的 幅度特性。此外, FIR 滤波器的单位抽样响应是有限长的,因而滤波器一定是稳定的。再有,只 要经过一定的延时,任何非因果有限长序列都能变成因果的有限长序列,因而总能用因果系统来 实现。最后, FIR 滤波器由于单位冲激响应是有限

10、长的,可以用快速傅立叶变换 (FFT算法来实现 过滤信号,从而可大大提高运算效率。但是,要取得很好的衰减特性, FIR 滤波器 H(z的阶次比 IIR 滤波器的要高。2.2线性相位 FIR 数字滤波器的特点2.2.1 单位冲激响应 h(n的特点FIR 滤波器的单位冲激响应 h(n是有限长 (0 n N-1, 其 Z 变换为:-= -= 1 0( N mm z n hzH在有限 Z 平面有 (N-1个零点 , 而它的 (N-1个极点均位于原点 z=0处。 2.2.2 线性相位的条件如果 FIR 滤波器的单位抽样响应 h(n为实数而且满足以下任一条件 : 偶对称:h(n=h(N-1-n奇对称:h(

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

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

13、数、偶对称时的相位函数、冲激响应、幅度函数波形图 图 2-2 长度 N 为奇数、偶对称时的相位函数、冲激响应、幅度函数波形图 图 2-3 长度 N 为偶数、奇对称时的相位函数、冲激响应、幅度函数波形图 图 2-4 长度 N 为奇数、奇对称时的相位函数、冲激响应、幅度函数波形图四种线性相位 FIR 滤波器的特性可以总结如下: 第一种情况,偶对称、奇数点,四种滤波器都可设计;第二种情况,偶对称、偶数点,可设计低、带通滤波器,不能设计高通和带阻; 第三种情况,奇对称、奇数点,只能设计带通滤波器,其它滤波器都不能设计; 第四种情况,奇对称、偶数点,可设计高、带通滤波器,不能设计低通和带阻。2.3 FI

14、R数字滤波器的设计原理一个截止频率为 c (rad/s的理想数字低通滤波器,其传递函数的表达式是:=-c c j j d e e H , 0, ( (式 2.3.1 由式 2.3.1可以看出,这个滤波器在物理上是不可实现的,因为冲激响应具有无限性和因果 性。为了产生有限长度的冲激响应函数,我们取样响应为 (n h ,长度为 N ,其系数函数为 (z H :-=-=1( (N n nzn h z H (式 2.3.2用 (n h 表示截取 (n h d 后冲激响应,即 ( ( (n h n n h d =,式子中 (n 为窗函数,长度为 N 。当 =(N-1/2时,截取的一段 (n h 对 (N

15、-1/2对称,可保证所设计的滤波器具有线性相位。一般来说, FIR 数字滤波器输出 (n y 的 Z 变换形式 Y(z与输入 (n x 的 Z 变换形式之间的关 系如下:( ( 1( 0( ( ( (1z X zn h zh h z X z H z Y n-+= (式 2.3.3从上面的 Z 变换和结构图可以很容易得出 FIR 滤波器的差分方程表示形式。 对式 2.3.3进行 反 Z 变换,可得:1( ( 1( 2( ( 1( (x n h n x h n x h n y +-+= (式 2.3.4 图 2-5卷积型滤波器式 (2.3.4为 FIR 数字滤波器的时域表示方法, 其中 (n x

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

17、往以频率 响应的幅度特性的允许误差来表征。以低通滤波器特性为例,频率响应有通带、过渡带及阻带三 个范围。在通带内: 1- AP (j e H 1 c c 在阻带中: (j e H st A st c 其中 c 为通带截止频率 , st 为阻带截止频率, Ap 为通带误差 , st A 为阻带误差。 图 2-6 低通滤波器的幅度特性与模拟滤波器类似,数字滤波器按频率特性划分为低通、高通、带通、带阻、全通等类型, 由于数字滤波器的频率响应是周期性的,周期为 2。由于频率响应的周期性, 频率变量以数字频率 来表示, 所以数字滤波器设计中必须给出抽 样频率。图 2-7为各种数字滤波器理想幅度,可以看出

18、:1、 一个高通滤波器相当于一个全通滤波器减去一个低通滤波器。 2、 一个带通滤波器相当于两个低通滤波器相减。3、 一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器。这里的相加相减都是相当于并联结构。图 2-7中所示的各种数字滤波器理想频率响应只表示了正频率部分,这样的理想频率响应是 不可能实现的,原因是频带之间幅度响应是突变的,因而其单位抽样响应是非因果的。因此要给 出实际逼近容限。数字滤波器的系统函数 (z H ,它在 z 平面单位圆上的值为滤波器频率响应(j eH ,表征数字滤波器频率响应特征的三个参量是幅度平方响应、相位响应和群延时响应。 图 2-7 各种理想数字滤波器的幅度频率响

19、应 .3窗函数设计法3.1窗函数设计原理分析设数字滤波器的传输函数为 (j e H , (n h d 是与其对应的单位脉冲响应 , (z H 为系统函数。-=-=1( (N n nj j en h eH (式 3.1.1d eeHn h nj j dd (21(-=(式 3.1.2低通高通带通带阻全通(j d eH(j d eH(j d eH(j d eH(j d eH-=-=1( (N n nzn h z H (式 3.1.3一般说来 , (n h d 是无限长的 , 需要求对 (j d e H 的一个逼近。采用窗函数设计法时 , 可通过 对理想滤波器的单位采样响应加窗设计滤波器( ( (n

20、 h n n h d = (式 3.1.4其中 , (n 是一个长度有限的窗 , 在区间 0 n N外值为 0 ,且关于中间点对称1( (n N n -= (式 3.1.5频率响应根据 (式 3.1.5 ,由卷积定理得出( (21 (j j d j eeH eH =(式 3.1.6理想的频率响应被窗函数的离散时间傅立叶变换 (j e “平滑”了。采用窗函数设计法设计出来的滤波器的频率响应对理想响应 (j d e H 的逼近程度 , 由两个因 素 决 定 : (j e 主 瓣 的 宽 度 ; (j e 旁 瓣 的 幅 度 大 小 。理想的情况是 (j e 主瓣的宽度窄 , 旁瓣的幅度小。但对于一

21、个长度固定的窗函数来说 , 这些不能 独立地达到最小。窗函数的一些通用性质为 :1、窗函数的长度 N 增加 , 主瓣的宽度减小 , 使得过渡带变小。关系为 :NB = C其中 :B是过渡带的宽度 ;C 是取决于窗函数的一个参数。如矩形窗为 4。调整 N 可以有效地控制过渡带的宽度 , 但 N 的改变不改变主瓣和旁瓣的相对比例。 随着 N 值增加 , 过渡带变窄 , 波动频率也随着增加 , 虽然总的幅度有所减少 , 但截止频率附近的肩峰并不减少 , 而只是随着 N 值的 增加 , 肩峰被抑制在愈来愈小的范围内 , 使肩峰宽度变窄。2、 窗函数的旁瓣的幅度大小取决于窗函数的选择。 选择恰当的窗函数

22、使主瓣包含更多的能量 , 相应旁瓣的幅度就减小。旁瓣幅度的减小 , 可以减少通带和阻带的波动 , 使通带尽可能趋近水平 , 阻 带尽可能达到最大衰减。但通常此时过渡带会变宽。3、取不同的窗函数对幅度特性的整形效果比单纯的增加窗口长度要强得多。 3.2设计方法这种方法也叫傅里叶级数法。 一般是先给出所要求的理想的滤波器的频率响应 (j d eH , 要求设计一个 FIR 滤波器频率响应 -=-=1( (N n nj j en h eH 来逼近 (j d eH 。 设计是在时域进行的,因而先由 (j d e H 的傅里叶反变换导出 (n h d ,即d eeH n h nj j d d (21 (

23、-=(式 3.2.1由于 (j d eH 是矩形频率响应特性,故 (n h d 一定是无限长序列,且是非因果的,而 FIR 滤波器的 (n h 必然是有限长的,所以要用有限长的 (n h 来逼近无限长的 (n h d ,最有效的方法是截断(n h d 或者说用一个有限长度的窗口函数序列 (n 来截取 (n h d ,即( ( (n h n n h d = (式 3.2.2因而窗函数序列的形状及长度的选择就是关键。我们以一个截止频率为 c 的线性相位的理想矩形幅度特性的低通滤波器为例来讨论。 设低通 特性的群延时为 ,即-=-cc cc j j d e eH , , 0, ( (式 3.2.3

24、这表明,在通带 c 范围内, (j d e H 的幅度是均匀的,其值为 1,相位是 -。利用 (1式可得( (sin 21 (-=-n n d een h c c c nj j d cc(式 3.2.4(n h d 是中心点在 的偶对称无限长非因果序列, 要得到有限长的 (n h , 一种最简单的方法就是取矩形窗 (n R N ,即( (n R n N =但是按照线形相位滤波器的约束, (n h 必须是偶对称的,对称中心应为长度的一半 (N-1/2,因 而必须 =(N-1/2,所以有=-=21-N , 010, ( ( ( (为其他 n N n n h n n h n h d d (式 3.2

25、.5 将 (式 3.2.4 代入 (式 3.25 ,可得10, , 021( 21(sin (-=N n n N n N n n h c cc为其他值 (式 3.2.6 此时,一定满足 1( (n N h n h -=这一线性相位的条件。下面求 (n h 的傅里叶变换, 也就是找出待求 FIR 滤波器的频率特性, 以便能看出加窗处理后 究竟对频率响应有何影响。按照复卷积公式,在时域是相乘、频域上是周期性卷积关系,即d eeH eH j j d j ( (21 (-=(式 3.2.7因而 (j e H 逼近 (j d e H 的好坏,完全取决于窗函数的频率特性 (j e W 。 窗函数 (n 的

26、频率特性 (j e W 为-=-=1( (N n nj j en eW (式 3.2.8对矩形窗 (n R N ,则有2sin(sin( (211N N eeeW N j N n nj j R -=-= (式 3.2.9 也可表示成幅度函数与相位函数 21( (-=N j R j N eW eW (式 3.2.10其中2sin(sin(N N W R =(式 3.2.11(j R eW 就是频域抽样内插函数, 其幅度函数 (R W 在 N /2±=之内为一个主瓣, 两侧形成许多衰减振荡的旁瓣,如果将理想频率响应也写成 21( (-=N j d j d eH e H (式 3.2.12

27、则其幅度函数为<=c cd H , 0, 1 ( (式 3.2.133.3窗函数介绍实际应用的窗函数,可分为以下主要类型:1、幂窗 -采用时间变量某种幂次的函数,如矩形、三角形、梯形或其它时间 (t的高次幂; 2、三角函数窗 -应用三角函数,即正弦或余弦函数等组合成复合函数,例如汉宁窗、海明窗等; 3、指数窗 -采用指数时间函数,如 ste-形式,例如高斯窗等。下面介绍几种常用窗函数的性质和特点。(1矩形窗矩形窗属于时间变量的零次幂窗,函数形式为 :=T t Tt T t , 0, 1( (式 3.3.1相应的窗谱为:TTW sin 2 (=(式 3.3.2矩形窗使用最多,习惯上不加窗就

28、是使信号通过了矩形窗。这种窗的优点是主瓣比较集中, 缺点是旁瓣较高,并有负旁瓣,导致变换中带进了高频干扰和泄漏,甚至出现负谱现象。 图 3-1矩形窗的时域及频域波形(2三角窗亦称费杰 (Fejer窗,是幂窗的一次方形式,其函数形式是:-=T t Tt TtT t , 0, 1(1 ( (式 3.3.3三角窗与矩形窗比较,主瓣宽约等于矩形窗的两倍,但旁瓣小,而且无负旁瓣,如图 3-2所 示。 图 3-2三角窗的时域及频域波形(3汉宁 (Hanning窗汉宁窗又称升余弦窗,其时域表达式为:+=T t T t T t T t , 0, cos 2121(1 ( (式 3.3.4相应的窗谱为:-+=T

29、 T T T TTW sin( sin(21sin ( (式 3.3.5由此式可以看出,汉宁窗可以看作是 3个矩形时间窗的频谱之和,或者说是 3个 sin(t型函 数之和,而括号中的两项相对于第一个谱窗向左、右各移动了 /T,从而使旁瓣互相抵消,消去 高频干扰和漏能。可以看出,汉宁窗主瓣加宽并降低,旁瓣则显著减小,从减小泄漏观点出发,汉 宁窗优于矩形窗。但汉宁窗主瓣加宽,相当于分析带宽加宽,频率分辨率下降。 (4海明 (Hamming窗海明窗也是余弦窗的一种,又称改进的升余弦窗,其时间函数表达式为:>+=T t T t T t T t , 0, cos 4. 054. 0(1( (式 3

30、.3.6其窗谱为:-+=T T T T TTW sin( sin(46. 0sin 08. 1 ( (式 3.3.7海明窗与汉宁窗都是余弦窗, 只是加权系数不同。 海明窗加权的系数能使旁瓣达到更小。 分析 表明,海明窗的第一旁瓣衰减为 -42dB 。海明窗的频谱也是由 3个矩形窗的频谱合成,但其旁瓣衰 减速度为 20dB /(10oct,这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。(5高斯窗高斯窗是一种指数窗。其时域函数为:>=-T t T t eTt at , 0, 1 (2 (式 3.3.8式中 a 为常数, 决定了函数曲线衰减的快慢。 a 值如果选取适当, 可以使截断点

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

32、布莱克曼窗主瓣宽,旁瓣小,频率识别精度最低, 但幅值识别精度最高。 图 3-3几种常用的窗函数的时域和频域波形对于窗函数,还有一些要求:1 3dB 带宽 B ,它是主瓣归一化的幅度下降到 -3dB 时的带宽。当数据长度为 N 时,矩形窗主 瓣两个过零点之间的宽度为 4/N。2最大边瓣峰值 A(dB。3边瓣谱峰渐进衰减速度 D(dB/cot。所以,理想的窗函数应当具有最小的 B 和 A ,和最大的 D 。3.4窗函数法设计步骤1、首先是给定所要求的频率响应函数 (j d e H ; 2、其次,求单位冲激响应 d eeH n h nj j d d (21 (-=;3、 再次, 有过渡带宽及阻带最小

33、衰减的要求, 查表选定窗函数及 N 的大小, 一般 N 的大小要 通过几次试探而后确定;4、求得所设计的 FIR 滤波器的单位冲激响应;( ( (n h n n h d =, n=0, 1, ,N-1;5、求 -=-=1( (N n nj j en h e H ,检验是否满足设计要求,如不满足,则需要重新设计。3.5设计实例线性相位 FIR 低通滤波器的设计 (用窗函数法 。指标要求 :通带截止频率 :0.2, 阻带起始频 率 :0.4, 阻带最小衰减 :-50dB。(1设 (j e H 为理想线性相位滤波器=-其他, 0 (, c j j e e H 由所需低通滤波器的过渡带求出理想低通滤波

34、器的截止数字频率 =0.3, 得出 :。为线性相位所需的移位 , 21, , ,(sin2121 (-=-=-N n n n n d e d e e n h c c n j n j j d cc (2由阻带衰减确定窗函数 , 由过渡带宽确定 N 值。阻带最小衰减 50dB, 比对 6种窗函数基本 参数选定窗函数为海明窗。 所要求的过渡带宽 : =0.4-0.2=0.2N=6.6/0.2=33,=(N-1/2=16(3由海明窗函数确定 FIR 滤波器的 h(n。21(1(sin (12cos(46. 054. 0n -=-=N n N n n h n R N n c d N (得出:(12cos

35、(46. 054. 016(16(3. 0sin ( ( (n R N n n n n h n h N d -=(4仿真检验各项指标 , 得出结论 :满足设计要求。 取 N=33,偶对称 ,得 :过渡带宽 :0.3476563, 第一通带波纹 :0.020837dB, 第一阻带最小衰减 :60.9159dB。 图 3-4 例中设计的线性相位 FIR 低通滤波器幅度响应曲线 (海明窗 3.6窗函数法计算中的主要问题首先当 (j d eH 很复杂或不能按 d eeHn h nj j dd (21 (-=是直接计算积分时, 则必须用求和代替积分,以便在计算机上计算,也就是要计算离散傅里叶反变换,一般

36、都采用 FFT 来计 算。将积分限分成 M 段,也就是令抽样频率为k Mk 2=,空, 2, M-1则有nk MjM k kMjd M eeH Mn h 212 (1 (-=频域的抽样造成时域序列的周期延拓,延拓周期是 M ,即-=+=r dM rM n hn h ( (由于 (n h d 有可能是无限长的序列,因而严格的说,必须 M 时 (n h M 才能等于 (n h d 而不 产生混叠现象,即(lim (n h n h M M d =实际上,由于 (n h d 随 n 的增加衰减很快,一般只要 M 足够大,即 M>>N,就足够了。其次,窗函数设计法的另一个困难就是需要先确定窗

37、函数的形状和窗函数需要的点数 N ,以满足给定的 频率响应指标。这一困难可利用计算机采用累试法加以解决。一般在设计凯泽窗时,则零阶变形贝塞尔函数可采用无穷级数来表达=+=12020 2(! 11 2(! 1 (k k k k x k x k x I这个无穷级数可用有限项级数去近似,项数多少有要求的精度来确定。窗函数法的优点是简单,有闭合形式的公式可循,因而很实用。其缺点是通带、阻带的截止 频率不易控制。这就需要对 FIR 滤波器进行优化设计,下面的内容是 FIR 滤波器的优化。4 FIR数字滤波器的优化4.1均方误差最小化准则这一准则是使误差能量最小,若用 (j d e H 表示要求的频率相应

38、,用 (j e H 表示实际得到 的滤波器频率响应,以 (j e E 表示频率相应误差,即( ( (j j d j e H eH eE -= (式 4.1.1则均方误差为:d eH eH d eE e j j d j 222( (21 (21-=(式 4.1.2设计的目的就是选择一组 ( (1j e H F n h -=使得 2e 最小。现将 (1式中的 (j d e H 和 (j eH 分别用他们的冲激响应表示,即nj n dj d en heH -=( ( (式 4.1.3nj N n j en h eH -=( (1(式 4.1.4由于用 FIR 数字滤波器来逼近,故 (n h 长度是有

39、限长的。将他们代入 (式 4.1.1 式,得nj dN n nj dj j d j en hen h n heH eH eE -=-+-=-= ( ( ( ( ( (1其他(式 4.1.5按照帕赛瓦公式有:+-=-=-其他21222( ( (21n h n h n h d eE e d N n d j (式 4.1.6由此式看出,等式右边第二个求和式之取决于给定的特性 (n h d ,它和设计值 (n h 无关,故是一 个常数,要是 2e 最小,就必须使第一个求和式最小,即希望0 ( (=-n h n h d , 10-N n在这一条件下,就有min(22e e =也就是说,要满足-=nN n

40、 n h n h d 其他 , 010, ( ( (式 4.1.7 此式恰好是矩形窗的结果,所以,矩形窗设计结果一定满足最小均方误差准则。矩形窗虽然过渡 带最窄, 但是由于吉布斯效应, 窗谱的肩峰过大, 造成所设计出的滤波器通带起伏不均匀且过大, 而阻带衰减则过小,不能满足要求。4.2切比雪夫最佳一致逼近定理切比雪夫最佳一致逼近的基本思想是,对于给定的区间 a, b上的连续函数 f(x,在所有 n次多项式的集合 n 中, 寻找一多项式 ( p x , 使它在 a, b上对 f(x的偏差和其他一切属于 ( n 的多项式 p(x对 f(x的偏差相比是最小的,即m ax |( ( |m inmax

41、|( ( |a x ba x bp x f x p x f x -=- (式 4.2.1切比雪夫逼近理论指出,这样的多项式是存在的,而且是唯一的,并指出了构造这种最佳一 致逼近多项式的方法,这就是有名的“交错点组定理” :设 f(x是定义在 a, b上的连续函数, p(x为 n 中一个一阶次不超过 n 的多项式,并令( max |( ( |a x bE x p x f x =-及 E(x=p(x-f(x, p(x是 f(x最佳一致逼近多项式的充要条件是,E(x在 a, b上至少存在 n+2个交错点 b x x x a n << 21,使得( i n E x E =±1(

42、( i i E X E x +=- (式 4.2.2这 n+2个点即是“交错点组” ,显然 x 1,x 2x n+2, 是 E(x的极值点。n 阶切比雪夫多项式( cos(arccos n C x n x = (式 4.2.3在区间 -1, 1上存在 n+1点cos(k x k n=, k=0,1,.,n, (式 4.2.4轮流使得 ( n C x 取得最大值 +1和最小值 -1。 ( n C x 是 x 的多项式, 且最高项 nx 的系数是 12n -,可以证明,在所有 n 阶多项式中,多项式1( 2n n C x -和 0的偏差为最小。这样,如果我们在寻找 p(x时,能使误差函数为某一个

43、( n C x ,那么,这样的 p(x将是对 f(x的最佳一致逼近。 4.3利用切比雪夫逼近理论设计 FIR 数字滤波器基于交错定理,最优线性相位 FIR 滤波器的设计步骤如下:1、 输入部分, 包括滤波器性能要求以及滤波器的类型, 前者指的是所需的频率响应的幅度函 数 d H ( , 加权函数 (W 和滤波器单位抽样响应长度 N , 后者是要指出所需要设计的是带通 (包括低通、带通、高通、带阻等 滤波器或是微分器或者是希尔伯特变换器。2、用公式表示逼近问题,也就是表示加权逼近误差 (E 。 3、用瑞米兹多次交换算法,求逼近问题的解。 4、计算滤波器的单位抽样响应。需要注意的是,在整个设计程序

44、中,瑞米兹算法是以子程序形式出现的。 最优线性相位 FIR 数字滤波器完成的以上四个步骤,更详细的解释:1、输入数据,滤波器性能要求以及滤波器的类型。2、根据滤波器的类型和单位抽样响应长度,确定逼近函数 cos( 的个数 r 。 3、在 从 0到 的频率区间,用密集的格点来表示离散频率,两格点的距离可以写成*r=总 格 点 数格 点 密 度将格点处频率用有下标的频率来表示并赋以标称频率值。4、用子程序 EFF 和 WATE 分别将格点频率上要求的函数值 (d H 和加权函数值 (W 。5、用公式表达逼近问题,将 (dH 和 (W 分别变成 (dH和 (W 。 6、瑞米兹交换算法求解逼近问题,这

45、包括设定 (r+1个极值频率的初始假设值 (猜想值 。 7、计算滤波器单位冲激响应。8、把最佳误差和单位抽样响应打印出来。 4.4瑞米兹算法单独使用切比雪夫逼近理论,在实际的应用上存在着一定的困难。原因有两个:一是交错点 组事先是不知道的,这样当然无法求解交错点组的系数矩阵。要确定一组交错点组并非易事,即 使对于较小的 M 也是如此。二是直接求解交错点组的系数矩阵比较困难。为此,我们利用数值分 析中的 Remez 算法,靠一次次的迭带来求得一组交错点组,而且在每一次迭带中都避免直接求解 交错点组的系数矩阵。瑞米兹算法的步骤为:第一步:首先在频率子集上等间隔地取 M+2个频率 0, 1, , M

46、 1+, 作为交错点组的初始 猜测位置,然后按照下式进行计算M 1kd k k 0M 1kk k k 0H (1/W (+=+=- (式 4.4.1式中M 1kk i 0,i ki k 1(1cos( cos(+=- (式 4.4.2把 0, 1, , M 1+代入上式, 可求出 , 它是相对第一次指定的交错点组所产生的偏差, 实际上即是 2。求出 以后,利用重心形式的拉格朗日插值公式,可以在不求出 0, M 的 情况下,得到一个 g H ( ,即Mkkk 0k g Mkk 0kCcos cos H ( cos cos =(-(-(式 4.4.3把 g H ( 代入 j j j j g d E

47、 (e W (e |H (e H (e |=-式 , 可求得误差函数 (E 。如果在子集 F 上,对所有的频率 ,都有 |E( |,这说明, (E 是纹波的极值,初 始猜测的 0, 1, M 1+恰是交错点组。这时,设计工作即可以结束。当然,对第一次猜测 的位置,不会恰好如此。一般,在某些频率处,总有 > (E ,这说明,需要交换上次猜测的 交错点组中的某些点,得到一组新的交错点组。第二步:对上次确定的交错点组 0, 1, M 1+中的每一个点,都在其附近检查时候在 某一个频率处有 > (E ,如若有,再在该点附近找出局部极值点,用这一局部极值点代替原 来的点。待这 M+2个点都检

48、查过后,便得到一组新的交错点组 0, 1, M 1+,再利用上式 是求出 , g H ( , (E ,这样就完成一次迭代,也就是完成了一次交错点组的交换。通过交 换算法,使得这一次的交错点组中的每一个 i 都是由上一次的交错点组所产生的 (E ,局部极 值频率点,因此,用这次的交错点组求出的 将增大。第三步:利用和第二步相同的方法,把在各频率处使 > (E 的点作为新的局部极值点, 从而又得到一组新的交错点组。重复上述步骤,因为新的交错点组的选择都是作为每一次求出的(E 的局部极值点,因此,在迭代中,每次的 |都是递增的。 最后收敛到自己的上限。也即 g H ( 最佳的一致逼近 d H ( 的解。 因此, 若再迭代一次, 新的误差曲线 (E 的峰值将不会 大于 |,这时迭代可以结束。由最后的交错点组可按Mkkk 0k g Mkk 0kCcos co