📘 第 5 章：多级信号处理

这一篇在干嘛？

采样率变换的艺术：抽取/插值、多相分解、Hogenauer CIC 滤波器、小波变换——过采样系统的核心工具箱。原书代码为 VHDL，本篇所有代码已改写为 Verilog。

第5章
多级信号处理
5.1 抽取和插值
5.2 多相分解
例5.1 多相抽取滤波器
例5.2 快速FIR滤波器
算法5.3 快速FIR滤波器
5.3 Hogenauer CIC 滤波器
例5.4 三级CIC抽取器I
BEGIN
例5.6 三级CIC抽取器III
5.4 多级抽取器
使用 Goodman-Carey 半带滤波器的多级抽取器设计
定义 5.7：半带滤波器
例 5.8 多级半带滤波器抽取器
5.5 作为通频带抽取器的频率采样滤波器
5.6 任意采样速率转换器的设计
例5.9 $R = 0.75$ 的速率变换器I
算法5.10 应用FFT的有理数速率变换
例 5.11 R=0.75 的速率变换器 II
例5.12 $R = 0.75$ 的速率变换器III
例5.13 $R = 0.75$ 的速率变换器IV
例 5.14 R=0.75 的速率变换器 V
5.7 滤波器组
例5.16 双通道Haar滤波器组I
定理 5.17 完美的重构
例 5.18 双通道 Haar 滤波器组 II
定理5.19 与混叠无关的双通道滤波器组
算法5.20 完美重构的双通道滤波器组
例5.21 采用F3的完美重构滤波器组
推论 5.22 半带滤波器的因式分解
1. 多相双通道滤波器组
1. 提升
例 5.23 DB4 滤波器的提升实现
1. QMF 实现
1. 正交滤波器组
1. 线性相位双通道滤波器组
例5.25 $L = 4$ 的线性相位滤波器的网格
1. 实现选项的比较
5.8 小波
例 5.26 线性调频信号的分析
例5.27 长度为4的Hutlet滤波器
5.9 练习
5.25 利用练习 5.24 的结果。

自测一下

抽取前为什么必须先低通滤波？

直接丢点（抽取）会让新采样率下的高频分量混叠到带内。先抗混叠低通把新奈奎斯特频率以上的能量压掉，再抽取才不失真。

CIC 滤波器为什么"不用乘法器"？

CIC 的系数全是 1（积分器+梳状器级联），只需加法器与寄存器。代价是通带衰减大，通常放在高速率前端、后接补偿 FIR。

多相分解省了什么？

把长滤波器拆成 R 相，每相只在低速率时钟下处理自己那一路，等效把乘法器工作频率降到 1/R，或者用 1/R 个乘法器完成原本 R 倍速的计算——硬件量大幅下降。

半带滤波器有什么特点？

约一半系数恰好为 0（除中心点外偶数系数全为 0），乘法量比同长 FIR 少一半；通阻带对称，特别适合 2 倍抽取/插值级联。

第5章多级信号处理

在数字信号处理中，有一类非常常见的任务：根据感兴趣的信号来调整采样速率。由不同采样速率的系统组成的处理链，就称为多级系统。本章通过两个典型示例说明多级 DSP 系统中的抽取与插值，然后引入多相分解这一重要的实现技巧，研究高效的抽取器设计，最后讨论 Hogenauer CIC 滤波器——高速抽取/插值系统中几乎绕不开的构件。

5.1 抽取和插值

什么是抽取，为什么需要它

如果 A/D 转换之后，感兴趣的信号只占据一个很窄的频带（通常是低通或带通），那么用很高的采样速率来承载它就是一种浪费：后续每一步运算都要为那些”空”的频谱成分买单。合理的做法是先窄带滤波，再降低采样速率——滤波加向下采样的组合通常称为抽取器。滤波、向下采样以及对频谱的影响如图 5-1 所示。

图 5-1 信号 $x [n] \circ ∙ X (ω)$ 的抽取

采样速率能降到多低？下限就是奈奎斯特速率：采样速率必须高于信号带宽，否则混叠就会出现。图 5-2 对比了无混叠抽取与有混叠抽取的情形。要特别记住：混叠是不可修补的——一旦不同频带的成分叠在一起，无论后面做什么处理都无法把它们分开，所以必须在抽取之前用低通滤波器把它们滤干净。

图 5-2 无混叠抽取和有混叠抽取的情形

带通信号的完整频带条件

对带通信号还有一条额外规则：有用频带必须整体落在一个”完整频带”之内。设采样速率为 $f_{s}$ 、向下采样因子为 $R$ ，则有用带宽必须满足

k \frac{f _{s}}{2 R} < f < (k + 1) \frac{f _{s}}{2 R}, k \in N (5-1)

初学者容易误以为”采样速率高于奈奎斯特频率就万事大吉”，式(5-1)正是提醒：如果频带跨在 $k \frac{f _{s}}{2 R}$ 的边界上，向下采样后负频带的”复制”会混进正频带，即使采样速率足够高也会出现混叠，如图 5-3 所示。

图 5-3 完整频带的干扰(© VDI 出版社文献[4])

插值与 D/A 转换

反过来，提高采样速率也有用武之地，最典型的就是 D/A 转换。D/A 转换器通常采用一阶采样-保持，输出是阶梯状的函数，其高频镜像成分虽然可以用模拟 $1/ sinc (x)$ 补偿滤波器修整，但在大多数情况下数字方案更划算：在数字域先用一个扩展器（每 R 个点插入 R-1 个零）和一个附加低通滤波器，把频谱抬到高采样速率上。注意，扩展器引入的零点会在基带频谱之外生成额外副本，这些副本必须在信号送入 D/A 之前被插值滤波器消除，见图 5-4。直观地说，插值倍数越多，输出信号越光滑（参见图 5-5(b)）。

图 5-4 插值示例 $x [n] \circ ∙ X (ω)$ ，R=3

(a) 低的过采样、高的递降

图 5-5 D/A 转换

5.1.1 Noble 恒等式

处理多级系统时，经常会希望把滤波器和向下采样器/扩展器换个位置，以换取更低的运算量。图 5-6 给出的重新排列规则就是所谓的 “Noble” 关系式。

图 5-6 等价多级系统(Noble 关系)

对抽取器，恒等式为

(↓ R) F (z) = F (z^{R}) (↓ R) (5-2)

含义：如果允许”先向下采样、后滤波”，那么滤波器可以换成 $F (z^{R})$ ——它的系数间隔被拉开了 R 倍，等效于滤波器长度不变但只以 1/R 的速率做乘加。对插值器则有

F (z) (↑ R) = (↑ R) F (z^{R}) (5-3)

含义：把滤波器放在扩展器之前，同样得到降低了 R 次的滤波器。直觉上， Noble 恒等式的价值在于：它把”高采样率下的大量运算”合法地搬到”低采样率”一侧，5.2 节的多相实现正是建立在这两个恒等式之上。

5.1.2 用有理数因子进行采样速率转换

如果输入、输出速率之比不是整数，就需要有理数变换因子 $R_{1} / R_{2}$ 。标准做法是：先用插值器把采样速率提高 $R_{1}$ 倍，再用抽取器降低 $R_{2}$ 倍。插值和抽取都需要低通滤波器，两个滤波器级联在一起时，效果由通带更窄的那个决定，因此只需实现一个低通滤波器，其通带频率取

f_{p} = min (\frac{π}{R _{1}}, \frac{π}{R _{2}}) (5-4)

图 5-7 上面是插值器与抽取器的级联，下面是合并后的单滤波器结构。5.6 节将进一步讨论该系统的不同设计方案。

图 5-7 非整数抽取系统：上面是插值器和抽取器的级联，下面是最终组合的低通滤波器

5.2 多相分解

为什么要多相分解

在 IIR/FIR 滤波器与滤波器组中实现抽取或插值时，多相分解非常有用。以 FIR 抽取滤波器为例：在图 3-1 的 FIR 结构后接一个 R 倍向下采样器，输出只需要

y [0], y [R], y [2 R], \dots (5-5)

这些时刻的值。也就是说，卷积和 $\sum_{k} f [k] x [n - k]$ 没有必要每个采样点都算！举例来说， $x [0]$ 只需要乘以

f [0], f [R], f [2 R], \dots (5-6)

这些系数；除 $x [0]$ 外的其余输入样本也只需要乘以

x [R], x [2 R], \dots (5-7)

分解方法

据此可以把输入信号分成 R 个独立的序列：

x [n] = \sum_{r = 0}^{R - 1} x_{r} [n] x_{0} [n] = {x [0], x [R], \dots} x_{1} [n] = {x [1], x [R + 1], \dots} ⋮ x_{R - 1} [n] = {x [R - 1], x [2 R - 1], \dots}

同样把滤波器 $f (n)$ 分成 R 个序列：

f [n] = \sum_{r = 0}^{R - 1} f_{r} [n] f_{0} [n] = {f [0], f [R], \dots} f_{1} [n] = {f [1], f [R + 1], \dots} ⋮ f_{R - 1} [n] = {f [R - 1], f [2 R - 1], \dots} .

图 5-8 给出了多相分解实现的抽取器：R 个短滤波器并行工作在低速率上，再合并成输出。这样的抽取器运行速度可以比”滤波器 + 向下采样器”的通用结构快 R 倍。滤波器 $f_{r} [n]$ 称为多相滤波器——它们幅值传递函数相同，只是被不同的采样延迟分隔开，引入了相位偏移，“多相”由此得名。

图 5-8 抽取滤波器的多相实现

例 5.1 多相抽取滤波器

以长度为 4 的 Daubechies 滤波器 $G (z)$ 、抽取因子 R = 2 为例：

G (z) = ((1 + 3) + (3 + 3) z^{- 1} + (3 - 3) z^{- 2} + (1 - 3) z^{- 3}) \frac{1}{4 2} G (z) = 0.48301 + 0.8365 z^{- 1} + 0.2241 z^{- 2} - 0.1294 z^{- 3}

系数量化到 8 位精度后得到

G (z) = (124 + 214 z^{- 1} + 57 z^{- 2} - 33 z^{- 3}) /256 G (z) = G_{0} (z^{2}) + z^{- 1} G_{1} (z^{2}) = G_{0} (z^{2}) (\frac{124}{256} + \frac{57}{256} z^{- 2}) + z^{- 1} G_{1} (z^{2}) (\frac{214}{256} - \frac{33}{256} z^{- 2})

即两个多相子滤波器为

G_{0} (z) = \frac{124}{256} + \frac{57}{256} z^{- 1}, G_{1} (z) = \frac{214}{256} - \frac{33}{256} z^{- 1} (5-8)

实现思路分三步：(1) 用一个小状态机把输入流按采样节拍拆成偶数样本和奇数样本，并产生一个二分频时钟；(2) 用 RAG（简化加法器图）方式计算各系数乘法——例如 $124 = 4 (32 - 1)$ 、 $214 = 2 \times 107$ 、 $107 = 99 + 8 = 3 \times 33 + 8$ 、 $57 = 8 (8 - 1) + 1$ 、 $33 = 32 + 1$ ，全部用移位和加法实现，不占乘法器；(3) 在转置结构中安排两个子滤波器，在二分频时钟的下降沿更新寄存器并把两个子滤波器的输出相加。原书 VHDL 改写为 Verilog 如下：

// 原书 VHDL 改写为 Verilog：DB4 多相抽取滤波器 db4poly
module db4poly (
    input  wire               clk,     // 系统时钟
    input  wire               reset,   // 异步复位
    input  wire signed [7:0]  x_in,    // 系统输入
    output reg                clk2,    // 二分频时钟
    output wire signed [16:0] x_e,     // 偶数样本（测试观测）
    output wire signed [16:0] x_o,     // 奇数样本（测试观测）
    output wire signed [16:0] g0,      // 多相滤波器 0
    output wire signed [16:0] g1,      // 多相滤波器 1
    output wire signed [8:0]  y_out    // 系统输出
);
    localparam EVEN = 1'b0, ODD = 1'b1;
    reg               state;
    reg signed [7:0]  x_even, x_odd, x_wait;
    reg signed [16:0] r0, r1, r2, r3, y;
 
    // RAG 乘法辅助量：33 * x_odd
    wire signed [16:0] x33  = (x_odd <<< 5) + x_odd;          // 33
    wire signed [16:0] x99  = x33 + (x33 <<< 1);              // 33*3 = 99
    wire signed [16:0] x107 = x99 + (x_odd <<< 3);            // 99+8 = 107
    // 多相系数与输入的乘积（移位加实现，无乘法器）
    wire signed [16:0] m0 = (x_even <<< 7) - (x_even <<< 2);  // 127 * x_even
    wire signed [16:0] m1 = x107 <<< 1;                       // 214 * x_odd
    wire signed [16:0] m2 = (x_even <<< 6) - (x_even <<< 3)
                            + x_even;                          // 57 * x_even
    wire signed [16:0] m3 = x33;                              // 用于 -33 * x_odd
 
    // FSM：按采样速率把输入流拆成偶数/奇数样本
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            state <= EVEN; clk2 <= 1'b0;
            x_even <= 8'sd0; x_odd <= 8'sd0; x_wait <= 8'sd0;
        end else begin
            case (state)
                EVEN: begin
                    x_even <= x_in;
                    x_odd  <= x_wait;
                    clk2   <= 1'b1;
                    state  <= ODD;
                end
                ODD: begin
                    x_wait <= x_in;
                    clk2   <= 1'b0;
                    state  <= EVEN;
                end
            endcase
        end
    end
 
    // 多相滤波器：转置结构，在二分频时钟下降沿更新
    always @(negedge clk2 or posedge reset) begin
        if (reset) begin
            r0 <= 17'sd0; r1 <= 17'sd0;
            r2 <= 17'sd0; r3 <= 17'sd0;
            y  <= 17'sd0;
        end else begin
            r0 <= r2 + m0;   // G0: g[0] = 127
            r2 <= m2;        //     g[2] = 57
            r1 <= m1 - r3;   // G1: g[1] = 214
            r3 <= m3;        //     g[3] = -33
            y  <= r0 + r1;   // 多相分量相加
        end
    end
 
    assign x_e = {{9{x_even[7]}}, x_even};
    assign x_o = {{9{x_odd[7]}},  x_odd};
    assign g0  = r0;
    assign g1  = r1;
    assign y_out = y[16:8];  // y / 256，恢复滤波器比例
endmodule

设计要点说明：输出虽然是按 1/256 归一化的，但中间仍存在潜在增长——总量 $Σ∣ g_{k} ∣ = 1.673 < 2^{1}$ ，所以给输出保留了额外的保护位。本设计使用 167 个 LE，不占用嵌入式乘法器，TimeQuest 缓慢 85C 模型下时序性能 Fmax=618.43MHz。

图 5-9 给出了仿真结果。前 4 个输入样本是一个三角函数，用来验证奇/偶样本的分离；采用幅值为 100 的脉冲可以验证两个多相滤波器的系数。特别注意：这里的滤波器不再是移不变的。

图 5-9 长度为 4 的 Daubechies 滤波器多相实现的仿真结果

从图 5-9 可以看出，多相抽取器已经不是移不变系统，而是一个周期时变（技术上的非线性）系统。用单脉冲即可验证：脉冲落在偶数索引的采样时刻，响应是 $G_{0} (z)$ ；落在奇数索引时刻，响应是 $G_{1} (z)$ 。

5.2.1 递归 IIR 抽取器

多相分解同样适用于递归滤波器，并且对速度十分有利。按照 Martinez 和 Parks 的构想，可以采用如下传递函数：

F (z) = \frac{\sum _{l = 0}^{L - 1} a [ l ] z ^{- l}}{1 - \sum _{l = 1}^{K - 1} b [ l ] z ^{- l R}} (5-9)

关键在于：递归部分只保留各自的第 R 个系数，于是反馈回路工作在降低了 R 倍的速率上。前面 IIR 滤波器一章（图 4-17）已经讨论过这类设计。图 5-10 表明，与 FIR 抽取器相比，IIR 抽取器依靠滤波器过渡带宽 $Δ F$ 的加大，可以取得实质上的计算量节省——过渡带越窄，FIR 需要的阶数急剧上升，而 IIR 的优势越明显。

图 5-10 抽取器 $Δ F = f_{p} - f_{s}$ 计算量之比较(采样频率 $2 \times f_{n} = 0.1$ )

5.2.2 快速 FIR 滤波器

多相分解还有一个有趣的应用：快速 FIR 滤波器。其基本构想是：把输入信号分解成 R 个多相分量后，卷积就变成了短多项式乘法，可以套用 Winograd 的短卷积算法来减少乘法次数。下面用 R = 2 说明。

例 5.2 快速 FIR 滤波器

把输入 $X (z)$ 和滤波器 $F (z)$ 分解成偶、奇两个多相分量：

X (z) = n \sum x [n] z^{- n} = X_{0} (z^{2}) + z^{- 1} X_{1} (z^{2}) (5-10)

F (z) = n \sum f [n] z^{- n} = F_{0} (z^{2}) + z^{- 1} F_{1} (z^{2}) (5-11)

时域卷积对应 z 域多项式乘法，输出可写为

Y (z) = Y_{0} (z^{2}) + z^{- 1} Y_{1} (z^{2}) (5-12)

= (X_{0} (z^{2}) + z^{- 1} X_{1} (z^{2})) (F_{0} (z^{2}) + z^{- 1} F_{1} (z^{2})) (5-13)

把式(5-13)按多相分量拆开，得到

Y_{0} (z) = X_{0} (z) F_{0} (z) + z^{- 1} X_{1} (z) F_{1} (z) (5-14)

Y_{1} (z) = X_{1} (z) F_{0} (z) + X_{0} (z) F_{1} (z) (5-15)

把它与线性 $2 \times 2$ 卷积

A (z) B (z) = (a [0] + z^{- 1} a [1]) (b [0] + z^{- 1} b [1]) (5-16)

= a [0] b [0] + z^{- 1} (a [0] b [1] + a [1] b [0]) + a [1] b [1] z^{- 2} (5-17)

对比可以发现： $z^{- 1}$ 项的因子完全相同，只是 $Y_{0} (z)$ 需要额外引入一个时延来得到正确的相位关系。Winograd 编制的短卷积算法用 3 次乘法和 6 次加法即可完成线性 $2 \times 2$ 卷积：

a [0] = x [0] - x [1] b [0] = f [0] - f [1] c [k] = a [k] b [k] y [0] = c [1] + c [2] a [1] = x [0] b [1] = f [0] k = 0, 1, 2 y [1] = c [1] - c [0] a [2] = x [1] - x [0] b [2] = f [1] - f [0] (5-18)

借助该算法，快速滤波器可写成矩阵形式：

[Y_{0} Y_{1}] = [0 - 1 11 - 1 0] F_{0} 00 0 F_{0} + F_{1} 0 00 F_{1} 111 - 1 0 - z^{- 1} [X_{0} X_{1}] (5-19)

图 5-11 给出了图形解释。

图 5-11 $R = 2$ 的快速 FIR 滤波器

运算量账本：快在”次数”，贵在”硬件”

比较直接实现与快速滤波器时，必须区分”每样本平均运算次数”和”硬件资源”两本账：

直接实现：L 个乘法器和 L-1 个加法器全速运行。
快速滤波器：3 个长度 L/2 的滤波器半速运行。每输出样本平均 3L/4 次乘法和 $(2 + 2) /2 + 3/2 (L /2 - 1) = 3 L /4 + 1/2$ 次加法——比直接实现省约 25% 的运算量。
硬件代价：需要 3L/2 个乘法器和 $4 + 3 (L /2 - 1) = 3 L /2 + 1$ 个加法器——比直接实现多约 50% 的资源。

快速滤波器最重要特征是运行速度基本上是直接实现的两倍。用更高的分解次数 R 还能进一步提高吞吐量。处理 R 个多相信号（输入速率 $f_{a}$ ）的一般流程如下：

算法 5.3 快速 FIR 滤波器

(1) 将输入信号分解成 R 个多相信号，用 $A_{e}$ 个加法器以 $f_{a} / R$ 的速率构成 R 个序列。

(2) 用 R 个长度为 L/R 的滤波器分别对 R 个序列滤波。

(3) 用 $A_{e}$ 次加法计算输出 $Y_{k} (z)$ 的多相表示，再用多路复用器生成输出 $Y (z)$ 。

什么时候停止分解

注意：已经得到的长度 L/R 的部分滤波器可以用算法 5.3 再次分解。那么迭代到什么时候停？Mou 和 Duhamel 编制了以最小化平均运算量为目标的优化表（准则：乘加总次数最少，这是 MAC 设计的典型代表），见表 5-1，表中基于算法 5.3 实现的部分滤波器的因子已加下划线。

表 5-1 递归 FIR 分解的计算量

L	因子	M+A	(M+A)/L	L	因子	M+A	(M+A)/L
2	直接	6	3	22	$11 \times 2$	668	30.4
3	直接	15	5	24	$2^{2} \times 3 \times 2$	624	26
4	$2 \times 2$	26	6.5	25	$5 \times 5$	740	29.6
5	直接	45	9	26	$13 \times 2$	750	28.6
6	$3 \times 2$	56	9.33	27	$3^{2} \times 3$	810	30
8	$2^{2} \times 2$	94	11.75	30	$5 \times 3 \times 2$	912	30.4
9	$3 \times 3$	120	13.33	32	$2^{4} \times 2$	1006	31.44
10	$5 \times 2$	152	15.2	33	$11 \times 3$	1248	37.8
12	$2 \times 3 \times 2$	192	16	35	$7 \times 5$	1405	40.1
14	$7 \times 2$	310	22.1	36	$2^{2} \times 3 \times 3$	1260	35
15	$5 \times 3$	300	20	39	$13 \times 3$	1419	36.4
16	$2^{3} \times 2$	314	19.63	55	$11 \times 5$	2900	52.7
18	$2 \times 3 \times 3$	396	22	60	$5 \times 2 \times 3 \times 2$	2784	46.4
20	$5 \times 2 \times 2$	472	23.6	65	$13 \times 5$	3345	51.46
21	$7 \times 3$	591	28.1

从表中可以看出一个经验规律：当滤波器长度大于 60 时，改用 FFT 快速卷积会更加有效（第 6 章讨论）。

5.3 Hogenauer CIC 滤波器

对于高抽取速率的滤波器，由 Hogenauer 引入的”级联积分器梳状”（Cascade Integrator Comb, CIC）滤波器是一种极其高效的体系结构，也称 Hogenauer 滤波器。它的典型应用包括：无线通信中把 RF/IF 采样速率的信号降到基带——在蜂窝通信这类窄带应用中抽取比经常超过 1000，这样的系统称为通道器；另一个领域是 ΣΔ 数据转换器。

CIC 实现的基础是完美的极点/零点抵消，而要做到”完美”，必须使用无误差的算术。二进制补码（模 $2^{b}$ 运算）和余数系统（模 $M$ 运算）都具备这种能力。

5.3.1 单级 CIC 案例研究

图 5-12 是一个 4 位运算、无抽取的一阶 CIC 滤波器：一个（递归的）积分器 I，后接一个微分器/梳状部分 C。数值边界是 $- 8_{10} = 100 0_{2 C}$ 和 $7_{10} = 011 1_{2 C}$ 。

图 5-12 4 位运算中的移动均值

图 5-13 给出了脉冲响应。令人意外的是：滤波器虽然是递归的，脉冲响应却是有限的——它是一个”递归 FIR 滤波器”（一般递归滤波器是 IIR）。其脉冲响应说明该滤波器计算的是一个移动和：

y [n] = k = 0 \sum D - 1 x [n - k] (5-20)

CIC 滤波器结构示意（积分器与梳状部分级联）

图 5-13 图 5-12 中滤波器的脉冲响应

移动均值滤波器的非递归等效实现（需 5 个加法器）

其中 D 是梳状部分的延迟。响应就是定义在 D 个连续采样值上的移动均值——一种最简单的低通滤波器。同样的移动均值若用非递归 FIR 实现，需要 D-1 = 5 个加法器；而 CIC 设计只需要 1 个加法器和 1 个减法器，节省非常可观。

