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)的结构它不再是简单的“零点块→极点块→零点块→极点块…”。它的计算流程是第一个二阶节只计算其零点部分前馈路径b11, b12。随后将第一个二阶节的极点部分反馈路径a11, a12与第二个二阶节的零点部分b21, b22合并计算。它们共享同一组延迟单元d(n-1), d(n-2)。以此类推中间的二阶节都将其极点与后一个二阶节的零点合并。最后一个二阶节只计算其极点部分。这样做的核心优势在于每个中间节点延迟单元d(n)的信号幅度得到了“再均衡”。在传统结构中信号先经过零点可能衰减再经过极点可能放大幅度变化剧烈。在优化结构中零点和极点的计算被交错安排使得信号在通过每个二阶节时其增益波动被平滑了。极点带来的放大效应可以部分抵消前级零点可能带来的衰减反之亦然。3.3 对定点实现的益处更温和的缩放需求这种平滑的增益分布直接降低了对输入缩放倍数的要求。在原文的对比实验中传统结构需要将输入左移9位缩小512倍来防止溢出而优化后的结构仅需左移5位缩小32倍。计算一下动态范围的改善假设输入是16位有符号整数范围-32768到32767。传统方法缩放后有效数据位只剩16-97位峰值约为±64。量化噪声相对增大。优化方法缩放后有效数据位有16-511位峰值约为±1024。后者比前者多保留了4个有效位相当于动态范围提升了约24dB输出信号幅度也大了16倍通常已能满足后级ADC或DAC的输入要求从而省去了那颗额外的运放。4. TMS320C54x上的高效代码实现与剖析理论需要代码落地。下面我们结合附录B的代码深入解析在C54x上实现这一优化结构的技巧。我将使用更清晰的注释和分段说明。4.1 内存与系数表规划C54x的汇编效率极高但需要精细规划。优化结构要求系数在内存中按特定顺序排列。; 假设我们实现一个4阶椭圆低通滤波器分解为2个二阶节N2 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关键点解析系数顺序B2, B1, B0, A2, A1/2。这种排列是为了配合后续的指令流实现流水线化计算。A1存储为一半因为C54x的乘法器在FRCT小数模式置1时会自动将乘积左移1位乘以2以修正两个Q15小数相乘产生的额外符号位。因此如果我们希望计算d(n-1)*a1可以预先将a1除以2存入然后利用自动左移恢复。这节省了一次显式的移位操作。状态缓冲区d长度为3*N。在计算过程中AR3指针会在这个缓冲区中滑动指向上一个二阶节计算出的d(n)它同时也是下一个二阶节计算所需的d(n-2)状态经过延迟操作后实现了状态的共享。4.2 主滤波循环指令级优化详解核心滤波循环是性能关键。C54x的RPTB块重复指令和MAC乘累加指令是主力。INLOOP: STM #d7, 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 ; 设置块重复计数器对于N2这里循环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, *AR30 ; 将结果饱和处理到高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 ; 处理下一个样本代码精要解析状态指针舞蹈AR3在延迟缓冲区d中的移动是算法核心。它精确地指向每个二阶节所需的d(n), d(n-1), d(n-2)并通过DELAY指令本质是数据移动来更新状态模拟了z^-1延迟操作。系数指针同步AR4与AR3完美配合顺序访问B2, B1, B0, A2, A1。MAC和MPY指令中的*AR4实现了系数的自动递增。累加器A的复用在整个计算过程中累加器A承载了中间结果。它先被初始化为第一个零点部分的输出然后在循环中不断累加上一个极点和下一个零点的贡献最后加上最后一个极点部分形成最终输出。这种设计最大限度地减少了中间结果的存储和重载。块重复循环对于高阶滤波器N2中间的二阶节计算被包装在RPTB循环中极大地提高了代码密度和执行效率。4.3 初始化与模式设置正确的初始化是稳定运行的前提。begin: STM #1111111110100000b, PMST ; 设置PMST例如IPTR指向0x80 MP/MC0微计算机模式 STM #0010001100000000b, ST1 ; 设置ST1 BRAF0块重复无效 CPL0DP直接寻址 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用于某些偏移寻址本例中未直接使用关键设置说明OVM1这是IIR滤波器在定点DSP上的“生命线”。没有它溢出会导致数据从正最大值跳变到负最大值产生灾难性失真。FRCT1DSP处理的是小数-1 ≤ x 1而非整数。此模式确保乘法结果正确。你的滤波器系数如0.591在程序中实际存储为0.591 * 32768 ≈ 19381Q15格式。状态清零滤波器初始状态必须为零否则会有瞬态响应。5. 实验配置、结果分析与性能对比理论分析和代码实现之后我们需要用实验数据说话。原文提供了一个基于特定指标的对比我们可以将其深化。5.1 滤波器指标与设计工具原文指定了一个低通滤波器通带截止频率200Hz阻带截止频率500Hz。假设采样频率为Fs2000Hz满足奈奎斯特定理。设计工具通常使用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 实验结果深度解读原文中的对比表格是核心成果。我们来逐一解读滤波器类型算法阶数输入缩放输出电平程序大小 (字)数据大小 (字)周期数/输出传统IIR椭圆49位左移 (1/512)-1000 ~ 1000 (小)56827优化IIR椭圆45位左移 (1/32)-16000 ~ 16000 (正常)61822FIR凯泽窗32无-16000 ~ 16000 (正常)653236对称FIR凯泽窗32无-16000 ~ 16000 (正常)533228分析输出电平优化IIR的输出达到了与FIR滤波器相当的满量程水平-16000~16000而传统IIR的输出被严重压缩。这直观证明了优化结构有效恢复了信号幅度。输入缩放缩放从1/512减少到1/32动态范围损失从9位减少到5位信噪比提升了约24dB20*log10(512/32) ≈ 24dB。这是一个质的飞跃。性能与资源周期数优化IIR22周期甚至比传统IIR27周期更快。这是因为优化结构减少了冗余的加载/存储操作计算流更紧凑。它比高阶FIR36周期快得多体现了IIR的效率优势。代码大小优化IIR61字比传统IIR56字略大因为循环控制和系数访问模式稍复杂但差异很小。数据内存两者相同8字都使用了共享延迟线的优化。与FIR对比要达到相似的滤波特性陡峭的过渡带FIR需要32阶是IIR的8倍。这导致其数据内存需求32字状态缓冲区和计算量36周期都显著更高。对称FIR利用系数对称性减少了乘法次数28周期但内存占用不变。5.3 实际测试波形观察在示波器或逻辑分析仪上观察时域输入一个满幅度的正弦波例如1kHz在阻带内。传统IIR输出幅度极小几乎淹没在噪声中优化IIR输出幅度清晰可见且波形光滑无明显失真。频域通过DAC输出接频谱分析仪观察滤波器的幅频响应。优化IIR和传统IIR在通带、阻带衰减上应基本一致但传统IIR的通带内信号底噪会明显更高这是过度缩放导致量化噪声增大的结果。6. 常见问题、调试技巧与进阶优化在实际工程中把代码跑起来只是第一步让它稳定、鲁棒地工作才是挑战。6.1 溢出与饱和调试即使采用了优化结构和缩放在极端输入或特定频率下溢出仍可能发生。调试方法监视OVC计数器C54x的ST0寄存器中有溢出计数器OVC。在调试阶段可以在滤波循环后检查并清零它。如果OVC持续增加说明发生了饱和可能需要进一步微调缩放因子或检查系数。使用仿真器观察累加器在CCS等IDE中单步执行并观察累加器A的值。关注它在关键计算节点如完成一个二阶节计算后是否接近饱和边界如±0.75 * 2^31。极限测试输入一个幅值等于最大正数0x7FFF的直流或低频信号观察输出是否稳定。这是最坏情况测试。6.2 系数量化与极限环定点系数量化可能改变滤波器的零极点位置轻微影响频响甚至可能将极点推到单位圆外导致不稳定尽管椭圆滤波器设计时本身是稳定的。预防措施系数缩放在量化前确保每个二阶节的极点部分分母系数满足|a1| |a2| 1对于直接I型这是一个保证该二阶节稳定的充分条件并非必要。如果量化后不满足可能需要微调滤波器设计如稍微增加通带波纹或使用更稳健的结构如耦合型。使用高精度累加C54x的32位累加器提供了足够的保护位。确保在关键求和步骤如MAC指令序列前累加器已被正确初始化或加载避免残留大值。注意零输入极限环对于极低电平或零输入由于舍入和溢出饱和IIR滤波器输出可能不会衰减到零而是在几个固定值间振荡。优化结构对此有一定改善但若应用对空闲噪声极度敏感可能需要考虑加入微小的输出抖动Dithering或使用更高位宽的处理器。6.3 性能优化进阶当处理速度成为瓶颈时可以考虑使用循环展开对于固定阶数如4阶的滤波器可以完全展开循环消除RPTB的开销。代码体积会增加但速度更快。利用双MAC单元某些C54x的衍生型号如C54xx具有双MAC单元。可以重新组织计算将部分并行的乘加运算安排到同一周期但汇编代码会非常复杂。数据放在片内DARAM确保延迟缓冲区d和系数表table位于零等待状态的片内RAM中。访问片外RAM会引入等待周期严重拖慢速度。使用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; // 第三步循环处理中间节本例N2无中间节略 // 第四步最后一个极点节 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的世界里对数据流动和范围保持敬畏是写出稳定、高效代码的不二法门。