FIR滤波器的设计.doc

基于窗函数法 FIR 数字滤波器的设计 陈敏敏 宁夏大学电子信息工程 专业 09 级 12009243656 摘要 简述了数字滤波器中的有限长单位冲激响应 FIR 滤波器的原理 对 FIR 滤波器的窗函数设 计方法进行了研究 窗函数法在 FIR 数字滤波器的设计中有着广泛的应用 通过 MATLAB 的仿真 实现 传统的数字滤波器设计方法繁琐且结果不直观 本文利用 MATLAB 具有强大的科学计算和图 形显示这一优点 与窗函数法设计理论相结合共同设计 FIR 数字滤波器 不但使设计结果更加直观 而且提高了滤波器的设计精度 从而更好地达到预期效果 关键字 FIR 数字滤波器 窗函数 MATLAB 仿真 目 录 1 引言 1 2 FIR 数字滤波器的介绍 2 2 1 FIR 数字滤波器的特点 2 2 2 线性相位 FIR 数字滤波器的特点 2 2 2 1 单位冲激响应 h n 的特点 2 2 2 2 线性相位的条件 2 2 2 3 线性相位特点和幅度函数的特点 2 2 3 FIR 数字滤波器的设计原理 4 2 4 数字滤波器的性能指标 5 3 窗函数设计法 6 3 1 窗函数设计原理分析 6 3 2 设计方法 7 3 3 窗函数介绍 9 3 4 窗函数法设计步骤 13 3 5 设计实例 13 3 6 窗函数法计算中的主要问题 14 4 MATLAB 简介与数字滤波器的 MATLAB 实现 19 4 1 MATLAB 简介 19 4 2 MATLAB 程序 20 结论 26 谢辞 26 参考文献 26 附录 27 1 引言 数字信号处理 DSP digital signal processing 是从 20 世纪 60 年代以来 随着信息 学科和计算机的高速发展而迅速发展起来的一门新兴学科 数字信号处理是把信号用数字 或符号表示的序列 通过计算机或通用 专用 信号处理设备 用数字的数值计算方法处理 例如滤波 变换 压缩 增强 估计 识别等 以达到提取有用信息便于应用的目的 数字滤波是数字信号处理的一部分 数字滤波器按照单位取样响应 h n 的时域特性可 以分为无限脉冲响应 IIR 系统和有限脉冲响应 FIR 系统 FIR 数字滤波器的优点在于它 可以做成具有严格线性相位 而同时可以具有任意的幅度特性 它的传递函数没有极点 这保 证了设计出的 FIR 数字滤波器一定是平稳的 所谓数字滤波器设计 简单地说 就是要找到一组能满足特定滤波要求的系数向量 a 和 b 而滤波器设计完成后还需要进一步考虑如何将其实现 即选择什么样的滤波器结构 来完成滤波运算 FIR 数字滤波器的设计方法很多 其中较为常用的是窗函数设计法 频 率采样设计法和最优化设计法 本文讨论利用窗函数法来实现各种 FIR 滤波器的设计 窗函数法设计的基本思想是把给定的频率响应通过 IDTFT Inverse Discrete Time Fourier Transform 求得脉冲响应 然后利用加窗函数对它进行截断和平滑 以实现一 个物理可实现且具有线性相位的 FIR 数字滤波器的设计目的 其核心是从给定的频率特性 通过加窗确定有限长单位取样响应 h n MATLAB 软件是由美国 Math works 公司推出的用于数值计算和图形处理的科学计算系 统环境 MATLAB 新的版本集中了日常数学处理中的各种功能 包括高效的数值计算 矩阵 运算 信号处理和图形生成等功能 在 MATLAB 环境下 用户可以集成地进行程序设计 数 值计算 图形绘制 输入输出 文件管理等各项操作 MATLAB 提供了一个人机交互的数 学系统环境 该系统的基本数据结构是矩阵 在生成矩陈对象时 不要求作明确的维数说 明 与利用 c 语言或 FORTRAN 语言作数值计算的程序设计相比 利用 MATLAB 可以节省大量 的编程时间 在工程技术界 MATLAB 被用来解决一些实际课题和数学模型问题 典型的应 用包括数值计算 算法预设计与验证 以及一些特殊的短阵计算应用 如自动控制理论 统计 数字信号处理 时间序列分拆 等 2 FIR 数字滤波器的介绍 2 1 FIR 数字滤波器的特点 数字信号处理主要是研究用数字或符号的序列来表示信号波形 并用数字的方式去处 理这些序列 把它们改变成在某种意义上更为希望的形式 以便估计信号的特征参量 或 削弱信号中的多余分量和增强信号中的有用分量 有限长单位冲激响应 FIR 数字滤波器可以做成具有严格的线性相位 同时又可以具有 任意的幅度特性 此外 FIR 滤波器的单位抽样响应是有限长的 因而滤波器一定是稳定 的 再有 只要经过一定的延时 任何非因果有限长序列都能变成因果的有限长序列 因 而总能用因果系统来实现 最后 FIR 滤波器由于单位冲激响应是有限长的 可以用快速 傅立叶变换 FFT 算法来实现过滤信号 从而可大大提高运算效率 但是 要取得很好的衰 减特性 FIR 滤波器 H z 的阶次比 IIR 滤波器的要高 2 2 线性相位 FIR 数字滤波器的特点 2 2 1 单位冲激响应 h n 的特点 FIR 滤波器的单位冲激响应 h n 是有限长 0 n N 1 其 Z 变换为 1 0 N m m znhzH 在有限 Z 平面有 N 1 个零点 而它的 N 1 个极点均位于原点 z 0 处 2 2 2 线性相位的条件 如果 FIR 滤波器的单位抽样响应 h n 为实数而且满足以下任一条件 偶对称 h n h N 1 n 奇对称 h n h N 1 n 其对称中心在 n N 1 2 处 则滤波器具有准确的线性相位 2 2 3 线性相位特点和幅度函数的特点 1 h n 偶对称 2 1 2 1 cos 1 0 N n N nhH k n 幅度函数 H 包括正负值 相位函数是严格线性相位 说明滤波器有 N 1 2 个抽样的 延时 它等于单位抽样响应 h n 长度的一半 图 2 1 中 线性相位无 90 附加相移 幅度函 数在 处存在零点 且对 呈奇对称 因此不适合作高通滤波器 图 2 2 所示 线性相 位无 90 附加相移 幅度函数对在 0 2 呈偶对称 因此适合作低通 高通滤波器 2 h n 奇对称 2 2 1 2 1 sin 1 0 N n N nhH 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 为偶数 偶对称时的相位函数 冲激响应 幅度函数波形图 图 2 2 长度 N 为奇数 偶对称时的相位函数 冲激响应 幅度函数波形图 图 2 3 长度 N 为偶数 奇对称时的相位函数 冲激响应 幅度函数波形图 图 2 4 长度 N 为奇数 奇对称时的相位函数 冲激响应 幅度函数波形图 四种线性相位 FIR 滤波器的特性可以总结如下 第一种情况 偶对称 奇数点 四种滤波器都可设计 第二种情况 偶对称 偶数点 可设计低 带通滤波器 不能设计高通和带阻 第三种情况 奇对称 奇数点 只能设计带通滤波器 其它滤波器都不能设计 第四种情况 奇对称 偶数点 可设计高 带通滤波器 不能设计低通和带阻 2 3 FIR 数字滤波器的设计原理 一个截止频率为 rad s 的理想数字低通滤波器 其传递函数的表达式是 c 式 2 3 1 c c j j d e eH 0 由式 2 3 1 可以看出 这个滤波器在物理上是不可实现的 因为冲激响应具有无限性 和因果性 为了产生有限长度的冲激响应函数 我们取样响应为 长度为 N 其系数 nh 函数为 zH 式 2 3 2 1 0 N n n znhzH 用表示截取后冲激响应 即 式子中为窗函数 nh nhd nhnnh d n 长度为 N 当 N 1 2 时 截取的一段对 N 1 2 对称 可保证所设计的滤波器具 nh 有线性相位 一般来说 FIR 数字滤波器输出的 Z 变换形式 Y z 与输入的 Z 变换形式之 ny nx 间的关系如下 式 2 3 3 1 0 1 zXznhzhhzXzHzY n 从上面的 Z 变换和结构图可以很容易得出 FIR 滤波器的差分方程表示形式 对式 2 3 3 进行反 Z 变换 可得 式 2 3 4 1 1 2 1 xnhnxhnxhny 图 2 5 卷积型滤波器 式 2 3 4 为 FIR 数字滤波器的时域表示方法 其中是在时间 n 的滤波器的输入 nx 抽样值 根据式 2 3 4 即可对滤波器进行设计 从上面的公式我们可以看出 在对滤波器 实际设计时 整个过程的运算量很大 设计完成后对已设计的滤波器的频率响应进行校核 运算量也很大 并且在数字滤波器设计的过程中 要根据设计要求和滤波效果不断地调整 以达到设计的最优化 在这种情况下 要进行大量复杂的运算 单纯靠公式计算和编制简 单的程序很难在短时间内完成 而利用 MATLAB 工具进行计算机辅助设计 则可以快速有效 地设计数字滤波器 大大的减少了计算量 2 4 数字滤波器的性能指标 我们在进行滤波器设计时 需要确定其性能指标 一般来说 滤波器的性能要求往往 以频率响应的幅度特性的允许误差来表征 以低通滤波器特性为例 频率响应有通带 过 渡带及阻带三个范围 在通带内 1 AP 1 j eH c c 在阻带中 j eH st A st c 其中为通带截止频率 为阻带截止频率 Ap 为通带误差 为阻带误差 c st st A 图 2 6 低通滤波器的幅度特性 与模拟滤波器类似 数字滤波器按频率特性划分为低通 高通 带通 带阻 全通等 类型 由于数字滤波器的频率响应是周期性的 周期为 2 由于频率响应的周期性 频率变量以数字频率来表示 所以数字滤波器设计中必须 给出抽样频率 图 2 7 为各种数字滤波器理想幅度 可以看出 1 一个高通滤波器相当于一个全通滤波器减去一个低通滤波器 2 一个带通滤波器相当于两个低通滤波器相减 3 一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器 这里的相加相减都是相当于并联结构 图 2 7 中所示的各种数字滤波器理想频率响应只表示了正频率部分 这样的理想频率 响应是不可能实现的 原因是频带之间幅度响应是突变的 因而其单位抽样响应是非因果 的 因此要给出实际逼近容限 数字滤波器的系统函数 它在 z 平面单位圆上的值 zH 为滤波器频率响应 表征数字滤波器频率响应特征的三个参量是幅度平方响应 j eH 相位响应和群延时响应 图 2 7 各种理想数字滤波器的幅度频率响应 3 窗函数设计法 3 1 窗函数设计原理分析 设数字滤波器的传输函数为 是与其对应的单位脉冲响应 为系 j eH nhd zH 低通 高通 带通 带阻 全通 j d eH j d eH j d eH j d eH j d eH 统函数 式3 1 1 1 0 N n njj enheH 式3 1 2 deeHnh njj dd 2 1 式3 1 3 1 0 N n n znhzH 一般说来 是无限长的 需要求对的一个逼近 采用窗函数设计法时 nhd j d eH 可通过对理想滤波器的单位采样响应加窗设计滤波器 式3 1 4 nhnnh d 其中 是一个长度有限的窗 在区间 0 n N 外值为 0 且关于中间点对称 n 式3 1 5 1 nNn 频率响应根据 式 3 1 5 由卷积定理得出 式3 1 6 2 1 jj d j eeHeH 理想的频率响应被窗函数的离散时间傅立叶变换 平滑 了 j e 采用窗函数设计法设计出来的滤波器的频率响应对理想响应的逼近程度 由 j d eH 两个因素决定 主瓣的宽度 旁瓣的幅度大小 j e j e 理想的情况是主瓣的宽度窄 旁瓣的幅度小 但对于一个长度固定的窗函数来说 这 j e 些不能独立地达到最小 窗函数的一些通用性质为 1 窗函数的长度N增加 主瓣的宽度减小 使得过渡带变小 关系为 NB C其中 B是过渡带的宽度 C是取决于窗函数的一个参数 如矩形窗为 4 调整N可以有效地控制过渡带的宽度 但N的改变不改变主瓣和旁瓣的相对比例 随着N 值增加 过渡带变窄 波动频率也随着增加 虽然总的幅度有所减少 但截止频率附近的肩峰 并不减少 而只是随着N值的增加 肩峰被抑制在愈来愈小的范围内 使肩峰宽度变窄 2 窗函数的旁瓣的幅度大小取决于窗函数的选择 选择恰当的窗函数使主瓣包含更多 的能量 相应旁瓣的幅度就减小 旁瓣幅度的减小 可以减少通带和阻带的波动 使通带尽可 能趋近水平 阻带尽可能达到最大衰减 但通常此时过渡带会变宽 3 取不同的窗函数对幅度特性的整形效果比单纯的增加窗口长度要强得多 3 2 设计方法 这种方法也叫傅里叶级数法 一般是先给出所要求的理想的滤波器的频率响应 要求设计一个 FIR 滤波器频率响应来逼近 j d eH 1 0 N n njj enheH j d eH 设计是在时域进行的 因而先由的傅里叶反变换导出 即 j d eH nhd 式 3 2 1 deeHnh njj dd 2 1 由于是矩形频率响应特性 故一定是无限长序列 且是非因果的 而 FIR j d eH nhd 滤波器的必然是有限长的 所以要用有限长的来逼近无限长的 最有效 nh nh nhd 的方法是截断或者说用一个有限长度的窗口函数序列来截取 即 nhd n nhd 式 3 2 2 nhnnh d 因而窗函数序列的形状及长度的选择就是关键 我们以一个截止频率为的线性相位的理想矩形幅度特性的低通滤波器为例来讨论 c 设低通特性的群延时为 即 式 3 2 3 cc cc j j d e eH 0 这表明 在通带 范围内 的幅度是均匀的 其值为 1 相位是 c j d eH 利用 1 式可得 式 3 2 4 sin 2 1 n n deenh c ccnjj d c c 是中心点在的偶对称无限长非因果序列 要得到有限长的 一种最简单的方 nhd nh 法就是取矩形窗 即 nRN nRn N 但是按照线形相位滤波器的约束 必须是偶对称的 对称中心应为长度的一半 N 1 nh 2 因而必须 N 1 2 所以有 式 3 2 5 2 1 N 0 10 为其他n Nnnh nnhnh d d 将 式 3 2 4 代入 式 3 25 可得 式 3 2 6 1 0 0 2 1 2 1 sin Nn n N n N n nh c c c 为其他值 此时 一定满足这一线性相位的条件 1 nNhnh 下面求的傅里叶变换 也就是找出待求 FIR 滤波器的频率特性 以便能看出加窗 nh 处理后究竟对频率响应有何影响 按照复卷积公式 在时域是相乘 频域上是周期性卷积关系 即 式 3 2 7 deeHeH jj d j 2 1 因而逼近的好坏 完全取决于窗函数的频率特性 j eH j d eH j eW 窗函数的频率特性为 n j eW 式 3 2 8 1 0 N n njj eneW 对矩形窗 则有 nRN 式 3 2 9 2 sin 2 sin 2 1 1 0 N N eeeW N j N n njj R 也可表示成幅度函数与相位函数 式 3 2 10 2 1 N j R j N eWeW 其中 式 3 2 11 2 sin 2 sin N N WR 就是频域抽样内插函数 其幅度函数在之内为一个主瓣 两 j R eW R W N 2 侧形成许多衰减振荡的旁瓣 如果将理想频率响应也写成 式 3 2 12 2 1 N j d j d eHeH 则其幅度函数为 式 3 2 13 c c d H 0 1 3 3 窗函数介绍 实际应用的窗函数 可分为以下主要类型 1 幂窗 采用时间变量某种幂次的函数 如矩形 三角形 梯形或其它时间 t 的高次幂 2 三角函数窗 应用三角函数 即正弦或余弦函数等组合成复合函数 例如汉宁窗 海明 窗等 3 指数窗 采用指数时间函数 如形式 例如高斯窗等 st e 下面介绍几种常用窗函数的性质和特点 1 矩形窗 矩形窗属于时间变量的零次幂窗 函数形式为 式 3 3 1 Tt Tt Tt 0 1 相应的窗谱为 式 3 3 2 T T W sin2 矩形窗使用最多 习惯上不加窗就是使信号通过了矩形窗 这种窗的优点是主瓣比较 集中 缺点是旁瓣较高 并有负旁瓣 导致变换中带进了高频干扰和泄漏 甚至出现负谱 现象 图 3 1 矩形窗的时域及频域波形 2 三角窗亦称费杰 Fejer 窗 是幂窗的一次方形式 其函数形式是 式 3 3 3 Tt Tt T t Tt 0 1 1 三角窗与矩形窗比较 主瓣宽约等于矩形窗的两倍 但旁瓣小 而且无负旁瓣 如图 3 2 所示 图 3 2 三角窗的时域及频域波形 3 汉宁 Hanning 窗 汉宁窗又称升余弦窗 其时域表达式为 式 3 3 4 Tt Tt T t Tt 0 cos 2 1 2 1 1 相应的窗谱为 式 3 3 5 T T T T T T W sin sin 2 1sin 由此式可以看出 汉宁窗可以看作是 3 个矩形时间窗的频谱之和 或者说是 3 个 sin t 型函数之和 而括号中的两项相对于第一个谱窗向左 右各移动了 T 从而使旁 瓣互相抵消 消去高频干扰和漏能 可以看出 汉宁窗主瓣加宽并降低 旁瓣则显著减小 从减小泄漏观点出发 汉宁窗优于矩形窗 但汉宁窗主瓣加宽 相当于分析带宽加宽 频 率分辨率下降 4 海明 Hamming 窗 海明窗也是余弦窗的一种 又称改进的升余弦窗 其时间函数表达式为 式 3 3 6 Tt Tt T t Tt 0 cos4 054 0 1 其窗谱为 式 3 3 7 T T T T T T W sin sin 46 0 sin 08 1 海明窗与汉宁窗都是余弦窗 只是加权系数不同 海明窗加权的系数能使旁瓣达到更 小 分析表明 海明窗的第一旁瓣衰减为 42dB 海明窗的频谱也是由 3 个矩形窗的频谱合 成 但其旁瓣衰减速度为 20dB 10oct 这比汉宁窗衰减速度慢 海明窗与汉宁窗都是 很有用的窗函数 5 高斯窗 高斯窗是一种指数窗 其时域函数为 式 3 3 8 Tt Tte Tt at 0 1 2 式中 a 为常数 决定了函数曲线衰减的快慢 a 值如果选取适当 可以使截断点 T 为 有限值 处的函数值比较小 则截断造成的影响就比较小 高斯窗谱无负的旁瓣 第一旁瓣 衰减达一 55 dB 高斯窗的主瓣较宽 故而频率分辨率低 高斯窗函数常被用来截断一些 非周期信号 如指数衰减信号等 不同的窗函数对信号频谱的影响是不一样的 这主要是 因为不同的窗函数 产生泄漏的大小不一样 频率分辨能力也不一样 信号的截断产生了 能量泄漏 而用 FFT 算法计算频谱又产生了栅栏效应 从原理上讲这两种误差都是不能消 除的 但是我们可以通过选择不同的窗函数对它们的影响进行抑制 图 3 3 是几种常用的 窗函数的时域和频域波形 其中矩形窗主瓣窄 旁瓣大 频率识别精度最高 幅值识别精 度最低 布莱克曼窗主瓣宽 旁瓣小 频率识别精度最低 但幅值识别精度最高 图 3 3 几种常用的窗函数的时域和频域波形 对于窗函数 还有一些要求 1 3dB 带宽 B 它是主瓣归一化的幅度下降到 3dB 时的带宽 当数据长度为 N 时 矩 形窗主瓣两个过零点之间的宽度为 4 N 2 最大边瓣峰值 A dB 3 边瓣谱峰渐进衰减速度 D dB cot 所以 理想的窗函数应当具有最小的 B 和 A 和最大的 D 3 4 窗函数法设计步骤 1 首先是给定所要求的频率响应函数 j d eH 2 其次 求单位冲激响应 deeHnh njj dd 2 1 3 再次 有过渡带宽及阻带最小衰减的要求 查表选定窗函数及 N 的大小 一般 N 的 大小要通过几次试探而后确定 4 求得所设计的 FIR 滤波器的单位冲激响应 n 0 1 N 1 nhnnh d 5 求 检验是否满足设计要求 如不满足 则需要重新设计 1 0 N n njj enheH 3 5 设计实例 线性相位 FIR 低通滤波器的设计 用窗函数法 指标要求 通带截止频率 0 2 阻带 起始频率 0 4 阻带最小衰减 50dB 1 设为理想线性相位滤波器 j eH 其他 0 c j j e eH 由所需低通滤波器的过渡带求出理想低通滤波器的截止数字频率 0 3 得出 为线性相位所需的移位 2 1 sin 2 1 2 1 N n nnn dedeenh c c njnjj d c c 2 由阻带衰减确定窗函数 由过渡带宽确定 N 值 阻带最小衰减 50dB 比对 6 种窗函 数基本参数选定窗函数为海明窗 所要求的过渡带宽 0 4 0 2 0 2 N 6 6 0 2 33 N 1 2 16 3 由海明窗函数确定 FIR 滤波器的 h n 2 1 2 1 sin 1 2 cos 46 0 54 0 n N n N n nh nR N n c d N 得出 1 2 cos 46 0 54 0 16 16 3 0sin nR N n n nnhnh Nd 4 仿真检验各项指标 得出结论 满足设计要求 取 N 33 偶对称 得 过渡带宽 0 3476563 第一通带波纹 0 020837dB 第一阻带最小衰减 60 9159dB 我们设计的低通 FIR 数字滤波器的理想性能指标是 通带截止频率 Wp 0 3 阻带截 止频率 Ws 0 5 阻带衰减 At 不小于 40dB 通带衰减不大于 3dB fs 4000 Matlab 程序 wp 0 3 pi ws 0 5 pi wdelta ws wp N ceil 8 pi wdelta Wn 0 3 0 5 pi 2 b fir1 N Wn pi hanning N 1 freqz b 1 512 N ceil 8 0 15 n 0 N 1 window hanning N h1 w freqz window 1 figu