【数字信号处理含matlab代码】第五篇：滤波器结构转换（一）——直接型与级联型的互转

2026/8/4 13:05:49

第五篇：滤波器结构转换（一）——直接型与级联型的互转

在滤波器设计（如fir1、ellip）中，MATLAB 默认返回的是直接型（Direct Form）的分子/分母系数b, a。这种结构简单直观，但对于高阶（如 > 20 阶）滤波器，直接型对系数量化误差极为敏感，极容易在极点靠近单位圆时产生不稳定。工程中更常用的是级联型（Cascade Form）——将高阶系统拆解为若干个二阶节（SOS）的串联，每个二阶节只处理一对共轭极点/零点，从而大幅降低数值敏感度。本篇将带领你深入dir2cas.m与cas2dir.m，探究零极点配对、增益归一化及二阶节合并的全部奥秘。

1. 为什么需要级联型？

一个 (N) 阶系统的传递函数为：

[
H(z) = \frac{b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}}{1 + a_1 z^{-1} + \cdots + a_N z^{-N}}
]

直接型（I 型或 II 型）将所有系数集中在一个差分方程中。当 (N) 较大时，(a_k) 和 (b_k) 的动态范围差异巨大，微小的量化误差可能导致极点移出单位圆。

级联型将 (H(z)) 分解为 (K) 个二阶节的乘积：

[
H(z) = b_0 \prod_{i=1}^{K} \frac{1 + B_{i1} z^{-1} + B_{i2} z^{-2}}{1 + A_{i1} z^{-1} + A_{i2} z^{-2}}
]

每个二阶节只处理一对共轭复数根（或两个实根）。
优点：极点/零点被“隔离”在各自的节中，舍入误差只影响局部，不会蔓延。
在 FPGA/DSP 硬件实现中，级联型是绝对主流。

我们的工具包提供了双向转换：

dir2cas：直接型 → 级联型（便于稳定实现）
cas2dir：级联型 → 直接型（便于理论分析或与教科书对照）

2. 深度剖析`dir2cas.m`

该函数接受直接型系数b, a，输出增益b0，以及 (K \times 3) 的矩阵B, A（每行对应一个二阶节的分子/分母系数）。

第一步：增益归一化与系数对齐

b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;

先将分子分母的首项系数提取出来，归一化为首项为 1 的多项式。
最后b0汇聚了总增益（分子首项除以分母首项）。例如 FIR 滤波器a=1时，b0 = b(1)。

M=length(b);N=length(a);ifN>M b=[bzeros(1,N-M)];elseifM>N a=[azeros(1,M-N)];N=M;end

强制使b和a等长，不足补零，以便后续求根时长度一致。

第二步：确定二阶节个数 (K)

K=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==N;b=[b0];a=[a0];end

(K = \lfloor N/2 \rfloor) 表示需要K个二阶节。
如果 (N) 为偶数（例如 4 阶），那正好组成 2 个二阶节。但为了后续配对方便，代码在末尾补了一个零，使得长度变为奇数（5），这样broots中一定会有 4 个根（偶数个），便于两两配对。这是一个非常巧妙的预处理技巧。

第三步：求根并配对

broots=cplxpair(roots(b));aroots=cplxpair(roots(a));

roots(b)求分子多项式的零点（Zeros），roots(a)求极点（Poles）。
cplxpair是核心：它将复数根按共轭配对排序，实根排在最后。保证前两个一定是一对共轭复数，便于组合成实系数的二阶节。

第四步：组合成二阶节

fori=1:2:2*K Brow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix((i+1)/2),:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix((i+1)/2),:)=Arow;end

从broots中每次取两个根，用poly还原成多项式系数（例如根为 (r_1, r_2)，则多项式为 ( (z - r_1)(z - r_2) = z^2 - (r_1+r_2)z + r_1r_2 )）。
real保证因为共轭配对，虚部恰好抵消，得到纯实数系数。
注意：poly返回的是 ( [1, -(r_1+r_2), r_1r_2] )，对应 ( 1 + B_1 z^{-1} + B_2 z^{-2} )。这里采用的正幂次多项式，但系数排列上正好对应filter(B,A,x)所需的顺序（因为 MATLAB 的filter默认按 (z^{-1}) 降幂）。

3. 深度剖析`cas2dir.m`

反向转换相对简单，就是将若干个二阶节的系数通过卷积乘起来。

function[b,a]=cas2dir(b0,B,A)[K,L]=size(B);b=[1];a=[1];fori=1:K b=conv(b,B(i,:));a=conv(a,A(i,:));endb=b*b0;

初始化为[1]（常数 1）。
遍历每个二阶节，将其系数与累积多项式进行卷积（即乘法）。
最后乘以总增益b0。
注意：conv的顺序无所谓，卷积满足交换律。

4. 实战验证：转换前后频率响应完全一致

我们来测试一个 6 阶椭圆低通滤波器，验证转换后是否与原滤波器等价的频率响应。

% 设计一个 6 阶椭圆低通[b,a]=ellip(6,0.5,40,0.4);% 转为级联型[b0,B,A]=dir2cas(b,a);% 再转回直接型（检验是否复原）[b2,a2]=cas2dir(b0,B,A);% 对比转换前后的频率响应[H1,w]=freqz(b,a,512);[H2,~]=freqz(b2,a2,512);% 绘制误差plot(w/pi,abs(H1-H2));title('转换前后幅频响应误差');grid on;

你会发现误差在 (10^{-15}) 量级，这是浮点舍入误差，说明转换完全可逆。

5. 特殊注意事项：文件`cas2dir(1).m`

在代码包中，文件名被命名为cas2dir(1).m，这是 Windows 系统下载时自动添加的副本标记。实际使用前，请务必将其重命名为cas2dir.m，否则 MATLAB 无法正常调用该函数。本系列博文中，我们统一称其为cas2dir.m。

6. 级联型排序的艺术（高阶话题）

dir2cas中的配对顺序是“按根在复平面上的位置依次配对”，但实际工程中，极点/零点的配对顺序会影响中间变量的动态范围。最佳实践是：

极点配对：将相距最近的极点配成一对（cplxpair默认按实部排序，往往也是距离最近的）。
零点配对：选择能够“抵消”该极点影响的零点（即距离最近的零点）进行组合，以减小每个节的峰值增益。
我们的基础版本没有实现优化配对，但对于教学和理解已经足够。如果你感兴趣，可以在dir2cas中增加基于距离的贪婪匹配算法。