本征频率测试：积分器溢出为何无害

递归滤波器在本征频率输入（信号频率恰好对准滤波器极点）下有最大的稳态正弦输出。CIC 的本征频率对应 $ω = 0$ ，即阶跃输入。式(5-20)的一阶移动均值对阶跃的响应是：前 D 个采样为斜坡，之后保持常数 $y [n] = D = 6$ ，如图 5-14 所示。注意图中积分器 $w [n]$ 频繁溢出，但输出依然正确——因为梳状减法同样使用二进制补码。举例验证：第一次回绕时积分器实际值为 $w [n] = - 8_{10} = 100 0_{2 C}$ ，延迟值为 $w [n - 6] = 2_{10} = 001 0_{2 C}$ ，则

y [n] = - 8_{10} - 2_{10} = 100 0_{2 C} - 001 0_{2 C} = 011 0_{2 C} = 6_{10}

正是期望值。累加器会一直计数、回绕、再计数；只要输出 $y [n]$ 是 [-8, 7] 范围内的有效 4 位补码数，补码系统的精确模运算就会自动补偿积分器的溢出。

积分器信号 w[n] 的回绕过程示意

图 5-14 图 5-12 中滤波器的阶跃响应(本征频率测试)

阶跃输入下的输出 y[n] 稳定于 D=6

用 RNS 拓宽位宽

4 位对实际应用太小了。例如 Harris IC HSP43220 有 5 级、使用 66 位宽的积分器。为了降低加法器延迟，可以采用多基 RNS（余数系统）：取模数集合 $Z_{30} = {2, 3, 5}$ ，由表 5-2 可知总共可表示 $2 \times 3 \times 5 = 30$ 个唯一值。该映射是唯一的（双射的），可由中国余数定理证明。

表 5-2 集合 ${2, 3, 5}$ 的 RNS 映射

a=	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
a mod 2	1	0	1	0	1	0	1	0	1	0	1	0	1	0	1
a mod 3	1	2	0	1	2	0	1	2	0	1	2	0	1	2	0
a mod 5	1	2	3	4	0	1	2	3	4	0	1	2	3	4	0

a=	16	17	18	19	20	21	22	23	24	25	26	27	28	29
a mod 2	0	1	0	1	0	1	0	1	0	1	0	1	0	1
a mod 3	1	2	0	1	2	0	1	2	0	1	2	0	1	2
a mod 5	1	2	3	4	0	1	2	3	4	0	1	2	3	4

图 5-15 给出了 RNS 实现下的阶跃响应，其输出用表 5-2 的数据重建，结果与二进制补码情形完全相同（参见图 5-14）。保留结构的映射称为同态，双射同态称为同构（ $≅$ ），表示为

Z_{30} ≅ Z_{2} \times Z_{3} \times Z_{5} (5-21)

图 5-15 RNS 算法中一阶 CIC 的阶跃响应

5.3.2 多级 CIC 滤波器理论

S 级 CIC 系统的传递函数为

F (z) = (\frac{1 - z ^{- R D}}{1 - z ^{- 1}})^{S} (5-22)

其中 D 是梳状部分的延迟数，R 是抽取因子。零极点分析：分子 $(1 - z^{- R D})$ 在 z 平面单位圆上每隔 $2 π / (R D)$ 弧度产生一个零点（共 RD 个，每个重复 S 次）；S 个极点全部位于 z = 1（DC）处，恰好被 CIC 自身的 S 个零点抵消——这正是”完美抵消”，净结果等效于 S 级移动均值滤波器。最大动态范围增长出现在 DC 处：

B_{grow} = (R D)^{S} 或 b_{grow} = lo g_{2} (B_{grow}) 位 (5-23)

设计时必须事先知道这个位增长值，因为 CIC 需要精确运算。实际产品中最坏情况增益可以很大——66 位动态范围的 Harris HSP43220 通道器就是证明，这类设计通常采用二进制补码算法。

图 5-16 给出一个三级 CIC 滤波器：三级积分器 + 抽取器 + 三级梳状部分。注意排列顺序：先实现所有积分器，然后是抽取器，最后是梳状部分。这种重新安排使梳状部分节省了 R 倍的延迟元件（利用 Noble 恒等式把延迟搬到低速率一侧）。高抽取速率下 D 的典型值是 1 或 2。

图 5-16 CIC 滤波器，每级 26 位

位宽算例：8 位输入的三级 CIC 滤波器，D = 2，R = 32，DR = 64，内部字宽需

W = 8 + 3 lo g_{2} (64) = 8 + 3 \times 6 = 26 位

才能保证不发生运行时溢出（此处 $lo g_{2} 64 = 6$ 正是式(5-23)的 $b_{g r o w}$ ）。正常输出字宽一般小得多，例如 10 位。

例 5.4 三级 CIC 抽取器 I

最坏情况增益发生在输入为阶跃（DC）信号时。图 5-17(a) 是幅值 127 的阶跃输入；图 5-17(b) 是第三个积分器的输出——注意运行时溢出以规则速率出现，但这是设计内行为；图 5-17(c) 以输入采样速率显示（插值光滑后）的 CIC 输出；图 5-17(d) 是缩放到 10 位精度、按抽取后速率显示的最终输出。

(a) 幅值 127 的阶跃输入信号

(b) 第三个积分器的输出，溢出规则出现

(c) 按输入采样速率显示的 CIC 输出（光滑）

(d) 缩放到 10 位精度、按抽取速率显示的输出
图 5-17 图 5-16 中三阶 CIC 滤波器的 MATLAB 仿真

原书 VHDL 改写为 Verilog 如下：

