定点DSP上IIR滤波器优化：规范型结构防溢出与高效实现

1. 项目概述：在定点DSP上驯服IIR滤波器

如果你在TMS320C54x这类经典的定点DSP上实现过IIR滤波器，大概率经历过这样的场景：精心设计的滤波器在仿真里表现完美，一旦上板运行，输出要么是满屏的噪声，要么信号微弱到几乎看不见。问题根源往往不是算法错了，而是定点DSP有限的动态范围（16位数据，32位累加器）与IIR滤波器反馈路径带来的潜在增益之间，发生了不可调和的矛盾——数据溢出了。

IIR滤波器因其反馈结构，能以比FIR滤波器更低的阶数实现更陡峭的滚降，在实时性要求高的DSP应用中极具吸引力。然而，这份“高效率”在定点世界是有代价的。传统的直接I型（Direct Form I）实现，为了绝对稳定而进行的激进输入缩放，常常把信号压得过小，导致输出信噪比恶化，最后不得不在硬件上用一颗运放把信号再“放大”回来，这多少有些讽刺。

本文要分享的，正是我在多个C54x项目中反复验证过的一套“组合拳”：通过一种优化的直接I型结构（常被称为规范型或Canonical Direct Form I），在软件层面巧妙地重组计算流，从而在保证稳定、避免溢出的前提下，最大化输出信号的幅度。这不仅仅是省掉一颗运放的成本，更是让算法在资源受限的嵌入式环境中，真正发挥出IIR应有的效能。接下来，我会拆解其背后的原理、手把手展示代码实现，并附上那些只有踩过坑才知道的调试心得。

2. IIR滤波器设计原理与定点实现的挑战

在浮点处理器上设计IIR滤波器，工程师可以暂时忽略数据范围的烦恼，专注于传递函数和频响曲线。但一旦切换到TMS320C54x这类定点DSP，设计思路必须转变：我们必须同时是算法设计师和“资源调度员”，时刻关注数据通路上每一个节点的动态范围。

2.1 IIR滤波器结构选型：为何是直接I型？

IIR滤波器有多种实现结构，如直接I型、直接II型（规范型）、级联型、并联型。在C54x上，级联的二阶级联（Biquad）结构是事实上的标准。它将一个高阶滤波器分解为多个二阶节的乘积，每个二阶节独立稳定，便于单独进行缩放和量化误差分析，鲁棒性远优于直接实现高阶差分方程。

那么，对于每个二阶节，是选直接I型还是直接II型？直接II型（一个延迟单元）内存占用最少，但它将零极点路径的量化误差耦合在一起，对系数量化误差更敏感。直接I型（两个独立的延迟线）虽然多用了一组延迟单元，但它将零点（前馈路径）和极点（反馈路径）的计算在结构上分离开。这种分离带来一个关键好处：我们可以更清晰、独立地分析并控制反馈环路带来的增益，从而为定点化缩放提供直观的依据。在C54x这种内存资源并非极度稀缺（片内DARAM通常有数K字）而计算稳定性优先的场合，直接I型是更稳妥的起点。

2.2 定点DSP的溢出难题：反馈路径是“罪魁祸首”

TMS320C54x的CPU核心是一个32位的累加器（ACC），它接收来自16位数据总线（来自内存或寄存器）和16位系数相乘的32位结果。虽然累加器是32位，但最终结果需要存回16位的数据内存。问题就出在反馈回路上。

IIR滤波器的极点（分母系数）决定了其频率响应的峰值。在某些频率点，滤波器的增益可能远大于1。当一个信号（即使是经过缩放的）进入反馈回路，被大于1的系数反复乘加，其数值很容易增长到超出32位累加器能表示的范围（特别是在进行多个二阶节级联时，增益是累积的）。一旦发生溢出，CPU的溢出模式（OVM）如果被使能，会将其饱和处理为最大正值或负值，这相当于向信号中注入了大量的非线性失真，严重时滤波器会进入持续饱和状态，完全失效。

2.3 传统应对策略及其弊端：过度缩放与硬件补偿

