📘 第 4 章：IIR 数字滤波器

这一篇在干嘛？

用反馈换取效率：IIR 滤波器的理论与系数计算、有限字长效应、快速 IIR 流水线、窄带与 lattice 结构设计。 原书代码为 VHDL，本篇所有代码已改写为 Verilog。

• 4.1 IIR 数字滤波器概述
• 例4.1 有损积分器I
• 4.2 IIR理论
• 4.3 IIR 系数的计算
• 重要 IIR 滤波器设计特性的总结
• 4.4 IIR 滤波器的实现
• • 速度
• 低：波形
• 例4.2 巴特沃思二阶系统
• 4.5 快速 IIR 滤波器
• 例4.3 有损积分器II
• 例4.4 群集方法
• 例4.5 分散预测方法
• 例4.6 有损积分器III
• BEGIN
• 4.6 窄带 IIR 滤波器
• 例4.7 窄带IIR滤波器的位宽要求
• BEGIN
• 例4.8 改进的并行双二阶IIR滤波器
• MATLAB 结果如下:
• 4.7 窄带 IIR 滤波器的全通滤波器设计
• 一种可行的实现：
• 4.8 练习

---

自测一下

IIR 为什么用比 FIR 少得多的阶数就能达标？

IIR 有反馈，极点在单位圆内可产生尖锐的谐振，传递函数在相同指标下阶数远低于 FIR；代价是相位非线性且必须警惕稳定性（极点必须都在单位圆内）。

有限字长效应在 IIR 里为什么比 FIR 严重？

IIR 有反馈回路，系数量化可能让极点移出单位圆（直接振荡），量化噪声还会被极点放大。FIR 无反馈，量化只影响频响精度，不涉及稳定性。

什么是"损耗积分器"变换？

把模拟原型中的积分器替换为带衰减因子的数字损耗积分器，得到可实现的数字 IIR；这是从模拟滤波器表格到数字实现的桥梁。

流水线 IIR 为什么不能像 FIR 那样随便插寄存器？

反馈回路里插寄存器会改变传递函数（极点位置变化）。必须先重分解（如状态空间分解/散射结构）再插流水，或选用 lattice 等天然适合流水的结构。

第4章 IIR 数字滤波器

4.1 IIR 数字滤波器概述

第 3 章我们学过 FIR 滤波器。在动手做设计之前，值得先把 FIR 的”长处”和”短板”都摆出来，这样才能理解为什么需要 IIR。

FIR 滤波器最有吸引力的地方在于：线性相位很容易实现；多带滤波器设计可行；用 Kaiser 窗函数法可以自由迭代设计；对于抽取器和插入器（见第 5 章），FIR 结构简单；非递归结构天然稳定，不会出现极限环；很容易得到高速流水线设计；系数和算法舍入误差的预算通常较低，量化噪声也有明确的定义。它的短板是：递归 FIR 滤波器由于极点/零点消除不彻底可能不稳定；设计极小极大滤波器必须借助复杂的 Parks-McClellan 算法；滤波器长度很长时实现代价大。

与 FIR 相比，IIR（Infinite Impulse Response，无限脉冲响应）滤波器引入了反馈机制，能够直接实现系统传递函数的极点和零点，而 FIR 是一种”全零”滤波器。因此在给定滤波器阶数时，IIR 往往能以更低的阶数达到同样的性能。IIR 设计的传统路线是把模拟滤波器的设计规范转换到数字域——这条路之所以合理，是因为模拟滤波器的设计技巧已经相当成熟，有大量标准表格可供查阅。本章将介绍四类最重要的模拟原型滤波器：巴特沃思滤波器、切比雪夫 I 型和 II 型滤波器以及椭圆滤波器。

IIR 的优点是：采用模拟原型滤波器的标准设计方法容易理解；低阶设计即可获得很好的滤波特性，并能以很高速度运行；设计时用查表甚至袖珍计算器就够了；在相同的容限裕度方案下，IIR 滤波器比 FIR 短得多；可以采用闭环设计算法。缺点是：相位响应通常是非线性的，想得到线性相位非常困难（用全通滤波器做相位补偿会使系统复杂性加倍）；整数实现时可能出现极限环；多带设计非常困难，一般只能设计低通、高通和带通滤波器；反馈会引入不稳定性（不过大多数情况下，把极点取为单位圆的镜像位置可以在保持同样幅值响应的同时得到稳定滤波器）；高速流水线设计比 FIR 难得多。

例 4.1 有损积分器 I

滤波器的一个基本功能是平滑噪声信号。假定信号 x[n] 中混有宽带零均值随机噪声，从数学上看，用积分器就能把噪声”平均掉”。如果输入信号的平均值在一个有限的时间间隔内能够保持住，就可以采用有损积分器来处理含噪信号。“有损”指的是反馈支路上乘了一个略小于 1 的系数 a，让历史输出逐渐”衰减”，否则积分会无界增长。

图 4-1 就是一个简单的一阶有损积分器，它满足离散时间差分方程（本例取 a = 3/4）：

$$y[n]=x[n]+ay[n-1]=x[n]+\frac{3}{4}y[n-1]\text{\hspace{1em}(4-1)}$$

图 4-1 用作有损积分器的一阶 IIR 滤波器

从图 4-2(a) 的脉冲响应可以看到，一个 15 抽头的 FIR 滤波器能够实现与一阶有损积分器同样的功能——也就是说，用一阶 IIR 就能顶替一个 15 抽头的 FIR，这正是 IIR 高效率的直观体现。有损积分器的阶跃响应如图 4-2(b) 所示。

(a) x[n] = 1000δ[n] 的脉冲响应

(b) x[n] = 100σ[n] 的阶跃响应
图 4-2 a=3/4 的有损积分器的仿真结果

下面给出该 IIR 滤波器的一种 FPGA 实现（原书 VHDL 改写为 Verilog）：

module iir (
    input  wire               clk,    // 系统时钟
    input  wire               reset,  // 异步复位
    input  wire signed [14:0] x_in,   // 系统输入
    output wire signed [14:0] y_out   // 输出结果
);
    // 用户自定义类型 S15：[-2^14, 2^14-1] 的 15 位有符号整数
    reg signed [14:0] x, y;
 
    // 输入与递归部分都用触发器实现
    always @(posedge clk or posedge reset) begin
        if (reset) begin            // 异步清零
            x <= 0;
            y <= 0;
        end else if (clk) begin     // 时钟上升沿
            x <= x_in;
            y <= x + (y >>> 2) + (y >>> 1);  // y = x + y/4 + y/2
        end
    end
 
    assign y_out = y;               // 把 y 连接到输出引脚
endmodule

原书在 PROCESS 模块内用寄存器实现递归部分，用 CSD 编码实现乘法和加法；本例中 y/4 与 y/2 只需算术右移即可完成。该设计使用了 62 个 LE，没有用到嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 $F_{max}=147.3\text{MHz}$。脉冲幅值为 1000 时的滤波器响应如图 4-3 所示，与图 4-2(a) 给出的 MATLAB 仿真结果相吻合。

图 4-3 有损积分器脉冲响应的 MODELSIM 仿真结果

另一种可选设计方案采用标准逻辑向量数据类型和 LPM_ADD_SUB 宏函数（练习 4.6 将讨论）。那种方法生成的代码更长，但优点是可以在位域上对符号扩展和乘法器进行直接控制。

4.2 IIR 理论

顾名思义，非递归滤波器（FIR）没有反馈，脉冲响应有限；递归滤波器（IIR）具有反馈，一般认为具有无限脉冲响应。图 4-4(a) 给出了递归部分和非递归部分相互独立的滤波器，把它们组合起来就得到图 4-4(b) 所示的规范（直接型）滤波器。图 4-4 中滤波器的传递函数可以写成：

$$F(z)=\frac{{\sum }_{l=0}^{L-1}b[l]z^{-l}}{1-{\sum }_{l=1}^{L-1}a[l]z^{-l}}\text{\hspace{1em}(4-2)}$$

对应的差分方程为：

$$y[n]=\sum _{l=0}^{L-1}b[l]x[n-l]+\sum _{l=1}^{L-1}a[l]y[n-l]\text{\hspace{1em}(4-3)}$$

图 4-4 带有反馈的滤波器

与 FIR 的差分方程 (3-2) 比较就会发现：递归系统不仅依赖输入 x[n] 的前 L 个值，还与输出 y[n] 的前 L-1 个值有关——这就是”反馈”的含义。

对 F(z) 求极点和零点：非递归部分（分子）给出零点 $p_{0l}$，分母给出极点 $p_{\infty l}$。利用极点/零点图可以直观地查出滤波器最重要的性质：在 z 域传递函数中用 $z=e^{j\omega T}$ 替换，就能用几何方法构造傅里叶传递函数：

$$F(\omega )=|F(\omega )|{\mathrm{e}}^{j\theta (\omega )}=\frac{{\prod }_{l=0}^{L-2}(p_{0l}-{\mathrm{e}}^{j\omega T})}{{\prod }_{l=0}^{L-2}(p_{\infty l}-{\mathrm{e}}^{j\omega T})}=\frac{\exp (\mathrm{j}{\sum }_l{\beta }_l){\prod }_{l=0}^{L-2}{v}_l}{\exp (\mathrm{j}{\sum }_l{\alpha }_l){\prod }_{l=0}^{L-2}{u}_l}\text{\hspace{1em}(4-4)}$$

图 4-5 给出了幅值（增益）与相位的几何解释：对于给定频率 ${\omega }_0$，增益是零点向量 $v_l$ 与极点向量 $u_l$ 的商。这些向量分别起始于各自的零点或极点，终止于感兴趣的频率点 $e^{j{\omega }_{0}T}$。图中示例的相位增益就是 $\theta ({\omega }_0)={\beta }_0+{\beta }_1-{\alpha }_0$。

图 4-5 用极点/零点分布图计算传递函数：幅值增益 = $u_0{u}_1/v_0$，相位增益 = ${\beta }_0+{\beta }_1-{\alpha }_0$

利用传递函数与极点/零点图之间的这种联系，可以推出以下几条性质：

1. 单位圆上的零点 $p_0=e^{j{\omega }_{0}T}$（无与之抵消的极点）生成傅里叶域传递函数在频率 ${\omega }_0$ 处的一个零点。
2. 单位圆上的极点 $p_{\infty }={\mathrm{e}}^{j{\omega }_{0}T}$（无与之抵消的零点）生成傅里叶域传递函数在频率 ${\omega }_0$ 处的无穷增益。
3. 所有极点都位于单位圆内部的稳定滤波器，可以接受任意类型的输入信号。
4. 实际滤波器的单个极点和零点位于实轴上，而复数极点和零点总是成对出现：如果 ${\alpha }_0+j{\alpha }_1$ 是一个极点或零点，则 ${\alpha }_0-j{\alpha }_1$ 也必定是。
5. 线性相位（恒时延）滤波器的所有极点和零点关于单位圆对称，或者位于 z = 0 处。

把性质 3) 与 5) 组合起来就会发现：对稳定的线性相位系统，所有零点必须关于单位圆对称，且极点只允许位于 z = 0 处。因此，IIR 滤波器（极点 $z\ne 0$）只能是近似线性相位。为实现近似，需要借用模拟滤波设计的一条著名准则：用具有单位增益、但引入非零相位增益的全通滤波器，在感兴趣的频率范围（通带）内实现相位线性化。

4.3 IIR 系数的计算

经典 IIR 滤波器设计中，先让数字滤波器的设计规范尽量贴近理想滤波器，再把理想数字滤波器模型规范在数学上转换成一组模拟滤波器模型规范，所用的工具就是双线性 z 变换：

$$s=\frac{z-1}{z+1}\text{\hspace{1em}(4-5)}$$

思路是：用双线性变换把数字规范映射为模拟规范，套用经典的模拟原型公式合成，再映射回数字 IIR 滤波器。下面逐一介绍四类原型。

巴特沃思滤波器的幅值平方频率响应为：

$$\mid F(\omega ){\mid }^2=\frac{1}{1+{\left(\frac{\omega }{{\omega }_s}\right)}^{2N}}\text{\hspace{1em}(4-6)}$$

$|F(\omega ){|}^2$ 的极点沿圆周分布，彼此相隔 $\pi /N$ 弧度；传递函数在 $\omega =0$ 处 N 次可微，说明它在 0 Hz 附近”局部最光滑”。图 4-6（上）给出了一个巴特沃思设计示例，其容差方案与 Kaiser 窗函数及图 3-7 的等波纹设计一致。

(a) 巴特沃思滤波器传递函数

(b) 巴特沃思滤波器通带组延迟

(c) 巴特沃思滤波器极点/零点分布图（×=极点，o=零点）

切比雪夫 I 型和 II 型滤波器按切比雪夫多项式 $V_N(\omega )=\cos (N\cos (\omega ))$ 定义，要求滤波器的极点驻留在一个椭圆上。I 型的幅值平方频率响应为：

$$\mid F(\omega ){\mid }^2=\frac{1}{1+{\epsilon }^2{V}_N^2\left(\frac{\omega }{{\omega }_s}\right)}\text{\hspace{1em}(4-7)}$$

典型 I 型滤波器的幅频与脉冲响应如图 4-7（上）所示：通带呈等波纹，阻带光滑。II 型的模型则是把波纹挪到阻带：

$$\mid F(\omega ){\mid }^2=\frac{1}{1+\left({\epsilon }^2{V}_N^2{\left(\frac{\omega }{{\omega }_s}\right)}^{-1}\right)}\text{\hspace{1em}(4-8)}$$

典型 II 型滤波器示例如图 4-7（下）所示：通带光滑，阻带呈波纹特性。

(a) 切比雪夫 I 型滤波器传递函数

(b) 切比雪夫 I 型滤波器通带组延迟

(c) 切比雪夫 I 型滤波器极点/零点分布图（×=极点，o=零点）

(a) 切比雪夫 II 型滤波器传递函数

(b) 切比雪夫 II 型滤波器通带组延迟

(c) 切比雪夫 II 型滤波器极点/零点分布图（×=极点，o=零点）
图 4-7 基于 MATLAB 工具箱的切比雪夫滤波器设计（上面是切比雪夫 I 型，下面是切比雪夫 II 型）

椭圆滤波器根据雅可比椭圆方程 $U_N(\omega )$ 的解定义，幅值平方频率响应的模型为：