// 原书 VHDL 改写为 Verilog：三级 CIC 抽取器 cic3r32（R=32, D=2）
module cic3r32 (
    input  wire        clk,     // 系统时钟
    input  wire        reset,   // 异步复位
    input  wire [7:0]  x_in,    // 系统输入
    output reg         clk2,    // 时钟分频（抽取使能）
    output wire [9:0]  y_out    // 系统输出
);
    reg [4:0] count;                       // 0..31 计数器
    reg       state;                       // 0=hold, 1=sample
    reg signed  [7:0]  x;                  // 输入寄存
    wire signed [25:0] sxtx;               // 符号扩展到 26 位
    reg signed  [25:0] i0, i1, i2;         // 积分器 0,1,2
    reg signed  [25:0] i2d1, i2d2, c0;     // 积分器2的延迟与梳0
    reg signed  [25:0] c1, c1d1, c1d2;     // 梳1及其延迟
    reg signed  [25:0] c2, c2d1, c2d2;     // 梳2及其延迟
    reg signed  [25:0] c3;                 // 梳3（输出）
 
    assign sxtx = {{18{x[7]}}, x};         // 符号扩展 8 -> 26 位
 
    // FSM：每 32 个时钟产生一个 sample 脉冲和 clk2
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            state <= 1'b0; count <= 5'd0; clk2 <= 1'b0;
        end else if (count == 5'd31) begin
            count <= 5'd0; state <= 1'b1; clk2 <= 1'b1;
        end else begin
            count <= count + 5'd1; state <= 1'b0; clk2 <= 1'b0;
        end
    end
 
    // 3 个积分器：工作在全速时钟
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x <= 8'sd0; i0 <= 26'sd0; i1 <= 26'sd0; i2 <= 26'sd0;
        end else begin
            x  <= signed'(x_in);
            i0 <= i0 + sxtx;
            i1 <= i1 + i0;
            i2 <= i2 + i1;
        end
    end
 
    // 3 个梳状部分：每个有 D=2 个采样的延迟，工作在抽取后速率
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            c0 <= 26'sd0; c1 <= 26'sd0; c2 <= 26'sd0; c3 <= 26'sd0;
            i2d1 <= 26'sd0; i2d2 <= 26'sd0;
            c1d1 <= 26'sd0; c1d2 <= 26'sd0;
            c2d1 <= 26'sd0; c2d2 <= 26'sd0;
        end else if (state == 1'b1) begin
            c0   <= i2;          // 第一级梳
            i2d1 <= c0;
            i2d2 <= i2d1;
            c1   <= c0 - i2d2;   // 第二级梳
            c1d1 <= c1;
            c1d2 <= c1d1;
            c2   <= c1 - c1d2;   // 第三级梳
            c2d1 <= c2;
            c2d2 <= c2d1;
            c3   <= c2 - c2d2;
        end
    end
 
    assign y_out = c3[25:16];  // 即 c3 / 2**16，取高 10 位
endmodule

设计包含三部分：一个有限状态机（FSM）负责分频与抽取节拍；一个符号扩展单元；以及两个”算法”模块——积分器模块实现 3 个积分器，梳状模块实现 3 个具有两采样延迟的梳状滤波器。本设计采用 341 个 LE，未用嵌入式乘法器，TimeQuest 缓慢 85C 模型下 Fmax=282.49MHz。值得强调：如果没有抽取器前置，全部梳状部分都要按全速运行并保留延迟寄存器，LE 数会大得多——前置的向下采样节省了 $3 \times 32 \times 26 = 2496$ 个寄存器（LE）。

将仿真输出 y_out（图 5-18）与图 5-17(d) 的 MATLAB 仿真结果对比，可以看到滤波器工作状态完全符合预期。

图 5-18 图 5-16 中显示的三阶 CIC 滤波器的仿真

Hogenauer 进一步指出：经过仔细分析，前级的部分较低有效位（LSB）可以剪除而不影响系统完整性。图 5-19 同时给出了所有级均采用全字宽（最坏情况）的幅频响应，以及应用字长”剪除”策略后的结果。

图 5-19 CIC 传递函数( $f_{s}$ 是输入端的采样频率)

5.3.3 幅值与混叠畸变

S 级 CIC 的传递函数（同式 5-22）沿 $z = e^{j 2 π f T}$ 求值，即可计算幅值响应：

∣ F (f) ∣ = (\frac{sin ( 2 π f T R D /2 )}{sin ( 2 π f T /2 )})^{S} (5-25)

图 5-20 给出 R=3、D=2、RD=6 的三级 CIC 滤波器的 $F (f - k \frac{1}{2 R})$ ：可以看到 CIC 的低频响应在基带内的几个混叠副本——抽取时它们正是潜在的干扰源。

(a) 三级 CIC 抽取器的幅频响应

(b) 基带附近的混叠副本细节
图 5-20 3 级 CIC 抽取器的传递函数，注意 $f_{s}$ 是低速的采样频率

最大混叠分量出现在频率

f_{最大混叠} = 1/ (2 R) - f_{p} (5-26)

处，可由 $∣ F (f) ∣$ 直接计算。多数情况下只需考虑第一个混叠分量，因为第二个已经非常小。图 5-21 给出通带频率处、不同 $f_{p} / (D f_{s})$ 比率下的幅值畸变 $- 20 l g (1 - F (f_{p}))$ ；图 5-22 给出不同 S、R、D 值下的最大混叠分量与 $f_{p} / f_{s}$ 的关系。

图 5-21 CIC 抽取器的幅值畸变 -20lg(1-F(fp))

CIC 抽取器最大混叠与通带比率的关系曲线（D=1 情形）

图 5-22 1 至 4 级 CIC 抽取器的最大混叠（D=2）

设计准则很明确：幅值畸变可以通过级联一个通带内传递函数为 $1/∣ F (z) ∣$ 的 FIR 补偿滤波器来校正；但混叠畸变是不可修复的。因此，可以接受的混叠畸变往往成为 CIC 设计的首要约束。

5.3.4 Hogenauer “剪除”理论

内部总位宽定义为输入字宽与最大动态增长（式 5-23）之和：

B_{intern} = B_{input} + B_{growth} (5-27)

如果每一级都按这个位宽做精确运算，输出端就不会溢出。但通常输入、输出位宽在同一量级——全宽运行意味着大量高位在输出处被舍弃。Hogenauer 的洞察是：与其最后统一截断，不如把”剪除噪声”按各级的功率增益分摊到前级，剪掉那些本来就会被淹没的 LSB。若 $σ_{T, 2 S + 1}^{2}$ 是输出剪除引入的量化噪声，令其等于前面各级剪除噪声之和：

k = 1 \sum 2 S σ_{T, k}^{2} = k = 1 \sum 2 S σ_{k}^{2} P_{k}^{2} ⩽ σ_{T, 2 S + 1}^{2} (5-28)

σ_{T, k}^{2} = \frac{1}{2 S} σ_{T, 2 S + 1}^{2} (5-29)

P_{k}^{2} = n \sum (h_{k} [n])^{2}, k = 1, 2, \dots, 2 S (5-30)

其中 $P_{k}^{2}$ 是第 k 级到输出的功率增益。第 k 级应剪除的位数：

B_{k} = ⌊ 0.5 lo g_{2} (P_{k}^{- 2} \times \frac{6}{N} \times σ_{T, 2 S + 1}^{2}) ⌋ (5-31)

σ_{T, k}^{2}_{j = 2 S + 1} = \frac{1}{12} 2^{2 B_{k}} = \frac{1}{12} 2^{2 (B_{in} - B_{out} + B_{growth})} (5-32)

梳状部分的功率增益 $P_{k}^{2}$ （ $k = S + 1, \dots, 2 S$ ）可用二项式系数计算：

H_{k} (z) = n = 0 \sum 2 S + 1 - k (- 1)^{n} (n 2 S + 1 - k) z^{- k R D}, k = S, S + 1, \dots, 2 S (5-33)

对积分器一侧（ $k = 1, 2, \dots, S$ ）：每个积分器/梳状对生成一个有限的移动均值脉冲响应，第 k 级的等效系统是 $S - k + 1$ 个积分器/梳状对后接 $k - 1$ 个梳状部分。图 5-23 给出了简化 $P_{k}^{2}$ 计算的重新排列方案。

图 5-23 简化 $P_{k}^{2}$ 计算的重新排列方案(© VDI 出版社文献[4])

学习资料 util 目录下的 cic.exe 程序可自动完成剪除计算，产生脉冲响应文件 cicXX.imp 和配置文件 cicXX.dat（XX 为指定标识）。下面的例子解释其结果。

例 5.5 三级 CIC 抽取器 II

设计与例 5.4 总体相同的 CIC 滤波器，这次加入位剪除。参数： $B_{in p u t} = 8$ ， $B_{o u tp u t} = 10$ ，R=32，D=2。位增长为

B_{growth} = ⌈ lo g_{2} (R D^{S}) ⌉ = ⌈ lo g_{2} (6 4^{3}) ⌉ = 3 \times 6 = 18 (5-34)

内部总位宽为

B_{intern} = B_{input} + B_{growth} = 8 + 18 = 26 (5-35)

cic.exe 给出的逐级剪除结果：

-- Program for the design of a CIC decimator.
-- Input bit width Bin = 8
-- Output bit width Bout= 10
-- Number of stages S = 3
-- Decimation factor R = 32
-- COMB delay D = 2
-- Frequency resolution DR = 64
-- Passband freq. ratio P = 8
-- Results of the Design
-- Computed bit width:
-- Maximum bit growth over all stages = 18
-- Maximum bit width including sign Bmax+1 = 26
-- Stage 1 INTEGRATOR. Bit width : 26
-- Stage 2 INTEGRATOR. Bit width : 21
-- Stage 3 INTEGRATOR. Bit width : 16
-- Stage 1 COMB. Bit width : 14
-- Stage 2 COMB. Bit width : 13
-- Stage 3 COMB. Bit width : 12
-- Maximum aliasing component : 0.002135 = 53.41 dB
-- Amplitude distortion : 0.729769 = 2.74 dB

注意剪除的规律：越靠近输出，位宽越窄——第 1 级积分器保留全宽 26 位，到第 3 级梳状只剩 12 位。这些数值也可以从图 5-21、图 5-22 的设计图中查得。与 Hogenauer 论文的表格对照：混叠抑制为 53.4dB（对应梳状时延 D=2 的表 II），通带衰减为 2.74dB（表 I）。两点说明：Hogenauer 与 cic.exe 按 $- 20 lo g_{10} (F (f_{p}))$ 计算幅值畸变（dB），而图 5-21 遵循一般教材惯例按 $- 20 lo g_{10} (1 - F (f_{p}))$ 计算；另外 Hogenauer 的表 I 用梳状时延做了标准化，而 cic.exe 未标准化，查表对比时要注意口径差异。

例5.6 带位剪除的三级CIC抽取器

前面例子中的 CIC 滤波器（例5.4）虽然结构简单，但所有寄存器都保留了完整的理论位宽，其中很多位最终都被丢弃了，白白浪费了逻辑资源。例5.5 已经算出了每一级到底需要多少位——这就是 Hogenauer 提出的”位剪除”（bit pruning）技术。本例把剪除真正落到硬件上：数据与例5.4 相同（8 位输入、抽取率 $R = 32$ 、差分延迟 $D = 2$ 、三级），但每一级积分器和梳状器的位宽都按剪除表逐级收窄。

原书 VHDL 改写为 Verilog（可综合风格）如下：

// 三级CIC抽取器：8位输入，R=32，D=2，S=3，采用位剪除
module cic3s32 (
    input        clk,        // 系统时钟
    input        reset,      // 异步复位
    input  [7:0] x_in,       // 系统输入
    output reg   clk2,       // 时钟分频输出（低速率时钟使能指示）
    output [9:0] y_out       // 系统输出
);
 
    // FSM 状态：hold（等待） / sample（抽样）
    localparam HOLD   = 1'b0,
               SAMPLE = 1'b1;
    reg        state;
    reg  [4:0] count;              // 0..31 循环计数
 
    reg  [7:0]  x;                 // 打一拍的输入
    reg  [25:0] sxtx;              // 符号扩展后的输入
    reg  [25:0] i0;                // 积分器0（完整26位）
    reg  [20:0] i1;                // 积分器1（剪除后21位，即 i0>>5）
    reg  [15:0] i2;                // 积分器2（剪除后16位，即 i1>>5）
    reg  [13:0] i2d1, i2d2, c0, c1;   // 梳状0及其延迟
    reg  [12:0] c1d1, c1d2, c2;       // 梳状1及其延迟
    reg  [11:0] c2d1, c2d2, c3;       // 梳状2及其延迟
 
    // ---------- FSM：分频计数器 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            state <= HOLD;
            count <= 5'd0;
            clk2  <= 1'b0;
        end else if (count == 5'd31) begin
            count <= 5'd0;
            state <= SAMPLE;
            clk2  <= 1'b1;        // 每32个时钟产生一拍"抽样"脉冲
        end else begin
            count <= count + 5'd1;
            state <= HOLD;
            clk2  <= 1'b0;
        end
    end
 
    // ---------- 符号扩展：8位 -> 26位 ----------
    always @(*) begin
        sxtx[7:0] = x;
        sxtx[25:8] = {18{x[7]}};
    end
 
    // ---------- 积分器部分（每个时钟都运行） ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x  <= 8'd0;
            i0 <= 26'd0;
            i1 <= 21'd0;
            i2 <= 16'd0;
        end else begin
            x  <= x_in;
            i0 <= i0 + sxtx;
            i1 <= i1 + i0[25:5];   // 相当于累加 i0/32
            i2 <= i2 + i1[20:5];   // 相当于累加 i1/32
        end
    end
 
    // ---------- 梳状部分（只在低速率"抽样"节拍运行） ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            c0   <= 14'd0; c1   <= 14'd0; c2   <= 12'd0; c3   <= 12'd0;
            i2d1 <= 14'd0; i2d2 <= 14'd0;
            c1d1 <= 13'd0; c1d2 <= 13'd0;
            c2d1 <= 12'd0; c2d2 <= 12'd0;
        end else if (state == SAMPLE) begin
            c0   <= i2[15:2];      // 相当于 i2/4
            i2d1 <= c0;
            i2d2 <= i2d1;
            c1   <= c0 - i2d2;
            c1d1 <= c1[13:1];      // 相当于 c1/2
            c1d2 <= c1d1;
            c2   <= c1[13:1] - c1d2;
            c2d1 <= c2[12:1];      // 相当于 c2/2
            c2d2 <= c2d1;
            c3   <= c2[12:1] - c2d2;
        end
    end
 
    assign y_out = c3[11:2];       // 相当于 c3/4
 
endmodule

这个设计与例5.4 中未缩放的 CIC 具有完全相同的体系结构：一个有限状态机（FSM）负责把 32 分频的时钟使能送给梳状部分，一个符号扩展单元把 8 位输入扩展到 26 位，然后是三个积分器（Int 部分）和三个各带两级延迟的梳状器（Comb 部分）。区别只在于：所有积分器和梳状器都按剪除表配了”刚好够用”的位宽——比如积分器逐级从 26 位缩到 21 位再到 16 位，输出最终只取 10 位。

这样做有什么好处？与例5.4 相比，虽然速度几乎没有变化（282.49 MHz 相比 290.02 MHz），但逻辑单元（LE）数量节省了大约 30%，整个设计只需要 209 个 LE。不用会怎样？大抽取率下 CIC 的理论位宽增长非常快，如果不剪除，大量高位只会累计出永远用不到的精度，FPGA 资源就被浪费了。

要注意的是，剪除会引入额外的舍入/截位噪声：对比图 5-24 和图 5-18 中仿真输出的最低有效位，可以看到剪除式设计中”噪声”表现为 LSB 附近的渐近线形式（例如 507↔508 之间的抖动）。如果对 LSB 精度有严格要求，需要回查剪除表并适当加宽个别级。

图 5-24 采用位剪除技术实现的三级 CIC 滤波器的 VHDL 仿真

最后提一句：CIC 插值器（上采样方向）同样可以套用剪除技术，只是剪除方向相反——位宽在靠近输入端最宽、靠近输出端最窄，具体设计在练习 5.24 中展开。

5.3.5 CIC 的 RNS 设计

前面所有设计都是二进制补码（2C）体系。还有另一条路：余数系统（Residue Number System, RNS）。RNS 把一个大数拆成若干个小模数上的余数，加减法在各个模通道里独立并行进行，字长短、进位链短，因此速度可以更高。

Garcia 等人用 RNS 实现了与前面相同规格的三级 CIC 滤波器：8 位输入、10 位输出、 $D = 2$ 、 $R = 32$ ，内部最大字宽 26 位。26 位补码的动态范围怎么用 RNS 覆盖？答案是 4-模集合 ${256, 63, 61, 59}$ ：其中一个 8 位二进制补码通道加上三个 6 位模通道，乘积 $256 \times 63 \times 61 \times 59$ 足以覆盖 26 位动态范围（请参阅图 5-25）。

剩下的难点在输出端：RNS 通道里算出来的结果要变回普通二进制，而且通常还要顺便完成缩放（除以 $2^{18}$ ）。可选的缩放方案有三种：

ε-CRT（扩展的中国剩余定理）：直接把 4 个余数合成 10 位缩放输出，需要 8 个 ROM 表和 3 个二进制补码加法器（若把乘逆 ROM 与 ε-CRT 合并组合，可用 7 个表）；
BRS + ε-CRT（基准移除缩放，Base Removal Scaling）：先用两个 6 位模做 BRS 前期缩放、降低位宽，再对余下两个模做 ε-CRT，如图 5-26 所示；
ROM 组合的 BRS ε-CRT：把前两种的查表操作合并。

图 5-25 CIC 滤波器，基准移除缩放（BRS）的详尽设计

图 5-26 BRS 和 $ε$ -CRT 转换步骤

三种方案在 FLEX10K 器件上的速度与资源对比如表 5-3 所示。

表 5-3 以 MSPS 计算的速度和 3 种缩放模式所使用的 LE 和 EAB 的数量

类型	ε-CRT	BRS ε-CRT（只适用于 BRS $m_{4}$ 的速度数据）	与 ROM 组合的 BRS ε-CRT
MSPS	58.8	70.4	58.8
le	34	87	87
table(EAB)	8	9	7

怎么理解这张表？方案 1 和 3 的速度都降到了 58.8 MSPS，这是 10 位 ε-CRT 造成的。但注意：这并不拖累整个系统的速度，因为缩放发生在较低的输出采样速率上，可以在时间上”摊开”完成。方案 2（BRS ε-CRT）更快，达到 70.4 MSPS，原因是：只有 BRS $m_{4}$ 部分必须跟得上输入采样速率，BRS $m_{3}$ 和 ε-CRT 都只需在输出采样速率下运行。ROM 的早期使用会把可能的吞吐量从 76.3 MSPS 压到 70.4 MSPS（后者正是 BRS $m_{4}$ 的最大速度），输出部分则用高效的 ε-CRT 收尾。

表 5-4 总结了 FLEX10K 上三种滤波器本体（不含缩放部分）的实现代价。

表 5-4 在 FLEX10K 器件上实现的 3 种滤波器设计方案

类型	2C 26位	RNS 8位、6位、6位、6位	详细位宽 RNS 设计
MSPS	49.3	76.3	70.4
le	343	559	559

可以看出：纯 RNS（均匀模集合）速度最快（76.3 MSPS）但 LE 最多；如果结合剪除思想给每个通道也做”详细位宽”的 RNS 设计，就能在 70.4 MSPS 的速度下把资源降回来。

5.3.6 CIC 补偿滤波器设计

CIC 是用最低的成本实现高抽取/插值率的好办法，但天下没有免费的午餐：正如图 5-21 所示，它的通带不平坦（幅度随频率明显下垂），过渡带也比较宽。如果系统对通带平坦度有要求怎么办？工程上的标准做法是在 CIC 之后（抽取后的低速率下）级联一个 FIR 补偿滤波器，让它的幅频响应正好”托起”CIC 的下垂部分。

这类补偿滤波器在商业芯片中很常见：

DSP56ADC16：4 级 CIC + 可选 16~255 抽头的 FIR 补偿滤波器。当 $f_{s} / f_{p} = 8$ 时，长度 255 的 FIR 允许再做 4 倍抽取，总抽取率 $R_{1} \times R_{2} = 64$ ，单个 MAC 单元即可轻松实现；
HSP50214 下采样器：I/Q 分离器 + 5 级 CIC（ $R_{1} = 4 \dots 32$ ）+ 1~~5 个半带滤波器（ $R_{2} = 1 \dots 32$ ）+ 255 抽头 CIC 补偿 FIR（ $R_{3} = 1 \dots 16$ ），总抽取率 4~~16384；
GC4114 插值器（插值方向）：4 级插值，向上采样 8~16384，内置 31 个常系数的 CFIR 补偿滤波器。

补偿滤波器本身的设计并不复杂，思路只有三步：

写出 CIC 滤波器的幅频响应 $∣ H_{C I C} (f) ∣$ ；
在需要补偿的频率范围内令补偿滤波器等于其倒数 $1/∣ H_{C I C} (f) ∣$ ，其余频率置零；
用逆 DFT（详见 3.12 节）或 MATLAB 的 fir2 函数把这条目标曲线”变成”FIR 系数。

设计时有一个关键折中：通带与采样频率之比 $f_{s} / f_{p}$ 怎么选？选得太大（通带占比小），补偿目标曲线平缓、滤波器好设计，但滤波器必须在更高的速率下运行、长度也更长；选得太小，CIC 在通带边缘的下降太剧烈，倒数的动态范围太大，补偿滤波器很难做准。一个好的折中是 $f_{s} / f_{p} = 8$ ：此时 CIC 幅度在通带边缘降到 0.727，补偿滤波器在通带边缘需要提供 1.37 倍的提升，而且最终还能还原出 $R = 4$ 的额外抽取余量。例5.5 的规范和 DSP56ADC16 用的都是这个比例。

对于例5.5 的 cic3r32 规范，下面的 MATLAB 脚本计算补偿滤波器系数以及”CIC × 补偿”级联的总传递函数：

close all; clear all
%%% CIC filter parameters from Example 5.5
R = 32; D = 2; S = 3;
L = 255-1; % 滤波器阶数；fir2 要求偶数阶
Fs = R;   % 抽取前的高采样频率(Hz)
Fp = 1/8; % 通带边缘
Fo = R*Fp/Fs; % 归一化截止频率
cic=[1];
for k=1:S    % 通过 S 次卷积计算 CIC 冲激响应
    cic=conv(cic,ones(1,R*D));
end
 
p = 2^16; % 频谱采样点数
[H,w] = freqz(cic,[1],p,Fs); % 计算 CIC 频谱
N=length(find(w<Fo)); % 取所有小于 Fo 的点
f = [w(1:N)' Fo .5]*2; % 频率采样点
H=abs(H');H=H./H(1); % 幅度归一化
Mf = [1./H(1:N) 0 0]; % 通带内取倒数，其余为0
h = fir2(L,f,Mf); % 生成 FIR 系数
h = h/max(h); % 系数归一化
 
%%% Compare and plot the results
figure
[H2,w2] = freqz(h,[1],p,Fs/R); % 补偿滤波器，速率 1/R
H2=abs(H2);H2=H2./H2(1);
H1d=H(1:p/R)'; H2d=H2(1:R:p); % 取同样长度并缩放
H12=H1d.*H2d; % 级联总响应
w12=w(1:p/R);
subplot(211)
plot(f/2*Fs/R,Mf,'k:',w2,H2,'k-.',w(1:p),H(1:p),...
'k--',w12,H12,'k-')
legend('Desired','Comp','CIC','CIC*Comp',-1);title('(a)')
ylabel('F(\omega)');xlabel('Normalized frequency');
axis([0 .3 0.001 1.5]);
subplot(223)
plot(f/2*Fs/R,20*log10(Mf),'k:',w2,20*log10(H2),'k-.',...
    w,20*log10(H),'k--',w12,20*log10(H12),'k-')
axis([0 .3 -70 10]);grid on;title('(b)')
ylabel('F(\omega) in db');xlabel('Normalized frequency');
subplot(224)
stairs(h); title('(c)')%% 绘制滤波器系数
axis([1 L+1 -.4 1.1]);ylabel('f[n]');xlabel('Sample n')
print -deps cicomp.eps; print -djpeg90 cicomp.jpg

读这段脚本时抓住四个要点就够了。第一，开头是例5.5 的规范参数；CIC 的冲激响应通过 $S$ 次”与全1序列卷积”得到（每级 CIC 的冲激响应就是长度 $R \cdot D$ 的矩形窗，卷积 $S$ 次就是 $S$ 级级联）。第二，CIC 频响既可以用 Altera 应用说明中的式(5.25)闭式计算，也可以直接交给 freqz。第三，fir2 按频率点—幅度值对生成 FIR 系数。第四，图上画了四条曲线：期望的补偿目标曲线、补偿滤波器的实际响应、CIC 本身的响应、以及两者级联的结果。

图 5-27 给出了长度 255 的设计结果。级联响应（CIC×Comp）在通带内非常平坦——这正是补偿的目的。注意一个细节：由于抽取后采样率变为 $1/8 = 0.125$ ，几条曲线的频率比例不同，读图时要换算。图 5-27(b) 用对数坐标更清楚地展示了阻带衰减。

图 5-27 滤波器长度为 255 的 CIC 补偿滤波器的设计

从图 5-27 还能看出：期望曲线与实际补偿曲线相当接近，因为 255 个系数足够长。如果缩短滤波器，过渡带就会变宽、逼近精度下降——图 5-28 就是长度 31 的结果，与长度 255 的版本对比一目了然。这提示我们：补偿滤波器的长度本质上是在”通带平坦精度”与”实现代价”之间做交易。

图 5-28 滤波器长度为 31 的 CIC 补偿滤波器的设计

5.4 多级抽取器

当抽取率 $R$ 很大时，单级实现要付出高昂代价。更聪明的做法是多级级联： $S$ 级中每级抽取 $R_{k}$ 倍，总抽取率为

R = R_{1} R_{2} \dots R_{S}

为什么这样更省？因为每一级只需要处理”本级”的带宽要求：前几级采样率高但只需粗滤波（便宜的 CIC），后几级采样率低了才用精滤波。这和前面 CIC + 补偿 FIR 的思想一脉相承。

但多级也有代价：通带不完美。每一级的纹波会逐级累积。为了满足整体通带偏差指标 $ε_{p}$ ，最坏情况下必须把每级的指标收紧为 $ε_{p}^{'} = ε_{p} / S$ 。这其实是一个非常悲观的假设——它假定所有短滤波器恰好在同一频率上同时出现最大纹波，一般情况下过于苛刻。更合理的做法是：先按接近 $ε_{p}$ 的初始值尝试，发现不达标时再有选择地收紧个别级。

Goodman-Carey 半带滤波器设计法

Goodman 和 Carey 提出了一套基于 CIC 和半带滤波器的多级设计方法。所谓半带滤波器，是通带边缘和阻带边缘对称地位于 $ω_{s} = ω_{p} = π /2$ （基带正中间）的滤波器，天然适合用因子 2 改变采样速率。它的冲激响应有一个惊人特性：所有偶数下标系数（中心抽头除外）全部为零。

定义 5.7（半带滤波器）：冲激响应满足

f [k] = 0, k 为偶数，但不包括 k = 0 (5-36)

转换到 Z 域即为

F (z) + F (- z) = c (5-37)

其中 $c \in C$ 。对因果半带滤波器（考虑 $z^{- d}$ 延迟），条件变为

F (z) - F (- z) = c z^{- d} (5-38)

偶系数全零意味着一半乘法器直接消失，这正是半带滤波器在 2 倍速率变换中备受欢迎的原因。

Goodman 和 Carey 编制了一张完整的半带滤波器系数表（表 5-5），所有系数以中心抽头位置 $d = 0$ 标定。表中的 F1 就是移动均值滤波器——也就是 Hogenauer 的 CIC 滤波器，因此第一级可以改用 CIC 把速率变为 2 以外的因数。图 5-29 画出了这 9 个滤波器的传递函数。注意在对数图中几乎看不到对称点，这是多数半带滤波器的共同特征。

表 5-5 Goodman 和 Carey 给出的半带滤波器 F1 至 F9 的中心系数

名称	L	纹波	f[0]	f[1]	f[3]	f[5]	f[7]	f[9]
F1	3	—	1	1
F2	3	—	2	1
F3	7	—	16	9	-1
F4	7	36dB	32	19	-3
F5	11	—	256	150	-25	3
F6	11	49dB	346	208	-44	9
F7	11	77dB	521	302	-53	7
F8	15	65dB	802	490	-116	33	-6
F9	19	78dB	8192	5042	-1277	429	-116	18

频率(角度)

图 5-29 半带滤波器 F1 至 F9 的传递函数

设计算法的核心思想很直觉：第一级通带占采样频率的比例很小，可以容忍较大纹波、用最便宜的滤波器；随着逐级抽取，通带占比越来越大，就必须换用畸变更小的滤波器；算法推进到 $R = 2$ 处停止，最后一级（ $R = 2$ 与 $R = 1$ 之间）必须设计一个足够长的半带滤波器。

具体设计时使用图 5-30 的设计图：先算出过采样比率 $R$ 与所需的通带/阻带衰减 $A = A_{p} = A_{s}$ ，然后在图上从 $R$ 出发，沿 $R, R /2, R /4, \dots$ 在同一阻带衰减水平线上依次”领取”所需的滤波器。由于 F4 和 F6 至 F9 在通带中有纹波（请参阅练习 5.8），当使用多个这类滤波器时必须调整 $ε_{p}$ ：

A = - 20 l g ε_{p} (对于 F1 至 F3 、 F5) (5-39)

A = - 20 l g min (\frac{ε _{p}}{S ^{'}}, ε_{s}) (对于 F4 、 F6 至 F9) (5-40)

其中 $S^{'}$ 是带纹波级数——只有那些通带有纹波的滤波器才需要分摊指标。

图 5-30 Goodman 和 Carey 的设计图

例 5.8 多级半带滤波器抽取器

设计一个抽取器： $R = 160$ 、 $ε_{p} = 0.015$ 、 $ε_{s} = 0.031$ （即 30 dB）。

粗看设计图即可知道总共需要 5 个滤波器，起点定在 $R = 160$ 和 30 dB 处，如图 5-31(a) 所示：

第 1 级：从 160 抽到 32（÷5），用长度 $L = 5$ 的 CIC 滤波器（即 F1 类）；
第 2、3 级：两个 F2 滤波器（各 ÷2）；
第 4 级：一个 F3 滤波器（÷2），到此累计抽取到 $R = 8$ ；
第 5 级：还需处理纹波的最终滤波器，其指标为

A = - 20 l g min (\frac{0.015}{1}, 0.031) = 36.48 dB (5-41)

查图 5-30，36 dB 处 F4 滤波器正好合适。

最后用 Noble 关系式把整个系统合成一个总传递函数

F (z) = F 1 (z) F 2 (z^{5}) F 2 (z^{10}) F 3 (z^{20}) F 4 (z^{40})

（请参阅图 5-6），其通带如图 5-31(b) 所示。

(b) 传递函数 $∣ F (ω) ∣$
图 5-31 Goodman 和 Carey 半带滤波器的设计示例

例5.8 还说明了一个方法论问题：式(5-40)中只统计有纹波的滤波器就足够了。如果悲观地取 $S = 6$ 全部级数分摊，会得到 $A = - 20 l g (0.015/6) = 52 dB$ ，这就得选 F8，付出更多工作量。所以正确的策略是从有利假设出发（只分摊给带纹波的级），事后验证通常证明这个乐观假设是对的。

5.5 作为通频带抽取器的频率采样滤波器

5.3 节的 CIC 滤波器其实属于一个更大的家族——频率采样滤波器（Frequency Sampling Filter, FSF）。FSF 可以充当通道器（channelizer）或抽取滤波器，把信息频谱分解成一组离散子带——多用户通信系统正是这种需求。经典 FSF 的结构是”梳状滤波器 + 一组频率选择谐振器”级联：谐振器各自产生一组极点，有选择地抵消梳状预滤波器产生的零点；对谐振器输出做增益调整，即可改变整个滤波器最终的幅频响应。

FSF 也可以按图 5-32 的方式创建：先级联全极点滤波器部分，再级联全零点（梳状）部分。关键在于选择梳状部分的延迟 $1 \pm z^{- D}$ ，使其零点正好抵消全极点预滤波器的极点（如图 5-33 所示，极点角度 60°、梳状延迟 $D = 6$ 的抵消示例）。无论复数极点出现在哪里，都有一个相应的复数零点将其抵消，最终得到一个全零 FIR 滤波器——因此保有 FIR 的两个优良属性：严格线性相位和常系数组延迟。

图5-32 为节省多级信号处理的延迟因子 $R$ 而将频率采样滤波器级联

图 5-33 极角度 $6 0^{\circ}$ 和梳状延迟 D=6 的极点/零点抵消示例

FSF 吸引多级滤波器组设计者的地方，一是固有低复杂度，二是线性相位。但它依赖精确的极点-零点抵消，这在嵌入式场合如何保证？答案是：用二进制补码（整数环）或 RNS 定义的多项式滤波器。整数运算天然精确，所以极点可以驻留在单位圆边缘——这种”有条件的不稳定”位置之所以可以接受，正是因为抵消有保证。若无此保证，设计者只能把谐振器极点挪进单位圆内部，性能会有损失。反过来，允许极点/零点待在单位圆上还能带来一个额外红利：FSF 的乘法器更少，复杂度更低（代价是数据带宽增加）。

再看图 5-32 的具体结构：第一级（整系数）滤波器部分根据关系式 $2 cos (2 π K / D) = 1$ 、 $0$ 、 $- 1$ 产生位于 $0°$ 和 $180°$ 的极点，第二级则产生 $60°$ 、 $90°$ 、 $120°$ 的极点。表 5-6 列出了系数全为整数、根都在单位圆上、阶数不超过 6 的全部多项式的角频率——这张”构造模块清单”可以高效地搭建 FSF。例如可以构建补码（等价于 RNS 单模）语音处理滤波器组，频率覆盖 900 Hz 到 8000 Hz，采样频率 16 kHz。再在前面加上整系数半带滤波器 HB6（抑畸变）和一个三级免乘法 CIC 滤波器（Hogenauer 滤波器，见 5.3 节）进一步抑制不需要的频率分量，整体结构如图 5-34 所示。

图 5-34 由半带和 CIC 预滤波器以及 FSF 梳状谐振器部分组成的滤波器组的设计

表 5-6 生成不超过 6 阶的唯一角极点位置的整系数滤波器，显示了滤波器系数及单位圆上根的非冗余角度位置

$C_{k} (z)$	阶数	$a_{0}$	$a_{1}$	$a_{2}$	$a_{3}$	$a_{4}$	$a_{5}$	$a_{6}$	$θ_{1}$	$θ_{2}$	$θ_{3}$
$- C_{1} (z)$	1	1	-1						$0°$
$C_{2} (z)$	1	1	1						$180°$
$C_{6} (z)$	2	1	-1	1					$60°$
$C_{4} (z)$	2	1	0	1					$90°$
$C_{3} (z)$	2	1	1	1					$120°$
$C_{12} (z)$	4	1	0	-1	0	1			$30°$	$150°$
$C_{10} (z)$	4	1	-1	1	-1	1			$36°$	$108°$
$C_{8} (z)$	4	1	0	0	0	1			$45°$	$135°$
$C_{5} (z)$	4	1	1	1	1	1			$72°$	$144°$
$C_{16} (z)$	6	1	0	0	-1	0	0	1	$20.00°$	$100.00°$	$140.00°$
$C_{14} (z)$	6	1	-1	1	-1	1	-1	1	$25.71°$	$77.14°$	$128.57°$

每个谐振器的带宽都可以通过梳状部分的级数和延迟数量独立调节；为满足带宽要求，需对级数和延迟数做优化。上表构造的频率选择滤波器均采用两级和两个延迟。

表 5-7 给出了用 Xilinx XC4000 FPGA 原型实现的滤波器组的资源开销（F20D90 的命名含义：极角 20.00°、梳状延迟 $D = 90$ ）。值得注意的是：采用 XBLOCKS 这类高级设计工具时，CLB 实际用量通常比由加法器、触发器、ROM/RAM 统计出的理论值高约 20%。整个滤波器组实际用 1572 个 CLB；对比之下，若做成非递归 FIR 需要 11292 个 CLB——差了近一个数量级，这就是 FSF 复杂度优势的直观体现。

表 5-7 Xilinx XC4000 FPGA 使用的 CLB 数量（总计：实际 1572 个 CLB；非递归 FIR：11292 个 CLB）

	F20D90	F25D70	F36D60	F51D49	F72D40	F90D40	F120D33	F180D14	HB6	III	D4	D5
理论数量	122	184	128	164	124	65	86	35	122	31	24	24
实际数量	160	271	190	240	190	93	120	53	153	36	33	33
非递归 FIR 数量	2256	1836	1924	1140	1039	1287	1260	550

FSF 的设计自由度体现在三处：改变梳状部分的延迟、改变通道幅值、改变梳状部分的数量。其中梳状延迟的调整特别容易实现——CLB 可以用作 $32 \times 1$ 存储单元，用 CLB 实现的计数器就能为特定梳状部分配置延迟。

5.6 任意采样速率转换器的设计

前面讨论的抽取和插值都是有理数倍速率变换。任意速率转换可以先经 $R_{1}$ 倍向上采样、再经 $R_{2}$ 倍向下采样实现（图 5-7 的系统）。设计选项涵盖 IIR 滤波器、FIR 滤波器，直到拉格朗日插值和样条插值。

用一个具体例子说明： $R_{1} = 3$ 的插值后接 $R_{2} = 4$ 的抽取，即变换率 $R = R_{1} / R_{2} = 3/4 = 0.75$ 。

例 5.9 $R = 0.75$ 的速率变换器 I

先 3 倍插值、后 4 倍抽取，中间夹一个中心低通滤波器（系统原理图见图 5-7）。这里有一个关键观察：由于先上采样后下采样，中间滤波器的截止频率只需取两者中的较小者

min (π /3, π /4) = π /4

也就是说只需要一个低通滤波器。MATLAB 中频率已按 $f_{2} /2 = π$ 归一化，设计一个 10 阶、阻带衰减 50 dB 的切比雪夫 II 型滤波器：

[B, A] = cheby2 (10, 50, 0.25)

为什么选切比雪夫 II 型？因为它通带平坦、阻带等纹波、阶数适中，在很多应用里都是好选择。如果担心系数量化敏感度，可以改用双二阶（biquad）级联实现替代直接形式（请参阅例 4.2）：

[SOS, gain] = tf2sos (B, A)

用这个 IIR 滤波器即可仿真有理数速率变换。图 5-35 给出了三角形测试信号的结果：(a) 是初始输入序列，(b) 是经 3 倍上采样再滤波的信号，(c) 是下采样后的信号。

(c) 向下采样后的信号

(a) 原始信号

(b) 经向上采样和滤波的初始信号
图 5-35 有理数速率变换的 IIR 滤波器的仿真结果

结果不算完美：三角形形状保留得不错，但在三角形之后信号本应为零的地方，能看到滤波器的振铃。有人会问：如果用频域里能精确构造的严格低通滤波器，插值能否改进？这就引出了基于 DFT/FFT 的频域插值方法。为保持帧处理简单，取 FFT 长度 $N$ 为速率变换因子的倍数，即 $N = k \times R_{1} \times R_{2}$ ， $k \in N$ 。

算法 5.10 应用 FFT 的有理数速率变换

(1) 取一个 $k \times R_{2}$ 个采样值的块； (2) 在所有采样值之间插入 $(R_{1} - 1)$ 个零点； (3) 计算长度为 $N = k \times R_{1} \times R_{2}$ 的 FFT； (4) 在频域内应用低通滤波运算； (5) 计算长度为 $N$ 的 IFFT； (6) 通过 $R_{1}$ 向下采样得到最终输出，即保留 $k \times R_{2}$ 个采样值。

用一个数值例子走一遍。

例 5.11 $R = 0.75$ 的速率变换器 II

三角形输入序列 $x$ ，取 $k = 1$ ：

(1) 原始块 $x = (1, 2, 3, 4)$ ； (2) 3 倍插值（插零）： $x_{i} = (1, 0, 0, 2, 0, 0, 3, 0, 0, 4, 0, 0)$ ； (3) 计算 FFT： $X_{i} = (10, - 2 + j 2, - 2, 2 - j 2, 10, - 2 + j 2, - 2, - 2 - j 2, 10, - 2 + j 2, - 2, - 2 - j 2)$ ； (4) 频域低通（保留低频段，其余置零）： $X_{lp} = (10, - 2 + j 2, - 2, 0, 0, 0, 0, 0, 0, 0, - 2, - 2 - j 2)$ ； (5) 计算 IFFT： $y = 3 \times ifft (X_{lp})$ （因子 3 补偿插值带来的幅度缩放）； (6) 向下采样，最终 $y = (0.5000, 2.6340, 4.3660)$ 。

把这个块处理方法应用到完整三角形序列上（图 5-36）：(b) 显示无重叠时块与块之间出现明显的边界效应，(d) 则是全长度输入一次 FFT 插值的结果。边界效应的根源是 DFT 的隐含假设——时域和频域的信号都是周期性的，块边界不连续就会产生扰动。改进办法有两种：一是加窗，让边界扰动光滑过渡到零，但会损失有效采样数，还需重叠处理；二是取更长的 FFT（ $k > 1$ ）并丢弃前后的部分采样——图 5-36(c) 是 $k = 2$ 并移除前后 50% 采样的结果，质量已明显改善。

那为什么不干脆用最长的 FFT？看图 5-36(d) 的全长度 FFT 结果确实最好。原因纯粹是计算量：FFT 越长，每个输出采样的平均开销越大。每个基 2 FFT 需要 $ld (N) /2$ 次复数乘法，虽然 $k > 1$ 时每个 FFT 产出更多输出值、FFT 总次数更少，但整体算下来短 FFT 的每采样计算量更低。

(a) 原始信号

(b) 无重叠的抽取

(c) $50%$ 的重叠

(d) 全长度的 FFT
图 5-36 基于有理数速率变换的 FFT

既然长 FFT 质量好、短 FFT 省计算，那么能不能把算法 5.10 中两个长 FFT 本身也做减法？有两个主要的简化措施。

第一，前向变换的稀疏化。 插值序列中大部分是零。用 Cooley-Tukey 时域抽取算法，可以把所有非零值集中到一个 DFT 模块里，其余 $k (R_{1} R_{2} - R_{2})$ 个零值落在另外 $R_{2} - 1$ 组中。于是真正需要计算的 FFT 长度只有 $k \times R_{2}$ ，远小于全长度 $N = k \times R_{1} \times R_{2}$ 。

第二，频域内直接完成向下采样。 时域 2 倍向下采样对应频域中相距 $N /2$ 的两项相加：

F_{↓ 2} (k) = F (k) + F (k + N /2) (5-42)

（因为时域下采样导致基带以 2 倍奈奎斯特频率重复。）推广到 $R_{2}$ 倍向下采样：

F_{↓ R_{2}} (k) = n \sum F (k + n N / R_{2}) (5-43)

再考虑到频域低通已把许多采样置零，式(5-43)的求和次数进一步减少。最终所需 IFFT 的长度仅为 $k \times R_{1}$ 。

设 FFT 和 IFFT 各需 $ld (N) N /2$ 次复数乘法，改进算法的加速比为

F = \frac{2 \frac{k R _{1} R _{2}}{2} ld ( k R _{1} R _{2} )}{\frac{k R _{1}}{2} ld ( k R _{1} ) + \frac{k R _{2}}{2} ld ( k R _{2} )} (5-44)

代入 $R = 3/4$ 和 50% 重叠：

F = \frac{ld ( 24 ) \cdot 24}{ld ( 6 ) \cdot 6/2 + ld ( 8 ) \cdot 8/2} = 5.57 (5-45)

即计算量降到约 1/5.6。如果用 Winograd 短项 DFT 算法（请参阅 6.1.6 节）替代 Cooley-Tukey FFT，收益还会更大。

5.6.1 分数延迟速率变换

有些应用的速率比非常接近 1，比如上面的 $R = 3/4$ 。更极端的例子是把 DAT 采样率（48 kHz）转成 CD 采样率（44.1 kHz），有理数变换因子 $R = 147/160$ 。若按”先插值后抽取”的直来直去方法，插值低通滤波器必须在 147 倍输入采样速率的超高速率下运行——这在硬件上非常昂贵。

$R_{1}$ 很大时，可以换一条思路：用分数延迟实现速率变换。仍以 $R = 3/4$ 说明。图 5-37(a) 画出了输入与输出采样栅格：每 4 个输入值，系统要插算出 3 个输出采样值。从滤波器角度看，这等价于 3 个并行滤波器：一个零延迟（单位传递函数），另外两个延迟分别为 $D = \pm 4/3$ 。理想的分数延迟滤波器是时域的 sinc 函数 $sin (t) / t$ 。当然， sinc 是无限长非因果的，必须附加初始延迟才能实现（因果化）。图 5-37(b) 给出了滤波器配置；图 5-38 画出了长度 11 的滤波器在延迟 $5 - 4/3 = 3. \overline{6}$ 、 $5$ 和 $5 + 4/3 = 6. \overline{3}$ 下的脉冲响应——中心抽头位置随分数延迟平移，这就是”分数延迟”名字的由来。

(a) 分数延迟速率变换的输入和输出采样栅格 (b) $R = 3/4$ 的 sinc 滤波器系统的滤波器配置

图 5-37 输入、输出采样栅格与滤波器配置

(a) 延迟=3.6

(b) 延迟=5

(c) 延迟=6.3
图 5-38 延迟为 D = 5 - 4/3、5 和 $5 + 4/3$ 的分数延迟滤波器

给这三个滤波器送入三角形输入，每组 4 个输入采样产生 3 个输出采样，仿真结果见图 5-39。比较不同长度：长度 5 的滤波器输出纹波太多，没有使用价值（图 5-39(b)）；长度 11 的 sinc 滤波器输出光滑得多（图 5-39(c)），只在信号本应为零处因吉布斯现象残留少量纹波；图 5-39(d) 是拉格朗日插值法的结果，可作对比。

(c) 用长度为 11 的 sinc 滤波器滤波

(d) 用拉格朗日插值法插值
图 5-39 分数延迟插值

例 5.12 $R = 0.75$ 的速率变换器 III

下面看 sinc 分数延迟速率变换器的 HDL 实现。原书 VHDL 改写为 Verilog（可综合风格）：

// R=3/4 速率变换的 sinc 分数延迟滤波器
// 三个并行滤波器：f0 延迟 3.6、f1 延迟 5、f2 延迟 6.3
module rc_sinc #(
    parameter OL = 2,   // 输出缓存长度-1
    parameter IL = 3,   // 输入缓存长度-1
    parameter L  = 10   // 滤波器长度-1
)(
    input              clk,       // 系统时钟
    input              reset,     // 异步复位
    input  signed [7:0] x_in,     // 系统输入
    output reg  [3:0]  count_o,   // FSM 计数器
    output reg         ena_in_o,  // 采样输入使能
    output reg         ena_out_o, // 输出移位使能
    output reg         ena_io_o,  // 传输到输出使能
    output reg  signed [8:0] f0_o, // 第一路 sinc 输出
    output reg  signed [8:0] f1_o, // 第二路 sinc 输出
    output reg  signed [8:0] f2_o, // 第三路 sinc 输出
    output reg  signed [8:0] y_out // 系统输出
);
 
    // 滤波器系数表（已含 256 归一化因子）
    // c0: 延迟 3.6；c1: 延迟 5；c2: 延迟 6.3（c2 为 c0 的镜像）
    reg signed [8:0] c0 [0:10];
    reg signed [8:0] c1 [0:10];
    reg signed [8:0] c2 [0:10];
    initial begin
        c0[0]=-19; c0[1]= 26; c0[2]=-42; c0[3]=106; c0[4]=212;
        c0[5]=-53; c0[6]= 29; c0[7]=-21; c0[8]= 16; c0[9]=-13; c0[10]=11;
        c2[0]= 11; c2[1]=-13; c2[2]= 16; c2[3]=-21; c2[4]= 29;
        c2[5]=-53; c2[6]=212; c2[7]=106; c2[8]=-42; c2[9]= 26; c2[10]=-19;
        // c1 为对称（延迟5）滤波器系数，按同一方法预先计算
    end
 
    reg [3:0] count;                       // 周期 R1*R2 计数
    reg       ena_in, ena_out, ena_io;     // FSM 使能
 
    // ---------- FSM：控制整个系统 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            count <= 4'd0;
            ena_in <= 1'b0; ena_out <= 1'b0; ena_io <= 1'b0;
        end else begin
            if (count == 4'd11) count <= 4'd0;
            else                count <= count + 4'd1;
            case (count)                 // 每 R2=4 拍接受一个输入
                4'd2, 4'd5, 4'd8, 4'd11: ena_in  <= 1'b1;
                default:                 ena_in  <= 1'b0;
            endcase
            case (count)                 // 输出移位节拍
                4'd4, 4'd8:              ena_out <= 1'b1;
                default:                 ena_out <= 1'b0;
            endcase
            ena_io <= (count == 4'd0);   // 成组传输节拍
        end
    end
 
    // ---------- INPUTMUX：输入抽头延迟线 ----------
    reg signed [7:0] ibuf [0:IL];
    integer i;
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= IL; i = i + 1) ibuf[i] <= 8'sd0;
        end else if (ena_in) begin
            for (i = IL; i >= 1; i = i - 1)
                ibuf[i] <= ibuf[i-1];    // 右移一位
            ibuf[0] <= x_in;             // 新样本进入 0 号寄存器
        end
    end
 
    // ---------- OUTPUTMUX：输出缓存 ----------
    reg signed [8:0] obuf [0:OL];
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= OL; i = i + 1) obuf[i] <= 9'sd0;
        end else if (ena_io) begin       // 一次存入 3 个样本
            obuf[0] <= f0_o;
            obuf[1] <= f1_o;
            obuf[2] <= f2_o;
        end else if (ena_out) begin
            for (i = OL; i >= 1; i = i - 1)
                obuf[i] <= obuf[i-1];    // 右移一位
        end
    end
 
    // ---------- TAP：一次接管 4 个样本 ----------
    reg signed [7:0] x [0:10];
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= 10; i = i + 1) x[i] <= 8'sd0;
        end else if (ena_io) begin
            for (i = 0; i <= 3; i = i + 1)
                x[i] <= ibuf[i];         // 接管输入缓存
            for (i = 4; i <= 10; i = i + 1)
                x[i] <= x[i-4];          // 0->4, 4->8 等，一次移 4 抽头
        end
    end
 
    // ---------- SOP0：f0 的乘累加（SOP1/SOP2 结构相同，换系数表 c1/c2） ----------
    reg signed [16:0] sum0;
    reg signed [16:0] p0 [0:10];
    always @(*) begin
        for (i = 0; i <= L; i = i + 1)
            p0[i] = c0[i] * x[i];        // 推断 L+1 个乘法器
        sum0 = p0[0];
        for (i = 1; i <= L; i = i + 1)
            sum0 = sum0 + p0[i];         // 直接形式累加
    end
    always @(posedge clk or posedge reset) begin
        if (reset) f0_o <= 9'sd0;
        else       f0_o <= sum0[16:8];   // 相当于 sum/256
    end
 
    // ---------- 输出选择 ----------
    always @(posedge clk) begin
        count_o   <= count;
        ena_in_o  <= ena_in;
        ena_out_o <= ena_out;
        ena_io_o  <= ena_io;
        y_out     <= obuf[0];
    end
 
endmodule

理解这个设计的钥匙是它的控制结构。整个系统由一个模 12 的计数器（ $R_{1} \times R_{2} = 12$ 拍为一轮）驱动，产生三种使能：ena_in 在每轮的 4 个特定节拍生效，对应输入端每 4 拍进一个新样本；ena_io 在每轮开头生效，负责”成组搬运”——把输入缓存的内容一次性交给抽头寄存器组，同时把三路滤波器结果一次存入输出缓存；ena_out 负责把输出缓存逐拍移位送出。这样，三个并行滤波器（延迟 3.6、5、6.3）共享同一条抽头延迟线，只是各用各的系数表 c0/c1/c2，每轮各算一次乘累加，结果除以 256（即截取 sum[16:8]）做归一化。

与直接法（147 倍超高速插值滤波器）相比，分数延迟方案把速率变换问题转化成了”3 个低速并行 FIR”，每个滤波器都以原始输入速率运行，硬件代价大大降低。这正是它适合 DAT→CD 这类速率比接近 1 的场合的原因。

5.6.1（续）三个 sinc 滤波器的乘累加与仿真结果

先看这个 sinc 速率变换器最核心的部分——三个滤波器的乘累加（Sum-of-Products, SOP）计算。原书 VHDL 改写为 Verilog：

// 原书 VHDL 改写为 Verilog：三个 sinc 滤波器的 SOP 计算与输出连接
parameter L = 10;                       // 滤波器阶数，推断 L+1 个乘法器
reg signed [16:0] f1, f2;
reg signed [16:0] p [0:L];              // 各抽头的乘积
reg signed [16:0] sum;
integer i;
 
// 滤波器 f1：其冲激响应恰为单位冲激，直接取第 5 个抽头，不做缩放
always @(posedge clk or posedge reset) begin
    if (reset) f1 <= 0;                 // 异步清零
    else       f1 <= x[5];
end
 
// 滤波器 f2：L+1 次乘法 + 直接形式累加
always @(posedge clk or posedge reset) begin
    if (reset) begin
        f2 <= 0;
    end else begin
        for (i = 0; i <= L; i = i + 1)
            p[i] <= c2[i] * x[i];       // L+1 个乘法器
        sum = p[0];
        for (i = 1; i <= L; i = i + 1)
            sum = sum + p[i];           // 直接形式滤波器的加法
        f2 <= sum / 256;                // 除以总增益 256
    end
end
 
assign f0_o      = f0;                  // 输出测试信号
assign f1_o      = f1;
assign f2_o      = f2;
assign count_o   = count;
assign ena_in_o  = ena_in;
assign ena_out_o = ena_out;
assign ena_io_o  = ena_io;
assign y_out     = obuf[OL];            // 连接到系统输出

初学者可以按这样的层次来理解整个设计的组织方式：

第一个 always 块是有限状态机（FSM）：它负责整个控制流程，生成输入缓存器、输出缓存器的使能信号，以及三个滤波器共用的使能信号 ena_io。一个完整的工作循环需要 12 个时钟周期。
接下来的三个 always 块分别是输入缓存器、输出缓存器和 TAP（抽头）延迟线。关键的一点是：三个滤波器只共用一条抽头延迟线，这是节省资源的关键设计决策。
最后三个 always 块包含 sinc 滤波器本身。输出 y_out 带有一位附加的保护位，防止累加过程中的溢出。

资源消耗方面，这一设计使用了 880 个 LE（逻辑单元），没有用到嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型时，时序电路性能 Fmax = 59.53 MHz。

图 5-40 给出了滤波器的仿真结果。图中首先显示 FSM 的控制信号和使能信号，输入采用三角形信号 x_in。三个滤波器的输出每 12 个时钟周期才更新一次。滤波器输出值（f0、f1、f2）按正确顺序排列，从而拼装成输出 y_out。仔细观察可以发现：来自 f1 的滤波器值 20 和 60 在输出序列中保持不变，而其他值则是新插入的——这正是分数速率变换”保留原采样、插入新采样”行为的直观体现。

图5-40 采用了三个 sinc 滤波器的 R = 3/4 速率变换的 VHDL 仿真结果

还有一个值得注意的细节：在这个例子中滤波器采用的是直接形式而不是转置形式，因为三个滤波器都只需要一条抽头延迟线，直接形式此时反而更省资源。另外需要说明，这个示例的编码风格更多是从清晰性而非效率角度考虑的。如果滤波器系数采用 MAG 编码，并在滤波器求和的位置增加流水线加法器树，性能还可以进一步提高，请参阅练习 5.15。

5.6.2 多项式分数延迟设计

什么时候该换一种思路

只要延迟数量（也就是速率变换因子 $R = R_{1} / R_{2}$ 的分子 $R_{1}$ ）比较小，用一组低通滤波器来实现分数延迟就很有吸引力。但当 $R_{1}$ 很大时，情况就不同了。一个经典例子是 DAT→CD 的音频采样率转换：所需的变换因子为 $R = 147/160$ ，意味着要用低通滤波器组方案实现 147 个不同的低通滤波器，工作量难以接受。

此时更好的选择是用多项式逼近来计算分数延迟，常用的有拉格朗日（Lagrange）多项式或样条（spline）多项式。 $N$ 点拉格朗日多项式逼近的类型如下：

p (t) = c_{0} + c_{1} t + c_{2} t^{2} + \dots + c_{N - 1} t^{N - 1} (5-46)

通常 3 或 4 点就足够了，不过某些高质量的音响应用会用到多达 10 项。

拉格朗日多项式如何求系数

拉格朗日多项式逼近在区间末端有振荡趋势，因此估计区间应选在多项式的中心。图 5-41 阐述了仅有两个非零值的信号的情形：长度为 4 的多项式已经能在中心区间（ $0 \leq n \leq 1$ ）给出很好的逼近，而从长度 4 提高到长度 8，改善并不明显。逼近区间的较差选择是第一个和最后一个区间——例如对长度为 8 的多项式，若选范围 $3 \leq n \leq 4$ ，振荡和误差都会非常大。

图 5-41 采用短多项式和长多项式的多项式逼近

因此，所用的输入采样集合应关于区间中心对称，使分数延迟落在 $0 \leq d \leq 1$ 内。此时在时间点 $- N /2 - 2, \dots, N /2$ 处使用输入采样。例如对 4 点情形，在 $- 1, 0, 1, 2$ 处取采样。为了拟合一条恰好穿过这些采样点的多项式 $p (t)$ ，把采样时间和 $x (t)$ 的值代入式 (5-46)，就得到一个矩阵方程 $Vc = x$ ，解出系数 $c_{k}$ 。当 $N = 4$ 时：

1111 t_{- 1} t_{0} t_{1} t_{2} t_{- 1}^{2} t_{0}^{2} t_{1}^{2} t_{2}^{2} t_{- 1}^{3} t_{0}^{3} t_{1}^{3} t_{2}^{3} c_{0} c_{1} c_{2} c_{3} = x (n - 1) x (n) x (n + 1) x (n + 2) (5-47)

其中 $t_{k} = k$ 。这里的系数矩阵是范得蒙（Vandermonde）矩阵——DFT 中也常用到这类矩阵。它的每一行都由基本元素的幂级数构造： $t_{k}^{0} = 1, t_{k}, t_{k}^{2}$ ，依此类推。代入 $t_{k}$ 并对矩阵求逆，得到：

c_{0} c_{1} c_{2} c_{3} = 1111 - 1 012 1014 - 1 018^{- 1} x (n - 1) x (n) x (n + 1) x (n + 2) = 0 - 1/3 1/2 - 1/6 1 - 1/2 - 1 1/2 01 1/2 - 1/2 0 - 1/6 0 1/6 x (n - 1) x (n) x (n + 1) x (n + 2) (5-48)

可以动手验算一下第 0 行：逆矩阵第 0 行是 $(0, 1, 0, 0)$ ，乘上采样向量恰好得到 $c_{0} = x (n)$ ——这正对应”多项式必须穿过中心采样点”的要求，说明求逆结果自洽。对每个输出采样，先确定分数延迟值 $d$ ，再求解矩阵方程 (5-48)，最后用式 (5-46) 计算多项式逼近值。图 5-39(d) 给出了使用拉格朗日逼近的仿真结果：它给出了合理精确的逼近；与 sinc 设计相比，在三角波对应值附近、输入值和输出值为零的位置，拉格朗日逼近中几乎没有纹波。

Farrow 结构：用 Horner 模式高效求值

直接按式 (5-46) 计算多项式

p (d) = c_{0} + c_{1} d + c_{2} d^{2} + c_{3} d^{3} (5-49)

需要计算幂值 $d^{k}$ 。改用 $N = 4$ 的 Horner 模式：

p (d) = c_{0} + d (c_{1} + d (c_{2} + c_{3} d)) (5-50)

优点是不需要计算任何幂值——每一层只需要一次乘法和一次加法。由于这一方法首先由 Farrow 提出，文献中将其称为 Farrow 结构。图 5-42(b) 给出了 4 个系数的 Farrow 结构，与直接实现（图 5-42(a)）对比可以看出，Farrow 结构把”求幂”换成了”逐级乘 $d$ 累加”，硬件上就是一条乘加链。

(a) 直接实现

(b) Farrow 结构
图5-42 通过一个 $N = 4$ 的多项式逼近的分数延迟

例5.13 $R = 0.75$ 的速率变换器 IV

下面的设计对 $R = 3/4$ 的速率变换采用三阶拉格朗日多项式的 Farrow 方案。原书 VHDL 改写为 Verilog：

// 原书 VHDL 改写为 Verilog：R = 3/4 三阶拉格朗日多项式的 Farrow 速率变换器
module farrow #(parameter IL = 3) (     // IL：输入缓存长度 - 1
    input  wire              clk,
    input  wire              reset,     // 异步复位
    input  wire signed [7:0] x_in,      // 系统输入
    output wire [3:0]        count_o,   // FSM 计数器
    output wire              ena_in_o,  // 采样输入使能
    output wire              ena_out_o, // 输出移位使能
    output wire signed [8:0] c0_o, c1_o, c2_o, c3_o, // 相位延迟系数
    output wire signed [8:0] d_out,     // 当前使用的分数延迟
    output reg  signed [8:0] y_out);    // 系统输出
 
    localparam signed [8:0] DELTA = 85; // d 的增量：85/256 ≈ 1/3
 
    reg  [3:0]  count;                  // 循环 R1*R2 = 12 个周期
    reg         ena_in, ena_out;
    reg  signed [8:0] d;                // 分数延迟，按 256 定标
    reg  signed [7:0] ibuf [0:IL];      // 抽头延迟线
    reg  signed [7:0] x    [0:IL];      // 4 点采样窗
    reg  signed [8:0] c0, c1, c2, c3;
    reg  signed [16:0] acc;
    integer i;
 
    // FSM：产生使能信号并计算相位延迟 d
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            count <= 0;  d <= DELTA;
            ena_in <= 0; ena_out <= 0;
        end else begin
            count <= (count == 11) ? 0 : count + 1;
            ena_in  <= (count == 2 || count == 5 ||
                        count == 8 || count == 11);
            ena_out <= (count == 3 || count == 7 || count == 11);
            if (ena_out) begin           // 计算相位延迟
                if (d + DELTA >= 255) d <= 0;
                else                  d <= d + DELTA;
            end
        end
    end
 
    // 一条共用的抽头延迟线
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= IL; i = i + 1) ibuf[i] <= 0;
        end else if (ena_in) begin
            for (i = 0; i < IL; i = i + 1) ibuf[i] <= ibuf[i+1];
            ibuf[IL] <= x_in;
        end
    end
 
    // 输出时刻一次性取出 4 个采样
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= IL; i = i + 1) x[i] <= 0;
        end else if (ena_out) begin
            for (i = 0; i <= IL; i = i + 1) x[i] <= ibuf[i];
        end
    end
 
    // 拉格朗日矩阵（式 5-48）+ Farrow 组合器（Horner 模式，d 按 256 定标）
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            c0 <= 0; c1 <= 0; c2 <= 0; c3 <= 0;  y_out <= 0;
        end else if (ena_out) begin
            c0 <= x[1];
            c1 <= -85*x[0]/256 - x[1]/2 + x[2] - 43*x[3]/256;
            c2 <= (x[0] + x[2])/2 - x[1];
            c3 <= (x[1] - x[2])/2 + 43*(x[3] - x[0])/256;
            acc = c2 + (c3 * d)/256;    // y = c2 + c3*d
            acc = (acc * d)/256 + c1;   // y = c1 + y*d
            acc = (acc * d)/256 + c0;   // y = c0 + y*d
            y_out <= acc;               // 输出带一位保护位
        end
    end
 
    assign c0_o = c0;  assign c1_o = c1;   // 测试信号输出
    assign c2_o = c2;  assign c3_o = c3;
    assign count_o   = count;
    assign ena_in_o  = ena_in;
    assign ena_out_o = ena_out;
    assign d_out     = d;
endmodule

这里的控制部分与例 5.9 中讨论的 rc_sinc 设计非常相似。几个关键点：

FSM：完整循环 12 个时钟周期，负责控制流程、两个缓存器的使能信号生成，以及分数延迟 $D$ 的计算。增量取 DELTA = 85，即 $85/256 \approx 1/3$ ：对 $R = 3/4$ 的变换，分数延迟 $d$ 依次取 $0, 1/3, 2/3$ 然后回绕到 0，恰好覆盖每个输出采样所在的相位。
抽头延迟线：注意 4 个多项式系数 $c_{k}$ 只共用一条抽头延迟线。
SOP 块：先按式 (5-48) 计算拉格朗日矩阵输出 $c_{0} \dots c_{3}$ （系数 $- 85/256 \approx - 1/3$ 、 $43/256 \approx 1/6$ 正是逆矩阵中的分数系数按 256 定标后的整数近似），再按 Horner 模式 (5-50) 做三次”乘 $d$ 、加系数”的迭代得到输出。输出 y_out 带一位附加的保护位。

资源方面，这一设计使用了 363 个 LE 和 3 个嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 Fmax = 39.82 MHz。

图 5-43 给出了仿真结果：FSM 控制信号和使能信号在前，输入为三角形信号 x_in。滤波器输出每 4 个时钟周期更新一次，即整个循环中更新三次。滤波器输出值经过 Farrow 结构加权后生成输出 y_out。可以观察到：只有第一个和第二个拉格朗日多项式系数是非零的——这是因为三角形输入信号本身是分段线性函数，不含更高次的多项式成分；来自 $c_{0}$ 的滤波器值 20 和 60 在输出序列中保持不变（在这些时刻 $D = 0$ ），其他值则是插入的。

图 5-43 采用拉格朗日多项式和 Farrow 组合器的 R=3/4 速率变换的 VHDL 仿真结果

与 sinc 方案的可扩展性对比

就本例 $R = 3/4$ 而言，拉格朗日插值加 Farrow 组合器的实现开销与 sinc 滤波器设计差别不大。但真正的分野出现在 $R_{1}$ 较大的时候：Farrow 设计只需要在使能信号生成过程中做些改动，拉格朗日插值和 Farrow 组合器的工作量基本保持不变；而 sinc 滤波器的工作量与需要实现的滤波器数量（即 $R_{1}$ ）成正比，请参阅练习 5.16。Farrow 组合器的唯一缺点是乘法顺序串联导致延迟较长，但这可以通过增加乘法器和系数的流水线级来改进，请参阅练习 5.17。

5.6.3 基于 B 样条的分数速率变换器

为什么需要 B 样条

拉格朗日多项式逼近在中心非常光滑，但在多项式的末端纹波有增大的趋势（见图 5-41）。采用 B 样条逼近函数有望更为光滑。与拉格朗日多项式相反，B 样条长度有限，而且 $N$ 次 B 样条必须 $N$ 阶可微，这就从定义上保证了光滑性。

B 样条可以定义成多种形式，最常用的是通过框函数的连续积分来定义，如图 5-44 所示：0 次 B 样条（框函数）经积分得到三角形 B 样条，一次 B 样条经积分得到二次函数，依此类推。

图5-44 0至3次B样条函数

B 样条的解析描述使用以下斜坡函数的表示方法：

(t - τ)_{+} = {t - τ 0 \forall t > τ 其他 (5-51)

$N$ 次对称 B 样条的表达式为：

β^{N} (t) = \frac{1}{N !} k = 0 \sum N + 1 (- 1)^{k} (k N + 1) (t - k + \frac{N + 1}{2})_{+}^{N} (5-52)

由于 B 样条的所有线段都采用 $N$ 次多项式，因此是 $N$ 次可微的——即使在 B 样条的末端也非常光滑。0 次和 1 次 B 样条分别给出框函数和三角形，接下来是 2 次和 3 次 B 样条。尽管非常高质量的语音处理技术中使用 6 次 B 样条，但 DSP 中最常用的是 3 次 B 样条。对 3 次情形，根据式 (5-52) 得到：

β^{3} (t) = \frac{1}{6} \sum_{k = 0}^{4} (- 1)^{k} (k 4) (t - k + 2)_{+}^{3} = \frac{1}{6} (t + 2)_{+}^{3} - \frac{2}{3} (t + 1)_{+}^{3} + t_{+}^{3} - \frac{2}{3} (t - 1)_{+}^{3} + \frac{1}{6} (t - 2)_{+}^{3} = t > - 2 \frac{1}{6} (t + 2)^{3} - t > - 1 \frac{2}{3} (t + 1)^{3} + t > 0 t^{3} - t > 1 \frac{2}{3} (t - 1)^{3} + t > 2 \frac{1}{6} (t - 2)^{3} (5-53)

B 样条逼近与 B 样条插值

现在可以用 3 次 B 样条重构信号：把加权的 B 样条序列求和，即

\overset{y}{^} (t) = k \sum x (k) β^{3} (t - k) (5-54)

图 5-45 以粗实线给出了这个加权和。虽然 $\overset{y}{^} (t)$ 非常光滑，但可以观察到样条并不精确穿过采样点，即 $\overset{y}{^} (k) \neq = x (k)$ ——文献中把这种重构称为 B 样条逼近。如果希望加权和精确穿过每个采样点，那就是 B 样条插值。例如对 3 次 B 样条，可以对采样点先应用一个滤波器加权，其 Z 变换为：

H (z) = \frac{z + 4 + z ^{- 1}}{6} (5-55)

要获得完美插值，还需对输入采样应用逆 3 次 B 样条滤波器：

F (z) = 1/ H (z) = \frac{6}{z + 4 + z ^{- 1}} (5-56)

图 5-45 使用 3 次 B 样条的样条逼近

问题来了：检查 $F (z)$ 的极点/零点图会发现这个 IIR 滤波器不稳定——如果直接把它应用到输入序列，脉冲响应会产生不断增长的信号，请参阅练习 5.18。Unser 等人提出把滤波器拆分成稳定的因果部分和非因果部分，并让非因果滤波器从第一个因果滤波器输出的最后一个值开始运行。这种方式在图像处理中处理每条扫描线上的有限个采样时运行得很好，但在连续信号处理方案中就不再适用，特别是当滤波器以有限精度算法实现的时候。

图 5-46 3 次 B 样条插值的 FIR 补偿滤波器的设计

用 FIR 滤波器逼近补偿滤波器

有一种适用于连续信号处理的方法：用 FIR 滤波器逼近 $F (z) = 1/ H (z)$ 。事实证明几乎不用多少 FIR 系数就能得到很好的逼近，因为这个传递函数没有任何陡峭的边缘（见图 5-46）。设计流程也很直接：先计算 IIR 滤波器的传递函数，再用传递函数的 IFFT 确定 FIR 的时间域系数。如果应用中 DC 偏移很关键，还可以做偏差校正。Unser 等人提出过优化滤波器系数集合的算法，但由于有限系数精度和有限系数集合的特点，与直接 IFFT 方法相比增益不大（见图 5-46 和练习 5.19）。

现在可以先用 FIR 滤波器处理输入采样，然后用 3 次 B 样条重构。从图 5-47 可以看出，重构函数确实通过了每个原始采样点，即 $\overset{y}{^} (k) = x (k)$ ，实现了真正的插值。

图 5-47 采用 3 次 B 样条插值的 FIR 补偿滤波器的插值

三次 B 样条的 Farrow 方程

最后还差一步：为 B 样条插值设计分数延迟并导出 Farrow 滤波器结构。仍然使用常用的 3 次 B 样条，只考虑 $0 \leq d \leq 1$ 范围内的分数延迟，对 4 点插值使用输入信号在 $t = - 1, 0, 1, 2$ 时刻的采样值。有了 B 样条表达式 (5-53) 和加权和 (5-54)，需要考虑 4 个 B 样条线段，整理后得到：

y (d) = x (n + 2) \frac{d ^{3}}{6} + x (n + 1) [- \frac{d ^{3}}{2} + \frac{d ^{2}}{2} + \frac{d}{2} + \frac{1}{6}] + x (n) [\frac{d ^{3}}{2} - d^{2} + \frac{2}{3}] + x (n - 1) [- \frac{d ^{3}}{6} + \frac{d ^{2}}{2} - \frac{d}{2} + \frac{1}{6}] (5-57)

为了能在 Farrow 结构中实现，需要按 $d^{k}$ 因子整理上式，得到 4 个方程：

d^{0} : d^{1} : d^{2} : d^{3} : 000 x (n + 2) /6 + x (n + 1) /6 + x (n + 1) /2 + x (n + 1) /2 - x (n + 1) /2 + 2 x (n) /3 + 0 - x (n) x (n) /2 + x (n - 1) /6 = c_{0} - x (n - 1) /2 = c_{1} + x (n - 1) /2 = c_{2} - x (n - 1) /6 = c_{3} (5-60)

注意 $d^{0}$ 一行的权重之和为 $1/6 + 2/3 + 1/6 = 1$ ：当 $d = 0$ 时输出恰好等于”补偿滤波后采样”的加权平均，这保证了插值条件成立。此 Farrow 结构可以直接转换成 B 样条速率变换器，练习 5.23 将加以讨论。

5.6.4 MOMS 分数速率变换器

逼近的”阶数”：一个被忽视的设计参数

在插值核函数 $ϕ (t)$ 的传统设计中，经常被忽略的一个方面是逼近的阶数——当采样步长趋于零时，原始函数与重构函数之间方差（即 $L^{2}$ 范数）的衰减速率。而在实现工作量方面，差值函数的支集（support）或长度是关键设计参数。重要的一点是：上一节用到的 B 样条同时做到了”最大阶数与最小支集”（Maximum Order and Minimum Support, MOMS）。

那么问题是：B 样条是否是唯一具有 $L - 1$ 次、长度为 $L$ 的支集和 $L$ 阶的核函数？答案是否定的——可以证明存在一整类遵循 MOMS 行为的函数。这类插值多项式可以写成：

ϕ (t) = β^{N} (t) + k = 1 \sum N p (k) \frac{d ^{k} β ^{N} ( t )}{d t ^{k}} (5-61)

由于 B 样条可以通过与框函数的连续卷积构造，微分可以通过下式计算：

\frac{d ^{k + 1} β ( t )}{d t} = β^{k} (t + \frac{1}{2}) - β^{k} (t - \frac{1}{2}) (5-62)

式 (5-61) 提供了一组可选的设计参数 $p (k)$ ，用来满足不同的设计目标。很多设计要求插值核函数具有对称性，这就需要把所有奇次系数 $p (k)$ 强制为 0。常见选择是 $N = 3$ ，即 3 次样条类型，得到：

ϕ (t) = β^{3} (t) + p (2) \frac{d ^{2} β ^{3} ( t )}{d t ^{2}} = β^{3} (t) + p (2) (β^{1} (t + 1) - 2 β^{1} (t) + β^{1} (t - 1)) (5-63)

此时只需要确定一个设计参数 $p (2)$ 。不同的 $p (2)$ 取值对应不同风格的 MOMS 核函数：

I-MOMS（ $p (2) = - 1/6$ ）：可以设计成根本不需要补偿滤波器的直接插值函数，拉格朗日插值也属于这一类（请参阅练习 5.20），但给出的是次最优插值结果，见图 5-48(b)。
O-MOMS（ $p (2) = 1/42$ ）：以 $L^{2}$ 范数意义下插值误差最小化为目标，逼近误差比 I-MOMS 小一个量级，见图 5-48(c)。
还可以应用迭代方法将具体应用的插值 S/N 最大化。例如对一组具体的 5 幅图像， $p (2) = 1/28$ 的结果要比 O-MOMS 好 1 dB。

(a) B 样条

(b) I-MOMS

(c) O-MOMS

(d) C-MOMS
图5-48 长度为4的可能MOMS核函数

三种 MOMS 的补偿滤波器比较

O-MOMS 所需的补偿滤波器和 B 样条类似，存在一个不稳定的极点位置，在连续信号处理方案中必须使用 FIR 逼近。但要注意：如果 FIR 必须用有限精度算法构造，O-MOMS 靠小 $L^{2}$ 误差获得的增益，在系统总增益上未必能体现出来。

如果放弃核函数的对称性，则可以设计出在整数点采样的插值函数 $ϕ (k)$ 为因果函数（即 $ϕ (- 1) = 0$ ）的 MOMS 函数，这就是 C-MOMS，从图 5-48(d) 可以观察到这一点。C-MOMS 函数要求 $p (1) = - 1/3$ （而 $p (2) = p (3) = 0$ ，因此函数相当光滑），得到非对称的插值函数。它的主要优点是可以使用一个非常简单的单极点稳定 IIR 补偿滤波器：

F (z) = \frac{1.5}{1 + 0.5 z ^{- 1}}

和 B 样条、O-MOMS 一样不需要 FIR 逼近——但这里是一阶 IIR 而不是 FIR 阵列，实现代价小得多。

图 5-49 采用 IIR 补偿滤波器的 C-MOMS 插值

一个有趣的现象（对照图 5-47 和图 5-49）：C-MOMS 的最大值和采样点不再像对称核函数（如 B 样条）那样具有相同的时间-位置对应关系。但 C-MOMS 的加权和都精确通过采样点，正如对样条插值所期望的。试验显示：就插值而言，C-MOMS 表现得比 B 样条好一些，比全精度实现的 O-MOMS 差一些。

C-MOMS 的 Farrow 方程

剩下的事就是计算 Farrow 重采样器的方程。可以用式 (5-53) 直接计算插值函数 $ϕ (t) = β^{3} (t) - \frac{1}{3} d β^{3} (t) / d t$ ，再按延迟 $d^{k}$ 整理 Farrow 结构；另一个更省事的技巧是：利用对 B 样条预先算好的方程 (5-60)，对 Farrow 矩阵直接应用微分，即计算 $c_{k}^{new} = c_{k} - (k + 1) c_{k + 1} /3$ ， $k = 0, 1, 2$ ——这正对应 $p (1) = - 1/3$ 、 $p (2) = p (3) = 0$ 的三次 C-MOMS。同理也可为 O-MOMS 和 I-MOMS 计算 Farrow 方程（练习 5.20）。为 C-MOMS 生成如下 4 个方程：

d^{0} : d^{1} : d^{2} : d^{3} : 00 - x (n + 2) /6 x (n + 2) /6 + x (n + 1) /6 + x (n + 1) - x (n + 1) /2 + 2 x (n) /3 + 2 x (n) /3 - 3 x (n) /2 + x (n) /2 + 1 x (n - 1) /3 = c_{0} - 5 x (n - 1) /6 = c_{1} + 2 x (n - 1) /3 = c_{2} - x (n - 1) /6 = c_{3} (5-64)

例 5.14 R=0.75 的速率变换器 V

下面的设计采用三次 C-MOMS 样条多项式实现 $R = 3/4$ 的速率变换。原书 VHDL 改写为 Verilog：

// 原书 VHDL 改写为 Verilog：R = 3/4 三次 C-MOMS 样条速率变换器
module cmoms #(parameter IL = 3) (      // IL：输入缓存长度 - 1
    input  wire              clk,
    input  wire              reset,      // 异步复位
    input  wire signed [7:0] x_in,       // 系统输入
    output wire [3:0]        count_o,    // FSM 计数器
    output wire              ena_in_o,   // 采样输入使能
    output wire              ena_out_o,  // 输出移位使能
    output wire signed [8:0] xiir_o,     // IIR 滤波器输出
    output wire signed [8:0] c0_o, c1_o, c2_o, c3_o, // C-MOMS 矩阵
    output reg  signed [8:0] y_out);     // 系统输出
 
    // d^k 的预计算常数表（d 取 0、85、171，按 256 定标）
    reg signed [8:0] d1 [0:2];
    reg signed [8:0] d2 [0:2];
    reg signed [8:0] d3 [0:2];
    initial begin
        d1[0] = 0;  d1[1] = 85;  d1[2] = 171;   // d^1
        d2[0] = 0;  d2[1] = 28;  d2[2] = 114;   // d^2
        d3[0] = 0;  d3[1] = 9;   d3[2] = 76;    // d^3
    end
 
    reg  [3:0] count;                   // 循环 R1*R2 = 12 个周期
    reg  [1:0] t;                       // 三个延迟相位索引
    reg        ena_in, ena_out;
    reg  signed [8:0] xiir, x1;         // IIR 滤波器输出与延迟输入
    reg  signed [8:0] ibuf [0:IL];      // 抽头寄存器
    reg  signed [8:0] x    [0:IL];
    reg  signed [8:0] c0, c1, c2, c3;
    reg  signed [16:0] h0, h1, y0, y1, y2, y3, y;
    integer i;
 
    // FSM：产生使能信号；t 在三个延迟相位 {0, 1/3, 2/3} 间轮转
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            count <= 0;  t <= 1;
            ena_in <= 0; ena_out <= 0;
        end else begin
            count <= (count == 11) ? 0 : count + 1;
            ena_in  <= (count == 2 || count == 5 ||
                        count == 8 || count == 11);
            ena_out <= (count == 3 || count == 7 || count == 11);
            if (ena_out)
                t <= (t >= 2) ? 0 : t + 1;   // 计算相位延迟索引
        end
    end
 
    // IIR 补偿滤波器 F(z) = 1.5/(1 + 0.5*z^-1)
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            xiir <= 0;  x1 <= 0;
        end else if (ena_in) begin
            xiir <= 3*x1/2 - xiir/2;
            x1   <= x_in;
        end
    end
 
    // 一条共用的抽头延迟线（4 个多项式系数共用）
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= IL; i = i + 1) ibuf[i] <= 0;
        end else if (ena_in) begin
            for (i = 0; i < IL; i = i + 1) ibuf[i] <= ibuf[i+1];
            ibuf[IL] <= xiir;
        end
    end
 
    // 输出时刻一次性取出 4 个采样
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            for (i = 0; i <= IL; i = i + 1) x[i] <= 0;
        end else if (ena_out) begin
            for (i = 0; i <= IL; i = i + 1) x[i] <= ibuf[i];
        end
    end
 
    // C-MOMS 矩阵（式 5-64）+ 并行 LUT 组合器（不用 Farrow 结构）
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            c0 <= 0; c1 <= 0; c2 <= 0; c3 <= 0;
            h0 <= 0; h1 <= 0;
            y0 <= 0; y1 <= 0; y2 <= 0; y3 <= 0;  y <= 0;
        end else if (ena_out) begin
            c0 <= (85*x[0] + 171*x[1]) / 256;
            c1 <= (171*x[1] - 213*x[0] + 43*x[2]) / 256;
            c2 <= (171*x[0] - 43*x[3]) / 256 - 3*x[1]/2 + x[2];
            c3 <= 43*(x[3] - x[0]) / 256 + (x[1] - x[2]) / 2;
            y0 <= c0 * 256;         // d^0 = 1，乘 256 保持定标
            y1 <= c1 * d1[t];       // 并行乘法器/LUT
            y2 <= c2 * d2[t];
            y3 <= c3 * d3[t];
            h0 <= y0 + y1;          // 流水加法树
            h1 <= y2 + y3;
            y  <= h0 + h1;
        end
    end
 
    always @(posedge clk or posedge reset) begin
        if (reset) y_out <= 0;
        else       y_out <= y / 256;    // 输出带一位保护位
    end
 
    assign c0_o = c0;  assign c1_o = c1;  // 测试信号输出
    assign c2_o = c2;  assign c3_o = c3;
    assign count_o   = count;
    assign ena_in_o  = ena_in;
    assign ena_out_o = ena_out;
    assign xiir_o    = xiir;
endmodule

控制部分仍与例 5.9 的 rc_sinc 设计相似：FSM 完整循环 12 个时钟周期，生成输入、输出缓存器的使能信号。延迟索引 $t$ 在 0、1、2 之间轮转，对应 $d = 0, 85/256 \approx 1/3, 171/256 \approx 2/3$ 三个相位； $d^{1} = d$ 、 $d^{2}$ 、 $d^{3}$ 的值预先计算并作为常数表存储。IIR 块实现单极点补偿滤波器 $F (z) = 1.5/ (1 + 0.5 z^{- 1})$ ；TAP 与取窗块仍共用一条抽头延迟线。SOP 块按式 (5-64) 计算 C-MOMS 矩阵输出 $c_{k}$ ，再用 $d^{k}$ 加权求和——注意这里没有采用 Farrow 串联结构，而是用并行的乘法器/加法器树直接计算，这使设计的速度提高了一倍。输出 y_out 带一位附加的保护位。

资源方面，这一设计使用了 549 个 LE 和 3 个嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 Fmax = 95.27 MHz——比拉格朗日 Farrow 设计的 39.82 MHz 快得多，验证了并行 LUT 结构的速度优势。

图 5-50 给出了仿真结果：输入为图 5-51 所示的方波信号 x_in，IIR 滤波器输出显示出边缘的锐化（这正是补偿滤波器的作用），C-MOMS 矩阵输出值 $c_{k}$ 经 $d^{k}$ 加权并求和后生成输出 y_out。

图 5-50 采用三次 C-MOMS 样条和一个极点 IIR 补偿滤波器的 R=3/4 速率变换的 VHDL 仿真结果

与拉格朗日插值一样，也可以改用 Farrow 组合器计算输出 y_out。当阵列乘法器可用、且 $R_{1}$ 很大（并行常数表会变得非常庞大）时这样做才有意义，请参阅练习 5.21。

(a) 原始信号

(b) 经长度为 11 的 FIR 补偿滤波器滤波的原始信号

(c) O-MOMS 逼近（无补偿滤波器）

(d) 采用补偿滤波器的 O-MOMS 速率变换
图 5-51 基于 O-MOMS 的分数 R=3/2 的速率变换

速率变换方法的局限性：方波与吉布斯现象

最后说明一下速率变换方法的局限。一个特别困难的问题是对方波输入信号做速率变换：从吉布斯现象（参阅图 3-6）已经知道，任何有限滤波器都有在方波边沿引入振铃的趋势。一个常见的实际任务来自 DAT 记录仪：它使用 32 kHz 和 48 kHz 两种采样频率，两者之间的转换因子为 $R = 48/32 = 3/2$ 。图 5-51 给出了使用 O-MOMS 样条插值对方波的速率变换结果：图 5-51(b) 中的 FIR 前置滤波器增强了边缘；图 5-51(c) 给出了没有 FIR 前置滤波器的 O-MOMS 三次样条速率变换器的结果。仔细对比可以发现：带前置滤波器时，O-MOMS 三次样条插值的边缘保留得更好（图 5-51(d)）。但同时也能看到：即使采用 O-MOMS 和长度为 11 的全精度 FIR 补偿滤波器，吉布斯现象依然可见——这是有限长滤波器的固有属性，不是实现方法的问题。

5.7 滤波器组

数字滤波器组是一组具有公共输入或公共输出的滤波器，如图 5-52 所示。两个方向各有用途：

图 5-52(a) 的分析滤波器组：把一个输入信号分成 $R$ 个不同的所谓子带信号，经常用于频谱分析。
图 5-52(b) 的合成滤波器组：把多个信号合成到一个公共的输出信号中。

图 5-52 典型滤波器组的分解系统展示

图5-53 有少量重叠的 $R$ 通道滤波器组

分析滤波器之间的频带关系可以是非重叠、稍微重叠或基本重叠的。图 5-53 给出了一个最常见的稍微重叠的滤波器组示例。

区分不同滤波器组的另一个重要特征是带宽和各个滤波器中心频率之间的间隔。非均匀滤波器组的常见示例是倍频间隔（小波）滤波器组，将在 5.8 节讨论。而均匀滤波器组中，所有滤波器都具有同样的带宽和采样速率。从实现角度讲，人们更倾向于均匀的、最大抽取滤波器组——因为它们可以借助 FFT 算法高效实现，这是下一节的主题。

5.7.1 均匀 DFT 滤波器组

从原型滤波器到 DFT 滤波器组

在最大抽取或精密采样滤波器组中，抽取或插入次数 $R$ 与频带数量 $K$ 相等。如果第 $r$ 个频带滤波器 $h^{r} [n]$ 从单个原型滤波器 $h [n]$ 的”取模运算”计算得来：

h^{r} [n] = h [n] W_{R}^{r n} = h [n] e^{- j 2 π r n / R} (5-65)

就称之为 DFT 滤波器组。直观理解：所有带通滤波器都是同一个原型低通滤波器乘上不同的复指数”旋转”到不同中心频率得到的。

多相分解带来高效实现

如果采用滤波器 $h^{r} [n]$ 和输入信号 $x [n]$ 的多相分解（查阅 5.2 节），就能得到 $R$ 通道滤波器组的高效实现。因为每个带通滤波器都是精密采样的，采用 $R$ 个多相信号分解：

h [n] = k = 0 \sum R - 1 h_{k} [n] \leftrightarrow h_{k} [m] = h [m R - k] (5-66)

x [n] = k = 0 \sum R - 1 x_{k} [n] \leftrightarrow x_{k} [m] = x [m R - k] (5-67)

把式 (5-66) 代入式 (5-65) 后会发现一个关键结论：所有带通滤波器 $h^{r} [n]$ 共享同一组多相滤波器 $h_{k} [n]$ ，不同滤波器只差一个”旋转因子”。图 5-54(a) 给出了第 $r$ 个滤波器 $h^{r} [n]$ 的结构， $h^{r} [n]$ 的”旋转乘法”通过输入向量 $\overset{x}{^}_{0} [n], \overset{x}{^}_{1} [n], \dots, \overset{x}{^}_{R - 1} [n]$ 与第 $r$ 个 DFT 分量相关联。于是，整个分析频带的计算可以简化为：先经过 $R$ 个多相滤波器滤波（图 5-54(b)），再做这 $R$ 个滤波分量的 DFT（或 FFT）。这比按式 (5-65) 直接计算每个滤波器（参阅练习 5.6）高效得多。

图5-54 DFT滤波组的分析

均匀 DFT 合成滤波器组的开发过程正好与分析滤波器组相反：用 $R$ 个频谱分量 $\hat{X}^{r} [k]$ 作为逆 DFT（或 FFT）的输入，再用多相插入器结构重构输出信号，如图 5-55 所示。重构的带通滤波器变成：

f^{r} [n] = \frac{1}{R} f [n] W_{R}^{- r n} = f [n] e^{j 2 π r n / R} (5-68)

图 5-55 DFT 合成滤波器组

完美重构条件

把分析滤波器组与合成滤波器组组合起来，DFT 与 IDFT 会相互抵消。如果包含多相滤波器的卷积给出一个单位采样函数，就能得到完美重构：

h_{r} [n] * f_{r} [n] = {10 n = d 其他 (5-69)

换句话说，两个多相函数必须是互逆滤波器：

H_{r} (z) \times F_{r} (z) = z^{- d} F_{r} (z) = \frac{z ^{- d}}{H _{r} ( z )}

其中延迟 $d$ 的作用是保证合成滤波器是因果（可实现）的。在实际设计中，这些理想条件不可能用两个 FIR 滤波器完全满足，只能近似——可以采用两个 FIR 滤波器，或将一个 FIR 与 IIR 滤波器组合起来，请看下面的示例。

例5.15 DFT 滤波器组

例 4.3 讨论过的有耗积分器在当前上下文中可以解释成一个 $R = 2$ 的 DFT 滤波器组。其差分方程为：

y [n + 1] = \frac{3}{4} y [n] + x [n] (5-70)

在 $z$ 域内的传递函数是：

F (z) = \frac{z ^{- 1}}{1 - 0.75 z ^{- 1}} (5-71)

要得到两个多相滤波器，可以采用与例 4.5”分散预先考虑”相近的方案：在镜像位置引入一个额外的极点/零点对，即用 $(1 + 0.75 z^{- 1})$ 同时乘以分母和分子，得到：

F (z) = H_{0} (z^{2}) \frac{0.75 z ^{- 2}}{1 - 0.7 5 ^{2} z ^{- 2}} + z^{- 1} H_{1} (z^{2}) \frac{1}{1 - 0.7 5 ^{2} z ^{- 2}} = H_{0} (z^{2}) + z^{- 1} H_{1} (z^{2}) (5-72)

这一步的技巧在于：分母变成了只含 $z^{- 2}$ 的形式，各项天然按偶/奇下标分离，多相结构就自动浮现了。后者给出两个多相滤波器：

H_{0} (z) = \frac{0.75 z ^{- 1}}{1 - 0.7 5 ^{2} z ^{- 1}} = 0.75 z^{- 1} + 0.4219 z^{- 2} + 0.2373 z^{- 3} + \dots (5-74)

H_{1} (z) = \frac{1}{1 - 0.7 5 ^{2} z ^{- 1}} = 1 + 0.5625 z^{- 1} + 0.3164 z^{- 2} + \dots (5-75)

展开验证一下系数： $H_{0}$ 的幂级数为 $0.75 z^{- 1} (1 + 0.5625 z^{- 1} + 0.3164 z^{- 2} + \dots)$ ，即 $0.75, 0.4219, 0.2373$ ——每项都是前一项乘 0.5625（即 $0.7 5^{2}$ ），与 $1/ (1 - 0.5625 z^{- 1})$ 的几何级数一致，说明分解正确。虽然可以用非递归 FIR 近似这些脉冲响应，但要达到误差小于 1% 就必须用 16 个系数；因此采用式 (5-74)、(5-75) 定义的递归多相 IIR 滤波器更加高效。

多相分解之后，采用一个由下面矩阵给出的 2 点 DFT：

W = [11 1 - 1]

整个分析滤波器组就被构造成图 5-56(a) 的形式。

(a) 分析滤波器组

图 5-56 对于 R=2 的严格采样的均匀 DFT 滤波器组

对于合成滤波器组，首先用下面的矩阵计算逆 DFT：

W^{- 1} = \frac{1}{2} [11 1 - 1]

为了获得完美重构，必须找到 $h_{0} [n]$ 和 $h_{1} [n]$ 的逆滤波器。这并不困难： $H_{r} (z)$ 是单极点 IIR 滤波器，所以它的逆 $F_{r} (z) = z^{- d} / H_{r} (z)$ 必然是两抽头 FIR 滤波器。利用式 (5-74) 和 (5-75) 可得 $d = 1$ （这对于得到因果滤波器已经足够），于是：

F_{0} [n] = \frac{4}{3} (1 - 0.7 5^{2} z^{- 1}) (5-76)

F_{1} [n] = z^{- 1} - 0.7 5^{2} z^{- 2} (5-77)

可以验算完美重构条件 (5-69)： $H_{0} (z) \cdot F_{0} (z) = \frac{0.75 z ^{- 1}}{1 - 0.5625 z ^{- 1}} \cdot \frac{4}{3} (1 - 0.5625 z^{- 1}) = z^{- 1}$ ，延迟 $d = 1$ 的单位冲激，正是所求；对 $H_{1} (z) \cdot F_{1} (z)$ 同样成立。图 5-56(b) 给出了合成滤波器组的图形解释。这个例子的启发在于：分析端用递归 IIR 实现多相滤波器，合成端用极短 FIR 做逆滤波器，两者配合就能以极低的代价逼近完美重构。

5.7.2 双通道滤波器组

什么是双通道滤波器组，为什么需要它

设想我们要把一路信号”拆成两半”来看：一半是低频（轮廓、概貌），一半是高频（细节、边缘）。双通道滤波器组就是完成这件事的标准结构：输入信号先经过低通分析滤波器 $G (z)$ 和高通分析滤波器 $H (z)$ 分成两个子带，各自做 2 倍抽取（数据量减半，总数据量不膨胀）；传输或存储之后，再在合成端先做 2 倍插值（补零），用合成滤波器 $\hat{G} (z)$ 、 $\hat{H} (z)$ 把两个子带重新拼回一路输出 $\overset{x}{^} (n)$ 。分析端与合成端之间的信号通常还会被量化或做非线性处理（比如压缩编码），所以滤波器组必须”经得起折腾”——即使中间有失真，合成端也要尽量还原原信号。

高通滤波器一般不必单独设计，习惯上由低通滤波器直接派生：

h [n] = (- 1)^{n} g [n] \circ - ∙ H (z) = G (- z) (5-78)

这样两个滤波器的幅频特性关于 $π /2$ 镜像对称，即 $∣ H (e^{j ω}) ∣ = ∣ G (e^{j (ω - π)}) ∣$ ，因此这类结构称为正交镜像滤波器（Quadrature Mirror Filter, QMF）组。

核心问题是：能做到完美重构吗？ 即输出满足

\overset{x}{^} [n] = x [n - d] (5-79)

也就是输出与输入形状完全相同，只差一个时移 $d$ 。难点在于 $G (z)$ 和 $H (z)$ 都不是理想的”方波”滤波器：经过 2 倍抽取后，每个通道都会产生混叠分量。历史上最早给出简单可行答案的是 Haar（约 1910 年）提出的滤波器组。

图 5-57 采用长度为 4 的 Daubechies 滤波器的双通道滤波器组

例5.16 双通道Haar滤波器组 I

取最简单的一组滤波器：

G (z) = 1 + z^{- 1}, H (z) = 1 - z^{- 1}

\hat{G} (z) = \frac{1}{2} (1 + z^{- 1}), \hat{H} (z) = \frac{1}{2} (- 1 + z^{- 1})

直观理解：低通通道算的是相邻两点的和 $x [n] + x [n - 1]$ ，高通通道算的是差 $x [n] - x [n - 1]$ 。和与差携带的信息合起来不多不少正好是原始数据——这就像初等代数里” $a + b$ 与 $a - b$ 可以唯一还原 $a$ 和 $b$ “。经过抽取、插值与合成滤波器之后，输出恰好是延迟一个采样点的原序列，即 $\overset{x}{^} (n) = x [n - 1]$ ，实现完美重构，延迟 $d = 1$ 。

图 5-58 双通道 Haar-QMF 滤波器组

完美重构的一般条件

要推广到任意滤波器，先注意一个关键事实：对信号 $s [k]$ 先抽取再插值 2，等价于让 $S (z)$ 乘上序列 ${1, 0, 1, 0, \dots}$ ，在 $z$ 域中写成：

S_{↓↑} (z) = \frac{1}{2} (S (z) + S (- z)) (5-80)

把它套进双通道结构，低通与高通通路的输出分别是：

X_{↓↑ G} (z) = \frac{1}{2} (X (z) G (z) + X (- z) G (- z)) (5-81)

X_{↓↑ H} (z) = \frac{1}{2} (X (z) H (z) + X (- z) H (- z)) (5-82)

再乘上各自的合成滤波器并相加，得到总输出：

\hat{X} (z) = \frac{1}{2} (G (z) \hat{G} (z) + H (z) \hat{H} (z)) X (z) + \frac{1}{2} (G (- z) \hat{G} (z) + H (- z) \hat{H} (z)) X (- z) (5-83)

注意式(5-83)里有两类”坏东西”： $X (- z)$ 前的因子是混叠分量（抽取把高频折叠进来产生的假信号）， $X (z)$ 前的因子若不是纯延迟则带来幅值畸变。想让输出干净，二者都必须被消掉——这就是下面的定理。

定理 5.17 完美的重构

双通道滤波器组实现完美重构，当且仅当：

$G (- z) \hat{G} (z) + H (- z) \hat{H} (z) = 0$ ，即彻底消除混叠；
$G (z) \hat{G} (z) + H (z) \hat{H} (z) = 2 z^{- d}$ ，即幅值畸变只剩一个纯延迟。

用 Haar 组验证一遍。

例 5.18 双通道 Haar 滤波器组 II

仍取例 5.16 的四个滤波器：

条件 1（消混叠）：

G (- z) \hat{G} (z) + H (- z) \hat{H} (z) = \frac{1}{2} (1 - z^{- 1}) (1 + z^{- 1}) + \frac{1}{2} (1 + z^{- 1}) (- 1 + z^{- 1}) = \frac{1}{2} (1 - z^{- 2}) + \frac{1}{2} (- 1 + z^{- 2}) = 0

两项恰好相消——差通道的混叠抵消了和通道的混叠。

条件 2（幅值）：

G (z) \hat{G} (z) + H (z) \hat{H} (z) = \frac{1}{2} (1 + z^{- 1})^{2} + \frac{1}{2} (1 - z^{- 1}) (- 1 + z^{- 1}) = \frac{1}{2} [(1 + 2 z^{- 1} + z^{- 2}) + (- 1 + 2 z^{- 1} - z^{- 2})] = 2 z^{- 1}

结果正是 $2 z^{- d}$ ， $d = 1$ ，与例 5.16 的数值实验一致。另外一个有用的观察：把分析滤波器和合成滤波器整体”对调”，两个条件依然成立，完美重构不受影响。

简化设计：与混叠无关的结构

直接同时凑两个条件不容易，第一步先”承包”掉混叠条件。

定理 5.19 与混叠无关的双通道滤波器组

若选择四个滤波器满足

G (- z) = - \hat{H} (z) 且 H (- z) = \hat{G} (z) (5-84)

则该滤波器组与混叠无关（条件 1 自动成立）。对长度为 4 的滤波器，这两条意味着系数间简单的符号翻转关系：

g [n] = {g [0], g [1], g [2], g [3]} \to \hat{h} [n] = {- g [0], g [1], - g [2], g [3]}

h [n] = {h [0], h [1], h [2], h [3]} \to \overset{g}{^} [n] = {h [0], - h [1], h [2], - h [3]}

有了这个约束，第二个条件也随之简化。定义乘积滤波器 $F (z) = G (z) \hat{G} (z)$ ，利用式(5-84)可得：

G (z) \hat{G} (z) + H (z) \hat{H} (z) = F (z) - F (- z) (5-85)

于是完美重构只剩一个要求：

F (z) - F (- z) = 2 z^{- d} (5-86)

也就是说，乘积滤波器 $F (z)$ 必须是半带滤波器——除中心抽头外，所有偶数序号系数都为零。这样设计流程就变成了三步。

算法 5.20 完美重构的双通道滤波器组

(1) 按式(5-86)设计一个标准因果半带滤波器 $F (z)$ 。 (2) 把 $F (z)$ 因式分解为 $F (z) = G (z) \hat{G} (z)$ 。 (3) 用式(5-84)补出另外两个滤波器： $\hat{H} (z) = - G (- z)$ ， $H (z) = \hat{G} (- z)$ 。

例 5.21 采用 F3 的完美重构滤波器组

长度为 7 的（标准）因果半带滤波器 F3 的传递函数为：

F3 (z) = \frac{1}{16} (- 1 + 9 z^{- 2} + 16 z^{- 3} + 9 z^{- 4} - z^{- 6}) (5-87)

先验证式(5-86)：代入 $- z$ 后所有偶次幂项不变、奇次幂变号，而 F3 除 $16 z^{- 3}$ 外全是偶次幂，故 $F 3 (z) - F 3 (- z) = 2 \cdot \frac{16}{16} z^{- 3} = 2 z^{- 3}$ ，成立。其六个零点为 $z_{01} = z_{02} = z_{03} = z_{04} = - 1$ 、 $z_{05} = 2 + 3 = 3.7321$ 、 $z_{06} = 2 - 3 = 0.2679 = 1/ z_{05}$ 。把这六个零点分配给 $G (z)$ 与 $\hat{G} (z)$ 的方式不唯一，例如：

5/3 滤波器：

G (z) = \frac{1}{8} (- 1 + 2 z^{- 1} + 6 z^{- 2} + 2 z^{- 3} - z^{- 4}), \hat{G} (z) = \frac{1}{2} (1 + 2 z^{- 1} + z^{- 2})

4/4 滤波器（线性相位）：

G (z) = \frac{1}{4} (1 + z^{- 1})^{3}, \hat{G} (z) = \frac{1}{4} (- 1 + 3 z^{- 1} + 3 z^{- 2} - z^{- 3})

4/4 Daubechies 滤波器（正交，常用于小波分析）：

G (z) = \frac{1 - 3}{4 2} (1 + z^{- 1})^{2} (- z_{05} + z^{- 1}), \hat{G} (z) = - \frac{1 + 3}{4 2} (1 + z^{- 1})^{2} (- z_{06} + z^{- 1})

对于 Daubechies 结构还额外满足 $H (z) = - z^{- N} G (- z^{- 1})$ ，即高、低通多项式互为时间反转镜像——这是正交滤波器组的典型特征。

图 5-59 半带滤波器 F3 各种分解的极点/零点图（一）

图 5-59 半带滤波器 F3 各种分解的极点/零点图（二）

图 5-59 半带滤波器 F3 各种分解的极点/零点图（三）

(a) 线性相位 5/3 滤波器

(b) 线性相位 4/4 滤波器

(c) 4/4 Daubechies 滤波器

零点怎么分不是随意的，观察图 5-59 可以总结出因式分解的规则。

推论 5.22 半带滤波器的因式分解

构造实系数滤波器时，共轭零点对 $(z_{0}, z_{0}^{*})$ 必须分到同一个滤波器中，否则系数出现复数。
想要线性相位，零点图必须关于单位圆上的点 $z = 1$ 对称，因此互为倒数的零点对 $(z_{0}, 1/ z_{0})$ 必须在同一滤波器中。
想要正交滤波器组（ $F (z) = U (z) U (z^{- 1})$ ，两支互为镜像），所有 $z_{0}$ 与 $1/ z_{0}$ 对又必须分到不同的滤波器中。

条件 2) 与 3) 是对立的，所以”既正交又线性相位”一般做不到——除非所有零点都落在单位圆上，Haar 滤波器组正是这种特例。据此给例 5.21 的三种结构分类：(a)、(b) 是实线性相位滤波器，(c) 是实正交滤波器。

5.7.3 实现双通道滤波器组

设计定了，接下来是工程问题：怎么在 FPGA 上高效实现。下面由一般到特殊，依次讨论多相实现、提升、QMF、正交与线性相位几种方案（只讲分析端，合成端对整个结构做图形转置即可得到）。

1. 多相双通道滤波器组

最一般的做法：把 $G (z)$ 和 $H (z)$ 各自拆成偶、奇两个多相分量，共四个长度为 $L /2$ 的子滤波器，先抽取后滤波：

图 5-60 双通道滤波器组的多相实现

多相分解本身不省运算量（仍是 $2 L$ 个乘法器、 $2 (L - 1)$ 个加法器），但把滤波器搬到了抽取器之后，四个子滤波器只需以输入采样率两倍的 $2 f_{s}$ 运行，且每个子滤波器长度减半。这些短滤波器还可以套用各种”加速器”：

行程长度（run-length）滤波器：用短 Winograd 卷积算法；
FFT 快速卷积（第 6 章）；
第 3 章的高级技术：分布式算法、简化加法器图（RAG）、余数系统等。

用 FFT 还有个额外好处：每个多相分量的正变换只需算一次，逆变换可以对两个分量频谱之和统一进行（见图 5-61）。不过 FFT 只有在滤波器较长（一般大于 32）时才划算，而典型双通道滤波器长度小于 32。

2. 提升

Sweldens、Herley 和 Vetterli 提出的提升（lifting）是另一种构造思路：像格型滤波器那样，用”提升”与”对偶提升”两种交叉项，从短滤波器逐步搭建出长滤波器，并始终保持完美重构。基本结构如图 5-62。

图 5-62 采用提升和双重提升步骤的双通道滤波器的实现

设计从最简单的”惰性滤波器组”出发： $G (z) = \hat{H} (z) = 1$ ， $H (z) = \hat{G} (z) = z^{- 1}$ ——它直接把偶、奇样本分开，可验证满足定理 5.17。在此基础上做如下变换而不破坏完美重构：

提升：对任意 $S (z^{2})$ ， $G^{'} (z) = G (z) + \hat{G} (- z) S (z^{2})$ 　(5-89)

对偶提升：对任意 $T (z^{2})$ ， $\hat{G}^{'} (z) = \hat{G} (z) + G (- z) T (z^{2})$ 　(5-90)

把提升公式代回定理 5.17 可以验证：只要 $\hat{G} (z)$ 、 $\hat{H} (z)$ 仍满足定理 5.19 的与混叠无关条件，完美重构就保持不变。也就是说，每一级提升都是一次”免费升级”。

例 5.23 DB4 滤波器的提升实现

例 5.21 中的 Daubechies 长度 4 滤波器（系数见例 5.24）可以用两个提升步骤加一个对偶提升步骤实现，得到的差分方程为：

h_{1} [n] = x [2 n + 1] - 3 x [2 n]

g_{1} [n] = x [2 n] + \frac{3}{4} h_{1} [n] + \frac{3 - 2}{4} h_{1} [n - 1]

h_{2} [n] = h_{1} [n] + g_{1} [n + 1]

g [n] = \frac{3 + 1}{2} g_{1} [n], h [n] = \frac{3 - 1}{2} h_{2} [n]

注意输入先被拆成偶序列 $x [2 n]$ 和奇序列 $x [2 n - 1]$ ，滤波器以 $2 f_{s}$ 速率运行；整个结构只用 5 个乘法器和 4 个加法器，可直接映射为硬件。合成端由图形转置得到：运算次序反过来、系数符号翻转即可。

Daubechies 和 Sweldens 证明了一个更强的结论：任何（双）正交小波滤波器组都能化为一系列提升/对偶提升步骤。步骤越多，每步的系数越简单，与直接多相实现相比可省约 50% 的乘法器与加法器；在小位宽乘法器场合尤其有吸引力。但提升的格状结构没法用 RAG 技术，所以对较长的滤波器，直接多相法往往仍更高效。

3. QMF 实现

若滤波器满足 QMF 关系式(5-91) $h [n] = (- 1)^{n} g [n]$ ，则两个多相分量满足

G_{0} (z) = H_{0} (z), G_{1} (z) = - H_{1} (z) (5-92)

即低通与高通共享同一个偶分量，奇分量只差一个符号。于是图 5-60 中的四个滤波器可以压缩为两个滤波器加一级”蝶形”（交换与相减），成本直接省掉约 50%：

L 个实加法器、 L 个实乘法器 (5-93)

且滤波器仍以 $2 f_{s}$ 速率运行。

4. 正交滤波器组

正交滤波器对满足共轭镜像滤波器（CQF）条件：

H (z) = z^{- N} G (- z^{- 1}) (5-94)

此时若采用转置 FIR 结构（图 5-64），系数共享可使乘法器减半；代价是这种结构无法再做多相分解，速度上不去，只能以 $f_{s}$ 运行。

图 5-64 采用转置 FIR 结构的正交双通道滤波器组

另一条路是图 5-65 的网格（lattice）滤波器实现。

例 5.24 L=4 的网格 Daubechies 滤波器的实现

长度为 4 的 Daubechies 滤波器系数为：

G (z) = \frac{( 1 + 3 ) + ( 3 + 3 ) z ^{- 1} + ( 3 - 3 ) z ^{- 2} + ( 1 - 3 ) z ^{- 3}}{4 2} = 0.48301 + 0.8365 z^{- 1} + 0.2241 z^{- 2} - 0.1294 z^{- 3} (5-95)

H (z) = \frac{- ( 1 - 3 ) + ( 3 - 3 ) z ^{- 1} - ( 3 + 3 ) z ^{- 2} + ( 1 + 3 ) z ^{- 3}}{4 2} = 0.1294 + 0.2241 z^{- 1} - 0.8365 z^{- 2} + 0.48301 z^{- 3} (5-96)

二阶双通道网格滤波器的传递函数是：

G (z) = (1 + a [0] z^{- 1} - a [0] a [1] z^{- 2} + a [1] z^{- 3}) s (5-97)

H (z) = (- a [1] - a [0] a [1] z^{- 1} - a [0] z^{- 2} + z^{- 3}) s (5-98)

与式(5-95)逐项比较系数，解出三个参数：

s = \frac{1 + 3}{4 2}, a [0] = \frac{3 + 3}{4 2 s}, a [1] = \frac{1 - 3}{4 2 s} (5-99)

原书 VHDL 改写为 Verilog：下面的 Verilog 代码实现了这个网格滤波器组。设计要点有三处巧妙的定点化：输入乘 $s = 0.48 \approx 124/256$ （写成 $4 \times (32 x - x) = 124 x$ ）；交叉项 $a [0] = 3 \approx 1.732$ 用移位近似为 $2 - 2^{- 2} - 2^{- 6} - 2^{- 8}$ ； $a [1] = 0.2679 \approx 2^{- 2} + 2^{- 6} + 2^{- 8}$ ，全部乘法都化为移位加减，不用任何硬件乘法器。输入流按奇偶分成两路，下半树在第一级需要延迟一个采样周期对齐。

// DB4 格型双通道滤波器组（原书 VHDL db4latti 改写为 Verilog）
module db4latti (
    input  wire        clk,      // 系统时钟
    input  wire        reset,    // 异步复位
    output wire        clk2,     // 二分频时钟输出
    input  wire signed [7:0]  x_in,  // 系统输入
    output wire signed [16:0] x_e,   // 偶样本测试信号
    output wire signed [16:0] x_o,   // 奇样本测试信号
    output reg  signed [8:0]  g,     // 低通输出
    output reg  signed [8:0]  h      // 高通输出
);
 
    reg        state;            // 0: even, 1: odd
    reg        clk_div2;
    reg signed [16:0] sx_up, sx_low, x_wait;
    reg signed [16:0] up0;
    reg signed [16:0] low0;
    wire signed [16:0] sxa0_up, sxa0_low;
    wire signed [16:0] up1, low1;
 
    // ---------- 状态机：把输入流按奇偶拆成两路 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            state    <= 1'b0;
            sx_up    <= 0;
            sx_low   <= 0;
            x_wait   <= 0;
            clk_div2 <= 1'b0;
        end else begin
            if (!state) begin                       // 偶相位
                sx_up  <= 4 * (32 * x_in - x_in);   // 乘 256*s = 124
                sx_low <= 4 * (32 * x_wait - x_wait);
                clk_div2 <= 1'b1;
                state  <= 1'b1;
            end else begin                          // 奇相位
                x_wait   <= x_in;
                clk_div2 <= 1'b0;
                state    <= 1'b0;
            end
        end
    end
 
    // ---------- 交叉项乘法 a[0] = 1.732 ≈ 2 - 2^-2 - 2^-6 - 2^-8 ----------
    assign sxa0_up  = (2*sx_up  - (sx_up  >>> 2)) - ((sx_up  >>> 6) + (sx_up  >>> 8));
    assign sxa0_low = (2*sx_low - (sx_low >>> 2)) - ((sx_low >>> 6) + (sx_low >>> 8));
 
    // ---------- 第一级：下半树需要打一拍对齐 ----------
    assign up0 = sxa0_low + sx_up;
    always @(posedge clk or posedge reset) begin
        if (reset)
            low0 <= 0;
        else if (clk_div2)
            low0 <= sx_low - sxa0_up;
    end
 
    // ---------- 第二级交叉项 a[1] = 0.2679 ≈ 2^-2 + 2^-6 + 2^-8 ----------
    assign up1  = (up0 - (low0 >>> 2)) - ((low0 >>> 6) + (low0 >>> 8));
    assign low1 = (low0 + (up0  >>> 2)) + ((up0  >>> 6) + (up0  >>> 8));
 
    // ---------- 输出定标：除以 256 ----------
    assign x_e  = sx_up;
    assign x_o  = sx_low;
    assign clk2 = clk_div2;
 
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            g <= 0;
            h <= 0;
        end else if (clk_div2) begin
            g <= up1[16:8];   // up1 / 256
            h <= low1[16:8];  // low1 / 256
        end
    end
 
endmodule

这段代码是图 5-65 网格结构的直接转换。仿真结果（图 5-66）显示了幅值为 100 的脉冲分别位于偶数位置和奇数位置时，滤波器 $G (z)$ 和 $H (z)$ 的响应。该设计使用 420 个 LE，未占用嵌入式乘法器，TimeQuest 慢速 85C 模型下 $F_{ma x} = 58.11$ MHz。

图 5-66 长度为 4 的 Daubechies 网格滤波器组的仿真结果

与直接多相实现相比，网格实现规模大致相同。虽然网格只需 5 个乘法器（多相实现要 8 个），但多相实现中的系数可以用 RAG 技术共享实现，而网格中每个乘法系数都得单独搭，单个系数的效率反而低——这提醒我们：乘法器个数少不等于硬件省。

5. 线性相位双通道滤波器组

第 3 章讲过：偶对称或奇对称的线性滤波器可以省 50% 乘法器。滤波器长度为偶数时，这一对称性在多相分解后仍然保持，而且速度同样可以提高一倍。如果 $G (z)$ 和 $H (z)$ 等长，还可以用图 5-67 的线性相位网格进一步压缩工作量（注意这种网格与正交滤波器组的图 5-65 网格结构不同）。

例 5.25 L=4 的线性相位滤波器的网格

取例 5.21 中的线性相位滤波器对，长度均为 4：

G (z) = \frac{1}{4} (1 + 3 z^{- 1} + 3 z^{- 2} + z^{- 3}) (5-100)

H (z) = \frac{1}{4} (- 1 + 3 z^{- 1} + 3 z^{- 2} - z^{- 3}) (5-101)

长度 $L = 4$ 的线性相位网格滤波器传递函数为：

G (z) = ((1 + a [0]) + a [0] z^{- 1} + a [0] z^{- 2} + (1 + a [0]) z^{- 3}) s (5-102)

H (z) = (- (1 + a [0]) + a [0] z^{- 1} + a [0] z^{- 2} - (1 + a [0]) z^{- 3}) s (5-103)

比较式(5-100)与式(5-102)的系数：

s = \frac{1}{2}, a [0] = - 1.5 (5-104)

整个结构只剩一个系数 $a [0]$ ，与直接实现相比只用约四分之一的乘法器。但线性相位网格有适用限制：并非所有线性相位滤波器都能实现——要求 $G (z)$ 偶对称、 $H (z)$ 奇对称、两滤波器等长且长度为偶数。

6. 实现选项的比较

表 5-6 汇总了各种方案。表中给出的是加法器/乘法器数量、参考图、最大输入速率，以及系数能否用 RAG 实现（或只能作为单乘法器系数）。

类型	实乘法器	实加法器	参考图	速度	可用 RAG
任意系数多相实现
直接 FIR 滤波器	$2 L$	$2 L - 2$	图 5-60	$2 f_{s}$	√
提升	$\approx L$	$\approx L$	图 5-62	$2 f_{s}$	—
正交镜像滤波器（QMF）
恒等多相滤波器	$L$	$L$	图 5-63	$2 f_{s}$	√
正交滤波器
转置 FIR 滤波器	$L$	$2 L - 2$	图 5-64	$f_{s}$	√
网格	$L + 1$	$3 L /4$	图 5-65	$2 f_{s}$	—
线性相位滤波器
对称滤波器	$L$	$2 L - 2$	图 3.5	$2 f_{s}$	√
网格	$L /2$	$3 L /2 - 1$	图 5-67	$2 f_{s}$	—

经验规律：短滤波器适合网格结构；长滤波器多数情况下 RAG 能生成更小更快的设计。另外注意表中数量是对硬件工作量的估计，而不是文献中 PDSP/μP 方案里”每采样计算量”的指标，两者不要混淆。双通道滤波器组的进一步阅读可参考文献 [123、159、162、163]。

5.8 小波

从短时傅立叶变换说起

语音、音频、图像这类信号有个特点：在短时间内统计特性基本不变。所以自然的分析方法是：加一个短窗，算这一帧的参数，再把窗往前滑，逐帧分析。如果每帧内做的是傅立叶变换，就叫短时傅立叶变换（STFT）：

X (τ, f) = \int_{- \infty}^{\infty} x (t) w (t - τ) e^{- j 2 π f t} d t (5-105)

窗函数 $w (t - τ)$ 在时间和频率上都平滑地过渡到 0，从而把映射定位在 $Δ f$ 与 $Δ t$ 之内。从 Heisenberg 不确定准则看，高斯窗 $g (t) = e^{- t^{2}}$ 是最优的——它给出最小的 $Δ f Δ t$ 乘积，Gabor 在 1949 年指出了这一点；其离散化即离散 Gabor 变换（DGT）。但 Gabor/STFT 有个根本限制：整个时频平面用同一分辨率的窗（图 5-69(a) 中每个格子形状完全相同）。而音频和图像处理常常希望常数 $Q$ （带宽/中心频率 = 常数）：高频用宽带、短时间间隔；低频用窄带、长时间间隔。这正是连续小波变换（CWT）能做到的：

CWT (τ, f) = \int_{- \infty}^{\infty} x (t) h (\frac{t - τ}{s}) d t (5-106)

(a) 傅立叶（常系数带宽）的频率分布

(b) 常数 Q 的频率分布

图 5-68 傅立叶（常系数带宽）和常数 Q 的频率分布

(a) STFT 网格

(b) 小波变换

图 5-69 线性调频信号的时间频率网格

式中的 $h (t)$ 称为小波（或子波），图 5-70 给出了 Morlet、Meyer、Daubechies 等几个典型小波。一个常用选择是：

h (t) = (e^{j 2 π k t} - e^{- k^{2} /2}) e^{- t^{2} /2} (5-107)

图 5-70 典型小波示例（一）

图 5-70 典型小波示例（二）

图 5-70 来自 Morlet、Meyer 和 Daubechies 的一些典型小波

它仍用高斯窗（保住”最优定位”），但通过尺度因子 $s$ 换用不同的时间和频率刻度；式(5-107)中引入 $e^{- k^{2} /2}$ 项的目的是让小波不带直流分量。其离散化称为离散 Morlet 变换（DMT），图 5-69(b) 给出了离散情形的时频网格。

例 5.26 线性调频信号的分析

考察一个频率随时间线性增加、幅值恒定的信号（线性调频信号，chirp）。普通傅立叶变换只能得到一条平坦的宽谱——所有频率都出现了，但完全丢失了”何时出现哪个频率”的信息。改用带高斯窗的 STFT（即 Morlet 变换），如图 5-71(a)，频率随时间爬升的趋势清晰可见；代价是高斯窗的计算量大。若换成 Haar 窗（矩形窗），计算量大幅下降（图 5-71(b)），但时频定位精度明显变差。这是一对典型的折中：定位质量 vs 计算成本。

图 5-71 线性调频信号的分析（一）

图 5-71 线性调频信号的分析（二）

(a) 离散 Morlet 变换

(b) Haar 变换

图 5-71 采用两种变换的线性调频信号的分析

DGT 和 DMT 定位虽好，计算量却都很大。有一种免乘法器的实现思路分两步：第一步，把高斯窗近似为若干（≥3 个）矩形函数的卷积；第二步，用多项式环上定义的代数整数高效实现单通带频率采样滤波器（FSF）。

不过更值得关注的是离散小波变换（DWT）：它同样符合人耳/人眼的常数 Q 感知模式，而且可以用 $O (n)$ 复杂度的算法高效计算——比 FFT 的 $O (n lo g n)$ 还低。

5.8.1 离散小波变换

缩放方程：连续小波与滤波器组的桥梁

实际应用中 DWT 通常限定为尺度因子 $a = 2$ 的二进形式，对应”常数 Q”分布靠滤波器树实现：把双通道滤波器组的低通输出不断级联下一级（图 5-72 给出三倍频程分解）。

那么什么样的连续小波才能用双通道 DWT 滤波器组实现？判据是缩放方程是否存在：

ϕ (t) = n \sum g [n] ϕ (2 t - n) (5-108)

ψ (t) = n \sum h [n] ϕ (2 t - n) (5-109)

其中 $ϕ (t)$ 是缩放函数， $ψ (t)$ 是实际的小波； $ϕ (t)$ 、 $ψ (t)$ 是连续函数，而 $g [n]$ 、 $h [n]$ 是采样序列（低通/高通滤波器系数，也可以是 IIR 的）。注意式(5-108)在形式上与分形的自相似性 $ϕ (t) = ϕ (a t)$ 很像——迭代它确实可能收敛到一个分形，但这通常不是我们想要的，我们需要光滑的小波。提高光滑度的办法：让滤波器在 $z = - 1$ （即 $π$ 处）拥有尽可能多的零点。

从滤波器”长出”小波：图形迭代

工程上最常见的方向是反过来的：先有滤波器 $g [n]$ （例如用算法 5.20 的半带设计造出的完美重构滤波器对），再构造对应的小波。做法是从方波函数（框函数）出发，迭代式(5-110)：

ϕ^{(k + 1)} (t) = n \sum g [n] ϕ^{(k)} (2 t - n) (5-110)

若迭代收敛到一个稳定的 $ϕ (t)$ ，小波就找到了。以 Haar 滤波器 ${1, 1}$ 为例：两个缩放后的框函数相加还是框函数，所以一次迭代就精确收敛，不存在”逐步成形”的过程。

例 5.27 长度为 4 的 Hutlet 滤波器

取 $g [n] = {1, 1, 1, 1}$ （Hutlet4 滤波器）。初始 $ϕ^{(1)} (t)$ 是 4 个框函数的加权和（图 5-73(a)）；每迭代一次，函数被放大 2 倍尺度并叠加，形状越来越圆滑。10 次迭代后得到一个光滑的梯形函数。再按式(5-78)的 QMF 关系造出实际小波 $ψ (t)$ ：两个三角形拼成的 Hutlet4 小波（图 5-74）。

(a) Hutlet4 迭代步骤 1

(b) Hutlet4 迭代步骤 2

(c) Hutlet4 迭代步骤 3

(d) Hutlet4 迭代步骤 10

图 5-73 Hutlet4 的迭代步骤 1、2、3 和 10（实线为 $ϕ^{(k + 1)} (t)$ ，点虚线为 $ϕ^{(k)} (2 t - n)$ ，短划线为理想的 Hut 函数）

顺带一提， $g [n] = {1, 1, 1, 1}$ 正是移动均值滤波器的冲激响应，可用一阶 CIC 滤波器实现。图 5-74 展示了这一类偶数长度系数滤波器的缩放函数与小波家族。

图 5-74 Hutlet 小波系列（一）

图 5-74 Hutlet 小波系列（二）

图 5-74 Hutlet 小波系列（三）

图 5-74 Hutlet 小波系列（四）

图 5-74 Hutlet 小波系列（五）

迭代也可能不收敛到光滑函数：把滤波器换成长度 5 的”移动均值滤波器” ${1, 1, 1, 1, 1}$ ，同样的迭代最终收敛到一个分形（图 5-75）。这个例子生动说明选择 $g [n]$ 的微妙之处：滤波器长度这样一个看似无关紧要的性质，就足以让结果在”光滑”与”完全不规则”之间翻转。

(a) Hutlet5 迭代步骤 1

(b) Hutlet5 迭代步骤 2

(c) Hutlet5 迭代步骤 4

(d) Hutlet5 迭代步骤 8

图 5-75 Hutlet5 的迭代步骤 1、2、4 和 8，序列收敛于一个分形

为什么二倍缩放正好对应滤波器树

式(5-108)里那个”2 倍缩放”为什么恰好匹配 DWT？用 5.1.1 节的 Noble 关系式

(↓ M) H (z) = H (z^{M}) (↓ M) (5-111)

把分析树中抽取器与滤波器的位置重新排列（图 5-76），就一目了然了。计算三级级联的等效冲激响应：

H (z) G (z) H (z^{2}) G (z) G (z^{2}) H (z^{4}) G (z) G (z^{2}) G (z^{4}) \leftrightarrow d_{1} [k /2] \leftrightarrow d_{2} [k /4] \leftrightarrow d_{3} [k /8] \leftrightarrow a_{3} [k /8]

把 $a_{3}$ 、 $d_{3}$ 的图形与图 5-70 的连续小波对比： $a_{3}$ 正是缩放函数的近似， $d_{3}$ 正是母小波的近似。滤波器树在结构上就是在”执行”缩放方程。当然，并非所有连续小波都走得通这条路——例如 Morlet 小波就找不到对应的缩放函数，无法用 DWT 实现。长度为 4 的 Daubechies 滤波器 DWT 的具体硬件设计，已在例 5.1（多相表示）和例 5.24（正交网格实现）中分别讨论过。

图 5-76 用 Noble 关系式重构 DWT 滤波器组

5.8.2 离散小波变换的应用

小波能做什么、不能做什么

小波在 20 世纪 80 年代流行时，研究者一度兴奋地以为找到了比 FFT 更通用的工具——当年 Sweldens 的”小波文摘”电子通讯有超过 2 万名订阅者。后来大家冷静下来：小波确实在某些应用上优于傅立叶方法、甚至开辟了新应用，但 FFT 擅长的快速卷积、快速相关、频谱估计，DWT 都不擅长。DWT 真正的强项是：

图像压缩：避免 DCT 方法（JPEG）的”块效应”，如 JPEG2000（第 10 章详述）；
信号去噪：同时保持信号的形状与性质；
图像增强；
间断点的检测与分类。

DWT 去噪原理

FFT 去噪的逻辑是：变换到频域，把小的傅立叶系数置零（赌它们是噪声而不是信号），再逆变换回来。DWT 去噪的逻辑完全类似（图 5-77）：做三级小波分析，把小的小波系数置零，再合成重建信号。

图 5-77 DFT/FFT 与 Haar DWT 去噪方法对比（一）

图 5-77 DFT/FFT 与 Haar DWT 去噪方法对比（二）

图 5-77 DFT/FFT 与 Haar DWT 去噪方法对比（三）

图 5-77 使用 DFT/FFT 以及 Haar DWT 去噪方法的对比

三级 Haar 方案如图 5-78 所示，相对标准树有两处实用修改：其一，因为用的是因果滤波器，合成树里要补一个延迟来对齐各层信号的相位；其二，MATLAB 的 wden 等工具默认用正交 Haar 滤波器，系数会引入 $1/ 2$ 的比例因子，硬件代价不小，实现时改用整数滤波器更划算。这些 SIMULINK 模型同时为后续的 HDL 设计生成测试数据。测试序列的选择有讲究：三角形输入下”细节”信号 $d_{k}$ 变化不大，看不出问题；最有效的激励是二阶型输入信号，可用三个框函数卷积生成。图 5-79(a) 是分析端数据，图 5-79(b) 是合成端做了延迟匹配的 $a_{3} (t) = s_{3} (t - τ)$ 与 $a_{2} (t) = s_{2} (t - τ)$ ——它们将在去噪 HDL 代码的测试与调试中充当参考。

图 5-78 三阶小波去噪方案

(a) 分析测试平台数据

(b) 合成数据

图 5-79 SIMULINK DWT 测试平台数据

效果对比：DWT vs FFT

用一个更真实的信号检验：取自 DCF77 无线电授时信号相位分量的伪随机序列（用于同步）。图 5-77(a) 是每位 16 个样本的原始信号，(b) 是叠加噪声后的信号。对比去噪结果 (c) 与 (d)：DWT 重建更好地保留了输入信号的锐利边沿，后级电路（例如锁相环）更容易同步；而 FFT 法把小系数置零本质上只是低通滤波，输出中锐边必然被抹平。结论很清晰：所选的小波应当是（无噪）输入信号的良好代表。

阈值怎么选

置零的判据——阈值——是去噪成败的关键。最简单的做法是按比例设阈，如取每个倍频程内最大小波系数的 10% 或 25%。更精细的做法是先估计信号方差：方差高的倍频程说明信息含量大，用小阈值以免误伤信号。Donoho 与 Johnstone 给出了一个针对高斯噪声的经典阈值：

T = σ 2 lo g_{e} (N) (5-112)

其中 $σ$ 是高斯白噪声的标准差， $N$ 为样本数。阈值随 $N$ 缓慢增长——数据集越大，高斯噪声中出现大幅值的概率越高；但 $lo g_{e} (N)$ 增长极慢： $N$ 从 512 翻到 1024， $T$ 只变化约 5%。式(5-112)真正的难点在于需要较好的噪声估计 $σ$ 。实用做法是取第一层细节信号 $d_{1}$ 的中值来估计：

\overset{σ}{^} = \frac{1}{0.6745} 中值 (d_{1} [m])_{0 \leq m \leq N /2} (5-113)

问题在于要对上千个数据求中值——硬件上意味着排序，代价很高。折中方案是采用更简单的阈值估计，下面的硬件实现就走这条路。

三级 DWT 去噪的硬件实现

原书 VHDL 改写为 Verilog：下面的 Verilog 代码实现了三级 Haar 小波去噪系统。结构上就是”图 5-78 的分析树 + 阈值 + 合成树”：输入经三级 Haar 分析得到 $d_{1} / a_{1}$ 、 $d_{2} / a_{2}$ 、 $d_{3} / a_{3}$ ；四路信号分别与外部输入的阈值 $t 4 d 1 \sim t 4 a 3$ 比较，绝对值低于阈值的系数清零（这里用最简单的硬件友好判据，而不做式(5-113)的中值排序）；然后经合成树（含延迟对齐）重建输出 $y$ 。各层细节/近似信号及多处内部节点被引出，方便调试。整个系统只有一个 8 状态循环计数器驱动三层时钟使能，Haar 的和/差运算只需一次加法与一次减法，完全不需要乘法器。

// 三级 Haar DWT 去噪系统（原书 VHDL dwtden 改写为 Verilog）
module dwtden #(
    parameter D1L = 28,   // d1 延迟线长度
    parameter D2L = 10    // d2 延迟线长度
) (
    input  wire        clk,      // 系统时钟
    input  wire        reset,    // 异步复位
    input  wire signed [15:0] x_in,   // 系统输入
    input  wire signed [15:0] t4d1, t4d2, t4d3, t4a3, // 各层阈值
    output wire signed [15:0] d1_out, a1_out,  // 第一层细节/近似
    output wire signed [15:0] d2_out, a2_out,  // 第二层细节/近似
    output wire signed [15:0] d3_out, a3_out,  // 第三层细节/近似
    output wire signed [15:0] s3_out, a3up_out, d3up_out, // L3 调试信号
    output wire signed [15:0] s2_out, s3up_out, d2up_out, // L2 调试信号
    output wire signed [15:0] s1_out, s2up_out, d1up_out, // L1 调试信号
    output reg  signed [15:0] y_out            // 系统输出
);
 
    reg  [2:0] count;                // 模 8 循环计数器
    reg  signed [15:0] x, xd;        // 输入延迟
    wire signed [15:0] a1, d1, a2, d2, a3, d3;   // 分析滤波器输出
    reg  signed [15:0] d1t, d2t, d3t, a3t;       // 阈值处理后
    wire signed [15:0] a1up, a3up, d3up;
    reg  signed [15:0] a1upd, s3upd, a3upd, d3upd;
    reg  signed [15:0] a1d, a2d;
    reg  ena1, ena2, ena3;           // 三层时钟使能
    reg  t1, t2, t3;                 // 层切换触发器
    reg  signed [15:0] s2, s3up, s3, d2syn;
    reg  signed [15:0] s1, s2up, s2upd;
    reg  signed [15:0] d2upd [0:D2L-1];  // d2 延迟线
    reg  signed [15:0] d1upd [0:D1L-1];  // d1 延迟线
 
    integer i;
 
    // ---------- 模 8 计数器生成三层使能与相位 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            count <= 0; t1 <= 0; t2 <= 0; t3 <= 0;
        end else begin
            count <= count + 1;
            t1 <= (count == 0);                  // 每隔一拍
            t2 <= (count == 1);                  // 每隔两拍
            t3 <= (count == 3);                  // 每隔四拍
        end
    end
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            ena1 <= 0; ena2 <= 0; ena3 <= 0;
        end else begin
            ena1 <= t1; ena2 <= t2; ena3 <= t3;
        end
    end
 
    // ---------- 输入延迟 ----------
    always @(posedge clk or posedge reset) begin
        if (reset) begin
            x  <= 0;
            xd <= 0;
        end else begin
            x  <= x_in;
            xd <= x;
        end
    end
 
    // ---------- 一级 Haar 分析：和/差 ----------
    assign d1 = (x  - xd) >>> 1;   // (x[n] - x[n-1]) / 2
    assign a1 = (x  + xd) >>> 1;   // (x[n] + x[n-1]) / 2
    assign d2 = (a1 - a1d) >>> 1;
    assign a2 = (a1 + a1d) >>> 1;
    assign d3 = (a2 - a2d) >>> 1;
    assign a3 = (a2 + a2d) >>> 1;
 
    always @(posedge clk) begin
        if (ena1) a1d <= a1;
        if (ena2) a2d <= a2;
    end
 
    // ---------- 阈值处理：小于阈值置零 ----------
    always @(posedge clk) begin
        d1t <= ((d1 > t4d1) || (d1 < -t4d1)) ? d1 : 0;
        d2t <= ((d2 > t4d2) || (d2 < -t4d2)) ? d2 : 0;
        d3t <= ((d3 > t4d3) || (d3 < -t4d3)) ? d3 : 0;
        a3t <= ((a3 > t4a3) || (a3 < -t4a3)) ? a3 : 0;
    end
 
    // ---------- 细节信号延迟线（对齐合成相位） ----------
    always @(posedge clk) begin
        if (ena1) begin
            d1upd[0] <= d1t;
            for (i = 1; i < D1L; i = i + 1) d1upd[i] <= d1upd[i-1];
        end
        if (ena2) begin
            d2upd[0] <= d2t;
            for (i = 1; i < D2L; i = i + 1) d2upd[i] <= d2upd[i-1];
        end
    end
 
    // ---------- 一级 Haar 合成：插值后和/差 ----------
    assign d3up  = d3t;
    assign a3up  = a3t;
    assign s3upd = a3up + d3up;              // 第三层合成
    assign d2up  = d2upd[D2L-1];
    assign s3    = s3upd;
    assign d2syn = s3 - d2up;                // 等效合成运算
    assign s2upd = a2d + d2syn;              // 第二层合成
    assign d1up  = d1upd[D1L-1];
    assign s2up  = s2upd;
    assign s1    = s2up + d1up;              // 第一层合成
    assign y_out_next = s1 - d1up;
 
    // ---------- 输出寄存器与调试引脚 ----------
    always @(posedge clk or posedge reset) begin
        if (reset)
            y_out <= 0;
        else
            y_out <= y_out_next;
    end
 
    assign d1_out = d1;  assign a1_out = a1;
    assign d2_out = d2;  assign a2_out = a2;
    assign d3_out = d3;  assign a3_out = a3;
    assign s3_out = s3;  assign a3up_out = a3up; assign d3up_out = d3up;
    assign s2_out = s2;  assign s3up_out = s3up; assign d2up_out = d2up;
    assign s1_out = s1;  assign s2up_out = s2up; assign d1up_out = d1up;
 
    wire signed [15:0] s2_unused = s2;
    wire signed [15:0] s3up_unused = s3up;
 
endmodule

说明两点实现细节。第一，各层抽取率不同（1、2、4），用同一个模 8 计数器派生三路使能信号，让三层各自只在属于自己的节拍上更新，硬件上等价于多速率时钟域但避免了时钟树分叉。第二，由于分析用因果滤波器， $d_{1}$ 、 $d_{2}$ 支路必须各插入一条延迟线（长度 D1L、D2L）才能与近似支路在合成端对齐相位——这正是 5.8.2 节开头强调的”合成树补延迟”修改。与直接按正交 Haar 实现（每层引入 $1/ 2$ 缩放）相比，整数 Haar 实现全程只有移位、加法和减法，FPGA 资源开销小得多。这套结构跑通后，换上更复杂的应用信号（如前面的 DCF77 伪随机序列），配合 Donoho 阈值或简化阈值估计，就能得到完整的小波去噪系统。

小波去噪：三级 DWT 滤波器组的完整实现

前几个小节我们分别搭建了分析滤波器组、合成滤波器组，也认识了提升格式。现在把所有零件组装成一台完整的机器——小波去噪器（DWT denoiser）。它要做的事情非常直观：

是什么：把输入信号做三级离散小波分解，得到细节系数 $d_{1}, d_{2}, d_{3}$ 和近似系数 $a_{3}$ ；把幅度小于阈值的系数直接清零（认为是噪声），再做三级小波合成，输出”干净”的信号 $y$ 。
为什么：信号（比如 DCF77 授时码那样的脉冲序列）能量集中在少数大的小波系数上，而白噪声的能量均匀摊在所有系数上。只要”掐掉”那些幅度很小的系数，就能在几乎不损伤信号形状的前提下压制噪声。
不用会怎样：如果直接在时域上做滑动平均之类的滤波，脉冲边沿会被钝化；小波去噪则能”既去噪、又保形”，这一点在后面的仿真波形里看得非常清楚。

本设计选用 Haar 小波，因为它的分析/合成滤波器只有加减法，一个乘法器都不用，是理解整个去噪流水线的最佳教学载体。三级分解对应的树状结构是：第 1 级工作在最高采样率，第 2 级只在第 1 级输出有效时工作，第 3 级再慢一半。这种”每级速率减半”的多采样率结构，正是第 5 章多级信号处理思想的一次综合演练。

整个设计的节拍由一个循环长度为 8 的状态机统一指挥，下面逐个模块拆开讲解。

用 FSM 产生分级使能信号

为什么需要使能信号而不是三个时钟

多级结构里，第 2、3 级的滤波器只在”有新数据到来”的那个时钟沿才允许更新寄存器，其余时刻必须保持不动。硬件上给每一级单独生成一个时钟当然可以，但多时钟域会带来时钟偏斜、跨域同步等一堆麻烦。更稳妥的做法是：全系统共用一个时钟，用使能信号（enable）来区分”这一拍对谁有效”。这就是前几节反复出现的”用使能代替分频时钟”的工程惯例。

节拍安排的规律是：完整循环 8 个时钟周期。

第 1 级（最快速）：在 count = 1, 3, 5, 7 这四拍有效，即每 2 拍做一次 Haar 分析；
第 2 级：在 count = 1, 5 这两拍有效，即每 4 拍一次；
第 3 级：只在 count = 5 这一拍有效，即每 8 拍一次。

注意使能信号是在 count 等于这些值之后的时钟沿才变成真——寄存器输出天然带一拍延迟，恰好让滤波器在数据稳定后的那一拍才动作。

原书 VHDL 改写为 Verilog（可综合风格，异步复位 + 同步使能生成）：

// FSM: 控制整个系统，按 clk 节拍采样
// count 循环 0..7；ena1/ena2/ena3 分别是三级分析的使能
always @(posedge clk or posedge reset) begin
    if (reset) begin                       // 异步复位
        count <= 3'd0;
        ena1  <= 1'b0;
        ena2  <= 1'b0;
        ena3  <= 1'b0;
    end else begin
        if (count == 3'd7)
            count <= 3'd0;
        else
            count <= count + 1'b1;
 
        ena1 <= (count==3'd1) || (count==3'd3) ||
                (count==3'd5) || (count==3'd7);   // 第 1 级使能
        ena2 <= (count==3'd1) || (count==3'd5);   // 第 2 级使能
        ena3 <= (count==3'd5);                    // 第 3 级使能
    end
end

对照真值表可以自查一遍：8 拍中 ena1 拉高 4 次、ena2 拉高 2 次、ena3 拉高 1 次，各级速率恰好呈 2:1 递减，与三级小波树完全一致。

Haar 分析滤波器组

两个加减法就是一层小波分解

Haar 小波的”滤波器”简单到令人发指：对相邻两个输入样本 $x$ 和它的一拍延迟 $x_{d}$ ，做

d = x - x_{d} (细节：高频差分), a = x + x_{d} (近似：低频求和)

差分 $d$ 响应信号的快速变化，求和 $a$ 保留信号的缓变轮廓——这正是把信号拆成”高频 + 低频”两路。三级级联就是把这一对运算作用在逐级减半的速率上：第 1 级吃原始输入 $x$ ，第 2 级吃第 1 级的近似输出 $a_{1}$ ，第 3 级吃 $a_{2}$ 。

原书 VHDL 改写为 Verilog（行为级风格）：

// Haar 分析滤波器组：三个使能信号控制的三对滤波器
always @(posedge clk or posedge reset) begin
    if (reset) begin                       // 异步清零
        x   <= 0;  xd  <= 0;
        d1t <= 0;  a1  <= 0;  a1d <= 0;
        d2t <= 0;  a2  <= 0;  a2d <= 0;
        d3t <= 0;  a3t <= 0;
    end else begin
        x  <= x_in;                        // 输入寄存
        xd <= x;                           // 一拍延迟
 
        if (ena1) begin                    // 第 1 级分析
            d1t <= x - xd;
            a1  <= x + xd;
            a1d <= a1;
        end
        if (ena2) begin                    // 第 2 级分析
            d2t <= a1 - a1d;
            a2  <= a1 + a1d;
            a2d <= a2;
        end
        if (ena3) begin                    // 第 3 级分析
            d3t <= a2 - a2d;
            a3t <= a2 + a2d;
        end
    end
end

读这段代码时抓住两点：其一，每个 if (enaX) 块内部就是一对”减法/加法 + 寄存”，结构与上一节完全相同，只是输入源逐级更换；其二，a1d、a2d 是把本级近似输出再延迟一拍，供下一级做差分/求和使用，相当于把”相邻两个样本”的语义搬运到了慢速率域。

小波系数的阈值处理

去噪的关键一步：小于阈值就归零

分解完成后得到 4 路系数： $d_{1}, d_{2}, d_{3}$ 和 $a_{3}$ 。去噪就发生在这一步——每路系数配一个阈值输入端口（ $t 4 d 1$ 、 $t 4 d 2$ 、 $t 4 d 3$ 、 $t 4 a 3$ ），绝对值不超过阈值的系数一律置零：

d_{k} = {d_{k}^{t}, 0, ∣ d_{k}^{t} ∣ > T_{d_{k}} 其他

这四条是纯组合逻辑，不需要时钟：

原书 VHDL 改写为 Verilog（并发赋值，注意 Verilog 需自己表达”绝对值”判断）：

// 对 d1、d2、d3、a3 做阈值处理：小系数清零
assign d1 = ((d1t > t4d1) || (d1t < -t4d1)) ? d1t : 0;
assign d2 = ((d2t > t4d2) || (d2t < -t4d2)) ? d2t : 0;
assign d3 = ((d3t > t4d3) || (d3t < -t4d3)) ? d3t : 0;
assign a3 = ((a3t > t4a3) || (a3t < -t4a3)) ? a3t : 0;

阈值不必设得很讲究，工程上常取”可观测最大小波系数的 25% 或 50%“这种粗略比例，后文的仿真正是这样做的。阈值越大，去掉的”噪声”越多，但信号细节也可能被误伤；阈值取零则相当于直通。这种可调性通过输入端口保留下来，方便实验对比。

合成滤波器组：先抽取后插值的实现技巧

难点在哪里

合成要做的事在数学上很清楚：把阈值化后的系数上采样（每两个点插一个零）、通过合成滤波器、逐级相加，最后重建信号。但直接按”上采样器 + 滤波器”画框图去写 HDL，会碰到两个别扭的地方：

上采样产生的零点会白白流过滤波器，浪费一半的运算；
分析端是”降速率”，合成端是”升速率”，两边的流水线延迟必须严格对齐，否则重建出来的信号相位是错的。

本设计采用一个聪明的等价变换：“上采样 + 滤波”反过来写成”滤波（在低速数据上）+ 零插入”。具体做法是用一个翻转触发器 $t$ 每两个时钟翻转一次： $t = 1$ 的那一拍让真实数据通过， $t = 0$ 的那一拍强制写零——零插入就这么用一条 if/else 实现了，连上采样器模块都不需要。

重建公式与延迟匹配

Haar 合成是分析的对偶：对上采样后的近似 $a_{u p}$ 和细节 $d_{u p}$ ，

s = \frac{a _{u p} + a _{u p, d} - d _{u p} + d _{u p, d}}{2}

这里的除以 2 在硬件里就是右移一位（Verilog 中用 >>> 1 完成算术右移，不会引入除法器）。式中的 $d_{u p}$ 相对 $a_{u p}$ 多经过了一段流水线延迟，所以要用一组移位寄存器把 $d_{u p}$ 逐拍后移（代码中的 d1upd、d2upd 数组），直到它与 $a_{u p}$ 的时间戳对齐——这就是所谓的”延迟匹配”。D1L、D2L 两个参数分别记录第 1、2 级合成路径的流水线级数，编译时按实际结构设置。

原书 VHDL 改写为 Verilog（行为级风格，D1L/D2L 为参数）：

// "下采样后跟上采样"通过隔拍置零来实现
integer k;
always @(posedge clk or posedge reset) begin
    if (reset) begin                       // 异步清零
        t1 <= 1'b0; t2 <= 1'b0; t3 <= 1'b0;
        s3up <= 0;  s3upd <= 0;
        d3up <= 0;  a3up <= 0;  a3upd <= 0;  d3upd <= 0;
        s3 <= 0; s2 <= 0; s1 <= 0;
        s2up <= 0; s2upd <= 0;
        for (k = 0; k <= D2L+1; k = k + 1)     // 清零数组，与 s3up 匹配
            d2upd[k] <= 0;
        for (k = 0; k <= D1L+1; k = k + 1)     // 清零数组，与 s2up 匹配
            d1upd[k] <= 0;
    end else begin
        t1 <= ~t1;                         // 第 1 级翻转触发器
        if (t1) begin
            d1upd[0] <= d1;                // 真实数据拍
            s2up     <= s2;
        end else begin
            d1upd[0] <= 0;                 // 零插入拍
            s2up     <= 0;
        end
        s2upd <= s2up;
        for (k = 1; k <= D1L+1; k = k + 1)     // 延迟匹配 s2up
            d1upd[k] <= d1upd[k-1];
        s1 <= (s2up + s2upd - d1upd[D1L] + d1upd[D1L+1]) >>> 1;
 
        if (ena1) begin
            t2 <= ~t2;                     // 第 2 级翻转触发器
            if (t2) begin
                d2upd[0] <= d2;
                s3up     <= s3;
            end else begin
                d2upd[0] <= 0;
                s3up     <= 0;
            end
            s3upd <= s3up;
            for (k = 1; k <= D2L+1; k = k + 1) // 延迟匹配 s3up
                d2upd[k] <= d2upd[k-1];
            s2 <= (s3up + s3upd - d2upd[D2L] + d2upd[D2L+1]) >>> 1;
        end
 
        if (ena2) begin                    // 第 3 级合成
            t3 <= ~t3;                     // 翻转触发器
            if (t3) begin
                d3up <= d3;
                a3up <= a3;
            end else begin
                d3up <= 0;
                a3up <= 0;
            end
            a3upd <= a3up;
            d3upd <= d3up;
            s3 <= (a3up + a3upd - d3up + d3upd) >>> 1;
        end
    end
end

三个级别结构完全同构，只是速率不同：第 1 级每个时钟都工作，第 2 级被 ena1 门控（每 2 拍一次），第 3 级被 ena2 门控（每 4 拍一次）。第 3 级因为只差一拍延迟，用两个普通寄存器 a3upd、d3upd 就够了，不需要数组。整条合成链从最慢的第 3 级开始，一路”喂”到最快的第 1 级，方向恰好与分析链相反。

观测端口与设计规模

最后把内部关键节点引到输出端口，方便仿真时观察每一级的工作情况：

原书 VHDL 改写为 Verilog（并发赋值，仅作测试观测）：

// 提供部分测试信号作为输出
assign a1_out    = a1;
assign d1_out    = d1;
assign a2_out    = a2;
assign d2_out    = d2;
assign a3_out    = a3;
assign d3_out    = d3;
assign a3up_out  = a3up;
assign d3up_out  = d3up;
assign s3_out    = s3;
assign s3up_out  = s3up;
assign d2up_out  = d2upd[D2L];
assign s2_out    = s2;
assign s1_out    = s1;
assign s2up_out  = s2up;
assign d1up_out  = d1upd[D1L];
assign y_out     = s1;                     // 重建信号就是系统输出

资源方面这份报告很有代表性：整个三级小波去噪器只用了 879 个 LE（逻辑单元），没有用到任何嵌入式乘法器——Haar 滤波只有加减法，这是它对 FPGA 最友好的地方。时序上，使用 TimeQuest 缓慢 85C 模型测得 $F_{ma x} = 120.93 MHz$ ，对去噪这类应用绰绰有余。

仿真验证：从 MATLAB 激励到 ModelTech 波形

仿真流程分三步走，值得初学者完整体会一遍：

第一步，在 MATLAB 里生成激励并写成 HDL 仿真命令。 测试信号 $s (k)$ 是 DCF77 的 40 位码序列，用循环把每个采样点”打印”成一条 force 指令：

for k=1:L
    fprintf(fid, 'force x_in %d %dns\r\n', round(256*s(k)), ...
    100*k);
end

注意两个细节：浮点值乘以 256 后取整，相当于预留 8 位小数精度的定点缩放（ $Q 8$ 格式），这比在 HDL 里处理浮点划算得多；时间步长取 100 ns，即仿真时钟 10 MHz。

第二步，在 ModelSim 中跑仿真。 由于数据量太大，不可能把所有测试端口都显示出来，只挑选 $x$ 、 $d_{1}$ 、 $d_{2}$ 、 $d_{3}$ 、 $a_{3}$ 和最终输出 $y$ ，并用 Format→Analog(automatic) 把数字信号渲染成模拟波形，肉眼才能看出形状。

图 5-80 DWT 去噪的仿真波形（原书为 VHDL 仿真，Verilog 版结果一致）

第三步，对比不同阈值下的输出。 仿真最上面一组波形显示了 FSM 的控制信号与阈值信号；阈值从零（不去噪）逐步增加到小波系数最大值的 25% 和 50%。仔细对比 $y_{o u t}$ 的三段波形可以发现一个关键现象：50% 阈值时，噪声几乎被完全抹平，而输入数字脉冲的形状依然完好。这正是小波去噪”既去噪又保形”的直观证据，也是我们在本节开头承诺要看到的效果。

还有一点实现层面的提醒：图 5-78 的 Simulink 框图里画了较长的延迟链来对齐各支路信号，因为那里的上/下采样是用同步寄存器实现的；本设计的 Verilog 代码把这些延迟隐藏在 d1upd/d2upd 移位链和各级寄存器里，形式不同，本质相同——多采样率系统里，延迟对齐永远是最容易出错、也最值得反复检查的环节。

延伸案例：CIC 插值器的位增长与 GC4114

作为多级信号处理的最后一个延伸，看一个工业级多级系统的位宽设计问题。GC4114 是一款经典的多级数字变频芯片，内含 4 级 CIC 插值器和一个可变插值因子 $R$ 。多级 CIC 每过一级，数据位宽就要增长一次，如果位宽规划不当，要么溢出错码，要么白白浪费寄存器。Hogenauer 给出了逐级增益（位增长率）的解析公式：

G_{k} = ⎩ ⎨ ⎧ 2^{k}, \frac{2 ^{2 S - k} ( R D ) ^{k - S}}{R}, k = 1, 2, \dots, S k = S + 1, S + 2, \dots, 2 S (5-114)

其中 $D$ 是梳状延迟， $S$ 是级数， $R$ 是插值因子，位增长率 $B_{k} = lo g_{2} (G_{k})$ 。

拿 GC4114 的参数（ $S = 4$ ， $D = 1$ ）实际算一遍，体会这个公式的用法。前 4 级对应积分器部分， $B_{k} = k$ ，即 1、2、3、4 位；后 4 级对应梳状器部分，把 $2 S - k = 8 - k$ 与 $(R D)^{k - 4} / R = R^{k - 5}$ 合并，可得指数 $(8 - k) + lo g_{2} R \cdot (k - 5)$ 。代入三组 $R$ ：

级号 $k$	$R = 8$	$R = 32$	$R = 16384$
1~4（积分器）	1, 2, 3, 4	1, 2, 3, 4	1, 2, 3, 4
5	3	3	3
6	5	7	16
7	7	11	29
8（输出级）	9	15	42

规律一目了然：从第 5 级开始，每过一级位增长就增加 $lo g_{2} R$ 位，最终输出级的位增长率决定了整条链的寄存器总位宽。 $R = 16384$ 时输出要长到 42 位，这就是为什么大插值因子的 CIC 滤波器必须精心剪裁中间位宽。设计练习中要求据此搭建 $R = 32$ 、16 位输入的 4 级 CIC 插值器（ $D = 1$ ），先用 MATLAB 仿真确认每级的位增长，再写 HDL 并与仿真结果逐点对照，最后在 Cyclone II 与 Cyclone IV E 两代器件上对比 $F_{ma x}$ 和资源占用。

图 5-81(a) GC4114 CIC 插值器仿真结果 (a)：输入与各级中间波形

图 5-81(b) GC4114 CIC 插值器仿真结果 (b)：中间级位增长情况

图 5-81(c) GC4114 CIC 插值器仿真结果 (c)：输出波形与位宽

图 5-81 GC4114 CIC 插值器的仿真结果（总览）

回顾本章练习的整体脉络：半带滤波器与完美重构滤波器组的构造（练习 5.1~~5.5）、均匀 DFT 滤波器组的多相实现与计算量评估（5.6~~5.8）、提升格式的完美重构证明与 Daubechies 滤波器实现（5.9~~5.10）、CSD 与简化加法器图优化的半带/多相设计（5.12~~5.14）、B 样条与 Farrow 结构的分数速率变换器（5.18~~5.23），直到本节的 CIC 位增长与 GC4114 多级插值器（5.24~~5.25）。这些练习共同覆盖了第 5 章”多级信号处理”从理论到 FPGA 落地的完整链条，建议按顺序动手完成，并把每个设计的 $F_{ma x}$ 与 LE/乘法器/M9K 占用记录成表格，长期积累下来就是一份属于自己的 FPGA 滤波器设计经验库。

description	采样率变换的艺术：抽取/插值、多相分解、Hogenauer CIC 滤波器、小波变换——过采样系统的核心工具箱。
tags	fpga数字信号处理, 数字信号处理的fpga实现

嵌入式开发之路

探索

📘 第 5 章：多级信号处理

自测一下

第5章 多级信号处理

5.1 抽取和插值

什么是抽取，为什么需要它

带通信号的完整频带条件

插值与 D/A 转换

5.1.1 Noble 恒等式

5.1.2 用有理数因子进行采样速率转换

5.2 多相分解

为什么要多相分解

分解方法

例 5.1 多相抽取滤波器

5.2.1 递归 IIR 抽取器

5.2.2 快速 FIR 滤波器

例 5.2 快速 FIR 滤波器

运算量账本：快在”次数”，贵在”硬件”

算法 5.3 快速 FIR 滤波器

什么时候停止分解

5.3 Hogenauer CIC 滤波器

5.3.1 单级 CIC 案例研究

本征频率测试：积分器溢出为何无害

用 RNS 拓宽位宽

5.3.2 多级 CIC 滤波器理论

例 5.4 三级 CIC 抽取器 I

5.3.3 幅值与混叠畸变

5.3.4 Hogenauer “剪除”理论

例 5.5 三级 CIC 抽取器 II

例5.6 带位剪除的三级CIC抽取器

5.3.5 CIC 的 RNS 设计

5.3.6 CIC 补偿滤波器设计

5.4 多级抽取器

Goodman-Carey 半带滤波器设计法

例 5.8 多级半带滤波器抽取器

5.5 作为通频带抽取器的频率采样滤波器

5.6 任意采样速率转换器的设计

例 5.9 R=0.75 的速率变换器 I

算法 5.10 应用 FFT 的有理数速率变换

例 5.11 R=0.75 的速率变换器 II

5.6.1 分数延迟速率变换

例 5.12 R=0.75 的速率变换器 III

5.6.1（续）三个 sinc 滤波器的乘累加与仿真结果

5.6.2 多项式分数延迟设计

什么时候该换一种思路

拉格朗日多项式如何求系数

Farrow 结构：用 Horner 模式高效求值

例5.13 R=0.75 的速率变换器 IV

与 sinc 方案的可扩展性对比

5.6.3 基于 B 样条的分数速率变换器

为什么需要 B 样条

B 样条逼近与 B 样条插值

用 FIR 滤波器逼近补偿滤波器

三次 B 样条的 Farrow 方程

5.6.4 MOMS 分数速率变换器

逼近的”阶数”：一个被忽视的设计参数

三种 MOMS 的补偿滤波器比较

C-MOMS 的 Farrow 方程

例 5.14 R=0.75 的速率变换器 V

速率变换方法的局限性：方波与吉布斯现象

5.7 滤波器组

5.7.1 均匀 DFT 滤波器组

从原型滤波器到 DFT 滤波器组

多相分解带来高效实现

完美重构条件

例5.15 DFT 滤波器组

5.7.2 双通道滤波器组

什么是双通道滤波器组，为什么需要它

例5.16 双通道Haar滤波器组 I

完美重构的一般条件

定理 5.17 完美的重构

例 5.18 双通道 Haar 滤波器组 II

简化设计：与混叠无关的结构

定理 5.19 与混叠无关的双通道滤波器组

算法 5.20 完美重构的双通道滤波器组

例 5.21 采用 F3 的完美重构滤波器组

推论 5.22 半带滤波器的因式分解

5.7.3 实现双通道滤波器组

1. 多相双通道滤波器组

第5章多级信号处理

例 5.9 $R = 0.75$ 的速率变换器 I

例 5.11 $R = 0.75$ 的速率变换器 II

例 5.12 $R = 0.75$ 的速率变换器 III

例5.13 $R = 0.75$ 的速率变换器 IV