最直接的防溢出方法是输入缩放（Input Scaling）。即在数据进入滤波器之前，先将其右移若干位（相当于除以2^N），强行降低信号幅度，为后续的增益留出“净空（Headroom）”。

这种方法简单粗暴，但副作用明显。假设为了防止最坏情况下的溢出，我们需要将输入右移9位（除以512）。这意味着输入信号的有效精度从16位骤降到7位（16-9），量化噪声大幅增加。更糟糕的是，经过整个滤波器链后，输出信号幅度可能只有满量程的几百分之一，信噪比（SNR）严重恶化。为了得到可用的输出电平，传统做法是在DAC之前增加一级模拟运算放大器进行硬件放大。这增加了BOM成本、板级面积和功耗，背离了使用高效IIR滤波器简化系统的初衷。

3. 优化策略：重组计算流的规范直接I型

我们需要一种更智能的方法，既能防止溢出，又能最大限度地保留信号的幅度和精度。这就是原文中提到的“Efficient IIR Filter Design”核心思想——通过改变二阶节的级联顺序和计算流程，实现动态范围的均衡分布。

3.1 从直接I型到规范直接I型：延迟单元的共享

首先看一个标准二阶直接I型（图1）。它有两组延迟单元（z^-1），分别用于输入序列和输出序列。根据线性时不变系统的性质，零点和极点的子系统是串联的，其顺序可以交换（图2）。交换后，两个子系统中间的信号m(n)和d(n)在数学上是等价的。

关键在于图3：当我们交换零极点后，发现前一个系统的输出延迟，正好是后一个系统所需要的输入延迟。因此，两组延迟单元可以合并为一组。这就是“规范型”（Canonical Form）的由来，它用最少的延迟元件实现了相同的系统函数。在硬件描述语言中，这节省了寄存器；在DSP软件中，这节省了数据存储单元（RAM）。

3.2 级联二阶节的优化：交叉组合计算

将规范型的思路扩展到多个二阶节级联的高阶滤波器（图4），就产生了关键的优化洞见。观察图4(b)的结构，它不再是简单的“零点块→极点块→零点块→极点块…”。

它的计算流程是：

1. 第一个二阶节只计算其零点部分（前馈路径，b11, b12）。
2. 随后，将第一个二阶节的极点部分（反馈路径，a11, a12）与第二个二阶节的零点部分（b21, b22）合并计算。它们共享同一组延迟单元d(n-1), d(n-2)。
3. 以此类推，中间的二阶节都将其极点与后一个二阶节的零点合并。
4. 最后一个二阶节只计算其极点部分。

这样做的核心优势在于：每个中间节点（延迟单元d(n)）的信号幅度得到了“再均衡”。在传统结构中，信号先经过零点（可能衰减），再经过极点（可能放大），幅度变化剧烈。在优化结构中，零点和极点的计算被交错安排，使得信号在通过每个二阶节时，其增益波动被平滑了。极点带来的放大效应，可以部分抵消前级零点可能带来的衰减，反之亦然。

3.3 对定点实现的益处：更温和的缩放需求

这种平滑的增益分布，直接降低了对输入缩放倍数的要求。在原文的对比实验中，传统结构需要将输入左移9位（缩小512倍）来防止溢出，而优化后的结构仅需左移5位（缩小32倍）。

计算一下动态范围的改善：假设输入是16位有符号整数（范围-32768到32767）。传统方法缩放后，有效数据位只剩16-9=7位，峰值约为±64。量化噪声相对增大。优化方法缩放后，有效数据位有16-5=11位，峰值约为±1024。后者比前者多保留了4个有效位，相当于动态范围提升了约24dB，输出信号幅度也大了16倍，通常已能满足后级ADC或DAC的输入要求，从而省去了那颗额外的运放。

4. TMS320C54x上的高效代码实现与剖析

理论需要代码落地。下面我们结合附录B的代码，深入解析在C54x上实现这一优化结构的技巧。我将使用更清晰的注释和分段说明。

4.1 内存与系数表规划

C54x的汇编效率极高，但需要精细规划。优化结构要求系数在内存中按特定顺序排列。