$${\left|F(\omega )\right|}^2=\frac{1}{1+{\epsilon }^2{U}_N^2{\left(\frac{\omega }{{\omega }_s}\right)}^{-1}}\text{\hspace{1em}(4-9)}$$

典型椭圆滤波器的响应如图 4-6（下）所示，其通带和阻带都呈现波纹。

(a) 椭圆滤波器传递函数

(b) 椭圆滤波器通带组延迟

(c) 椭圆滤波器极点/零点分布图（×=极点，o=零点）
图 4-6 基于 MATLAB 工具箱的滤波器设计（上面是巴特沃思滤波器，下面是椭圆滤波器）

在图 3-8 相同的容差设计方案下比较这四类实现：巴特沃思滤波器需要 19 阶，切比雪夫滤波器 8 阶，椭圆滤波器只要 6 阶。同时观察图 4-6 和图 4-7 可以发现：波纹随阶数的减少而增加，群延迟也变得高度非线性。大多数情况下，切比雪夫 II 型是良好的折中——阶数适中、通带平坦、群延迟可以容忍。

重要 IIR 滤波器设计特性的总结

每种原型都为设计者提供了一种折中方案，其特性总结如下：

• 巴特沃思：最平的通带、平坦的阻带、宽的过渡带
• 切比雪夫 I 型：等波纹的通带、平坦的阻带、适中的过渡带
• 切比雪夫 II 型：平坦的通带、等波纹的阻带、适中的过渡带
• 椭圆型：等波纹的通带、等波纹的阻带、狭窄的过渡带

给定设计要求后，一般可以把握下列观察特征：

• 滤波器阶数：最低——椭圆型；中间——切比雪夫 I 型或 II 型；最高——巴特沃思。
• 通带特征：等波纹——椭圆型、切比雪夫 I 型；平坦——巴特沃思、切比雪夫 II 型。
• 阻带特征：等波纹——椭圆型、切比雪夫 II 型；平坦——巴特沃思、切比雪夫 I 型。
• 过渡带特征：最窄——椭圆型；中间——切比雪夫 I 型和 II 型；最宽——巴特沃思。

4.4 IIR 滤波器的实现

借助 MATLAB 之类的软件，获得 IIR 滤波器的传递函数只是简单的练习；真正的问题是在哪种体系结构下实现它。最重要的结构有：

• 直接 I 型（图 4-8）
• 直接 II 型（图 4-9）
• 一阶或二阶系统的级联（图 4-10(a)）
• 一阶或二阶系统的并联实现（图 4-10(b)）
• 基本级联或并联设计中典型二阶部分的双二阶实现（图 4-11）
• 正交形式：一阶或二阶状态变量系统的级联（图 4-10(a)）
• 并行正交：并联的一阶或二阶状态变量系统（图 4-10(b)）
• 连续的分数结构
• 网格滤波器（Gray-Markel 结构，图 4-12）
• 波形的数字实现（Fettweis 方法）
• 一般状态空间滤波器

图 4-8 采用乘法器模块的转置直接 I 型 IIR 滤波器

图 4-9 采用乘法器模块的转置直接 II 型 IIR 滤波器

图 4-10 F(z) 的级联和并联实现

图 4-11 传递函数 $F(z)=(b[0]+b[1]z^{-1}+b[2]z^{-2})/(a[0]+a[1]z^{-1}+a[2]z^{-2})$ 的双二阶实现：传递函数中的双二阶包含两个二次方程式

图 4-12 网格 IIR 滤波器

每种体系结构都有各自擅长的场合，通用的选择规律如下：

• 速度：高——直接 I 型和 II 型；低——波形结构。
• 定点算法舍入误差的灵敏度：高——直接 I 型和 II 型；低——正交、网格。
• 定点系数舍入误差的灵敏度：高——直接 I 型和 II 型；低——并联、波形。
• 特殊性质：正交加权输出——网格；最佳二阶部分——正交；任意 IIR 技术规范——状态变量。

借助 MATLAB 之类的软件工具，可以很容易地把系数从一种体系结构转换到另一种。下面举例说明。

例 4.2 巴特沃思二阶系统

假定要设计一种用二阶系统（双二阶）级联实现的巴特沃思滤波器（阶数 N=10，通带 Fp=0.3Fs），可以用下面的 MATLAB 代码生成系数：

N=10; Fp=0.3;
[B, A]=butter(N, Fp)
[sos, gain]=tf2sos(B,A)

先用 butter() 计算巴特沃思滤波器系数，再用传递函数到双二阶部分的转换函数 tf2sos() 把它拆成级联的二阶部分。得到的二阶部分系数如表 4-1 所示，总增益为 $4.961410\times 10^{-5}$。

表 4-1 二阶部分的系数

b[0,i]	b[1,i]	b[2,i]	a[0,i]	a[1,i]	a[2,i]
1.0000	2.1181	1.1220	1.0000	-0.6534	0.1117
1.0000	2.0703	1.0741	1.0000	-0.6831	0.1622
1.0000	1.9967	1.0004	1.0000	-0.7478	0.2722
1.0000	1.9277	0.9312	1.0000	-0.8598	0.4628
1.0000	1.8872	0.8907	1.0000	-1.0435	0.7753

图 4-13 给出了该滤波器的传递函数、群延迟和极点/零点图。注意所有零点都位于 $z_{0i}=-1$ 附近——这一点从分子系数也能看出：理想情况下应有 $b[1,i]=2$、$b[0,i]=b[2,i]=1$（对应 $(1+z^{-1}{)}^2$），表中出现的正是这种理想值的舍入误差。

(a) 幅值

(b) 群延迟响应

(c) 极点/零点分布图（×=极点，o=零点）
图 4-13 10 阶巴特沃思滤波器的仿真结果

4.4.1 有限字长效应

Crochiere 和 Oppenheim 已经证明：数字滤波器所需的系数字长与系数灵敏度密切相关，因此同一个 IIR 滤波器在不同结构下所需的字长范围差别很大。以他们分析过的 8 阶椭圆滤波器为例，把该传递函数分别用波形、级联、并联、网格、直接 I 型和 II 型以及连续分数体系结构实现，表 4-2 第 2 列给出了满足具体最大通带误差准则所预计的系数字长的保守估计。可以看出：直接形式需要的字长明显大于波形或并联结构；而就位宽 W 与乘法器数 M 的乘积而言，波形结构给出了最佳的复杂度 M×W（见第 6 列）。

表 4-2 Crochiere 和 Oppenheim 给出的 8 阶椭圆滤波器的数据（根据成本乘积 M×W 分类）

类型	字长 W	乘积 M	加法	延迟	成本 M×W
波形	11.35	12	31	10	136
级联	11.33	13	16	8	147
并联	10.12	18	16	8	182
网格	13.97	17	32	8	238
直接 I 型	20.86	16	16	16	334
直接 II 型	20.86	16	16	8	334
连续分数结构	22.61	18	16	8	408

再从 FIR 滤波器（第 3 章）的角度看：为了简化包含若干乘法器的模块设计，需要引入简化加法器图（RAG）技术。Dempster 和 Macleod 从 RAG 乘法器实现策略的角度重新计算了该 8 阶椭圆滤波器，结果见表 4-3。第 2 列是乘法器模块的规模：直接 II 型需要两个规模分别为 9 和 7 的模块；波形结构因为没有两个系数共用输入，无法设计乘法器模块，需要实现 11 个单独乘法器。第 3 列是规范有符号位（CSD）设计需要的加法器/减法器数量 B，第 4 列是单个优化的乘法器加法器图（MAG）的相同结果，第 5 列是 RAG 的结果，第 6 列是 RAG 设计的整体加法器/字宽乘积。表 4-3 表明：采用 RAG 算法时，乘法器模块的规模是一条基本准则，级联和并联形式可以得到与波形数字滤波器相当或更好的结果。对 FPGA 设计而言延迟不是问题，因为所有逻辑单元都自带触发器。

表 4-3 采用 CSD、MAG 和 RAG 策略实现 8 阶椭圆滤波器的数据

类型	模块规模	CSD B	MAG B	RAG B	W(B+A)
级联	4×3、2×1	26	26	24	453
并联	11×9、4×2、1×1	31	30	29	455
波形	11×1	58	63	22	602
网格	1×9、8×1	33	31	29	852
直接 I 型	1×16	103	83	36	1085
直接 II 型	1×9、1×7	103	83	41	1189
连续分数结构	18×1	118	117	88	2351

4.4.2 滤波器增益系数的优化

一般情况下，IIR 整数系数是从浮点系数导出的：先把最大系数规格化为 1，再乘以增益因数（通常是位宽 $2^{\mathrm{round}(W)}$）。但大多数情况下，在 $2^{\lfloor W\rfloor }\dots 2^{\lceil W\rceil }$ 范围内选择增益系数更为有效，而且传递函数基本不发生变化——因为在乘以增益系数之后系数就被舍入了。例如，在 $2^{\lfloor 11.33\rfloor }\dots 2^{\lceil 11.33\rceil }$ 范围内寻找 Crochiere-Oppenheim 设计示例中级联滤波器（表 4-3 给出的增益 $2^{11.33-1}=1024$）的最优系数，结果如表 4-4 所示。

表 4-4 最小化级联滤波器复杂性的增益因子的变化

不同增益	CSD	MAG	RAG
最优增益	1122	1121	1121
最优增益的加法器	23	21	18
增益=1024 的加法器	26	26	24
提高	12%	19%	25%

比较可见，实现乘法器所需加法器的数量有实质性改进。本例中 MAG 和 RAG 的最佳增益系数恰好相同，但两者并不必然一致。

4.5 快速 IIR 滤波器

第 3 章里，FIR 滤波器的时序性能几乎可以”零成本”地靠流水线技术提高。但流水线 IIR 滤波器就复杂了：仅为反馈回路上的加法器插入流水线寄存器，就相当于移动极点位置，从而改变了传递函数。不过，文献中已有若干在不改变传递函数的前提下提高吞吐量的策略，包括：

• 在时域中预测交叉
• 成群（群集）的预测极点/零点设置
• 分散预测极点/零点的设置
• IIR 抽取滤波器的设计
• 并行处理
• RNS 实现

前 5 种方法基于滤波器体系结构或信号流技术，最后一种建立在计算机算法的基础上（第 2 章）。为了简化每个示例的代码，下面都只考虑一阶 IIR 滤波器，但同样的思路对高阶 IIR 滤波器也适用。

4.5.1 时域交叉

先研究一阶 IIR 系统的差分方程：

$$y[n+1]=ay[n]+bx[n]\text{\hspace{1em}(4-10)}$$

输出 $y[n+1]$ 依赖上一步的 $y[n]$，反馈环路限制了速度。预测算法的做法是：把 $y[n+1]$ 代入 $y[n+2]$ 的差分方程，“向前多算一步”：

$$y[n+2]=ay[n+1]+bx[n+1]=a^{2}y[n]+abx[n]+bx[n+1]\text{\hspace{1em}(4-11)}$$

相应的系统如图 4-14 所示：反馈环路上现在有两个延迟，乘法系数变成了 $a^2$。

图 4-14 采用预测算法的有损积分器

把预测推广到 $(S-1)$ 步，结果为：

$$y[n+S]=a^{S}y[n]+\underset{(\eta )}{\underbrace{\sum _{k=0}^{S-1}{a}^{k}bx[n+S-1-k]}}\text{\hspace{1em}(4-12)}$$

从式 (4-12) 可以看到：表达式 $(\eta )$ 定义了一个系数为 $\{b,ab,a^{2}b,\dots ,a^{S-1}b\}$ 的 FIR 滤波器，可以用第 3 章的流水线技术（流水线乘法器和流水线加法器树）实现；递归部分则可以用系数为 $a^S$ 的 S 级流水线乘法器实现。下面用一个示例说明预测设计算法。

例 4.3 有损积分器 II

回到例 4.1 的有损积分器（a = 3/4），这次加上两步预测：

$$\begin{array}{rl}y[n+2]& =\frac{3}{4}y[n+1]+x[n+1]=\frac{3}{4}\left(\frac{3}{4}y[n]+x[n]\right)+x[n+1]\\ & =\frac{9}{16}y[n]+\frac{3}{4}x[n]+x[n+1]\\ y[n]& =\frac{9}{16}y[n-2]+\frac{3}{4}x[n-2]+x[n-1]\end{array}\text{\hspace{1em}(4-13)}$$

预测结构的有损积分器由一个非递归部分（对 x 的 FIR 滤波器）和一个延迟为 2、系数为 9/16 的递归部分构成。实现该预测形式 IIR 滤波器的代码如下（原书 VHDL 改写为 Verilog）：

module iir_pipe (
    input  wire               clk,    // 系统时钟
    input  wire               reset,  // 异步复位
    input  wire signed [14:0] x_in,   // 系统输入
    output wire signed [14:0] y_out   // 系统输出
);
    reg signed [14:0] x, x3, sx, y, y9;
 
    // 输入、输出及各级流水线均用触发器实现
    always @(posedge clk or posedge reset) begin
        if (reset) begin            // 异步清零
            x  <= 0; x3 <= 0; sx <= 0;
            y9 <= 0; y  <= 0;
        end else begin
            x  <= x_in;
            x3 <= (x >>> 1) + (x >>> 2);   // 计算 x*3/4
            sx <= x + x3;                  // x 各项之和 = FIR 部分
            y9 <= (y >>> 1) + (y >>> 4);   // 计算 y*9/16
            y  <= sx + y9;                 // 计算输出
        end
    end
 
    assign y_out = y;               // 把寄存器 y 连接到输出引脚
endmodule

这个示例的流水线乘法器和加法器分两步实现：第一步计算 $\frac{9}{16}y[n]$，第二步计算 $\frac{3}{4}x[n]+x[n+1]$，再与 $\frac{9}{16}y[n]$ 相加。该设计使用了 123 个 LE，未使用嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 $F_{\max }=215.05\text{MHz}$。滤波器对幅值为 1000 的脉冲的响应如图 4-15 所示。

图 4-15 预测算法有损积分器的脉冲响应的 VHDL 仿真

与例 4.1 的 62 个 LE、147.3 MHz 方案相比：预测算法的流水线需要更多资源，但速度提升约 30%；对比图 4-3 和图 4-15 可以看到，预测设计方案有额外的总延迟；两种方法的量化效果差值为 ±2。练习 4.7 讨论了一种采用标准逻辑向量数据类型和 LPM_ADD_SUB 宏函数的可选设计方案：代码更长，但可以在位域上对符号扩展和乘法器进行直接控制。

