IIR 滤波器从零到设计Ⅷ:IIR 滤波器工程实现,从数字零极点到可靠的 SOS

IIR 滤波器从零到设计Ⅷ:IIR 滤波器工程实现,从数字零极点到可靠的 SOS 前面的设计流程已经得到数字滤波器的零点、极点和整体增益理论上可以直接把所有零极点展开成一个高阶传递函数但在实际工程中尤其是高阶、窄带、低采样率或定点滤波器中一般不会直接使用高阶多项式而是组织成二阶节级联每个二阶节写成本文解决四个实际问题如何判断一个二阶节是否容易出现较大的内部响应零点和极点应当怎样配成 SOS多个 SOS 应当按照什么顺序级联整体增益应当怎样分配避免中间节点过大。一、为什么不直接使用高阶多项式将所有零极点展开成高阶多项式理论响应不会改变但系数可能对舍入误差非常敏感。特别是当多个极点靠得很近时微小的系数量化误差都可能明显改变极点位置严重时甚至会把极点推到单位圆外使滤波器不稳定。因此工程中更常采用零极点形式 ----- 二阶节配对 ----- SOS 排序 ----- 节间缩放这样并不是不用多项式而是只把零极点展开为多个二阶多项式不再展开成一个高阶多项式。二、Q 表示什么考虑模拟二阶系统其中表示这一对极点对应的固有频率Q 表示阻尼强弱和响应集中程度。若极点为则并且因此Q 大 -----小 ----- 极点靠近虚轴极点越靠近稳定边界自然响应衰减得越慢所以对于数字二阶节若极点为其分母为此时 r 越接近 1极点越靠近单位圆系统记忆越长。数字滤波器中不一定需要强行计算一个统一的数字 Q。工程实现时通常直接检查极点半径 r每个二阶节的峰值实际运行时的内部状态范围。所以Q 可以理解为一个风险提示高 Q 二阶节通常需要更谨慎地配对、排序和缩放。Q1: 怎么推导得到和Q已知共轭极点1. 根据极点写出分母二阶系统的分母为代入极点利用得到所以2. 与标准二阶形式比较为了方便描述二阶系统人们把分母统一写成将两个表达式放在一起对应项的系数必须相等。常数项相等所以再比较一次项系数整理得到Q2: 为什么模拟系统极点的实部决定自然响应的衰减速度假设一对模拟极点为对应的特征多项式是展开因此对应的齐次微分方程为为什么解与极点有关对于这种常系数微分方程假设解为它的导数为代入微分方程因为所以解这个方程正好得到所以齐次解是也就是把复指数变成正弦和余弦先把衰减项提出来根据欧拉公式对于实系数系统两个系数也会形成共轭关系因此最终可以整理成实数形式其中 B 和 C 由初始状态决定。利用正弦、余弦合成关系因此它在时域中对应的自然响应大致是可以把它拆成两部分理解振幅包络决定衰减速度决定振荡速度极点位置和衰减速度极点的实部是所以极点越靠近虚轴就说明越小而指数包络中的越小下降得越慢。例如比较两个极点对应的振幅包络分别为和在时所以第一个响应几乎已经完全消失第二个响应还剩下约36.8%。这就是“极点越靠近虚轴衰减越慢”的直接原因。三、零极点怎样组成二阶节1. 先将共轭极点配对实系数滤波器中的复极点一定成共轭对出现这一对极点组成实系数二阶分母展开得到若则2. 共轭零点同样配对若零点为则对应的二阶分子为这样便能保证 SOS 系数为实数。3. 实数根如何处理实数零点或极点可以两两组成二阶节。如果滤波器阶数为奇数最后可能剩下一个一阶节。工程中仍可将其放进 SOS 表中只需把缺少的二阶系数补零。四、零点对与极点对怎样匹配将共轭根分别组成零点对和极点对以后还需要决定哪一对零点应当和哪一对极点放在同一个 SOS 中不同配对方式的最终传递函数相同但单个二阶节的峰值可能完全不同。对于候选二阶节定义它的峰值为实际配对时可以遵循以下思路优先处理最靠近单位圆的极点对将它与不同的剩余零点对分别组合计算每一种组合的选择峰值较小的组合重复处理剩余零极点。直观上靠近该极点频率的零点通常能够抑制二阶节的局部峰值。但真正实现时最好直接计算频率响应峰值不要只根据零极点距离判断。这是一种贪心配对方法。它不一定得到全局最优结果但计算量远小于枚举所有组合而且在实际滤波器中通常已经足够有效。五、零极点配对和 SOS 排序不是一回事这两个步骤容易混淆。零极点配对决定的是每一个 SOS 内部放哪些零点和极点SOS 排序决定的是已经组成的 SOS 按什么顺序串联例如已经得到理论上所以顺序不会改变理想传递函数。但是实际计算中每经过一节都会产生一个中间信号因此不同顺序会产生不同的中间信号范围。六、为什么 SOS 顺序会影响实现假设级联顺序为经过前 m 节后的累计响应为对应的累计峰值为虽然最终输出响应固定但中间节点的最大增益会随着顺序变化。如果某个顺序使中间累计增益达到 30而另一个顺序只达到 2那么前一种顺序更容易在定点运算中溢出在浮点运算中损失有效精度产生很大的 DF-II Transposed 内部状态放大系数量化和舍入误差。七、SOS 实际怎样排序一种简单经验是低峰值、低风险的节放前面高 Q 节放后面但这不是绝对规则因为前后各节的响应可能在不同频率上互相补偿。更可靠的方法是计算累计峰值。假设当前累计响应为对于每个尚未选择的 SOS计算选择使最小的 SOS 作为下一节然后更新累计响应。这就是贪心排序如果共有 M 个二阶节贪心算法需要比较计算量约为而枚举全部顺序需要检查种排列。因此贪心算法减少的是相对于全排列搜索的计算量。它不能保证全局最优但很适合节数较多的工程滤波器。对于节数很少的滤波器可以枚举所有顺序对于节数较多的滤波器可以先贪心排序再尝试交换相邻 SOS 继续优化。八、为什么还需要增益分配假设完成配对和排序后其中是暂时没有分配整体增益的二阶节K 是滤波器的整体增益。如果直接把全部 K 放到第一节第一节输出可能过大如果全部放到最后一节最后一节的内部状态又可能过大。因此把整体增益拆成每个二阶节变为只要满足最终传递函数就不会改变。九、怎样计算节间缩放量假设希望每个中间节点的频率响应峰值不超过 L。L 是人为设定的目标范围。例如定点系统中可以留出一定余量而不是直接顶到满量程。第一节计算第一节未缩放时的峰值选择加入后第一节输出端的峰值便被调整到 L。第二节此时第一节已经带有第二节还没有分配。计算然后选择加入后前两节输出端的累计峰值就被调整到 L。后续各节同样地选择实际设计中一般只对前 M-1 个中间节点这样缩放。最后一节负责恢复整体增益于是最终传递函数保持不变。十、一个简单的缩放例子假设第一节未缩放峰值为因此加入后计算第一节与第二节的累计峰值所以此时前两节的累计增益系数为如果整体增益为并且一共只有三个 SOS则最后一节需要最终仍然满足这个例子同时说明缩放不能凭空消除增益只能重新分配增益。如果最后一节得到的增益过大就需要重新调整SOS 顺序中间节点目标值 L整体增益的分配方式零极点配对结果。然后重新检查各节的内部状态。十一、频率响应峰值不等于内部状态峰值前面的累计峰值方法控制的是每个 SOS 输出端的频率响应它非常适合进行第一轮排序和缩放但不能完全代表 DF-II Transposed 内部状态变量的最大值。原因是二阶节内部状态具有自己的传递函数。即使某一节的输出没有超过范围内部状态仍可能更大。因此最终工程验证还必须使用实际实现结构记录其中、是第 k 个 DF-II Transposed 二阶节的两个状态变量。如果内部状态仍然过大就继续调整配对、顺序或缩放。十二、频率峰值怎样实际计算对于任意累计响应在范围内建立频率网格逐点计算然后取最大值需要注意不能用各节峰值直接相乘不同 SOS 的峰值可能出现在不同频率高 Q 滤波器的峰值可能很窄需要使用更密的频率网格找到大致峰值后可以在附近继续细化搜索。十三、从数字零极点到 SOS 的完整工程流程整个过程可以整理为十四、总结从数学上看只要零点、极点和整体增益不变零极点配对、SOS 顺序和增益位置都不会改变理想传递函数。但在真实计算机中理想响应相同 ≠ 数值表现相同其中Q 和极点半径用于判断二阶节的数值风险零极点配对决定单个 SOS 是否容易产生大峰值SOS 排序决定级联中间节点的信号范围增益分配和缩放用于控制中间信号最终仍需在真实 DF-II Transposed 结构中检查内部状态。至此完整路线已经形成后续我会介绍下窗函数设计FIR滤波器同时使用代码实现上述滤波器。