; 假设我们实现一个4阶椭圆低通滤波器，分解为2个二阶节（N=2） N .set 2 ; 二阶节的数量 .bss d, 3*N ; 延迟缓冲区：每个二阶节需要d(n), d(n-1), d(n-2)三个状态 ; 优化结构下，相邻节共享状态，所以总共需要3*N个单元，而非4*N .bss X, 1 ; 输入缓冲区 .bss Y, 1 ; 输出缓冲区 .data ; 系数表排列顺序是此算法的关键！ ; 格式为: [B2, B1, B0, A2, A1] 对于每个二阶节。 ; 注意：这里存储的A1是实际A1值的一半，因为后续计算利用了MAC指令的自动左移。 table: ; 第1个二阶节 (Section #1) .word 19381 ; B2 .word -23184 ; B1 .word 19381 ; B0 .word -26778 ; A2 (实际存储 -a2) .word 29529 ; A1/2 (实际存储 a1/2，符号为负？这里原文是正，但计算时按负处理，需结合代码逻辑) ; 第2个二阶节 (Section #2) .word 11363 ; B2 .word -20735 ; B1 .word 11363 ; B0 .word -30497 ; A2 .word 31131 ; A1/2

关键点解析：

1. 系数顺序：B2, B1, B0, A2, A1/2。这种排列是为了配合后续的指令流，实现流水线化计算。
2. A1存储为一半：因为C54x的乘法器在FRCT（小数模式）置1时，会自动将乘积左移1位（乘以2），以修正两个Q15小数相乘产生的额外符号位。因此，如果我们希望计算d(n-1)*a1，可以预先将a1除以2存入，然后利用自动左移恢复。这节省了一次显式的移位操作。
3. 状态缓冲区d：长度为3*N。在计算过程中，AR3指针会在这个缓冲区中滑动，指向上一个二阶节计算出的d(n)，它同时也是下一个二阶节计算所需的d(n-2)状态（经过延迟操作后），实现了状态的共享。

4.2 主滤波循环：指令级优化详解

核心滤波循环是性能关键。C54x的RPTB（块重复）指令和MAC（乘累加）指令是主力。

INLOOP: STM #d+7, AR3 ; AR3指向延迟缓冲区末端（具体初始位置需根据节数调整） STM #table, AR4 ; AR4指向系数表开头 ; 第一步：计算第一个二阶节的零点部分（只有前馈路径） MPY *AR4+, *AR3-, A ; A = d(n-2)_old * B2_s1 MAC *AR4+, *AR3, A ; A += d(n-1)_old * B1_s1 DELAY *AR3- ; 内存移动：d(n-2) = d(n-1)，为下一节准备 MAC *AR4+, *AR3, A ; A += d(n)_old * B0_s1 DELAY *AR3 ; 内存移动：d(n-1) = d(n) ; 此时，A中包含了经过第一个零点处理后的中间结果 ; AR3指向了最新的d(n)位置（即将被新输入覆盖） ; 第二步：读入新输入并缩放 PORTR 100H, *AR3 ; 从端口读入新样本x(n) LD *AR3, B ; 加载到B寄存器 STH B, 11, *AR3- ; 左移5位（缩放1/32）后存回d(n)，并修改AR3 ; 注意：这里左移11位是指将高16位（STH）存储时，源操作数B左移11位后的高16位。 ; 更常见的做法是使用`LD *AR3, 5, A`或`ST #0, ASM`结合移位。此处代码是一种特定写法。 ; 第三步：循环处理中间的二阶节（合并零极点计算） STM #N-2, BRC ; 设置块重复计数器，对于N=2，这里循环0次 RPTB ELOOP-1 ; 开始块重复 LOOP: ; 对于多于2个二阶节的情况，此循环处理中间节 ; 计算当前节的极点部分（反馈）并更新状态d(n) MAC *AR4+, *AR3-, A ; A += d(n-2) * (-A2_current) （累加上一节的零点结果） MAC *AR4, *AR3, A ; A += d(n-1) * (-A1_current) （注意AR4未更新，指向A1/2） MAC *AR4+, *AR3-, A ; 完成A1计算并更新AR4。此时A = 零点结果 + 反馈项 STH A, *AR3+0 ; 将结果饱和处理到高16位，存入作为新的d(n) ; 计算当前节的零点部分，为下一节或输出准备 MPY *AR4+, *AR3-, A ; A = d(n-2) * B2_next MAC *AR4+, *AR3, A ; A += d(n-1) * B1_next DELAY *AR3- ; d(n-2) = d(n-1) MAC *AR4+, *AR3, A ; A += d(n) * B0_next DELAY *AR3- ; d(n-1) = d(n) ELOOP: ; 第四步：处理最后一个二阶节的极点部分 MAC *AR4+, *AR3-, A ; A += d(n-2) * (-A2_last) MAC *AR4, *AR3, A ; A += d(n-1) * (-A1_last) MAC *AR4+, *AR3, A ; 完成最后一个反馈项累加 DELAY *AR3 ; 最后的延迟移位（可选，为下一次迭代准备） STH A, *AR3 ; 最终的滤波结果y(n) PORTW *AR3, 200h ; 输出结果 B INLOOP ; 处理下一个样本

代码精要解析：

1. 状态指针舞蹈：AR3在延迟缓冲区d中的移动是算法核心。它精确地指向每个二阶节所需的d(n), d(n-1), d(n-2)，并通过DELAY指令（本质是数据移动）来更新状态，模拟了z^-1延迟操作。
2. 系数指针同步：AR4与AR3完美配合，顺序访问B2, B1, B0, A2, A1。MAC和MPY指令中的*AR4+实现了系数的自动递增。
3. 累加器A的复用：在整个计算过程中，累加器A承载了中间结果。它先被初始化为第一个零点部分的输出，然后在循环中不断累加上一个极点和下一个零点的贡献，最后加上最后一个极点部分，形成最终输出。这种设计最大限度地减少了中间结果的存储和重载。
4. 块重复循环：对于高阶滤波器（N>2），中间的二阶节计算被包装在RPTB循环中，极大地提高了代码密度和执行效率。

4.3 初始化与模式设置

正确的初始化是稳定运行的前提。

begin: STM #1111111110100000b, PMST ; 设置PMST，例如IPTR指向0x80， MP/MC=0（微计算机模式） STM #0010001100000000b, ST1 ; 设置ST1， BRAF=0（块重复无效）， CPL=0（DP直接寻址）， etc. STM #0, SWWSR ; 零等待状态（根据实际存储器速度调整） SSBX OVM ; 至关重要！使能溢出饱和模式。当溢出发生时，结果饱和为最大正值(7FFF FFFF)或负值(8000 0000)，而不是绕回。 SSBX FRCT ; 使能小数模式。两个Q15数相乘后自动左移1位，结果保持在Q31格式。 SSBX SXM ; 使能符号扩展模式。数据加载到累加器时进行符号扩展。 STM #d, AR3 RPTZ A, #7 ; 将累加器A清零，并重复下一条指令8次 STL A, *AR3+ ; 将延迟缓冲区d全部初始化为0 STM #2, AR0 ; 设置AR0为2，用于某些偏移寻址（本例中未直接使用）

关键设置说明：

• OVM=1：这是IIR滤波器在定点DSP上的“生命线”。没有它，溢出会导致数据从正最大值跳变到负最大值，产生灾难性失真。
• FRCT=1：DSP处理的是小数（-1 ≤ x < 1），而非整数。此模式确保乘法结果正确。你的滤波器系数（如0.591）在程序中实际存储为0.591 * 32768 ≈ 19381（Q15格式）。
• 状态清零：滤波器初始状态必须为零，否则会有瞬态响应。

5. 实验配置、结果分析与性能对比

理论分析和代码实现之后，我们需要用实验数据说话。原文提供了一个基于特定指标的对比，我们可以将其深化。

5.1 滤波器指标与设计工具

原文指定了一个低通滤波器：通带截止频率200Hz，阻带截止频率500Hz。假设采样频率为Fs=2000Hz（满足奈奎斯特定理）。

• 设计工具：通常使用MATLAB的ellip函数（椭圆滤波器）或butter（巴特沃斯）等。例如，在MATLAB中设计一个4阶椭圆低通滤波器：

  Fs = 2000; % 采样率 Fpass = 200; % 通带截止 Fstop = 500; % 阻带截止 Apass = 1; % 通带衰减，单位dB Astop = 40; % 阻带衰减，单位dB [N, Wn] = ellipord(Fpass/(Fs/2), Fstop/(Fs/2), Apass, Astop); [b, a] = ellip(N, Apass, Astop, Wn); [sos, g] = tf2sos(b, a); % 转换为二阶节形式

  得到二阶节系数后，需要将其量化为Q15格式：coeff_q15 = round(sos * 32768)。注意，a系数（极点部分）需要取负并存储为-a，因为我们的差分方程实现是y(n) = b0*x(n)+... - a1*y(n-1)-a2*y(n-2)。

5.2 实验结果深度解读

原文中的对比表格是核心成果。我们来逐一解读：

分析：

1. 输出电平：优化IIR的输出达到了与FIR滤波器相当的满量程水平（-16000~16000），而传统IIR的输出被严重压缩。这直观证明了优化结构有效恢复了信号幅度。
2. 输入缩放：缩放从1/512减少到1/32，动态范围损失从9位减少到5位，信噪比提升了约24dB（20*log10(512/32) ≈ 24dB）。这是一个质的飞跃。
3. 性能与资源：
   • 周期数：优化IIR（22周期）甚至比传统IIR（27周期）更快。这是因为优化结构减少了冗余的加载/存储操作，计算流更紧凑。它比高阶FIR（36周期）快得多，体现了IIR的效率优势。
   • 代码大小：优化IIR（61字）比传统IIR（56字）略大，因为循环控制和系数访问模式稍复杂，但差异很小。
   • 数据内存：两者相同（8字），都使用了共享延迟线的优化。
4. 与FIR对比：要达到相似的滤波特性（陡峭的过渡带），FIR需要32阶，是IIR的8倍。这导致其数据内存需求（32字状态缓冲区）和计算量（36周期）都显著更高。对称FIR利用系数对称性减少了乘法次数（28周期），但内存占用不变。

5.3 实际测试波形观察

在示波器或逻辑分析仪上观察：

1. 时域：输入一个满幅度的正弦波（例如1kHz，在阻带内）。传统IIR输出幅度极小，几乎淹没在噪声中；优化IIR输出幅度清晰可见，且波形光滑，无明显失真。
2. 频域（通过DAC输出接频谱分析仪）：观察滤波器的幅频响应。优化IIR和传统IIR在通带、阻带衰减上应基本一致，但传统IIR的通带内信号底噪会明显更高，这是过度缩放导致量化噪声增大的结果。

6. 常见问题、调试技巧与进阶优化

在实际工程中，把代码跑起来只是第一步，让它稳定、鲁棒地工作才是挑战。

6.1 溢出与饱和调试

即使采用了优化结构和缩放，在极端输入或特定频率下，溢出仍可能发生。

• 调试方法：
  1. 监视OVC计数器：C54x的ST0寄存器中有溢出计数器（OVC）。在调试阶段，可以在滤波循环后检查并清零它。如果OVC持续增加，说明发生了饱和，可能需要进一步微调缩放因子或检查系数。
  2. 使用仿真器观察累加器：在CCS等IDE中，单步执行并观察累加器A的值。关注它在关键计算节点（如完成一个二阶节计算后）是否接近饱和边界（如±0.75 * 2^31）。
  3. 极限测试：输入一个幅值等于最大正数（0x7FFF）的直流或低频信号，观察输出是否稳定。这是最坏情况测试。

6.2 系数量化与极限环

定点系数量化可能改变滤波器的零极点位置，轻微影响频响，甚至可能将极点推到单位圆外导致不稳定（尽管椭圆滤波器设计时本身是稳定的）。

• 预防措施：
  1. 系数缩放：在量化前，确保每个二阶节的极点部分（分母系数）满足|a1| + |a2| < 1（对于直接I型），这是一个保证该二阶节稳定的充分条件（并非必要）。如果量化后不满足，可能需要微调滤波器设计（如稍微增加通带波纹）或使用更稳健的结构（如耦合型）。
  2. 使用高精度累加：C54x的32位累加器提供了足够的保护位。确保在关键求和步骤（如MAC指令序列）前，累加器已被正确初始化或加载，避免残留大值。
  3. 注意零输入极限环：对于极低电平或零输入，由于舍入和溢出饱和，IIR滤波器输出可能不会衰减到零，而是在几个固定值间振荡。优化结构对此有一定改善，但若应用对空闲噪声极度敏感，可能需要考虑加入微小的输出抖动（Dithering）或使用更高位宽的处理器。

6.3 性能优化进阶

当处理速度成为瓶颈时，可以考虑：

1. 使用循环展开：对于固定阶数（如4阶）的滤波器，可以完全展开循环，消除RPTB的开销。代码体积会增加，但速度更快。
2. 利用双MAC单元：某些C54x的衍生型号（如C54xx）具有双MAC单元。可以重新组织计算，将部分并行的乘加运算安排到同一周期，但汇编代码会非常复杂。
3. 数据放在片内DARAM：确保延迟缓冲区d和系数表table位于零等待状态的片内RAM中。访问片外RAM会引入等待周期，严重拖慢速度。
4. 使用C语言内联汇编：对于复杂项目，可以用C语言编写框架，将核心滤波循环用asm()语句嵌入，在可读性和性能间取得平衡。

6.4 从汇编到C：可移植性考虑

虽然汇编效率最高，但现代开发更注重可维护性和可移植性。你可以用C54x的C编译器实现同样的算法：

#pragma DATA_SECTION(d, ".bss:d_buffer") #pragma DATA_SECTION(coeff, ".const:coeff_table") static short d[6]; // 3 * N const short coeff[10] = {19381, -23184, 19381, -26778, 29529, 11363, -20735, 11363, -30497, 31131}; // B2,B1,B0,A2,A1/2 for each section short iir_filter(short input) { long acc; short *p_d = &d[5]; // 指向缓冲区末端，模拟汇编指针初始化 const short *p_c = coeff; short i; // 第一步：第一个零点节 acc = (long)(*p_d--) * (*p_c++); // d(n-2)*B2 acc += (long)(*p_d) * (*p_c++); // d(n-1)*B1 p_d--; // 模拟DELAY: d(n-2) = d(n-1) (实际需要数据移动) // 这里需要手动移动数据，简化起见，先忽略... acc += (long)(*p_d) * (*p_c++); // d(n)*B0 p_d--; // 模拟DELAY // 第二步：缩放输入并存入 input >>= 5; // 缩放1/32 (实际是算术右移，需注意符号) *p_d = input; // 第三步：循环处理中间节（本例N=2，无中间节，略） // 第四步：最后一个极点节 acc += (long)(*p_d--) * (*p_c++); // + d(n-2)*(-A2) acc += (long)(*p_d) * (*p_c); // + d(n-1)*(-A1) (注意A1存储为一半) acc <<= 1; // 补偿A1/2的存储，因为C编译器不会自动像FRCT那样左移 acc += (long)(*p_d) * (*p_c++); // 完成累加 (原代码中合并了) // 更新延迟线（需要循环移动数据，此处简化） // ... 实际需要将d数组元素向后移动 // 饱和处理 if (acc > 0x7FFF0000L) acc = 0x7FFF0000L; else if (acc < (long)0x80000000) acc = (long)0x80000000; return (short)(acc >> 16); // 取高16位作为输出 }

C代码清晰，但编译器生成的效率通常低于手写汇编。关键循环仍需用汇编优化。此外，在C中实现DELAY操作（数据移动）需要仔细处理，避免使用memcpy造成低效。

最后，分享一个我调试时的实用技巧：在项目初期，可以先用MATLAB或Python生成一个浮点参考模型，并导出测试向量。然后在CCS中，将同样的测试向量加载到DSP的输入缓冲区，运行你的汇编滤波器，比较输出结果。这能快速定位是算法逻辑错误、系数量化问题还是溢出问题。记住，在定点DSP的世界里，对数据流动和范围保持敬畏，是写出稳定、高效代码的不二法门。