4.5.2 群集和分散预测的流水线技术

群集和分散预测方案通过为设计添加”自消除”的极点/零点对，来简化滤波器递归部分的流水线化。群集方法在传递函数的分母中引入额外的极点/零点，使得 $z^{-1}$、$z^{-2}$、…、$z^{-(S-1)}$ 的系数为 0。

例 4.4 群集方法

假定二阶传递函数的两个极点分别位于 1/2 和 3/4 处：

$$F(z)=\frac{1}{1-1.25z^{-1}+0.375z^{-2}}=\frac{1}{(1-0.5z^{-1})(1-0.75z^{-1})}\text{\hspace{1em}(4-16)}$$

在 $z=-1.25$ 处增加一个相互抵消的极点/零点对，得到新的传递函数：

$$F(z)=\frac{1+1.25z^{-1}}{1-1.1875z^{-2}+0.4688z^{-3}}\text{\hspace{1em}(4-17)}$$

这样分母中 $z^{-1}$ 项消失了，滤波器的递归部分就可以用额外的流水线级来实现。但群集方法有个隐患：被抵消的极点/零点对可能位于单位圆之外，如果抵消不完全，就会引入不稳定性。一般地，极点在 $r_1$、$r_2$ 处并拥有一个额外抵消对的二阶系统，必然有一个极点位于 $-(r_1+r_2)$，而 $|r_1+r_2|>1$，即位于单位圆外部。Soderstrand 等人给出了一种稳定的群集方法：引入多个抵消极点/零点对。

分散预测方法则没有稳定性问题。对于在 $p$ 点有极点的原始滤波器，引入 $(S-1)$ 个位于 $z_k=p{e}^{j\pi k/S}$ 处的极点/零点抵消对，结果是在传递函数的分母中只有 $z^0$、$z^S$、$z^{-2S}$ 等项非零。

例 4.5 分散预测方法

仍考虑极点位于 $z_{\infty 1}=0.5$ 和 $z_{\infty 2}=0.75$ 的二阶系统，要增加两个额外的流水线级。先注意一个一般结论：位于 $p$ 和 $p^{\ast }$ 的极点/零点对导出

$$(1-p{z}^{-1})(1-p^{\ast }{z}^{-1})=1-(p+p^{\ast })z^{-1}+r{r}^{\ast }{z}^{-2}$$

特别是当 $p=r\times \exp (j2\pi /3)$ 时：

$$(1-p{z}^{-1})(1-p^{\ast }{z}^{-1})=1-2r\cos (2\pi /3)z^{-1}+r^2{z}^{-2}=1+r{z}^{-1}+r^2{z}^{-2}$$

也就是说，在 $0.5e^{\pm j2\pi /3}$ 和 $0.75e^{\pm j2\pi /3}$ 处添加相互抵消的极点/零点对后：

$$\begin{array}{rl}F(z)& =\frac{1}{1-1.25z^{-1}+0.375z^{-2}}\times \frac{(1+0.5z^{-1}+0.25z^{-2})(1+0.75z^{-1}+0.5625z^{-2})}{(1+0.5z^{-1}+0.25z^{-2})(1+0.75z^{-1}+0.5625z^{-2})}\\ & =\frac{1+1.25z^{-1}+1.1875z^{-2}+0.4687z^{-3}+0.1406z^{-4}}{1-0.5469z^{-3}+0.0527z^{-6}}\\ & =\frac{512+640z^{-1}+608z^{-2}+240z^{-3}+72z^{-4}}{512-280z^{-3}+27z^{-6}}\end{array}$$

分母中只剩 $z^{-3}$ 与 $z^{-6}$ 项，递归部分就可以用两个额外的流水线级来实现。

值得注意的是：对于一阶 IIR 系统，群集方法和分散预测方法引入的是同样的相互抵消的极点/零点对——位置处于以原点为圆心的圆上，相邻角度相差 $2\pi /S$。此时非递归部分可以按”2 的幂分解”形式实现：

$$(1+a{z}^{-1})(1+a^2{z}^{-2})(1+a^4{z}^{-4})\dots \text{\hspace{1em}(4-19)}$$

图 4-16 给出了一阶部分极点/零点对的表示方法，其中递归部分可以通过 4 个流水线级实现。

(a) $F_1(z)=(1+\alpha z^{-1})$

(b) $F_2(z)=1+a^2{z}^{-2}$

(c) $F_3(z)=1/(1-a^4{z}^{-4})$

(d) $F(z)={\prod }_k{F}_k(z)=(1+a{z}^{-1})(1+a^2{z}^{-2})/(1-a^4{z}^{-4})=1/(1-a{z}^{-1})$
图 4-16 分散预测的一阶 IIR 滤波器的极点/零点图

4.5.3 IIR 抽取器设计

在抽取滤波器的基础上，Martinez 和 Parks 提出了一种基于极小极大方法的滤波器设计算法（参阅第 5 章）。导出的传递函数满足：

$$F(z)=\frac{{\sum }_{l=0}^{L}b[l]z^{-l}}{1-{\sum }_{n=0}^{N/S}a[n]z^{-nS}}\text{\hspace{1em}(4-20)}$$

也就是说，分母中只有每隔 S 个的系数非零，因此递归部分（分母）天然可以采用 S 级流水线。观察最终的极点/零点分布图（图 4-17(b)）：所有零点都位于单位圆上（与常见的椭圆滤波器一样），而极点位于一个圆上，主轴彼此相差 $2\pi /S$ 的角度。

(a) S=5 的 37 阶 Martinez-Parks IIR 滤波器传递函数

(b) 极点/零点分布图
图 4-17 S=5 的 37 阶 Martinez-Parks IIR 滤波器的传递函数及极点/零点分布图

4.5.4 并行处理

并行处理实现中，形成 P 条并联的 IIR 通路，每条通路都以 1/P 的输入采样速率运行，最后在输出端通过多路复用器组合在一起，如图 4-18 所示。这种方法速度快的原因有二：一般情况下多路复用器比乘法器和/或加法器快；而且每条通路都有更多时间（P 倍）来计算自己负责的输出。

图 4-18 并联 IIR 滤波器的实现，抽头延迟线（TDL）以 1/P 的输入采样速率运行

为了便于说明，考虑 P=2 的一阶系统，预测结构与式 (4-11) 相同：

$$y[n+2]=ay[n+1]+x[n+1]=a^{2}y[n]+ax[n]+x[n+1]\text{\hspace{1em}(4-21)}$$

把输出序列按偶数 n=2k 和奇数 n=2k-1 分成两路，得到：

