有限脉冲响应数字滤波器的设计.ppt

第7章有限脉冲响应数字滤波器的设计 7 1线性相位FIR数字滤波器的条件和特点7 2利用窗函数法设计FIR滤波器7 3利用频率采样法设计FIR滤波器7 4利用等波纹最佳逼近法设计FIR数字滤波器7 5IIR和FIR数字滤波器的比较 有限脉冲响应 FIR 滤波器在保证幅度特性满足技术指标的同时 很容易做到有严格的线性相位特性 用N表示FIR滤波器单位脉冲响应h n 的长度 其系统函数H z 为H z 是z 1的N 1次多项式 有N 1个零点 在原点z 0处有一个N 1重极点 因此 H z 永远稳定 稳定和线性相位特性是FIR滤波器最突出的优点 FIR滤波器设计任务 选择有限长度的h n 使频率响应函数H ej 满足技术指标要求 7 1线性相位FIR数字滤波器的条件和特点 1 线性相位FIR数字滤波器对于长度为N的h n 传输函数为 式中 Hg 称为幅度特性 称为相位特性 注意 这里Hg 不同于 H ej Hg 为 的实函数 可能取负值 而 H ej 总是正值 7 1 1 7 1 2 H ej 线性相位是指 是 的线性函数 即 为常数 7 1 3 如果 满足下式 0 0是起始相位 7 1 4 严格地说 此时 不具有线性相位 但以上两种情况都满足群时延 相位特性曲线的斜率 是一个常数 即也称这种情况为线性相位 满足 7 1 3 式是第一类线性相位 满足 7 1 4 式是第二类线性相位 0 2是第二类线性相位特性常用的情况 2 线性相位FIR的时域约束条件线性相位FIR滤波器的时域约束条件是指满足线性相位时 对h n 的约束条件 1 第一类线性相位对h n 的约束条件 第一类线性相位FIR数字滤波器的相位函数 可得 7 1 5 由上式得到 7 1 6 将 7 1 6 式中两式相除得到 即 移项并用三角公式化简得到 7 1 7 满足 7 1 7 式的一组解是 函数h n sin n 关于求和区间的中心 N 1 2奇对称 因为sin n 关于n 奇对称 如果取 N 1 2 则要求h n 关于 N 1 2偶对称 所以要求 和h n 满足如下条件 7 1 8 即 如果要求单位脉冲响应为h n 长度为N的FIR数字滤波器具有第一类线性相位特性 严格线性相位特性 则h n 应当关于n N 1 2点偶对称 当N确定时 FIR数字滤波器的相位特性是一个确知的线性函数 即 N 1 2 N为奇数和偶数时 h n 的对称情况如表7 1 1中的情况1和情况2所示 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 2 第二类线性相位对h n 的约束条件 第二类线性相位FIR数字滤波器的相位函数 2 可得 满足式 7 1 9 的一组解是 函数h n cos n 关于求和区间的中心 N 1 2奇对称 因为cos n 关于n 偶对称 如果取 N 1 2 则要求h n 关于 N 1 2奇对称 所以要求 和h n 满足如下条件 7 1 9 7 1 10 即 如果要求单位脉冲响应为h n 长度为N的FIR数字滤波器具有第二类线性相位特性 则h n 应当关于n N 1 2点奇对称 当N确定时 FIR数字滤波器的相位特性是一个确知的线性函数 即 2 N 1 2 N为奇数和偶数时h n 的对称情况如表7 1 1中情况3和情况4 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 2 线性相位FIR滤波器幅度特性Hg 的特点 幅度特性的特点就是线性相位FIR滤波器的频域约束条件 将时域约束条件h n h N n 1 代入下式设h n 为实序列 即可推导出线性相位条件对FIR数字滤波器的幅度特性Hg 的约束条件 当N取奇数和偶数时对Hg 的约束不同 因此 对于两类线性相位特性 下面分四种情况讨论其幅度特性的特点 这些特点对正确设计线性相位FIR数字滤波器具有重要的指导作用 式中 表示取不大于 N 1 2的最大整数 显然 仅当N为奇数时 M N 1 2 情况1 h n h N n 1 N为奇数 将时域约束条件h n h N n 1 和 代入式 7 1 1 和 7 1 2 得到 为了推导方便 引入两个参数符号 所以 7 1 11 幅度特性Hg 分析 因为cos n 关于 0 2 三点偶对称 所以Hg 关于 0 2 三点偶对称 因此情况1可以实现各种 低通 高通 带通 带阻 滤波器 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 情况2 h n h N n 1 N为偶数 仿照情况1的推导方法得到 7 1 12 式中 因为 是偶数 所以当时 幅度特性Hg 分析 因为cos n 关于过零点奇对称 关于 0和2 偶对称 所以Hg 0 Hg 关于 奇对称 关于 0和2 偶对称 因此 情况2不能实现高通和带阻滤波器 对N 12的低通情况 Hg 如表7 1 1中情况2所示 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 情况3 h n h N n 1 N为奇数 将时域约束条件 h n h N n 1 和 2 代入式并考虑得到 式中 N是奇数 N 1 2是整数 所以 当 0 2 时 sin n 0 而且sin n 关于过零点奇对称 因此Hg 关于 0 2 三点奇对称 由此可见 情况3只能实现带通滤波器 对N 13的带通滤波器举例 Hg 如表7 1 1中情况3所示 幅度特性Hg 分析 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 情况4 h n h N n 1 N为偶数 用情况3的推导过程可以得到 幅度特性Hg 分析 式中 N是偶数 N 1 2 N 2 1 2 当 0 2 时 sin n 0 当 时 sin n 1 n N 2 为峰值点 而sin n 关于过零点 0和2 两点奇对称 关于峰值点 偶对称 因此Hg 关于 0和2 两点奇对称 关于 偶对称 由此可见 情况4不能实现低通和带阻滤波器 对N 12的高通滤波器举例 Hg 如表7 1 1中情况4所示 表7 1 1线性相位FIR数字滤波器的时域和频域特性一览 应当注意 对每一种情况表7 1 1仅画出满足幅度特性要求的一种例图 例如 情况1仅以低通的幅度特性曲线为例 当然也可以画出满足情况1的幅度约束条件 Hg 关于 0 2 三点偶对称 的高通 带通和带阻滤波器的幅度特性曲线 所以 仅从表7 1 1就认为情况1只能设计低通滤波器是错误的 3 线性相位FIR数字滤波器的零点分布特点将h n h N 1 n 代入上式 得到 7 1 14 由上式可以看出 如z zi是H z 的零点 其倒数也必然是其零点 又因为h n 是实序列 H z 的零点必定共轭成对 因此也是其零点 这样 线性相位FIR滤波器零点必定是互为倒数的共轭对 确定其中一个 另外三个零点也就确定了 如图7 1 1中 当然 也有一些特殊情况 如图7 1 1中z1 z2和z4情况 图7 1 1线性相位FIR数字滤波器的零点分布 4 线性相位FIR滤波器网络结构设N为偶数 则有 令m N n 1 则有 因为 如果N为奇数 则将中间项h N 1 2 单独列出 和直接型结构比较 N取偶数 直接型需要N个乘法器 而线性相位结构减少到N 2个乘法器 节约了一半的乘法器 N取奇数 乘法器减少到 N 1 2个 也近似节约了近一半的乘法器 FIR滤波器直接型结构 图5 5 1第一类线性相位网络结构流图 图5 5 2第二类线性相位网络结构流图 7 2利用窗函数法设计FIR滤波器 7 2 1窗函数法设计原理 设希望逼近的滤波器频率响应函数为Hd ej 其单位脉冲响应是hd n 通常以理想滤波器作为Hd ej 其幅度特性逐段恒定 在边界频率处有不连续点 因而hd n 是无限时宽的 且是非因果序列 例如 线性相位理想低通滤波器的频率响应函数Hd ej 为 其单位脉冲响应hd n 为 7 2 2 7 2 1 hd n 是无限长非因果序列 hd n 的波形如图7 2 1 a 所示 图7 2 1窗函数设计法的时域波形 矩形窗 N 30 为了构造一个长度为N的第一类线性相位FIR滤波器 只有将hd n 截取一段 并保证截取的一段关于n N 1 2偶对称 设截取的一段用h n 表示 即 RN n 是一个矩形序列 长度为N 波形如图7 2 1 b 所示 当 取值为 N 1 2时 截取的一段h n 关于n N 1 2偶对称 保证所设计的滤波器具有线性相位 h n 就是所设计滤波器的单位脉冲响应 长度为N 其系统函数为H z 7 2 3 图7 2 1窗函数设计法的时域波形 矩形窗 N 30 RN n 矩形序列 就是起对无限长序列的截断作用 可以形象地把RN n 看做一个窗口 h n 则是从窗口看到的一段hd n 序列 所以称h n hd n RN n 为用矩形窗对hd n 进行加窗处理 用有限长的序列h n 去代替hd n 肯定会引起误差 在频域会出现吉布斯 Gibbs 效应 1 引起过渡带加宽2 通带和阻带内出现波动 尤其使阻带的衰减小 吉布斯效应是由于将hd n 直接截断引起的 也称截断效应 图7 2 2吉普斯效应 Hd ej 是一个以2 为周期的函数 可以展为傅里叶级数 即 傅里叶级数的系数为hd n 就是Hd ej 对应的单位脉冲响应 设计FIR滤波器就是 根据要求找到N个傅里叶级数系数h n n 0 1 2 N 1 以N项傅氏级数近似代替无限项傅氏级数 这样在一些频率不连续点附近会引起较大误差 这种误差就是截断效应 窗函数法也称为傅氏级数法 根据傅里叶变换的频域卷积定理 式的傅里叶变换 Hd ej 和WR ej 分别是hd n 和RN n 的傅里叶变换 分析用矩形窗截断的影响和改进的措施 WRg 称为矩形窗的幅度函数 如图7 2 3 b WRg 的主瓣 图中 2 N 2 N 区间上的一段波形WRg 的旁瓣 其余较小的波动 按照式 Hd ej Hdg e j 理想低通滤波器的幅度特性函数 图7 2 3 a 为 图7 2 3矩形窗加窗效应 将Hd ej 和WR ej 代入式 7 2 6 则 得到 加窗后的滤波器的幅度特性等于理想低通滤波器的幅度特性Hdg 与矩形窗幅度特性WRg 的卷积 图7 2 3矩形窗加窗效应 图7 2 2 f 表示Hdg 与WRg 卷积形成的Hg 波形当 0时 Hg 0 等于图7 2 2 a 与 b 两波形乘积的积分 相当于对WRg 在 c之间一段波形的积分 当 c 2 N时 近似为 之间波形的积分 将H 0 值归一化到1 当 c时 情况如图7 2 2 c 所示 当 c 2 N时 积分近似为WRg 一半波形的积分 对Hg 0 归一化后的值近似为1 2 3 当 c 2 N时 情况如图7 2 2 d 所示 WR 主瓣完全在区间 c c 之内 而最大的一个负旁瓣移到区间 c c 之外 因此Hg c 2 N 有一个最大的正峰 当 c 2 N时 情况如图7 2 2 e 所示 WRg 主瓣完全移到积分区间外边 由于最大的一个负旁瓣完全在区间 c c 内 因此Hg c 2 N 形成最大的负峰 图7 2 2表明 Hg 最大的正峰与最大的负峰对应的频率相距4 N 通过以上分析可知 对hd n 加矩形窗处理后 Hg 与原理想低通Hdg 的差别有以下两点 对hd n 加矩形窗处理后 Hg 与Hdg 的差别 1 在理想特性不连续点 c附近形成过渡带 过渡带的宽度近似等于WRg 主瓣宽度4 N 2 通带内产生了波纹 最大的峰值在 c 2 N处 阻带内产生了余振 最大的负峰在 c 2 N处 通带与阻带中波纹的情况与窗函数的幅度谱有关 WRg 旁瓣幅度的大小直接影响Hg 波纹幅度的大小 以上两点就是对hd n 用矩形窗截断后的吉布斯效应 通带内的波纹影响滤波器通带的平稳性 阻带内的波纹影响阻带内的衰减 可能使最小衰减不满足技术指标要求 如何减少吉布斯效应的影响 增加矩形窗的长度就可减少吉布斯效应的影响 分析一下N加大时WRg 的变化 在主瓣附近 按照式 7 2 5 WRg 可近似为该函数的性质 N加大时 主瓣幅度加高 同时旁瓣也加高 保持主瓣和旁瓣幅度相对值不变 N加大时 WRg 的主瓣和旁瓣宽度变窄 波动的频率加快 三种不同长度的矩形窗函数的幅度特性WRg 曲线如图7 2 4 a b c 所示 用这三种窗函数设计的FIR滤波器的幅度特性Hg 曲线如图7 2 4 d e f 所示 因此 当N加大时 Hg 的波动幅度没有多大改善 带内最大肩峰比H 0 高8 95 阻带最大负峰值为H 0 的8 95 使阻带最小衰减只有21dB 加大N只能使Hg 过渡带变窄 过渡带近似为主瓣宽度4 N 因此加大N 并不是减小吉布斯效应的有效方法 图7 2 4矩形窗函数长度的影响 分析表明 调整窗口长度N只能有效地控制过渡带的宽度 而要减少带内波动以及增大阻带衰减 只能从窗函数的形状上找解决问题的方法 构造新的窗函数形状 使其谱函数的主瓣包含更多的能量 相应旁瓣幅度更小 旁瓣的减小可使通带 阻带波动减小 从而加大阻带衰减 但这样总是以加宽过渡带为代价的 7 2 2典型窗函数介绍 介绍几种常用窗函数的时域表达式 时域波形 幅度特性函数 衰减用dB计量 曲线 以及用该窗函数设计的FIR数字滤波器的单位脉冲响应和损耗函数曲线 Hd ej 取理想低通 c 2 窗函数长度N 31 窗函数的几个参数 旁瓣峰值 n 窗函数的幅频函数 Wg 的最大旁瓣的最大值相对主瓣最大值的衰减值 dB 过渡带宽度Bg 用该窗函数设计的FIR数字滤波器 FIRDF 的过渡带宽度 阻带最小衰减 s 用该窗函数设计的FIRDF的阻带最小衰减 图7 2 4所示的矩形窗的参数为 n 13dB Bg 4 N s 21dB 1 矩形窗 RectangleWindow wR n RN n 其幅度函数为 2 三角形窗 BartlettWindow 7 2 8 其频谱函数为 其幅度函数为 三角窗的四种波形如图7 2 5所示 参数为 n 25dB Bg 8 N s 25dB 图7 2 5三角窗的四种波形 3 汉宁 Hanning 窗 升余弦窗 7 2 11 当N 1时 N 1 N汉宁窗的幅度函数WHng 由三部分相加 旁瓣互相对消 使能量更集中在主瓣中 汉宁窗的四种波形如图7 2 6所示 参数为 n 31dB Bg 8 N s 44dB 图7 2 6汉宁窗的四种波形 4 哈明 Hamming 窗 改进的升余弦窗其频谱函数WHm ej 为 其幅度函数WHmg 为当N 时 其可近似表示为 这种改进的升余弦窗 能量更加集中在主瓣中 主瓣的能量约占99 963 旁瓣峰值幅度为40dB 但其主瓣宽度和汉宁窗的相同 仍为8 N 哈明窗是一种高效窗函数 所以MATLAB窗函数设计函数的默认窗函数就是哈明窗 哈明窗的四种波形如图7 2 7所示 参数为 n 41dB Bg 8 N s 53dB 图7 2 7哈明窗的四种波形 7 2 13 5 布莱克曼 Blackman 窗 其频谱函数为 其幅度函数为其幅度函数由五部分组成 它们都是移位不同 且幅度也不同的WRg 函数 使旁瓣再进一步抵消 旁瓣峰值幅度进一步增加 其幅度谱主瓣宽度是矩形窗的3倍 布莱克曼窗的四种波形如图7 2 8所示 参数为 n 57dB B 12 N s 74dB 图7 2 8布莱克曼窗的四种波形 6 凯塞 贝塞尔窗 Kaiser BaselWindow 以上五种窗函数都称为参数固定窗函数 每种窗函数的旁瓣幅度都是固定的 凯塞 贝塞尔窗是一种参数可调的窗函数 是一种最优窗函数 7 2 15 式中 I0 是零阶第一类修正贝塞尔函数 可用下面级数计算 一般I0 取15 25项 便可以满足精度要求 参数可以控制窗的形状 一般 加大 主瓣加宽 旁瓣幅度减小 典型数据为4 9 当 5 44时 窗函数接近哈明窗 7 865时 窗函数接近布莱克曼窗 在设计指标给定时 可以调整 值 使滤波器阶数最低 所以其性能最优 凯塞 Kaiser 给出的估算 和滤波器阶数N的公式如下 7 2 17 式中 Bt s p 是数字滤波器过渡带宽度 应当注意 因为式 7 2 17 为阶数估算 所以必须对设计结果进行检验 另外 凯塞窗函数没有独立控制通带波纹幅度 实际中通带波纹幅度近似等于阻带波纹幅度 凯塞窗的幅度函数为 7 2 16 7 2 18 对 的8种典型值 将凯塞窗函数的性能列于表7 2 1中 供设计者参考 由表可见 当 5 568 时 各项指标都好于哈明窗 6种典型窗函数基本参数归纳在表7 2 2中 可供设计时参考 表7 2 1凯塞窗参数对滤波器的性能影响 表7 2 26种窗函数的基本参数 表中过渡带宽和阻带最小衰减是用对应的窗函数设计的FIR数字滤波器的频率响应指标 图7 2 4常用的窗函数 随着数字信号处理的不断发展 学者们提出的窗函数已多达几十种 除了上述6种窗函数外 比较有名的还有Chebyshev窗 Gaussian窗 5 6 MATLAB信号处理工具箱提供了14种窗函数的产生函数 下面列出上述6种窗函数的产生函数及其调用格式 wn boxcar N 列向量wn中返回长度为N的矩形窗函数w n wn bartlett N 列向量wn中返回长度为N的三角窗函数w n wn hanning N 列向量wn中返回长度为N的汉宁窗函数w n wn hamming N 列向量wn中返回长度为N的哈明窗函数w n wn blackman N 列向量wn中返回长度为N的布莱克曼窗函数w n wn kaiser N beta 列向量wn中返回长度为N的凯塞 贝塞尔窗函数w n 7 2 3用窗函数法设计FIR滤波器的步骤用窗函数法设计FIR滤波器的步骤如下 1 根据对过渡带及阻带衰减的指标要求 选择窗函数的类型 并估计窗口长度N a 窗函数类型选择 按照阻带衰减选择 原则 在保证阻带衰减满足要求的情况下 尽量选择主瓣窄的窗函数 b 窗口长度N估计 根据过渡带宽度 2 构造希望逼近的频率响应函数Hd ej 即 对所谓的 标准窗函数法 就是选择Hd ej 为线性相位理想滤波器 理想低通 理想高通 理想带通 理想带阻 理想滤波器的截止频率 c近似位于最终设计的 FIRDF的过渡带的中心频率点 幅度函数衰减一半 约 6dB 所以如果设计指标给定通带边界频率和阻带边界频率 p和 s 一般取 3 计算hd n 如果给出待求滤波器的频响函数为Hd ej 那么单位脉冲响应用下式求出 如果Hd ej 较复杂 或者不能用封闭公式表示 则不能用上式求出hd n 可以采用频域采样法求取 4 加窗得到设计结果 h n hd n w n 5 验算技术指标是否满足要求 设计出的滤波器频率响应用下式计算 例7 2 1 用窗函数法设计线性相位高通FIRDF 要求通带截止频率 p 2rad 阻带截止频率 s 4rad 通带最大衰减 p 1dB 阻带最小衰减 s 40dB 解 1 选择窗函数w n 计算窗函数长度N 已知阻带最小衰减 s 40dB 由表 7 2 2 可知汉宁窗和哈明窗均满足要求 选择汉宁窗 过渡带宽度Bt p s 4 汉宁窗的精确过渡带宽度Bt 6 2 N 所以要求Bt 6 2 N 4 解之得N 24 8 对高通滤波器N必须取奇数 取N 25 有 2 构造Hd ej 式中 3 求出hd n 将 12代入得 n 12 对应全通滤波器 是截止频率为3 8的理想低通滤波器的单位脉冲响应 二者之差就是理想高通滤波器的单位脉冲响应 4 加窗 7 2 4窗函数法的MATLAB设计函数简介实际设计时可调用MATLAB工具箱函数fir1实现窗函数法设计步骤 2 4 的解题过程 1 fir1用窗函数法设计线性相位FIR数字滤波器的工具箱函数 以实现线性相位FIR数字滤波器的标准窗函数法设计 标准 是指在设计低通 高通 带通和带阻FIR滤波器时 Hd ej 分别表示相应的线性相位理想低通 高通 带通和带阻滤波器的频率响应函数 因而将所设计的滤波器的频率响应称为标准频率响应 fir1的调用格式及功能 hn fir1 M wc 返回6dB截止频率为wc的M阶FIR低通滤波器系数向量hn 默认选用哈明窗 滤波器单位脉冲响应h n 与向量hn的关系为 h n hn n 1 n 0 1 2 M而且满足线性相位条件 h n h N 1 n 其中wc为对 归一化的数字频率 0 wc 1 当wc wcl wcu 时 得到的是带通滤波器 其 6dB通带为wcl wcu hn fir1 M wc ftype 可设计高通和带阻FIR滤波器 当ftype high时 设计高通FIR滤波器 当ftype stop时 且wc wcl wcu 时 设计带阻FIR滤波器 应当注意 在设计高通和带阻FIR滤波器时 阶数M只能取偶数 h n 长度N M 1为奇数 不过 当用户将M设置为奇数时 fir1会自动对M加1 hn fir1 M wc window 可以指定窗函数向量window 如果缺省window参数 则fir1默认为哈明窗 例如 hn fir1 M wc bartlett M 1 使用Bartlett窗设计 hn fir1 M wc blackman M 1 使用blackman窗设计 hn fir1 M wc ftype window 通过选择wc ftype和window参数 含义同上 可以设计各种加窗滤波器 2 fir2为任意形状幅度特性的窗函数法设计函数 用fir2设计时 可以指定任意形状的Hd ej 它实质是一种频率采样法与窗函数法的综合设计函数 主要用于设计幅度特性形状特殊的滤波器 如数字微分器和多带滤波器等 例7 2 1的设计程序如下 例7 2 1用窗函数法设计线性相位高通FIR数字滤波器wp pi 2 ws pi 4 Bt wp ws 计算过渡带宽度N0 ceil 6 2 pi Bt 根据表7 2 2汉宁窗计算所需 h n 长度N0 ceil x 取大于等 于x的最小整数N N0 mod N0 1 2 确保h n 长度N是奇数wc wp ws 2 pi 计算理想高通滤波器通带截止 频率 关于 归一化 hn fir1 N 1 wc high hanning N 调用fir1计算高通FIR数字滤波 器的h n M 1024 hk fft hn M n 0 N 1 subplot 2 2 1 stem n hn line 0 30 0 0 xlabel n ylabel h n k 1 M 2 w 2 0 M 2 1 M subplot 2 2 2 plot w 20 log10 abs hk k axis 0 1 80 5 xlabel ylabel 20lg Hg gridon 运行程序得到h n 的25个值 h n 0 0004 0 00060 00280 0071 0 0000 0 0185 0 02100 01650 06240 03550 1061 0 28980 6249 0 2898 0 10610 03550 06240 0165 0 02100 0185 0 00000 00710 0028 0 0006 0 0004 高通FIR数字滤波器的h n 及损耗函数如图7 2 9所示 图7 2 9高通FIR数字滤波器的h n 波形及损耗函数曲线 例7 2 2 对模拟信号进行低通滤波处理 要求通带0 f 1 5kHz内衰减小于1dB 阻带2 5kHz f 上衰减大于40dB 希望对模拟信号采样后用线性相位FIR数字滤波器实现上述滤波 采样频率Fs 10kHz 用窗函数法设计满足要求的FIR数字低通滤波器 求出h n 并画出损耗函数曲线 为了降低运算量 希望滤波器阶数尽量低 解 1 确定相应的数字滤波器指标 通带截止频率为 阻带截止频率为阻带最小衰减为 s 40dB 2 用窗函数法设计FIR数字低通滤波器 为了降低阶数选择凯塞窗 根据式 7 2 16 计算凯塞窗的控制参数为 指标要求过渡带宽度Bt s p 0 2 根据式 7 2 17 计算滤波器阶数为 取满足要求的最小整数M 23 所以h n 长度为N M 1 24 理想低通滤波器的通带截止频率 c s p 2 0 4 所以由式 7 2 2 和式 7 2 3 得到 式中 w n 是长度为24 3 395 的凯塞窗函数 实现本例设计的MATLAB程序 用凯塞窗函数设计线性相位低通FIR数字滤波器fp 1500 fs 2500 rs 40 wp 2 pi fp Fs ws 2 pi fs Fs Bt ws wp 计算过渡带宽度 alph 0 5842 rs 21 0 4 0 07886 rs 21 根据 7 2 16 式计算kaiser窗的控制参数 N ceil rs 8 2 285 Bt 计算kaiser窗所需阶数Nwc wp ws 2 pi 计算理想高通滤波器通带截止频率 关于 归一化 hn fir1 N wc kaiser N 1 alph 调用kaiser计算低通FIRDF的h n M 1024 hk fft hn M n 0 N subplot 2 2 1 stem n hn line 0 30 0 0 xlabel n ylabel h n k 1 M 2 w 2 0 M 2 1 M subplot 2 2 2 plot w 20 log10 abs hk k axis 0 1 80 5 xlabel ylabel 20lg Hg gridon 运行程序得到h n 的24个值 h n 0 00390 0041 0 0062 0 01470 00000 02860 0242 0 0332 0 07550 00000 19660 37240 37240 1966 0 0000 0 0755 0 03320 02420 02860 0000 0 0147 0 00620 00410 0039 低通FIR数字滤波器的h n 波形和损耗函数曲线如图7 2 10所示 图7 2 10低通FIR数字滤波器的h n 波形及损耗函数曲线 例7 2 3 窗函数法设计一个线性相位FIR带阻滤波器 要求通带下截止频率 lp 0 2 阻带下截止频率 ls 0 35 阻通带上截止频率 us 0 65 通带上截止频率 up 0 8 通带最大衰减 p 1dB 阻带最小衰减 s 60dB 解本例直接调用fir1函数设计 因为阻带最小衰减 s 60dB 所以选择布莱克曼窗 再根据过渡带宽度选择滤波器长度N 布莱克曼窗的过渡带宽度Bt 12 N 所以 解之得N 80 调用参数 设计程序为ep723 m 参数计算也由程序完成 ep723 m 例7 2 3用窗函数法设计线性相位带阻FIR数字滤波器 wlp 0 2 pi wls 0 35 pi wus 0 65 pi wup 0 8 pi 设计指标参数赋值 B wls wlp 过渡带宽度 N ceil 12 pi B 计算阶数N ceil x 为大于等于x的最小整数 wp wls wlp 2 pi wus wup 2 pi 设置理想带通截止频率 hn fir1 N wp stop blackman N 1 带阻滤波器要求h n 长度为奇数 所以取N 1 省略绘图部分程序运行结果 N 81 由于h n 数据量太大 因而仅给出h n 的波形及损耗函数曲线 如图7 2 11所示 图7 2 11带阻FIR数字滤波器的h n 波形及损耗函数曲线 窗函数设计法简单方便 易于实现 但存在以下缺点 滤波器边界频率不易精确控制 窗函数设计法总使通带和阻带波纹幅度相等 不能分别控制通带和阻带波纹幅度 但是工程上对二者的要求是不同的 希望能分别控制 所设计的滤波器在阻带边界频率附近的衰减最小 距阻带边界频率越远 衰减越大 所以 如果在阻带边界频率附近的衰减刚好达到设计指标要求 则阻带中其他频段的衰减就有很大富余量 说明这种设计法存在较大的资源浪费 或者说所设计滤波器的性能价格比低 7 4等波纹最佳逼近法设计FIR数字滤波器等波纹最佳逼近法是一种优化设计法 它克服了窗函数设计法和频率采样法的缺点 使最大误差 波纹的峰值 最小化 并在整个频段上均匀分布 等波纹 用等波纹最佳逼近法设计的FIR数字滤波器的幅频响应在通带和阻带都是等波纹的 而且可以分别控制通带和阻带波纹幅度 最佳逼近 在滤波器长度给定条件下 使加权误差波纹幅度最小化 本方法设计的滤波器性价比最高 阶数相同时 最大逼近误差最小 即通带最大衰减最小 阻带最小衰减最大指标相同时 阶数最低 7 4等波纹最佳逼近法设计FIR数字滤波器7 4 1等波纹最佳逼近法的基本思想用Hd 表示希望逼近的幅度特性函数 要求设计线性相位FIR数字滤波器时 Hd 必须满足线性相位约束条件 用Hg 表示实际设计的滤波器幅度特性函数 定义加权误差函数E 为 7 4 1 W 误差加权函数 用来控制不同频段 一般指通带和阻带 的逼近精度 等波纹最佳逼近基于切比雪夫逼近 在通带和阻带以 E 的最大值最小化为准则 采用Remez多重交换迭代算法求解滤波器系数h n 所以W 取值越大的频段 逼近精度越高 开始设计时应根据逼近精度要求确定W 在Remez多重交换迭代过程中W 是确知函数 等波纹最佳逼近设计中 把数字频段分为 逼近 或研究 区域 和 无关区域 逼近区域一般指通带和阻带 而无关区域一般指过渡带 设计过程中只考虑对逼近区域的最佳逼近 应当注意 无关区宽度不能为零 即Hd 不能是理想滤波特性 利用等波纹最佳逼近准则设计线性相位FIR数字滤波器数学模型的建立及其求解算法的推导复杂 求解计算必须借助计算机 幸好滤波器设计专家已经开发出MATLAB信号处理工具箱函数remezord和remez 只要简单地调用这两个函数就可以完成线性相位FIR数字滤波器的等波纹最佳逼近设计 在介绍MATLAB工具箱函数remezord和remez之前 先介绍等波纹滤波器的技术指标及其描述参数 图7 4 1给出了等波纹滤波器技术指标的两种描述参数 图7 4 1 a 用损耗函数描述 即 p 2 p 2dB s 11 20 s 20dB 这是工程实际中常用的指标描述方法 但是 用等波纹最佳逼近设计法求滤波器阶数N和误差加权函数W 时 要求给出滤波器通带和阻带的振荡波纹幅度 1和 2 图7 4 1 b 给出了用通带和阻带的振荡波纹幅度 1和 2描述的技术指标 显然 两种描述参数之间可以换算 如果设计指标以 p和 s给出 为了调用MATLAB工具箱函数remezord和remez进行设计 就必须由 p和 s换算出通带和阻带的振荡波纹幅度 1和 2 对比图7 4 2 a 和 b 得出关系式 7 4 2 7 4 3 由式 7 4 2 和 7 4 3 得到 7 4 4 7 4 5 图7 4 1等波纹滤波器的幅频特性函数曲线及指标参数 误差加权函数W 的作用 以及滤波器阶数N和波纹幅度 1和 2的制约关系 设期望逼近的通带和阻带分别为 0 4 和 5 16 对下面四种不同的控制参数 等波纹最佳逼近的损耗函数曲线分别如图7 4 2 a b c 和 d 所示 图中 W w1 w2 表示第一个逼近区 0 4 上的误差加权函数W w1 第二个逼近区 5 16 上的误差加权函数W w2 图7 4 2 a 中 通带频段 上的W 1 阻带频段 5 16 上的W 10 图7 4 2误差加权函数W 和滤波器阶数N对逼近精度的影响 比较图7 4 2 a b c 和 d 可以得出结论 1 当N一定时 误差加权函数W 较大的频带逼近精度较高 W 较小的频带逼近精度较低 2 如果改变W 使通 阻 带逼近精度提高 则必然使阻 通 带逼近精度降低 3 滤波器阶数N增大才能使通带和阻带逼近精度同时提高 所以 W 和N由滤波器设计指标 即 p和 s以及过渡带宽度 确定 1 根据给定的逼近指标估算滤波器阶数N和误差加权函数W 2 采用remez算法得到滤波器单位脉冲响应h n MATLAB工具箱函数remezord和remez就是完成以上2个设计步骤的有效函数 用等波纹最佳逼近法设计FIR数字滤波器的过程是 7 4 2remez和remezord函数及滤波器设计指标1 remez和remezord函数 1 remez 采用remez算法可实现线性相位FIR数字滤波器的等波纹最佳逼近设计 其调用格式为 hn remez M f m w 调用结果返回单位脉冲响应向量hn remez函数的调用参数 M f m w 一般通过调用remezord函数来计算 调用参数含义如下 M为FIR数字滤波器阶数 hn长度N M 1 f和m给出希望逼近的幅度特性 f为边界频率向量 0 f 1 要求f为单调增向量 即f k f k 1 k 1 2 而且从0开始 以1结束 1对应数字频率 模拟频率Fs 2 Fs表示时域采样频率 m是与f对应的幅度向量 m与f长度相等 m k 表示频点f k 的幅度响应值 如果用命令Plot f m 画出幅频响应曲线 则k为奇数时 频段 f k f k 1 上的幅频响应就是期望逼近的幅频响应值 频段 f k 1 f k 2 为无关区 简言之 Plot f m 命令画出的幅频响应曲线中 起始频段为第一段 奇数频段为逼近区 偶数频段为无关区 例如 对图7 4 2 f 0 1 4 5 16 1 m 1 1 0 0 plot f m 画出的幅度特性曲线如图7 4 3所示 图中奇数段 第一 三段 的水平幅度为希望逼近的幅度特性 偶数段 第二段 的下降斜线为无关部分 逼近时形成过渡带 并不考虑该频段的幅频响应形状 w为误差加权向量 其长度为f的一半 w i 表示对m中第i个逼近频段的误差加权值 图7 4 2 a 中 w 1 10 缺省w时 默认w为全1 即每个逼近频段的误差加权值相同 除了设计选频FIR数字滤波器 remez函数还可以设计两种特殊滤波器 希尔伯特变换器和数字微分器 调用格式分别为 hn remez M f m w hilbert hn remez M f m w defferentiator 希尔伯特变换器和数字微分器设计和应用的详细内容请参考文献 10 19 2 remezord 采用remezord函数 可根据逼近指标估算等波纹最佳逼近FIR数字滤波器的最低阶数M 误差加权向量w和归一化边界频率向量f 使滤波器在满足指标的前提下造价最低 其返回参数作为remez函数的调用参数 其调用格式为 M fo mo w remezord f m rip Fs 参数说明 f与remez中的类似 这里f可以是模拟频率 单位为Hz 或归一化数字频率 但必须从0开始 到Fs 2 用归一化频率时对应1 结束 而且其中省略了0和Fs 2两个频点 Fs为采样频率 缺省时默认Fs 2Hz 但是这里f的长度 包括省略的0和Fs 2两个频点 是m的两倍 即m中的每个元素表示f给定的一个逼近频段上希望逼近的幅度值 例如 对图7 4 3 f 1 4 5 16 m 1 0 图7 4 3希望逼近的幅度特性曲线 注意 省略Fs时 f中必须为归一化频率 有时估算的阶数M略小 使设计结果达不到指标要求 这时要取M 1或M 2 必须注意对滤波器长度N M 1的奇偶性要求 所以必须检验设计结果 如果无关区 过渡带 太窄 或截止频率太接近零频率和Fs 2时 设计结果可能不正确 rip表示f和m描述的各逼近频段允许的波纹幅度 幅频响应最大偏差 f的长度是rip的两倍 一般以 N fo mo w remezord f m rip Fs 返回的参数作为remez的调用参数 计算单位脉冲响应 hn remez N fo mo w 对比前面介绍的remez调用参数 可清楚地看出remezord返回参数N fo mo和w的含义 综上所述 调用remez和remezord函数设计线性相位FIR数字滤波器 关键是根据设计指标求出remezord函数的调用参数f m rip和Fs 其中Fs一般是题目给定的 或根据实际信号处理要求 按照采样定理 确定 下面给出由给定的各种滤波器设计指标确定remezord调用参数f m和rip的公式 编程时直接套用即可 2 滤波器设计指标1 低通滤波器设计指标 逼近通带 0 p 通带最大衰减 p dB 逼近阻带 s 阻带最小衰减 s dB remezord调用参数 7 4 6 其中 f向量省去了起点频率0和终点频率1 1和 2分别为通带和阻带波纹幅度 由式 7 4 4 和 7 4 5 计算得到 下面相同 2 高通滤波器设计指标 逼近通带 p 通带最大衰减 pdB 逼近阻带 0 s 阻带最小衰减 s dB remezord调用参数 7 4 7 3 带通滤波器设计指标 逼近通带 pl pu 通带最大衰减 pdB 逼近阻带 0 sl su 阻带最小衰减 sdB remezord调用参数 7 4 8 4 带阻滤波器设计指标 逼近阻带 sl su 阻带最大衰减 sdB 逼近通带 0 pl pu 通带最小衰减 p dB remezord调用参数 7 4 9 工程实际中常常给出对模拟信号的滤波指标要求 设计数字滤波器 对输入模拟信号采样后进行数字滤波 这时 调用参数f可以用模拟频率表示 但是 调用remezord时一定要加入采样频率参数Fs 这种情况的调用参数及调用格式见例7 4 3 例7 4 1 利用等波纹最佳逼近法重新设计FIR带阻通滤波器 指标与例7 2 3相同 即 逼近通带 0 0 2 0 8 通带最大衰减 p 1dB 逼近阻带 0 35 0 65 阻带最小衰减 s 60dB 解调用remezord和remez函数求解 由调用格式知道 首先要根据设计指标确定remezord函数的调用参数 再直接编写程序调用remezord和remez函数设计得到h n 将设计指标带入式 7 4 9 即可得到remezord函数的调用参数f m和rip 本例设计程序为ep741 m 即 ep741 m 例7 4 1用remez函数设计带阻滤波器 f 0 2 0 35 0 65 0 8 省略了0和1 m 1 0 1 rp 1 rs 60 由式 7 4 4 和式 7 4 5 求通带和阻带波纹幅度dat1 dat2和rip dat1 10 rp 20 1 10 rp 20 1 dat2 10 rs 20 rip dat1 dat2 dat1 M fo mo w remezord f m rip hn remez M fo mo w 以下绘图检验部分略去程序运行结果 M 28 即h n 的长度N 29 h n 及其损耗函数曲线如图7 4 4所示 例7 2 3中N 80 由此例可见 等波纹最佳逼近设计方法可以使滤波器阶数大大降低 图7 4 4调用remez函数设计的带阻FIR数字滤波器的h n 及损耗函数曲线 注意 设计结果应当是单位脉冲响应h n 的数据序列 但是一般N较大 列出h n 的全部数据序列所占篇幅太大 而且从一大堆数据中看不出h n 的变化规律 所以每道例题只给出其波形 读者运行程序就可以得到h n 的数据 例7 4 2 利用等波纹最佳逼近法设计FIR数字低通滤波器 实现对模拟信号的数字滤波处理 要求与例7 2 2相同 求出h n 并画出损耗函数曲线 解将例7 2 2中所给指标重和相频特性曲线写如下 通带截止频率 fp 1500Hz 通带最大衰减 p 1dB 阻带截止频率 fs 2500Hz 阻带最小衰减 s 40dB 对模拟信号的采样频率 Fs 10kHz 调用remezord和remez函数设计的程序为ep742 m 即 ep742 m 例7 4 2用remez函数设计低通滤波器 Fs 10000 对模拟信号采样的频率为1kHz f 1500 2500 边界频率为模拟频率 Hz m 1 0 rp 1 rs 40 dat1 10 rp 20 1 10 rp 20 1 dat2 10 rs 20 rip dat1 dat2 M fo mo w remezord f m rip Fs 边界频率为模拟频率 Hz 时必须加入采样频率Fs M M 1 估算的M值达不到要求 加1后满足要求 hn remez M fo mo w 以下绘图检验部分省略运行结果 滤波器阶数M 15 h n 及损耗函数曲线如图7 4 5所示 请读者注意 例7 2 2中用窗函数设计的滤波器阶数为23 图7 4 5调用remez函数设计的低通FIR数字滤波器的h n 及损耗函数曲线 7 5IIR和FIR数字滤波器的比较 IIR与FIR滤波器各有所长 在实际应用中应该考虑多方面因素加以选择 同时也要兼顾经济方面的要求以及计算工具的实际条件