$$y[n+2]=\{\begin{array}{l}y[2k+2]=a^{2}y[2k]+ax[2k]+x[2k+1]\\ y[2k+1]=a^{2}y[2k-1]+ax[2k-1]+x[2k]\end{array}\text{\hspace{1em}(4-22)}$$

其中 $n,k\in Z$。这两个方程就是并联 IIR 滤波器 FPGA 实现的基础。

例 4.6 有损积分器 III

作为例 4.1 和例 4.3 的扩展，考虑 a = 3/4 的并联有损积分器。如图 4-19 所示，双通道并联有损积分器是两个非递归部分（对 x 的 FIR 滤波器）和两个延迟为 2、系数为 9/16 的递归部分的组合。设计代码如下（原书 VHDL 改写为 Verilog）：

module iir_par (
    input  wire               clk,    // 系统时钟
    input  wire               reset,  // 异步复位
    input  wire signed [14:0] x_in,   // 系统输入
    output wire signed [14:0] x_e,    // 偶路输入（测试观察）
    output wire signed [14:0] x_o,    // 奇路输入（测试观察）
    output wire signed [14:0] y_e,    // 偶路输出（测试观察）
    output wire signed [14:0] y_o,    // 奇路输出（测试观察）
    output wire               clk2,   // 二分频时钟
    output wire signed [14:0] y_out   // 系统输出
);
    localparam EVEN = 1'b0, ODD = 1'b1;
    reg              state;         // 状态机：even / odd
    reg signed [14:0] x_even, xd_even, sum_x_even, y_even;
    reg signed [14:0] x_odd,  xd_odd,  sum_x_odd,  y_odd;
    reg signed [14:0] x_wait, y_wait, y;
    reg               clk_div2;
 
    // 复用进程：把 x 拆分成偶/奇样点，并以 clk 速率重组输出 y
    always @(posedge clk or posedge reset) begin
        if (reset) begin            // 异步复位
            state    <= EVEN;
            x_even   <= 0; x_odd  <= 0;
            y        <= 0;
            x_wait   <= 0; y_wait <= 0;
            clk_div2 <= 1'b0;
        end else begin
            if (state == EVEN) begin
                x_even   <= x_in;
                x_odd    <= x_wait;
                clk_div2 <= 1'b1;
                y        <= y_wait;
                state    <= ODD;
            end else begin
                x_wait   <= x_in;
                y        <= y_odd;
                y_wait   <= y_even;
                clk_div2 <= 1'b0;
                state    <= EVEN;
            end
        end
    end
 
    assign y_out = y;
    assign clk2  = clk_div2;
    assign x_e   = x_even;          // 输出一些额外的测试信号
    assign x_o   = x_odd;
    assign y_e   = y_even;
    assign y_o   = y_odd;
 
    // 算术进程：两条通路各自按式(4-22)运算，用 clk_div2 作时钟
    always @(negedge clk_div2 or posedge reset) begin
        if (reset) begin            // 异步清零
            xd_even    <= 0; sum_x_even <= 0; y_even <= 0;
            xd_odd     <= 0; sum_x_odd  <= 0; y_odd  <= 0;
        end else begin
            xd_even    <= x_even;
            sum_x_even <= (((xd_even <<< 1) + xd_even) >>> 2) + x_odd;
            y_even     <= (((y_even  <<< 3) + y_even ) >>> 4) + sum_x_even;
            xd_odd     <= x_odd;
            sum_x_odd  <= (((xd_odd  <<< 1) + xd_odd ) >>> 2) + xd_even;
            y_odd      <= (((y_odd   <<< 3) + y_odd  ) >>> 4) + sum_x_odd;
        end
    end
endmodule

设计用两个 always 块实现：第一个（复用）把 x 拆分成偶数和奇数索引部分，输出 y 以 clk 速率重组，同时还生成了速率降为 clk/2 的第二时钟信号；第二个模块按式 (4-22) 实现滤波器算法（其中 $(8y+y)/16=9y/16$、$(2x+x)/4=3x/4$ 用移位与加法完成）。该设计使用了 236 个 LE，没有用到嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 $F_{max}=479.39\text{MHz}$。图 4-20 给出了仿真结果。

图 4-19 双通道并联 IIR 滤波器的实现

图 4-20 并联 IIR 滤波器对脉冲 1000 的响应的 VHDL 仿真

与前面几种方法相比，并联实现的缺点是实现成本较高——本设计用了 236 个 LE（速度则达到 479 MHz）。

4.5.5 采用 RNS 的 IIR 设计

余数系统（RNS）使用固有的短字长进行并行运算，因此是实现快速（递归）IIR 滤波器的优秀候选方案。典型的 IIR-RNS 设计中，系统被实现为递归系统和非递归系统的集合，每个系统都按 FIR 结构定义（图 4-21）。每个 FIR 都可以在 RNS-DA 中用四分之一平方乘法器实现，也可以在第 2 章开发的索引域中实现。

图 4-21 采用两个 FIR 部分和缩放的 IIR 滤波器的 RNS 实现

对于稳定滤波器，递归部分必须做缩放运算，用来控制动态范围的增长率。缩放可以用混合基数转换、中国余数定理（CRT）或 $ϵ$-CRT 方法实现。对于高速设计，推荐在群集或分散预测流水线技术的基础上增加一个额外的流水线延迟。5.3 节将详细研究 RNS 递归滤波器的设计，可以看到 RNS 设计能把速度提升 40% 以上。

4.6 窄带 IIR 滤波器

关于不同 IIR 架构所需的（小数）位宽，文献中有很多研究。4.4.1 节讨论的 Crochiere 和 Oppenheim 的结果后来又被 Dempster 和 Macleod 基于简化加法器图技术重新审视过。结论是：二阶部分的串联或并联配置优势明显，波形数字滤波器（WDF）和网格滤波器分列第三、四位。

Crochiere 和 Oppenheim 研究的基础是带通滤波器，但很多实际 IIR 设计更倾向于低通滤波器。由于假定滤波器（单位）增益为 1，不同滤波器结构对全系统整数位的要求是另一个不常被提及的重要因素。

此外，对窄带 IIR 滤波器的最新研究表明：由于振幅的增益较小，并联结构的全通滤波器往往被采用。最初针对 WDF 设计的准则也可以用于二阶类型的全通滤波器，或与网格滤波器结合使用。这样，最具潜力的结构包括：

• 转置的直接 I 型方案
• 直接型二阶部分级联或并联方案
• 每部分拥有 1 或 2 个乘法器的网格滤波器
• 波形数字滤波器梯型滤波器
• 全通转置的直接 I 型方案
• 全通直接型二阶部分级联或并联方案
• 每部分拥有 1 或 2 个乘法器的全通网格滤波器
• 全通波形数字滤波器，又叫作网格波形数字滤波器

为了更好地理解这些设计选择，下面从一个典型的窄带低通滤波器设计入手。

4.6.1 窄带设计示例

在相同设计指标下，如果 FIR 滤波器太长、需要过多硬件资源并引入过长延迟，就应使用 IIR 滤波器。本例的指标是：采样频率 $60\times 128=7680\text{Hz}$，通带 040 Hz，阻带 403840 Hz，通带波纹 1 dB，阻带波纹 40~50 dB。50 dB 的阻带衰减稍高于 8 位系统所需的 $6\times 8=48\text{dB}$，以便为量化噪声留出空间。电力系统中常用这种滤波器来计算谐波畸变。它的过采样率很高，而通带和阻带适中，大约只需 8 位精度。

椭圆（或考尔）滤波器阶数最低，可作起点。它是一个 5 阶滤波器，但从极点-零点分布图（图 4-22(e)）可以推断，需要较高的精度才能匹配传递函数的误差方案。先用 MATLAB 的 ellip 函数计算滤波器系数：

$$[B,A]=\text{ ellip }(5,1,50,40/7680\ast 2)$$

参数依次是滤波器阶数、通带波纹（dB）、阻带波纹（dB）以及截止频率，采样频率归一化到 2。前向滤波器系数置于数组 A 中，反馈系数置于数组 B 中。通过 MATLAB 对直接型滤波器实现的仿真结果如下。

图 4-22 给出了传递函数、完整的极点/零点分布图以及忽略 z = -1 处零点的”局部放大”极点分布图。作为椭圆滤波器的典型特征，所有零点都位于单位圆上。

(a) 5 阶椭圆滤波器传递函数（幅值）

(b) 极点/零点分布图

(c) 通带

(d) 极点/零点分布图（放大）

(e) 放大的极点/零点分布图
图 4-22 5 阶椭圆滤波器：(a-c) 幅值，(d-e) 极点/零点分布图

由于滤波器频带很窄，脉冲响应非常长。要知道任意滤波器的最大增益，可以把特征频率作为输入信号：对于特征频率为 0 的低通滤波器（即直流），输入是阶跃函数，输出就是阶跃响应。图 4-23 给出了脉冲响应和阶跃响应。

(a) 5 阶椭圆滤波器脉冲响应

(b) 阶跃响应
图 4-23 5 阶椭圆滤波器

利用极点/零点值可以选择结构，并在 SIMULINK 中对滤波器仿真以确保运行正确。在 SIMULINK 中可以清楚地计算所需资源、确定最长路径，并增加额外的流水线寄存器来提升速度。这些额外的流水线寄存器会增加总延迟，但只要其他极点/零点的位置不变，传递函数就不会发生变化。

要在定点算法中实现滤波器，需要确定小数部分、每个加法器的整数位以及系数乘法器所需的位数。要把基于嵌入式乘法器和使用 LE 的设计区分开：采用 LE 时，可以使用三值（-1、0、1）CSD 编码，并把乘法器系数整合到精简加法器图（RAG-n）中。两种情况下，都需要先估计所需的小数位宽和整数位宽，然后（在基于 LE 实现时）用 CSD 或 RAG 编码精心调整系数。

估计小数部分精度的方法是：将系数量化为 N 位小数——乘以 $2^N$、四舍五入、再除以 $2^N$；然后用量化后的系数 $(A_q,B_q)$ 计算滤波器传递函数，检查是否达到误差方案要求。达到要求则结束；否则令 N 加 1（$N\to N+1$）再次尝试。

整数精度的估计则可以不用循环：只需确定滤波器的”特征频率”并将该频率的信号加到滤波器输入上。对于窄带低通滤波器，合理的假设是特征频率 $f=0$，即用阶跃函数作为输入；接下来对每一个加法器的最大幅值执行 ${\log }_2$ 运算，由此确定每一个节点所需的整数位数。

窄带 IIR 滤波器的 FPGA 实现结构与位宽设计

本单元围绕同一个 5 阶椭圆窄带 IIR 滤波器，比较三种实现结构——直接型、级联二阶（双二阶）型、并联二阶型，外加一种思路完全不同的网格（Lattice）滤波器。核心矛盾只有一个：窄带意味着极点离单位圆极近，反馈回路增益极大，稍不留神位宽就爆炸。读完你就能回答三个问题：每种结构需要多少位？为什么差这么多？工程上怎么选？

例4.7 窄带 IIR 滤波器的位宽要求（直接型的困境）

为什么窄带滤波器对位宽这么挑剔

先说”是什么”：这个 5 阶椭圆滤波器的通带极窄，反馈系数 A 中出现了 -4.97、+9.87、-9.82、+4.88、-0.97 这样一组数。它们不是随便来的——极点几乎贴着单位圆，导致分母多项式系数在数值上相互抵消得很厉害。说”为什么重要”：系数量化稍差一点，极点位置就跑偏，阻带衰减直接不达标；说”不用会怎样”：如果按普通滤波器的习惯给个 10~12 位小数，量化后的滤波器可能根本不是你要的那个滤波器。

小数位：30 位是怎么算出来的

方法很直接：把系数逐次增减小数位并量化，画出量化后的传递函数（幅频、相频、群延迟、衰减），与理想传递函数比较，直到误差满足设计规范为止。

图 4-24(a) 量化系数后直接型滤波器的幅频响应

图 4-24(b) 量化系数后直接型滤波器的相频响应

图 4-24(c) 量化系数后直接型滤波器的群延迟

图 4-24(d) 量化系数后直接型滤波器的衰减特性

图 4-24(e) 量化系数后的极点/零点分布

图 4-24 直接型滤波器量化系数的传递函数

由图 4-24 的结果可以看到：小数精度必须达到 30 位，阻带衰减才勉强满足误差方案。这比一般 IIR 滤波器高出一个数量级，根源就是前面说的极点贴圆问题。

整数位：33 位又是怎么来的

小数位解决”精度”，整数位解决”溢出”。做法：在 SIMULINK（或 MATLAB）里搭建滤波器，输入阶跃信号，测量每一个内部加法器输出的最大幅值。

图 4-25 SIMULINK 中转置形式的直接 I 型设计

图 4-26(a) 阶跃响应测量结果

图 4-26(b) 反馈增益的测量结果

(c) 当前结构中所有 10 个加法器幅值信号的 ${\log }_2$ 值 图 4-26 运算结果

图 4-26(c) 把每个加法器输出的最大幅值做了 ${\log }_2$ 运算，最坏情况的增益达到约 $6\times 10^9$，因此内部整数精度需要

$$\left|{\log }_2(6\times 10^9)\right|=\lceil 32.48\rceil =33\text{ 位}$$

这个数字触目惊心：一个 5 阶滤波器的内部中间值竟然涨到 33 位整数——反馈回路里信号被放大了近 60 亿倍。这正是窄带直接型的典型症状。

64 位数据通路：原书 VHDL 改写为 Verilog

小数 30 位 + 整数 33 位 + 1 位符号，数据通路需要

$$1+30+33=64\text{ 位}$$

64 位远超 VHDL INTEGER 类型的 32 位上限，原书因此借助 VHDL-2008 的 sfixed 定点类型来实现。下面的 Verilog 版本用”有符号 64 位定点数（Q34.30，即 34 位整数 30 位小数）“表达同样的设计思想：所有系数在编译期量化为 $2^{30}$ 的整数倍，每次算术运算后做饱和/截位处理，防止位增长。

// 5 阶 IIR 直接型实现（原书 VHDL 改写为 Verilog）
// 前馈系数 B = 0.000304 -0.000909 0.000605 0.000605 -0.000909 0.000304
// 反馈系数 A = 1.000000 -4.968203 9.874754 -9.815007 4.878564 -0.970108
module iir5sfix (
    input  wire         clk,      // 系统时钟
    input  wire         reset,    // 系统复位
    input  wire         swtch,    // 反馈开关（调试用）
    input  wire signed [63:0] x_in, // 系统输入，Q34.30
    output wire signed [39:0] t_out, // 反馈信号，Q24.16
    output wire signed [39:0] y_out  // 系统输出，Q24.16
);
    // 系数统一量化为 Q34.30（数值 = 真实系数 * 2^30 取整）
    localparam signed [63:0] A2 = -64'sd5334566906; // -4.9682025852
    localparam signed [63:0] A3 =  64'sd10602936011; //  9.8747536754
    localparam signed [63:0] A4 = -64'sd10538783414; // -9.8150069021
    localparam signed [63:0] A5 =  64'sd5238318091;  //  4.8785639415
    localparam signed [63:0] A6 = -64'sd1041645681;  // -0.9701081227
    localparam signed [63:0] B1 =  64'sd325960;      //  0.0003035737
    localparam signed [63:0] B2 = -64'sd975522;      // -0.0009085259
    localparam signed [63:0] B3 =  64'sd649566;      //  0.0006049556
    localparam signed [63:0] B4 =  64'sd649566;      //  0.0006049556
    localparam signed [63:0] B5 = -64'sd975522;      // -0.0009085259
    localparam signed [63:0] B6 =  64'sd325960;      //  0.0003035737
 
    reg signed [63:0] h, s1, s2, s3, s4, s5;
    reg signed [63:0] x, y, t, r2, r3, r4, r5;
 
    // 饱和处理函数：把 128 位乘加结果饱和到 Q34.30
    function signed [63:0] sat64;
        input signed [127:0] v;
        begin
            if (v > 128'sd9223372036854775807)      sat64 = 64'sd9223372036854775807;
            else if (v < -128'sd9223372036854775808) sat64 = -64'sd9223372036854775808;
            else sat64 = v[63:0];
        end
    endfunction
 
    // 输入重定为 Q34.30 有符号数
    always @(posedge clk) begin
        x <= x_in;
        // 反馈开关：断开时 h=x（可逐个观察脉冲响应中的反馈系数），闭合时 h=x-t
        h <= (swtch == 1'b0) ? x : sat64({{64{x[63]}}, x} - {{64{t[63]}}, t});
 
        if (reset) begin  // 异步复位全部寄存器
            t <= 0; y <= 0; r2 <= 0; r3 <= 0; r4 <= 0; r5 <= 0;
            s1 <= 0; s2 <= 0; s3 <= 0; s4 <= 0; s5 <= 0;
        end else begin    // IIR 直接型：反馈 5 项 + 前馈 6 项，逐级累加
            r5 <= sat64(                    A6 * h);
            r4 <= sat64({{64{r5[63]}}, r5} + A5 * h);
            r3 <= sat64({{64{r4[63]}}, r4} + A4 * h);
            r2 <= sat64({{64{r3[63]}}, r3} + A3 * h);
            t  <= sat64({{64{r2[63]}}, r2} + A2 * h);
            s5 <= sat64(                    B6 * h);
            s4 <= sat64({{64{s5[63]}}, s5} + B5 * h);
            s3 <= sat64({{64{s4[63]}}, s4} + B4 * h);
            s2 <= sat64({{64{s3[63]}}, s3} + B3 * h);
            s1 <= sat64({{64{s2[63]}}, s2} + B2 * h);
            y  <= sat64({{64{s1[63]}}, s1} + B1 * h);
        end
    end
 
    // Q34.30 截位为 Q24.16 输出：保留小数位 bit29..bit14
    assign t_out = t[43:14];
    assign y_out = y[43:14];
endmodule

这段代码里有几个初学者容易忽略的细节：

• 反馈开关 switch 的作用是调试。开关断开时，反馈回路被切断，t_out 输出的脉冲响应每次会显现一个反馈系数，方便逐项核对。
• 每个中间寄存器对应一次乘加，反馈支路 r5→r4→r3→r2→t 与前馈支路 s5→s4→s3→s2→s1→y 是两条独立的累加链。
• 每次运算后必须处理位宽。Verilog 版本用饱和函数代替原书的 resize（fixed_wrap, fixed_truncate）参数，目的相同：不让乘加结果的位宽无限制增长。

实现结果：该设计占用 2474 个 LE、128 个嵌入式乘法器，TimeQuest 缓慢 85C 模型下 $F_{max}=46.99\text{MHz}$。对于 34.30 定点格式（幅值 $1.0=40000000_{16}$），脉冲响应的仿真如图 4-27，与 MATLAB 结果一致。100 ns 系统时钟下，60 μs 的仿真对应 600 个时钟周期，最大脉冲出现在约第 133 个周期处。

图 4-27 直接型 5 阶椭圆滤波器脉冲响应的 MODELSIM 仿真结果

不做位宽处理会付出什么代价

对比三个实验非常有说服力：如果不指定 wrap/truncate 之类的位宽处理参数，综合器只能保守地扩展位宽，设计变大（6279 个 LE）且变慢（29.94 MHz）；如果干脆改用 64 位浮点，代码确实更简单（不用再做截位），但规模膨胀到 42682 个 LE、180 个乘法器，速度掉到 6.5 MHz。可见定点 + 显式位宽管理才是 FPGA 上性价比最高的路线。

此外，基于 LE（查找表）实现乘法器时还有一些流行的系数优化手段，可单独或组合使用：

1. 符号 2 的幂分解：把乘法器的二进制系数拆成若干个带符号的 2 的幂项之和，先每系数用 2 项，逐步增加直到满足误差指标。
2. 系数微调 + 非线性优化：先取整，再对每个系数做 ±1 的扰动，看对传递函数是否有利。可用整数线性规划（ILP）、遗传算法（GA）等实现，MATLAB 优化工具箱的 fminimax 就能做这件事。
3. 简化加法器图（RAG）：只有当多个系数共享相同输入时才有效。注意 IIR 的反馈回路里不能插入额外寄存器，所以要用加法器深度最小的 RAG；且 RAG 中比例因子固定，不必再精心调整系数。

要提醒的是：这些非线性优化不能保证找到全局最优，与结构选择带来的收益相比往往只是锦上添花——这就引出了下一节：换结构。

4.6.2 级联二阶系统窄带滤波器设计

为什么级联结构是首选

把 5 阶滤波器拆成若干一阶/二阶节（双二阶）级联，有三个直接的好处：

• 每个节内极点/零点可以就近配对，各节增益可控，内部信号不会像直接型那样爆炸，量化噪声也小；
• 椭圆滤波器零点在单位圆上，二阶零点部分形如 $1+b(k,1)z^{-1}+z^{-2}$，首尾系数都是 1，能省乘法器；
• 滤波器阶数已经是最短的（椭圆滤波器的优势），级联不会引入额外阶数。

极点/零点怎么配对：Jackson 规则

配对与排序有条理可循：距离单位圆最近的极点增益最大、波纹最大，应先与最近的零点配对以压低增益；排序时最靠近单位圆的节放最后，即先实现离单位圆最远的极点。对照极点/零点分布图，实轴上的一阶极点离单位圆最远，所以级联顺序是：一阶节 → 二阶节 → 二阶节。MATLAB 的 tf2sos 函数自动给出这种配对：

$$[\text{sos},\text{gain}]=\text{tf2sos(B, A)}$$

对 5 阶椭圆滤波器运行结果如下（每行一个二阶节，b 为前馈、a 为反馈系数）：

b[i,1]	b[i,2]	b[i,3]	a[i,1]	a[i,2]	a[i,3]
1.0000	1.0000	0	1.0000	-0.9887	0
1.0000	-1.9950	1.0000	1.0000	-0.9847	0.9852
1.0000	-1.9977	1.0000	1.0000	-0.9948	0.9959

注意 tf2sos 只给出一个整体增益 0.000304。为避免某一节内量化噪声过大，要把这个增益按比例分配到三节：让每节输出的幅值归一化（低通情形可按阶跃响应输出归一化），得到 $s(1)=0.0057$、$s(2)=0.1147$、$s(3)=0.4677$，且 $s(1)s(2)s(3)=0.000304$，与整体增益一致。以直接 II 型双二阶加独立比例因子搭建，得到图 4-28 的结构。

图 4-28 基于直接 II 型双二阶部分和单独比例的 5 阶椭圆滤波器的实现

级联结构的位宽收益

小数位仍用老办法：量化所有系数、比较传递函数误差。

图 4-29(a) 二阶节系数量化后的幅频响应

图 4-29(b) 二阶节系数量化后的相频响应

图 4-29(c) 二阶节系数量化后的群延迟

图 4-29(d) 二阶节系数量化后的衰减特性

图 4-29(e) 二阶节量化后的极点/零点分布

图 4-29(f) 不同量化位数的对比

图 4-29 二阶部分设计中量化系数的传递函数和极点零点图

由图 4-29：只需 14 位小数就满足误差方案，而直接型要 30 位——直接省下 16 位。整数位同样用 SIMULINK 测阶跃响应：

图 4-30(a) 级联设计的 SIMULINK 阶跃响应

图 4-30(b) 反馈增益测量结果

(c) 图 4-30(c) 各加法器最大幅值的 ${\log }_2$ 结果

最坏情况下需要 $\lceil {\log }_2(1700)\rceil =\lceil 10.7\rceil =11$ 位内部整数精度（对比直接型的 33 位）。总位宽从 64 位降到约 $1+14+11=26$ 位量级，已经落回普通整数类型能处理的范围。

改进型双二阶：先零点、后极点

还可以再抠出几位。改进的核心思想：让零点部分先于极点部分实现，先做前馈再做反馈，反馈支路里的信号幅值会更小。对图 4-31 的结构，记反馈部分输出寄存器为 $w[n]$，在 z 域：

$$W(z)=X(z)z^{-1}(b[1]+b[2]z^{-1}+b[3]z^{-2})\text{\hspace{1em}(4-23)}$$

借助 $W(z)$ 得到输出：

$$Y(z)=W(z)z^{-2}-Y(z)a[2]z^{-1}-Y(z)a[3]z^{-2}=\frac{W(z)z^{-2}}{1+a[2]z^{-1}+a[3]z^{-2}}\text{\hspace{1em}(4-24)}$$

于是改进双二阶的传递函数为：

$$F(z)=\frac{Y(z)}{X(z)}=\frac{b[1]+b[2]z^{-1}+b[3]z^{-2}}{1+a[2]z^{-1}+a[3]z^{-2}}{z}^{-3}\text{\hspace{1em}(4-25)}$$

也就是说，它与常规双二阶完全等价，只是多了 3 个采样周期的延迟——用延迟换低增益，对多数非实时严格的场景非常划算。

图 4-31 降低整数增益的改进的双二阶系统

另一个附带的好处是：最坏情况下，信号路径从一个”乘法器 + 加法器”缩短为”加法器 → 乘法器 → 加法器”中的更短组合，关键路径变短。重新测量整数位宽：

图 4-32(a) 改进双二阶的阶跃响应测量

图 4-32(b) 反馈增益测量结果

图 4-32 改进的双二阶的二阶设计下整数位宽的计算结果

由图 4-32 可见整数位宽进一步缩短。注意代价：(a) 不能再对所有双二阶系数统一做单一 RAG 优化；(b) 系统多了额外延迟。

4.6.3 并联二阶系统窄带滤波器设计

并联结构与部分分式展开

并联实现把传递函数做部分分式展开，各二阶节并行工作、输出相加。优点是不用操心各节的先后次序（并联没有顺序），且各节误差互不传播；缺点是部分分式展开会让系数变得更”碎”，直接形式中等于 0 或 ±1 的漂亮系数消失了。

在 MATLAB 中用 residue 函数做展开：

$$[R,P,D]=\text{residue}(B,A)$$

向量 R 存留数（分子项），P 存极点，标量 D 存直接项，即：

$$F(z)=D+\sum _{k=1}^N\frac{R(k)}{z-P(k)}\text{\hspace{1em}(4-26)}$$

对 5 阶椭圆滤波器得到：

$$\begin{array}{rl}& \mathrm{R}(1,2)=0.0015\pm j0.0007\\ & \mathrm{R}(3,4)=-0.0073\pm j0.0026\\ & \mathrm{R}(5)=0.0122\\ & \mathrm{P}(1,2)=0.9974\pm 0.0325\\ & \mathrm{P}(3,4)=0.9923\pm 0.0226\\ & \mathrm{P}(5)=0.9887\end{array}$$

直接项 $D=0.00030357$。FPGA 里不想处理复数，所以把共轭复极点对合并成实系数二阶节。对第一对：

$$\frac{R(1)}{z-P(1)}+\frac{R(2)}{z-P(2)}=\frac{(R(1)+R(2))z+R(1)P(2)+R(2)P(1)}{z^2+2z\Re (P(1))+|P(1){|}^2}\text{\hspace{1em}(4-27)}$$

合并工作同样可以交给 residue（输入三个参数即从留数/极点重建传递函数）：

[B1, A1] = residue([R(1) R(2)], [P(1) P(2)], 0)

得到第一、二个二阶节的系数：

$$\begin{array}{rl}B1& =0.0031-0.00320\\ A1& =1.0000-1.99480.9959\end{array}\begin{array}{rl}B2& =-0.01460.01460\\ A2& =1.0000-1.99470.9952\end{array}$$

一个值得注意的细节：共轭对合并后，分子退化为一阶多项式（常数 × 一次项），不再是与分母同阶的二阶——实现零点时省下 1 个乘法器和 1 个加法器。

并联结构的位宽计算

用 SIMULINK 搭好模型后，老三样：量化小数位、测阶跃响应定整数位。

图 4-33(a) 17 位量化下的传递函数

图 4-33(b) 17 位量化的局部放大

图 4-33(c) 18 位量化下的传递函数

图 4-33(d) 18 位量化的局部放大

图 4-33 并行二阶部分设计下不同量化系数的传递函数图 图 4-33(e) 量化对比汇总

图 4-33(f) 误差方案的判定

(c) 图 4-34 并行双二阶的二阶系统设计的整数位宽的计算结果

结论：17 位量化时阻带不达标，18 位小数即满足误差方案；标准并联双二阶实现的最大整数增益约 11 位。与级联型的 14 位相比略差一点小数位，但如前所述换来了”不用排序”的自由。

改进型并联结构

图 4-31 的”先零点后极点”技巧同样适用于并联各节，可以再省几位。改进后各节的整数位宽计算结果如图 4-35。

图 4-35(a) 改进并联结构的阶跃响应

图 4-35(b) 反馈增益测量结果

图 4-35 改进的并行双二阶的二阶设计的整数位宽的计算结果

由图 4-35 可以看出一个惊喜：使用改进双二阶后，内部不需要整数位（信号始终小于 1）。但延迟要补：由式 (4-25)，每个二阶节引入 3 个采样周期延迟，一阶节引入 1 个采样周期延迟。各支路延迟不同会导致输出求和时信号没对齐——补救办法是给短延迟支路补上相应的延迟线，并可通过测量 4 个子系统的脉冲响应与标准双二阶对比来验证。完整系统如图 4-36。

图 4-36 基于 SIMULINK 的改进的并行双二阶部分的设计

例4.8 改进的并行双二阶 IIR 滤波器

系数与精度选择

采用 residue 算出的系数。内部定点格式选 2.19（2 位整数、19 位小数）。定小数位的经验法则：每个十进制小数数字约值 3.32 位二进制，19 位需要 $\lceil 19/3.32\rceil =6$ 位十进制数字。注意 MATLAB 默认只显示 4 位，要用 sprintf 之类把系数的完整精度打印出来再量化——直接复制界面上的 4 位数会引入可观误差。

下面是并行双二阶 5 阶 IIR 的 Verilog 实现（原书 VHDL 改写为 Verilog）。输入输出为 16.16 格式，内部为 2.19 格式，系数为 1.18 格式。

// 5 阶 IIR 并行（改进双二阶）实现（原书 VHDL 改写为 Verilog）
// 系数：D = 0.00030357
//       B1 = 0.0031  -0.0032  0     A1 = 1.0000 -1.9948 0.9959
//       B2 = -0.0146  0.0146  0     A2 = 1.0000 -1.9847 0.9852
//       B3 = 0.0122                A3 = 0.9887（一阶节）
module iir5para (
    input  wire         clk,
    input  wire         reset,
    input  wire signed [31:0] x_in,    // 输入，Q16.16
    output wire signed [31:0] y_Dout,  // 直接（0 阶）部分
    output wire signed [31:0] y_1out,  // 一阶部分
    output wire signed [31:0] y_21out, // 二阶部分 1
    output wire signed [31:0] y_22out, // 二阶部分 2
    output wire signed [31:0] y_out    // 系统输出
);
    // 系数量化为 Q1.18（数值 = 真实系数 * 2^18 取整）
    localparam signed [19:0] A12 = -20'sd522937; // -1.99484680
    localparam signed [19:0] A13 =  20'sd261072; //  0.99591112
    localparam signed [19:0] B11 =  20'sd805;    //  0.00307256
    localparam signed [19:0] B12 = -20'sd829;    // -0.00316061
    localparam signed [19:0] A22 = -20'sd520271; // -1.98467605
    localparam signed [19:0] A23 =  20'sd258276; //  0.98524428
    localparam signed [19:0] B21 = -20'sd3838;   // -0.01464265
    localparam signed [19:0] B22 =  20'sd3839;   //  0.01464684
    localparam signed [19:0] A32 =  20'sd259175; //  0.98867974（一阶节）
    localparam signed [19:0] B31 =  20'sd3190;   //  0.012170
    localparam signed [19:0] D   =  20'sd80;     //  0.000304（直接部分）
 
    reg signed [19:0] x, y;
    reg signed [19:0] s11, s12, s21, s22, s31;
    reg signed [19:0] r12, r13, r22, r23, r32, r41, r42, r43;
 
    always @(posedge clk) begin
        if (reset) begin   // 异步复位所有寄存器
            y <= 0; r12 <= 0; r13 <= 0; r22 <= 0; r23 <= 0; r32 <= 0;
            r41 <= 0; r42 <= 0; r43 <= 0;
            s11 <= 0; s12 <= 0; s21 <= 0; s22 <= 0; s31 <= 0;
        end else begin
            x <= x_in >>> 3;            // Q16.16 -> Q2.19
            // 第 1 节：二阶（先零点后极点的改进双二阶）
            s12 <= (B12 * x) >>> 0;
            s11 <= s12 + (B11 * x);
            r13 <= s11 - (A13 * r12);
            r12 <= r13 - (A12 * r12);
            // 第 2 节：二阶
            s22 <= (B22 * x);
            s21 <= s22 + (B21 * x);
            r23 <= s21 - (A23 * r22);
            r22 <= r23 - (A22 * r22);
            // 第 3 节：一阶
            s31 <= (B31 * x);
            r32 <= s31 + (A32 * r32);
            // 第 4 节：常数（直接部分）
            r41 <= (D * x);
            // 输出加法器树
            r42 <= r41;
            r43 <= r42 + r32;
            y   <= r12 + r22 + r43;
        end
    end
 
    // Q2.19 -> Q16.16 输出（算术右移 3 位，符号自动扩展）
    assign y_out   = y   >>> 3;
    assign y_Dout  = r42 >>> 3;
    assign y_1out  = r32 >>> 3;
    assign y_21out = r22 >>> 3;
    assign y_22out = r12 >>> 3;
endmodule

解读几个要点：

• 端口上除了总输出，还引出 4 个并行子系统的输出（直接部分、一阶部分、两个二阶部分），调试时可以逐节观察各自的贡献。
• 每个二阶节只有 4 次乘法、3 次加法：因为合并共轭对后分子是一阶的（$b[i,3]=0$），比标准双二阶省了一次乘加。
• 一阶节 r32 <= s31 + A32*r32 是标准的单极点递推；直接部分 r41 <= D*x 只是一个常数乘法。

实现结果：624 个 LE、51 个嵌入式乘法器，$F_{max}=87.69\text{MHz}$（TimeQuest 缓慢 85C 模型）。对于 16.16 格式（幅值 $1.0=65536_{10}$），阶跃响应如图 4-37，与 MATLAB 仿真一致。100 ns 时钟下 60 μs 仿真对应 600 个周期，最大阶跃响应出现在约第 215 个周期处（各支路补齐延迟后的对齐时刻）。

图 4-37 基于改进的双二阶的 5 阶并行 IIR 滤波器阶跃响应的 MODELSIM 仿真结果

与直接型对比，结论一目了然：并行实现速度翻倍，LE 只需四分之一，乘法器减半。这是结构选择的胜利，不是什么优化技巧的胜利。

4.6.4 窄带 IIR 滤波器的网格滤波器设计

网格滤波器是什么、为什么要学它

直接型、级联型、并联型之间可以方便地互相转换（差分方程、极点/零点图、脉冲响应一一对应），而网格（Lattice）滤波器与它们的映射要绕一些路。即便如此它仍值得学，因为有两个独特的优点：

• 系数量化敏感性低；
• 稳定性一眼可验：所有网格系数幅值小于 1 即稳定。这在自适应滤波等系数经常变化的场合极其方便。

第一步：从极点求网格系数

设计通常分两步：先由极点算出网格反馈系数 $k(n)$，再把不同中间输出组合成前馈系数。以三阶为例，观察图 4-38(a) 的网格单元，每个单元的级联矩阵为：

$$\left[\begin{array}{c}{X}_{n+1}\\ Y_{n+1}\end{array}\right]=\left[\begin{array}{cc}1& k(n)z^{-1}\\ k(n)& z^{-1}\end{array}\right]\left[\begin{array}{c}{Y}_n\\ X_n\end{array}\right]\text{\hspace{1em}(4-29)}$$

图 4-38 网格滤波器（(a) 单元信号流，(b) 三阶级联）

三阶就是把三个单元级联（$X_1=Y_1$）。全系统传函取决于输出取哪一点：取 $Y_1$ 得全极点滤波器 $H(z)=Y_1/X_4$，取 $Y_4$ 得全通滤波器 $H(z)=Y_4/X_4$。用 Symbolic 工具箱可以验证全极点传函：

syms k1 k2 k3 z
disp('Two-pair matrix description:')
A1=[1 k1*z^-1;k1 z^-1];
A2=[1 k2*z^-1;k2 z^-1];
A3=[1 k3*z^-1;k3 z^-1];
disp('Third order transfer X(z)/Y(x):')
P=A3*A2*A1
S=P(1,1)+P(1,2);
collect(S,z)

脚本先定义符号，创建三个单元矩阵并算级联乘积 P。因为 $X_1=Y_1$，$P(1,1)+P(1,2)$ 之和是 $X_4/Y_1$，其倒数正是 $H(z)=Y_1/X_4$。运行结果（用 $z^{-n}$ 表示）：

$$H(z)=\frac{1}{1+(k_3{k}_2+k_1+k_2{k}_1)z^{-1}+(k_3{k}_2{k}_1+k_2+k_3{k}_1)z^{-2}+k_3{z}^{-3}}\text{\hspace{1em}(4-30)}$$

用 freqz 把它转成标准全极点多项式，可得网格系数与分母系数的关系：

$$\begin{array}{rl}{a}_4& =k_3\\ a_3& =k_3{k}_2{k}_1+k_2+k_3{k}_1\\ a_2& =k_3{k}_2+k_1+k_2{k}_1\end{array}$$

再算 $\mathrm{P}(2,1)+\mathrm{P}(2,2)$，得到一个互反多项式（分子分母互为镜像，$\mathrm{Num}(z^{-1})=\mathrm{Den}(z)z^L$），极点/零点关于单位圆镜像对称——这正是全通滤波器的特征。

第二步：Gray-Markel 组合出完整 IIR

全极点/全通本身有用（语音处理等），但通用 IIR 还需把各中间抽头加权求和：

$$Y(z)=\sum _{m=1}^{M+1}v(m)Y_m(z)X(z)\text{\hspace{1em}(4-31)}$$

各输出 $Y_m(z)$ 与权重 $v(m)$ 可递归计算。这套流程已封装在 MATLAB 函数 $[k,v]=\text{tf2latc}$ 中，生成网格参数 k 与前馈系数 v。对前面那个 5 阶椭圆滤波器：

k(1) = -0.999690100
k(2) =  0.999875949
k(3) = -0.999815189
k(4) =  0.999660865
k(5) = -0.970108307
v(1) = 0.000000004445
v(2) = 0.000000009318
v(3) = 0.000003644655
v(4) = 0.000005053596
v(5) = 0.000599705718
v(6) = 0.000303581790

从参数就能读出滤波器性格：所有 $|k|\approx 1$ 说明极点贴着单位圆（稳定，但勉强）；v 系数小到 $10^{-9}$ 量级，预示着输出是深度凹陷的窄带振幅。完整 5 阶结构如图 4-39。

图 4-39 5 阶 IIR 网格滤波器

网格滤波器的位宽

小数位：同样量化所有系数测传递函数。

图 4-40 5 阶网格滤波器量化系数的传递函数和极点零点图

由图 4-40：13 位不够，14 位小数才满足误差方案。比预期的多，原因在于反馈系数的跨度太大：$\max (v)/\min (v)=1.3\times 10^5$，即 17.1 位——最小和最大系数相差五个数量级，小数位必须照顾到最小那个系数。

整数位：SIMULINK 测阶跃响应。

图 4-41(a) 网格滤波器阶跃响应（全通输出良好）

图 4-41(b) 第一级输出（全极点输出）的增益测量

(c) 图 4-41(c) 各加法器最大幅值的 ${\log }_2$ 结果

全通输出和最终输出都很正常，但最坏情况出现在第一级的全极点输出：需要 $\lceil {\log }_2(2\times 10^8)\rceil =\lceil 27.6\rceil =28$ 位整数精度。合计 $14+28=42$ 位——比级联/并联差很多，但比直接型（64 位）好。

单乘法器改进结构与增益调节

Gray-Markel 提出的单乘法器结构可以改善整数增益，见图 4-42。它还有一个额外好处：每节只需一个乘法器，且两种版本可以在节间调节增益。先推级联矩阵：左上角加法器输出为 $X_2-X_1{z}^{-1}$，两个输出信号为：

$$Y_1=X_2+k_1(X_2-X_1{z}^{-1})\text{\hspace{1em}(4-32)}$$ $$Y_2=X_1{z}^{-1}+k_1(X_2-X_1{z}^{-1})\text{\hspace{1em}(4-33)}$$

把第二个索引变量移到左侧、第一个移到右侧，由式 (4-32) 解出：

$$X_2=\frac{1}{1+k_1}(Y_1+k_1{z}^{-1}{X}_1)\text{\hspace{1em}(4-34)}$$

代入式 (4-33)：

$$Y_2=z^{-1}{X}_1(1-k_1)+k_1\frac{1}{1+k_1}(Y_1+k_1{z}^{-1}{X}_1)=\frac{1}{1+k_1}(k_1{Y}_1+z^{-1}{X}_1)\text{\hspace{1em}(4-35)}$$

于是得到生成矩阵：

$$\left[\begin{array}{c}{X}_{n+1}\\ Y_{n+1}\end{array}\right]=\frac{1}{1+k(n)}\left[\begin{array}{cc}1& k(n)z^{-1}\\ k(n)& z^{-1}\end{array}\right]\left[\begin{array}{c}{Y}_n\\ X_n\end{array}\right]\text{\hspace{1em}(4-36)}$$

对比式 (4-29)：唯一区别是比例因子 $1/(1+k(n))$。若对第二种乘法器类型（图 4-42(b)）重复推导，比例因子变成 $1/(1-k(n))$。这样我们就掌握了控制每节增益的旋钮：要压低增益就令 $1/|1\pm k(n)|>1$，要抬升就令其小于 1。这是个折中游戏——增益不能太大（否则位宽又涨回去），也不能太小（否则量化噪声相对变大）。

(c) 图 4-42 网格滤波器（(a) 单乘法器结构之一，(b) 第二种乘法器类型）

图 4-42 改进网格滤波器结构

图 4-43(a) $e(n)=[-1,1,-1,1,-1]$ 配置下的阶跃响应

图 4-43(b) 全极点信号显著增大

图 4-43 基于相同操作 $e(n)=\text{sign}(k(n))=[-1,1,-1,1,-1]$（与因子 k 的符号匹配）的改进网格滤波器整数位宽的计算结果

暴力法是把 5 阶滤波器的全部 $2^5=32$ 种配置都测一遍；聪明法是用 Gray-Markel 算法，它保证在增益不超过上限 1 的前提下使增益尽可能大。对本滤波器所有 $k(n)\approx \pm 1$，若取 $e(n)=\text{sign}(k(n))$，各节因子 $|1\pm k|$ 趋近 0，增益暴涨——图 4-43(b) 里全极点信号变得很大，位宽反而比双乘法器结构涨得更多。

另一个极端是取 $e(n)=-\text{sign}(k(n))$，所有节的因子都大于 2，增益全小于 1：

图 4-44(a) $e(n)=[1,-1,1,-1,1]$ 配置下的阶跃响应

图 4-44(b) 全极点信号变得很小

图 4-44 基于操作 $e(n)=-\text{sign}(k(n))=[1,-1,1,-1,1]$（与因子 k 的符号相反）的改进网格滤波器整数位宽的计算结果

此时只需要小数位。而 Gray-Markel 原始方法取 $e(n)=1$，增益大于/小于 1 交替出现，各信号（包括全极点输出）都被”合理”地调整：

图 4-45(a) $e(n)=[1,1,1,1,1]$ 配置下的阶跃响应

(b) 图 4-45(b) 全极点输出被合理调整

(c) 图 4-45 基于 $e(n)=[1,1,1,1,1]$ 符号操作的改进网格滤波器整数位宽的计算结果

每节增益变了，输出权重也要相应调整——本质上是比例因子的连乘 $\prod (1\pm k(m))$：

$$\begin{array}{rl}u(m)& =v(m)/\pi (m)\\ \pi (m)& =\{\begin{array}{ll}1& m=M+1\\ {\prod }_{n=m}^M(1+e(n)k(n))& m=1,2,\dots ,M\end{array}\end{array}\text{\hspace{1em}(4-37)}$$

改进结构中系数 k 不变，只有权重 v→u。对 $e(n)=[1,1,1,1,1]$ 的 5 阶滤波器，一小段 MATLAB 脚本即可算出：

$\mathrm{u}(1) = 0.649195217878$ $\mathrm{u}(2) = 0.000421764763$ $\mathrm{u}(3) = 0.329929109321$ $\mathrm{u}(4) = 0.000084546124$ $\mathrm{u}(5) = 0.020062621320$ $\mathrm{u}(6) = 0.000303581790$

效果立竿见影：系数跨度从双乘法器结构的 17.1 位降到 ${\log }_2(\max (u)/\min (u))=12.9$ 位。

网格滤波器的固有代价：长延迟路径

位宽问题解决了，但还有一个网格滤波器绕不开的硬伤：从输入到输出存在一条很长的无寄存器组合路径。标准网格结构是 5 个加法器 + 1 个乘法器的延迟；Gray-Markel 结构更糟，达 10 个加法器 + 5 个乘法器。这条路径直接压低 $F_{max}$，使网格滤波器在 HDL/FPGA 设计中的吸引力大打折扣。一句话总结本单元：窄带 IIR 在 FPGA 上的最优解通常不是最”聪明”的结构（网格），而是位宽、速度、复杂度均衡得最好的改进型并联双二阶。

4.6.5 窄带 IIR 滤波器的波形数字滤波器设计

WDF 是什么：从模拟梯形电路”整体搬家”到数字域

前面几节介绍的都是”从传递函数出发”的设计方法：先算出分子分母系数，再套进直接型、级联型或并行型结构。波形数字滤波器（Wave Digital Filter，WDF）走的是另一条路——它不拆散电路，而是把一个现成的 R/L/C 模拟梯形滤波器整体变换到数字域。典型原型是一类 5 阶模拟考尔滤波器（又称椭圆滤波器），如图 4-46 所示。

图 4-46 WDF 模拟 5 阶原型滤波器

为什么要这样做？因为梯形模拟电路的灵敏度特性非常好：元件值的微小偏差对通带/阻带指标影响很小。WDF 通过双线性变换和分量变换，把这种低灵敏度”继承”到数字域——系数量化后传递函数几乎不动，这正是它对 FPGA 实现有吸引力的原因。

变换规则很简单：模拟域的电容变成单位延迟 $z^{-1}$，电感变成负增益 $-z^{-1}$。源阻抗和负载阻抗通常归一化为 1，在 z 域中可以直接忽略。变换后每个模块的输入”波形”记为 $A_k$，输出”波形”记为 $B_k$（注意这是数字信号，只是沿用了”波形”这个名字）。剩下的问题是：各个 C、L 模块的值不同，怎么把它们连起来？答案是适配器（Adaptor）——它负责处理串行/并行两种连接类型，并调整各端口的分量值。对照图 4-46 的梯形结构：“竖直”的并联支路需要并行适配器，“水平”的串联支路需要串行适配器。

并行三端口适配器

图 4-47(b) 是并行三端口适配器的电路符号，图 4-47(a) 是 SIMULINK 实现。提醒一点：后面的硬件实现里，普通乘法可以换成常系数乘法器以省资源，但设计阶段先用通用乘法器更方便——便于试探系数量化到多少位才够。

图 4-47 WDF 并行三端口适配器

适配器的输入输出关系由基尔霍夫定律导出（节点电流之和为零、回路电压之和为零）。本书符号与 Fettweis 1986 年的经典文献一致。对并行三端口适配器可得：

$$\left[\begin{array}{l}{B}_1\\ B_2\\ B_3\end{array}\right]=\left[\begin{array}{ccc}{c}_1& 2-c_1-c_3& c_3\\ c_1& 1-c_1-c_3& c_3\\ c_1& 2-c_1-c_3& c_3-1\end{array}\right]\left[\begin{array}{l}{A}_1\\ A_2\\ A_3\end{array}\right]\text{\hspace{1em}(4-39)}$$

无延迟循环：为什么必须”匹配”一个端口

数字滤波器里最怕出现无延迟环路——一个组合环里没有任何寄存器，硬件根本无法计算（每个输出都依赖自己，且”立刻”依赖）。所以 WDF 规定：三端口适配器必须有一个端口是”无映射”的，即 $B_k\ne f(A_k)$，这样的端口称为匹配端口。

不用会怎样？如果不匹配，例如第 3 端口的 $B_3$ 依赖 $A_3$，而 $A_3$ 又来自另一个适配器的输出 $B_3$，两个模块互相等待，环里没有延迟元件，电路就是死的。由式 (4-39) 可见，要使 $B_3$ 不含 $A_3$，必须令 $c_3=1$。代回得匹配的并行三端口适配器：

$$\left[\begin{array}{l}{B}_1\\ B_2\\ B_3\end{array}\right]=\left[\begin{array}{ccc}{c}_1& 1-c_1& 1\\ c_1& -c_1& 1\\ c_1& 1-c_1& 0\end{array}\right]\left[\begin{array}{l}{A}_1\\ A_2\\ A_3\end{array}\right]\text{\hspace{1em}(4-40)}$$

注意矩阵里出现了两个 1 和一个 0——匹配不仅避免了无延迟环，还白送一个”乘 1”和”乘 0”，乘法器数量减一。简化的 SIMULINK 设计和电路符号见图 4-48。

图 4-48 WDF 并行匹配的三端口适配器

串行三端口适配器

串行适配器推导方式相同，公式和符号流图见图 4-49，矩阵方程为：

$$\left[\begin{array}{l}{B}_1\\ B_2\\ B_3\end{array}\right]=\left[\begin{array}{ccc}1-c_1& -c_1& -c_1\\ -c_3& c_1+c_3-1& c_1+c_3-2\\ c_1+c_3-2& -c_3& 1-c_3\end{array}\right]\left[\begin{array}{l}{A}_1\\ A_2\\ A_3\end{array}\right]\text{\hspace{1em}(4-41)}$$

图 4-49 WDF 串行三端口适配器

同样要求第三个端口无映射（$B_3\ne f(A_3)$），令 $c_3=1$，式 (4-41) 简化为：

$$\left[\begin{array}{l}{B}_1\\ B_2\\ B_3\end{array}\right]=\left[\begin{array}{ccc}1-c_1& -c_1& -c_1\\ -1& c_1& c_1-1\\ c_1-1& -1& 0\end{array}\right]\left[\begin{array}{l}{A}_1\\ A_2\\ A_3\end{array}\right]\text{\hspace{1em}(4-42)}$$

这个无延迟版本对硬件同样友好：矩阵中的常数 1、0、-1 都不用乘法器，比未匹配版本少 1 个乘法器。匹配的串行三端口适配器见图 4-50。

图 4-50 WDF 串行匹配的三口适配器

组装 5 阶 WDF 并生成系数

有了并行/串行适配器及其匹配版本，就可以把图 4-46 的模拟梯形电路完整搬到数字域，得到图 4-51 的 5 阶 WDF 数字滤波器。

图 4-51 5 阶 WDF 数字滤波器

还缺系数。用代尔夫特理工大学 H. Lincklaen Arriens 开发的公开 WDF 工具箱，几行 MATLAB 就能生成全部系数与结构。目标指标：$F_{pass}=40\text{Hz}$、$F_{stop}=60\text{Hz}$、$F_s=7680\text{Hz}$、$A_{pass}=1\text{dB}$、$A_{stop}=48\text{dB}$。

先生成归一化低通梯形网络：

%% Generate the filter coefficients for system:
%% X = NLP_LADDER('cauer', filterOrder, passBandRipple_dB, ...
%%     stopBandRipple_dB, skwirNorm, freqNormMode)
Ladder = nlp_ladder('cauer', 5, 1, 48, 'a', 1);

工具箱给出 7 个适配器的系数：

Adaptor 1: 3p parallel, p3 matched : alpha1 = 0.00815
Adaptor 2: 3p parallel, p3 matched : alpha1 = 0.99882
Adaptor 3: 3p serial,   p3 matched : alpha1 = 0.10329
Adaptor 4: 3p parallel, p3 matched : alpha1 = 0.07774
Adaptor 5: 3p parallel, p3 matched : alpha1 = 0.99946
Adaptor 6: 3p serial,   p3 matched : alpha1 = 0.19518
Adaptor 7: 3p parallel             : alpha1 = 0.46866
                                     alpha3 = 0.01472

再做去归一化（截止频率换算到 40 Hz）并计算 WDF 结构与脉冲响应：

%% Denormalize the ladder filter
%% LpLadder = NLADDER2LP(NlpLadder, cutOffFrequency)
dLadder = nladder2lp(Ladder, fz2fs(40/7680));
 
%% Compute the WDF filter data and impulse responses:
%% [WDF, fwdB, revB, allB] = LADDER2WDF(Ladder, wdfType, ...
%%     impulseResponseLength, figNo)
[WDF, fwdB, revB, allB] = ladder2WDF(dLadder, '3p', 8*512, 3);

工具箱同时给出模拟梯形的元件值（两种镜像配置）：

Configuration 1:
Rs 1.00000 Ohm
C01 1.99211 F in shunt arm
L02 0.97767 H, parallel with
C02 0.23108 F in series arm
C03 2.46127 F in shunt arm
L04 0.76116 H, parallel with
C04 0.64686 F in series arm
C05 1.68564 F in shunt arm
RL 1.00000 Ohm

Configuration 2:
Rs 1.00000 Ohm
C01 1.68564 F in shunt arm
L02 0.76116 H, parallel with
C02 0.64686 F in series arm
C03 2.46127 F in shunt arm
L04 0.97767 H, parallel with
C04 0.23108 F in series arm
C05 1.99211 F in shunt arm
RL 1.00000 Ohm

位宽怎么定：整数位与小数位分开看

工具箱通过逆 FFT 估计每个端口变量的最坏情况幅值增长，并给出脉冲响应，如图 4-52 所示。

Brev needs negation

图 4-52 5 阶滤波器所有端口的 WDF 最大幅值（图片来自 H. Lincklaen Arriens 的 WDF 工具箱）

由图可见，只有匹配的并行适配器 2 和 5 需要可观的整数位宽：幅值接近 100，取 $\lceil {\log }_{2}100\rceil =7$ 位整数即可；其余适配器只需 1~2 位整数（保险起见）。整数位和小数位分开估计是 WDF 实用化的关键——若笼统给全部信号同一位宽，会浪费大量资源。

小数位的确定用系数量化试验：置 $c_q(i)=\text{round}(c(i)\times 2^N)/2^N$（先把系数乘 $2^N$ 取整，再除回 $2^N$），然后从最小的 N 试起，直到传递函数重新满足误差指标。图 4-53 的对比显示：14 位系数量化已足以满足通带与阻带要求——这正是 WDF 低灵敏度特性的直接体现。

图 4-53 波形数字滤波器量化系数的传递函数图，上图——使用浮点精度（FP），下图——使用 14 位系数量化

WDF 的优点与”阿喀琉斯之踵”

整体上 WDF 整数增益行为好、分量数少。但它有一个硬伤：从输入到输出存在一条组合直通路径。对前向输出，“波形”要穿过从输入到输出的一串适配器；更麻烦的是，即便不需要反向输出，为了维持所有延迟元件的状态，下一个采样到来前也必须算完整个反向路径。前向、反向两条长组合链，直接限制了最高时钟频率。

补救办法是用”对称”WDF（非匹配适配器放在中间而非末端，只适用奇数阶），可显著缩短前向/反向路径——文献[117]给出过中间放置非匹配适配器的 7 阶 WDF。但这类结构极难流水线化，因为插入流水线寄存器会改变环内延迟数，极点/零点位置随之改变，传递函数就不再保真。

作为过渡，这里先按资源、整数/小数位、最长路径三个维度对几种结构做个粗评（无 HDL 实现时数据是估计值，但足以选出靠谱结构）：

• 直接型：算术运算数量看着不差，但位宽要求很高，浮点下要双精度（64 位）。
• 级联（椭圆滤波器）：零点在单位圆上，乘法器很省；但恢复块增益为 1 需要额外乘法器；用转置直接 I 型（双二阶）可把最坏路径做短。
• 并行双二阶：加法器更少、无级联的节顺序问题、天然并行；配合改进双二阶可做到 1 加法器 + 1 乘法器的最短路径，是很有竞争力的方案。
• 网格滤波器：结构复杂些，小数位宽不错，但整数位宽太大，反而不及并行/级联型。
• WDF：乘法器最省，但加法器多，且长延迟路径让直接实现的 WDF 不具吸引力——除非改用下面要讲的网格 WDF（LWDF）。

4.7 窄带 IIR 滤波器的全通滤波器设计

用两个全通滤波器”做减法”造窄带

全通滤波器的定义：所有频率下幅值响应恒为 1，只改变相位，通常用来做相位延迟。读者自然会问：幅值处处为 1，怎么可能构成窄带滤波器？诀窍是同时运行两个全通滤波器，把它们的输出相加（或相减）。

原理藏在相位里：设两个全通滤波器在通带内相位延迟都为零，但在阻带内二者的相位差为 $180^{\circ }$。由于 $\sin (\omega t+180^{\circ })=-\sin (\omega t)$，把两路信号相加，阻带频率分量正好抵消——通带照常通过，一个窄带滤波器就成了。图 4-54 给出了典型全通设计中的相位关系。此外，因为 $|F(\omega )|=1$，合成滤波器的整数增益很小，硬件上又是一重实惠。

(a) 传递函数

(b) 相位 ${\mathrm{A}}_1(\mathrm{z})$、${\mathrm{A}}_2(\mathrm{z})$

(c) 极点/零点 ${\mathrm{A}}_1(z)$

(d) 传递函数

(e) 相位差

(f) 极点/零点 $A_2(z)$
图 4-54 并行全通滤波器设计的传递函数、相位和极点零点图

这种结构有约束：低通/高通型全通分解只能取奇数阶（偶数阶会导致复数系数，实系数实现代价更高[120]）；而实系数带通/带阻全通分解则要求总阶数为偶数。

求和/求差结构与极点分配

两个全通支路求和或求差的结构，画出来和单节网格滤波器很像（参见图 4-38）。这一思路最早用于组合 WDF，所以文献里叫它网格 WDF（LWDF）——但方法本身适用于一切全通滤波器：级联/并行双二阶、Gray-Markel 网格都可以。带单乘法器的网格滤波器尤其高效，因为它天然就是全通结构。

输出由下式合成：

$$H_{1,2}=\frac{1}{2}\left(A_1(z)\pm A_2(z)\right)\text{\hspace{1em}(4-43)}$$

设计时要让全通滤波器的极点与原始滤波器的极点匹配。Gazsi 指出：两个全通滤波器的极点相对单位圆交替分布；且巴特沃思、切比雪夫、椭圆等主流滤波器（都是经双线性变换设计的）都可以转换为并行全通形式[113]。对本 5 阶系统：最靠近单位圆的极点（含实极点）归入 $A_1(z)$，其余极点对归入 $A_2(z)$——即一个三阶加一个二阶的并行系统，见图 4-54(c) 与 (f)。

MATLAB 的 tf2cl（传递函数到全通网格对）直接支持这种设计：

$$[\text{latt1},\text{latt2}]=\text{tf2cl(B, A)}$$

对 5 阶系统运行，再把网格结构转回传递函数与直接形式以便画图：

[num1, den1] = latc2tf(latt1, 'allpass')
[num2, den2] = latc2tf(latt2, 'allpass')

得到（式 4-44）：

num1 = 0.9852  -1.9847   1.0000
den1 = 1.0000  -1.9847   0.9852
num2 = -0.9846  2.9682  -2.9835   1.0000
den2 = 1.0000  -2.9835   2.9682  -0.9846

latt1 = -0.9997   0.9852
latt2 = -0.9996   0.9998  -0.9846

留意分子分母的镜像关系：num1 = [0.9852, -1.9847, 1] 与 den1 = [1, -1.9847, 0.9852] 系数首尾对调——这正是全通滤波器”分子是分母镜像多项式”的特征，后面 4.7.3 节会利用它省一半乘法器。

图 4-54(a)(d) 是两个全通的传递函数；(b) 是各自相位，(e) 是相位差——阻带处清晰可见 $180^{\circ }$ 相位差；极点交替分布见图 4-54(c)(f)。

在这一全通框架下，再次从整数/小数位角度比较四种实现（级联/并行双二阶、WDF、网格滤波器），结论放在 4.7.5 之后的表 4-5 中统一汇总。这里先说前提：并行和级联用改进双二阶，网格用 Gray-Markel 单乘法器型，全通 WDF 用三端口适配器（性能最佳）。

4.7.1 窄带 IIR 滤波器的全通波形数字滤波器设计

用工具箱直接生成 LWDF

仍用同一工具箱，LWDF 的 MATLAB 脚本与标准 WDF 类似：

% 2- and 3-port circuit design as in Anderson et al (1995)
% Generate the filter coefficients for system:
% [Hs,wp] = HS_CAUER(filterOrder, passBandRipple_dB, ...
%     stopBandRipple_dB, skwirMode, cutOffFrequency, freqNormMode)
Hs = Hs_cauer(5, 1, 50, 'a', 1, 1);
 
% NLP2LP Normalized lowpass to lowpass transformation
dHs = nlp2lp(Hs, fz2fs(40.75/7680));
 
% Calculates the coefficients of a Lattice WDF
% [LWDF, Hz, Messages] = HS2LWDF(Hs, figNo)
[LWDF, Hz, Messages] = Hs2LWDF(dHs, 3);
% showLWDF
Hz2 = LWDF2Hz(LWDF); plotHz(Hz2, 1);

注意两处指标微调：通带 $40\text{Hz}\to 40.75\text{Hz}$、阻带波纹 $48\text{dB}\to 50\text{dB}$。原因：LWDF 设计成低通/高通功率互补、交汇于 3 dB 点；在 ±1 dB 波纹规范下，若恰好指定 40 Hz 通带，量化后通带会过早跌落，所以预留一点余量，确保量化后仍有 40 Hz 的实际通带[118]。脚本绘出的两端口基本结构见图 4-55。

图 4-55 使用两端口适配器的 LWDF 5 阶数字滤波器

从两端口系数换算三端口系数

先提取各节乘法因子：

%% Upper 1. order
beta(1) = LWDF.gamma(1,1,1);
%% Upper 2. order
beta(2) = LWDF.gamma(1,2,1);
beta(3) = LWDF.gamma(2,2,1);
%% Lower 2. order
beta(4) = LWDF.gamma(1,1,2);
beta(5) = LWDF.gamma(2,1,2);

两端口系数值为：

$$0.9887-0.99590.9995-0.98530.9997\text{\hspace{1em}(4-45)}$$

工具箱没有三端口实现，但按 Anderson/Sunnerfield/Lawson[121]的方法换算并不费事。为避免混淆，把两端口系数记为 ${\beta }_k$、三端口系数记为 ${\gamma }_k$。两端口适配器的传递函数为：

$$h_2(z)=\frac{-{\beta }_1-{\beta }_2(1-{\beta }_1)z^{-1}+z^{-2}}{1-{\beta }_2(1-{\beta }_1)z^{-1}-{\beta }_1{z}^{-2}}\text{\hspace{1em}(4-46)}$$

三端口设计的 z 域传递函数为：

$$h_3(z)=\frac{-{\gamma }_1-{\gamma }_2-1+({\gamma }_2-{\gamma }_1)z^{-1}+z^{-2}}{1+({\gamma }_2-{\gamma }_1)z^{-1}+(-{\gamma }_1-{\gamma }_2-1)z^{-2}}\text{\hspace{1em}(4-47)}$$

令 $h_2(z)=h_3(z)$，对比系数可得换算公式：

$${\gamma }_1=-\frac{1}{2}(1-{\beta }_1)(1-{\beta }_2)\text{\hspace{1em}(4-48)}$$ $${\gamma }_2=-\frac{1}{2}(1-{\beta }_1)(1+{\beta }_2)\text{\hspace{1em}(4-49)}$$

展开算一遍第一组（${\beta }_1=0.9987$ 对应一阶节直接沿用，${\beta }_2=-0.9959$、${\beta }_3=0.9995$ 对应二阶节）：$(1-{\beta }_2)=1.9959$，$(1-{\beta }_3)=0.0005$，故 ${\gamma }_2=-0.5\times 1.9959\times 0.0005\approx -0.000528$；$(1+{\beta }_3)=1.9995$，故 ${\gamma }_3=-0.5\times 1.9959\times 1.9995\approx -1.995400$。一阶节无需换算，${\gamma }_1={\beta }_1$。对应代码：

%% Upper 1. order
gamma(1) = beta(1);
%% Upper 2. order
gamma(2) = -0.5*(1-beta(2))*(1-beta(3));
gamma(3) = -0.5*(1-beta(2))*(1+beta(3));
%% Lower 2. order
gamma(4) = -0.5*(1-beta(4))*(1-beta(5));
gamma(5) = -0.5*(1-beta(4))*(1+beta(5));

最终三端口配置的 $\gamma$ 值为：

$$0.988727-0.000528-1.995400-0.000282-1.985024\text{\hspace{1em}(4-50)}$$

三端口 LWDF 结构如图 4-56 所示。注意 ${\gamma }_2$、${\gamma }_4$ 都接近 0，${\gamma }_3$、${\gamma }_5$ 接近 $-2$——系数越”整”，乘法器实现越省。

图 4-56 使用三端口适配器的 LWDF 5 阶数字滤波器

位宽评估：15 位小数、6 位整数

与 4.6.5 节相同，实现前先定位宽。对两端口和三端口 LWDF，15 个小数位即可满足指标，量化前后传递函数对比如图 4-57。

图 4-57 使用三端口适配器的网格波形数字滤波器的量化系数的传递函数（TF）图：上图——使用浮点（FP）精度，下图——使用 15 位系数量化

再逐个记录加法器输出的（绝对）最大值以确定整数位，结果见图 4-58：两种配置下最大值都是 64，6 位整数精度足够。

图 4-58 波形数字滤波器最大幅值的计算结果：上图——使用两端口配置，下图——使用三端口配置

资源与延迟：LWDF 的真正杀手锏

与标准 WDF 相比，LWDF 的资源和延迟估计都大幅改善：

• 两端口设计：16 个加法器、5 个乘法器（含低通网格输出用加法器；若还要高通网格输出，再加 1 个加法器）。
• 三端口设计：乘法器相同（5 个），加法器降到 12 个。
• 对比标准 WDF 的 29 加法器 + 8 乘法器，LWDF 节省约 50% 硬件。

最大优势在延迟：LWDF 没有反馈环穿多个节——每个两端口/三端口节后面都可以插流水线寄存器，把关键路径压缩到”两端口：2 加法器 + 1 乘法器；三端口：3 加法器 + 1 乘法器”。唯一要注意的是上、下两支路的流水线延迟必须相等（保证输出对齐），为此在两个分支里各补 2 个寄存器即可。回忆 4.6.5 节说的”WDF 难流水线化”，LWDF 因为结构规整、环短，恰好绕开了这个问题。

4.7.2 窄带滤波器的全通网格设计

4.6.4 节介绍过基本网格滤波器：它同时具备全极点输出和全通输出（图 4-38）。正因为网格滤波器天生带全通输出，不需要通用 IIR 设计里额外的输出求和网络，做全通分解几乎是”免费”的：两个并行全通（即式 4-43）的整体运算量相比标准单网格设计有本质改进。

窄带情形有个坑：每节两个乘法器的网格结构，整数增益往往偏大。改用 Gray-Markel 每节单乘法器结构后，增益可在网格架构内精细调节，整数位要求从双乘法器配置的 16 位降到 6 位；最长延迟路径也改进为 3 加法器 + 1 乘法器。系数就是 4.7 节由 tf2cl 算出的式 (4-44)。对单乘法器结构，Gray-Markel 优化方法给出二阶节符号 $(-,-)$、三阶节符号 $(+,+,+)$。

4.7.3 窄带滤波器的全通直接型设计

直接型的机会来自全通滤波器的数学性质：分子是分母的镜像多项式，极点/零点相对单位圆对称——极点在圆内（稳定性要求），零点在圆外的镜像位置。于是转置直接 I 型里一半的系数乘法器可以共享，乘法器数量减半。这类结构运算量最小：5 阶椭圆滤波器采用 16 位即可。

代价有两点：一是为了利用系数相似性，无法采用改进双二阶那种”零点优先”的缩放设计（参见图 4-31），结果是整体整数位宽高达 12 位（各方案中最高），小数位也要 22 位；二是位宽总量虽从单个直接型的 63 位降到全通两支路的 34 位，但仍明显高于其他结构。

4.7.4 窄带滤波器的全通级联双二阶设计

改用双二阶配置可对直接型小幅改进。三阶的下半支路（上半支路是二阶，已无可拆）可以分解为一个双二阶加一个一阶系统，从而把小数位降到 16 位、整体位宽降到 28 位。加法器/乘法器数量和延迟等其他参数不变。

4.7.5 窄带滤波器的全通并行双二阶设计

非全通情形下我们偏爱并行形式，但全通设计中并行配置有个致命弱点：并行分解会破坏系数对称性。以上半二阶支路无需处理、看下半三阶支路为例，原始系数（式 4-51）：

$$\begin{array}{rcl}\text{num2}& =& -0.98462.9682-2.98351.0000\\ \text{den2}& =& 1.0000-2.98352.9682-0.9846\end{array}\text{\hspace{1em}(4-51)}$$

分子分母镜像对称。先用部分分式展开把它拆成单极点之和：

$$[\text{R},\text{P},\text{D}]=\text{residuez(num2, den2)}$$

把两个复极点合并成一个实二阶节，单实极点写成常见分子/分母形式：

[b1, a1] = residuez([R(1) R(2)], [P(1) P(2)], 0);
[b2, a2] = residuez(R(3), P(3), 0);

得到的并行全通双二阶系数为：

D  = -1.0156
b1 = 0.0063  -0.0065   0
a1 = 1.0000  -1.9948   0.9959
b2 = 0.0246   0
a2 = 1.0000  -0.9887

可见并行滤波器的系数已完全丧失与式 (4-51) 的对称关系（b1、b2 变成很小的杂乱数值，还多出一个直接项 D）。后果是运算量增加：全通并行实现需 19 次操作，而级联双二阶只要 16 次；延迟与整数/小数位要求与级联相近。

可以在 MATLAB 里验证重建是否正确（把并联各节重新合并，应还原原系统）：

num2 = conv(b1, a2) + conv(b2, a1)
den2 = conv(a1, a2)
num2 = den2.*D + num2

全通设计总比较

表 4-5 汇总了各全通设计选项。从位宽看 LWDF 最好（整数 6 位 + 小数 15 位 = 21 位），且资源和延迟与其他方案相当；三端口 LWDF 的操作次数（17 次）与直接型/级联型（16 次）接近，延迟路径最好的是直接型、级联型和两端口 LWDF（2+1×）。但直接型和级联型需要超过两倍的位宽——综合权衡，LWDF 是胜出者，后续示例也将实现 LWDF。

表 4-5 窄带 IIR 滤波器网格全通设计选项的比较

结构		操作			延迟	位
		加法器	乘法器	$\Sigma$		整数	小数	$\Sigma$
直接		11	5	16	2+1×	12	22	34
级联	双二阶	11	5	16	2+1×	12	16	28
并行	双二阶	10	9	19	2+1×	12	16	28
网格	2乘法器	11	10	21	3+1×	16	14	30
网格	1乘法器	16	5	21	6+3×	7	14	21
LWDF	2端口	16	5	21	2+1×	6	15	21
LWDF	3端口	12	5	17	3+1×	6	15	21

例 4.9 全通三端口波形数字滤波器的 HDL 设计

下面的设计实现了图 4-56 的三端口 LWDF 结构 5 阶 IIR 滤波器。原书代码是 VHDL，此处已改写为等价的可综合 Verilog HDL：系数 $\gamma$ 定点化为 Q6.15 格式（即 round(gamma*2^15)），内部数据 22 位有符号，输入输出为 16.16 格式的 32 位向量；VHDL 中 resize(..., fixed_wrap, fixed_truncate) 的”截断+回绕”语义用 Verilog 的乘法后右移取低位实现。

// 原书 VHDL 改写为 Verilog
// 描述: 5 阶网格波形数字滤波器 (全通三端口 LWDF, 图 4-56)
// gamma 系数: 0.988727 -0.000528 -1.995400 -0.000282 -1.985024
module iir5lwdf (
    input  wire        clk,      // 系统时钟
    input  wire        reset,    // 异步高有效复位
    input  wire [31:0] x_in,     // 系统输入, 16.16 有符号定点
    output wire [31:0] y_ap1out, // 全通支路 1 输出
    output wire [31:0] y_ap2out, // 全通支路 2 输出
    output wire [31:0] y_ap3out, // 全通支路 3 输出
    output wire [31:0] y_out     // 系统输出
);
    localparam W    = 22;  // 内部 6.15 有符号定点, 共 22 位
    localparam FRAC = 15;  // 小数位数
 
    // gamma 系数, 定点值 = round(gamma * 2^15)
    localparam signed [W-1:0] G1 =  22'sd32405;  //  0.988727
    localparam signed [W-1:0] G2 = -22'sd17;     // -0.000528
    localparam signed [W-1:0] G3 = -22'sd65377;  // -1.995400
    localparam signed [W-1:0] G4 = -22'sd9;      // -0.000282
    localparam signed [W-1:0] G5 = -22'sd65037;  // -1.985024
 
    // 常系数乘法: 积取 44 位, 右移 FRAC 后取低 W 位
    // (等价于 VHDL 的 fixed_truncate + fixed_wrap)
    function signed [W-1:0] mul;
        input signed [W-1:0] a, b;
        reg   signed [2*W-1:0] p;
        begin
            p   = a * b;
            mul = p[W+FRAC-1:FRAC];
        end
    endfunction
 
    // 输入 16.16 -> 内部 6.15 (算术右移 1 位, wrap 截断)
    wire signed [W-1:0] x = $signed(x_in) >>> 1;
 
    // 状态寄存器 (每次赋值对应一个寄存器)
    reg signed [W-1:0] c1, c2, c3, l2, l3;
    reg signed [W-1:0] ap1, ap2, ap3, ap3r, y;
    // 组合中间变量 (不 infer 寄存器)
    reg signed [W-1:0] p1, a4, a5, a6, a8, a9, a10;
 
    always @(posedge clk or posedge reset) begin
        if (reset) begin                      // 复位全部寄存器
            y <= 0;    c1 <= 0;  ap1 <= 0;
            c2 <= 0;   l2 <= 0;  ap2 <= 0;
            c3 <= 0;   l3 <= 0;  ap3 <= 0;  ap3r <= 0;
        end else begin
            // ---- 1. 全通节: 一阶 ----
            p1  = mul(G1, c1 - x);
            c1  <= x + p1;
            ap1 <= c1 + p1;
            // ---- 2. 全通节: 二阶 (上支路) ----
            a4  = ap1 - l2 + c2;
            a5  = mul(G2, a4) + c2;
            a6  = mul(G3, a4) - l2;
            c2  <= a5;
            l2  <= a6;
            ap2 <= -a5 - a6 - a4;
            // ---- 3. 全通节: 二阶 (下支路) ----
            a8   = x - l3 + c3;
            a9   = mul(G4, a8) + c3;
            a10  = mul(G5, a8) - l3;
            c3   <= a9;
            l3   <= a10;
            ap3  <= -a9 - a10 - a8;
            ap3r <= ap3;      // 为对齐 AP1 增加的额外寄存器
            // ---- 输出加法器 ----
            y <= ap3r + ap2;
        end
    end
 
    // 内部 6.15 -> 输出 16.16 (左移 1 位并符号扩展)
    function [31:0] to_1616;
        input signed [W-1:0] v;
        to_1616 = {{9{v[W-1]}}, v, 1'b0};
    endfunction
 
    assign y_out    = to_1616(y);
    assign y_ap1out = to_1616(ap1);
    assign y_ap2out = to_1616(ap2);
    assign y_ap3out = to_1616(ap3);
 
endmodule

几点设计要点，对应原 VHDL 的写法逐条说明：

• 接口：除系统输入输出外，还引出三个全通支路的输出 $y_{ap1},y_{ap2},y_{ap3}$，方便在仿真中逐段核对（图 4-56 上半支路为一阶 + 二阶，下半支路为二阶）。
• 系数：5 个 $\gamma$ 用 localparam 直接按 Q6.15 定点常数给出，省去运行时配置。
• 一个 always 块覆盖三个全通节：复位时清零全部寄存器；posedge clk 分支里按”一阶节 → 上二阶节 → 下二阶节 → 输出加法器”的顺序书写，每个非阻塞赋值综合为一个寄存器。
• 位宽控制：所有乘法都在 mul 函数里右移 FRAC 并截取低 22 位，模拟原书 resize(..., fixed_wrap, fixed_truncate) 的回绕+截断语义。原书特别指出：若不用 wrap（改用饱和），LE 用量会翻倍以上——因为饱和比较逻辑远贵于自然截断。
• ap3r 寄存器：下半支路比上半支路少一级延迟，额外加一级寄存器使两支路输出对齐后才能相加。

该设计在 Cyclone IV E（EP4CE115F29C7）上使用 764 个 LE、12 个 9×9 嵌入式乘法器，TimeQuest 缓慢 85C 模型下 $F_{max}=55.97\text{MHz}$。滤波器对幅值 $1.0=65536_{10}$（16.16 格式）阶跃输入的响应仿真如图 4-59 所示，与 MATLAB 仿真结果一致。100 ns 系统时钟下，60 μs 仿真对应 600 个时钟周期；最大阶跃响应出现在约 215 个时钟周期处。

图 4-59 5 阶全通波形数字网格滤波器的阶跃响应的 MODELSIM 仿真结果

与例 4.7 中最好的非全通设计相比，LWDF 在嵌入式乘法器用量上优势明显（省 50%），LE 数量和速度也稍优于改进并行双二阶设计。表 4-6 汇总了本章窄带滤波器各设计的综合结果。

表 4-6 窄带 IIR 滤波器的 HDL 合成结果

结构	封装	LE	乘法器9×9	$F_{max}$ (MHz)
直接	—	6279	128	29.96
直接	√	2474	128	46.99
并行双二阶	—	1874	51	54.00
并行双二阶	√	624	51	87.69
全通三端口 LWDF	—	1465	12	33.83
全通三端口 LWDF	√	764	12	55.97

最后补充几张原书练习部分的配图。练习 4.10 要求设计 8 位输入的 10 阶巴特沃思 IIR 滤波器（内部采用 14.12 位格式），仿真结果可与图 4-60 给出的测试平台（其中 t_out 是递归部分的输出）相对照；练习 4.11 的 PREP 基准 5 是一个 4×4 无符号阵列乘法器加 8 位累加器（图 4-61）；练习 4.12 的 PREP 基准 6 则是一个带异步复位的 16 位累加器（图 4-62）。

图 4-60 练习 4.10 的 IIR 巴特沃思测试平台

图 4-61 PREP 基准 5

图 4-62 PREP 基准 6