IIR滤波器与相关函数的多采样并行优化：从算法原理到DSP汇编实现
1. 多采样编程：从单核到多核的思维跃迁
在数字信号处理（DSP）的世界里，性能就是生命线。无论是手机通话的降噪、音乐播放器的均衡器，还是雷达信号的分析，背后都是滤波器在实时处理着海量的数据流。传统的单样本处理模式，就像一个熟练的工人，一次只能加工一个零件，虽然稳定，但面对现代应用对实时性和能效的苛刻要求，就显得力不从心了。多采样编程技术，正是为了解决这个瓶颈而生。它不再是一个个地处理数据，而是像一条现代化的流水线，一次抓取、处理、输出多个样本，将硬件的并行计算能力压榨到极致。这不仅仅是代码的改写，更是一种设计思维的彻底转变——从串行思维转向并行思维，从关注单个操作的延迟转向关注整体数据流的吞吐量。我经历过从单ALU代码到为四ALU甚至更多核心架构重写算法的全过程，这种思维转变带来的性能提升是颠覆性的，但其中的设计门道，也远比想象中复杂。
2. 核心原理：为什么反馈让IIR滤波器的并行化如此不同？
要理解多采样编程的优化，必须先吃透IIR滤波器的本质。它与我们更熟悉的FIR（有限脉冲响应）滤波器有根本区别。FIR滤波器就像一个没有记忆的加权求和器，当前输出只依赖于当前及过去的一系列输入。它的计算是纯粹的前馈，结构上天然适合展开和并行。
而IIR滤波器引入了“反馈”。它的当前输出y(n)，不仅依赖于当前输入x(n)和过去的输入，更关键的是，它还依赖于过去的输出y(n-1),y(n-2)... 这个反馈回路带来了两个核心特征：一是可以用更少的阶数实现更陡峭的滤波特性，效率更高；二是给并行计算带来了巨大的挑战。
为什么是挑战？想象一下，你要同时计算四个输出：y(n),y(n+1),y(n+2),y(n+3)。对于y(n+1)，它的计算公式里包含y(n)。但y(n)本身也正在被计算中！这就形成了一个“先有鸡还是先有蛋”的数据依赖循环。在严格的串行计算中，这不成问题，因为我们是按顺序算完y(n)再算y(n+1)。但在并行计算中，四个输出我们希望同时算，这种直接的依赖会形成硬件难以处理的数据冲突，导致流水线停滞。
因此，多采样IIR滤波器的设计核心，就在于打破这个即时的数据依赖，通过数学上的展开和重构，将原本交织在一起的串行计算，转化为一批可以同时进行的独立乘加运算。这需要巧妙地重新组织计算公式，让四个并行ALU在计算时，所需要的“过去输出”数据，都是已经可用的、来自更早时间点的值，而不是本次循环中正在计算的其他结果。这个过程，就是算法从串行到并行的“重构”。
3. 内核结构设计：从通用内核到基础内核的演化之路
多采样优化的精髓，集中体现在“内核”的设计上。所谓内核，就是被循环执行的核心计算片段。我们的目标是用最少的指令，喂饱所有的ALU。
3.1 理想与现实的差距：通用内核的困境
最初，我们会根据并行计算公式，画出一个理想的“通用内核”。以四样本并行（Quad-Sample）IIR滤波器为例，我们同时计算y(n)到y(n+3)。观察计算过程会发现，在计算某个乘积项时，例如涉及系数C1的项，它会被用于多个输出的计算中。一个最直接的想法是：加载一次系数C1和一个延迟值D，然后让四个ALU分别完成C1*D对四个不同输出的累加。
这就是“通用内核”的雏形：一次加载两个操作数（一个系数，一个延迟值），然后执行四次乘累加（MAC）操作。理论上，这完美利用了并行性。但这里隐藏着一个关键操作：寄存器复制。为了让C1在后续的计算中继续被使用（比如下一轮计算C1会变成C2的角色），我们需要在加载新系数时，把旧系数在寄存器间进行拷贝移动（例如，将C3复制到C4，C2复制到C3，C1复制到C2，然后加载新的C1）。
问题来了：在像StarCore SC140这样的DSP架构中，单条指令能执行的数据移动（Load/Move）操作数量是有限的。一个通用内核需要2次加载（新系数、新延迟值）加上3次寄存器拷贝，总共5次数据移动，但硬件可能只支持2-3次。这就造成了“数据搬运带宽”不足，内核无法在单周期内完成，理想结构在现实中碰壁。
3.2 折叠与复用：基础内核的诞生
既然硬件不允许我们在一个内核里做那么多事，我们就扩大内核的视野。解决方案是将多个通用内核组合成一个更大的“基础内核”。
我们不再试图在一个小循环内完成所有数据的复用和搬运，而是把眼光放长。假设滤波器有8个系数（C1到C8）。我们设计一个基础内核，它连续执行4个通用内核的计算量。在这个更大的范围内，我们重新安排加载和计算顺序。
核心技巧是“折叠”加载操作。我们提前为未来几步计算加载好需要的数据。例如，在基础内核的开始，我们就一次性加载多个将来需要的系数和延迟值到不同的寄存器中。在后续的计算中，我们不再进行寄存器间的拷贝，而是直接通过改变寄存器在计算公式中的引用顺序来实现数据的“复用”。
具体来说，我们准备四组寄存器来存放系数（C1,C2,C3,C4）和延迟值。在基础内核的四行核心计算中，我们让同一组数据在不同的计算行中扮演不同的角色。比如，第一行我们用C4和D计算四个输出，第二行我们用C3和新的D计算，但此时C4寄存器里的数据在计算中被引用的位置发生了变化（可能从计算y(n)变为计算y(n+1)）。通过精心设计这个“旋转引用”的模式，我们完全消除了耗时的寄存器拷贝操作。
最终得到的基础内核，每行只需要2次加载操作（加载一个新的延迟值和一个新的系数），然后执行4次MAC。这完美匹配了SC140架构单指令4MAC和2次数据移动的能力，使得每个基础内核都能在单周期内执行，实现了理论上的峰值性能。这个从“通用内核”到“基础内核”的转化过程，是多采样编程优化中最具艺术性的部分，它要求开发者对数据流有极深的洞察力。
4. IIR滤波器多采样实现全解析
理论需要代码来落地。我们以一个8阶IIR滤波器、一次处理4个样本为例，拆解其完整的C语言仿真实现。这段代码清晰地展示了基础内核的运作机制。
4.1 初始化与指针管理
首先，我们定义滤波器的参数和状态。IirSize=8表示8个系数和8个延迟单元。DataBlockSize=40是要处理的数据块大小。延迟线数组Delay初始化全零（或历史值），它存储的是过去的输出y(n-1)到y(n-8)。
指针管理是关键。我们使用CoefPtr和DelayPtr两个指针，以模寻址的方式在系数数组和延迟线环形缓冲区中循环。初始化时，它们指向数组末尾（IirSize-1）。DecMod函数实现指针的递减和模环绕，确保指针在数组范围内循环。
#define DataBlockSize 40 #define IirSize 8 double Coef[IirSize] = {0.4, -0.3, 0.25, -0.20, -0.15, 0.10, -0.10, 0.05}; double Delay[IirSize] = {0}; // 初始化为0 int CoefPtr = IirSize - 1; int DelayPtr = IirSize - 1;4.2 四样本处理主循环
主循环以4为步进遍历输入数据。每次循环处理DataIn[i],DataIn[i+1],DataIn[i+2],DataIn[i+3]四个样本。
for (i = 0; i < DataBlockSize; i += 4) { sum1 = DataIn[i]; // y(n) 的初始值为 x(n) sum2 = DataIn[i+1]; // y(n+1) 的初始值为 x(n+1) sum3 = DataIn[i+2]; // y(n+2) 的初始值为 x(n+2) sum4 = DataIn[i+3]; // y(n+3) 的初始值为 x(n+3) // ... 后续计算 }4.3 核心计算：展开的基础内核
这是算法的心脏。为了消除反馈依赖，计算被分为两部分：前馈部分（使用过去的延迟值）和反馈部分（处理本次并行计算产生的内部依赖）。
第一步：前馈计算（处理历史延迟值）我们从延迟线中依次取出y(n-8),y(n-7)... 等历史值，并与对应的系数相乘累加到四个输出和上。这个过程完全并行，因为所有用到的历史值都是已知的。
// 以第一组计算为例，对应系数 C8, C7, C6, C5 C4 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); sum1 += C4 * D; // 只有 sum1 需要 C8*y(n-8) C3 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); sum1 += C3 * D; // sum1: C7*y(n-7) sum2 += C4 * D; // sum2: C8*y(n-7) 注意这里C4就是上一步的C8，实现了复用！ // 以此类推，加载 C2/C1，并完成 sum1, sum2, sum3 对历史延迟值的累加。第二步：基础内核循环（处理主要系数）完成历史值计算后，系数指针指回C4（实际上是原系数数组的C4），开始核心循环。这个循环实现了我们设计的基础内核。
for (j = 0; j < IirSize/4 - 1; j++) { // 对于8阶，循环 (8/4 - 1)=1 次 // 行1: 使用 C1, C2, C3, C4 和当前 D sum1 += C1 * D; sum2 += C2 * D; sum3 += C3 * D; sum4 += C4 * D; D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); C4 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); // 行2: 使用 C4, C1, C2, C3 和新的 D。注意寄存器角色的“旋转”！ sum1 += C4 * D; sum2 += C1 * D; sum3 += C2 * D; sum4 += C3 * D; D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); C3 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); // 行3: 使用 C3, C4, C1, C2 sum1 += C3 * D; sum2 += C4 * D; sum3 += C1 * D; sum4 += C2 * D; D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); C2 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); // 行4: 使用 C2, C3, C4, C1 sum1 += C2 * D; sum2 += C3 * D; sum3 += C4 * D; sum4 += C1 * D; D = Delay[DelayPtr]; DelayPtr = DecMod(DelayPtr, IirSize); C1 = Coef[CoefPtr]; CoefPtr = DecMod(CoefPtr, IirSize); }仔细观察这个循环：每行只有2次加载（一个D，一个Cx），却完成了4次MAC。C1-C4在四行中像齿轮一样旋转被使用，完全避免了寄存器拷贝。这正是“基础内核”的精髓。
第三步：处理内部反馈经过前两步，sum1到sum4已经包含了输入和所有历史延迟值与系数的乘积。但IIR的反馈还包括本次计算的输出之间的依赖（如y(n+1)依赖于y(n)）。这部分在循环外以级联方式处理：
// sum1 (即 y(n)) 现已计算完成 sum2 += C1 * sum1; // y(n+1) 加上 C1*y(n) sum3 += C2 * sum1; // y(n+2) 加上 C2*y(n) sum4 += C3 * sum1; // y(n+3) 加上 C3*y(n) // sum2 (即 y(n+1)) 现已计算完成 sum3 += C1 * sum2; // y(n+2) 加上 C1*y(n+1) sum4 += C2 * sum2; // y(n+3) 加上 C2*y(n+1) // sum3 (即 y(n+2)) 现已计算完成 sum4 += C1 * sum3; // y(n+3) 加上 C1*y(n+2) // sum4 (即 y(n+3)) 计算完成第四步：更新延迟线最后，将新计算出的四个输出sum1-sum4写回延迟线，覆盖最旧的四个延迟值，为下一组数据块的处理做好准备。
Delay[DelayPtr] = sum1; DelayPtr = DecMod(DelayPtr, IirSize); Delay[DelayPtr] = sum2; DelayPtr = DecMod(DelayPtr, IirSize); Delay[DelayPtr] = sum3; DelayPtr = DecMod(DelayPtr, IirSize); Delay[DelayPtr] = sum4; DelayPtr = DecMod(DelayPtr, IirSize);4.4 性能分析
经过这样的优化，性能提升是立竿见影的。对于一个N阶IIR滤波器：
- 单样本处理：每个输出样本需要大约N次乘加操作。
- 四样本并行处理：每4个样本需要执行N次乘加操作（因为基础内核每行做4个MAC，共N/4行，但每行处理4个样本的贡献，总MAC次数约为 N * 4 / 4？需要厘清）。更准确的分析来自指令周期：每个样本的指令周期数 = (4条指令/基础内核) * (N/4 个基础内核) / 4个样本 = N/4。这意味着计算速度提升到近4倍。
- 内存访问：每个样本的内存移动次数也减少到 N/2，降低了对内存带宽的压力。
关键心得：这个优化成功的核心在于将算法固有的串行反馈依赖，通过数学展开和巧妙的寄存器调度，转化为了可并行的计算模式。它牺牲了一些代码的直观性（循环结构变得复杂），换来了硬件利用率的极大提升。在实际编码时，务必用图表画出数据流和寄存器在每个时钟周期的状态，这是理解和解耦依赖关系的唯一方法。
5. 相关函数的多采样并行优化
相关函数是信号处理中另一个基础且计算密集的操作，用于衡量信号与自身时移版本之间的相似性。其公式为：R(n) = Σ (x(i) * x(i+n))，其中求和范围i从0到WindowSize-1。
与IIR滤波器不同，相关计算没有反馈，但它的挑战在于数据访问模式。计算不同滞后量n的相关值R(n),R(n+1)...时，需要重复访问输入数据序列的不同偏移位置。优化目标同样是利用多ALU同时计算多个滞后量（例如4个）的相关值。
5.1 数据流分析与内核构建
假设我们同时计算R(n),R(n+1),R(n+2),R(n+3)。观察计算公式可以发现，计算R(n)时需要x(0)*x(n),x(1)*x(n+1)...；计算R(n+1)时需要x(0)*x(n+1),x(1)*x(n+2)...。这里存在大量的数据复用机会。
我们定义两组指针：BasePtr指向基础序列x(i)，OffsetPtr指向偏移序列x(i+n)。理想的内核希望一次加载一个基础数据xb和四个偏移数据xd1,xd2,xd3,xd4（分别对应x(i+n),x(i+n+1),x(i+n+2),x(i+n+3)），然后并行计算四个乘积并累加到四个相关结果中：Cor1 += xb * xd1; Cor2 += xb * xd2; ...。
然而，和IIR优化类似，直接实现这个“通用内核”会遇到数据移动瓶颈：需要加载xb和xd1，同时还要将xd1,xd2,xd3在寄存器间移位（为下一次迭代准备）。这又超出了单指令的数据移动能力。
5.2 基础内核与旋转引用模式
解决方案与IIR滤波器的思路同源：构建一个更大的基础内核，并通过旋转引用避免寄存器拷贝。我们不再试图在一个小循环内维护完整的xd1-xd4滑动窗口，而是接受在基础内核的多次计算中，xb和xd的配对关系是变化的。
设计一个由4个“通用内核”组成的基础内核。在基础内核内部，我们让xd1-xd4的引用顺序每行旋转一次：
- 第1行：
Cor1 += xb * xd1;Cor2 += xb * xd2;Cor3 += xb * xd3;Cor4 += xb * xd4; 然后加载新的xb和xd1。 - 第2行：
Cor1 += xb * xd2;Cor2 += xb * xd3;Cor3 += xb * xd4;Cor4 += xb * xd1; 然后加载新的xb和xd2。 - 第3行：
Cor1 += xb * xd3; ...Cor4 += xb * xd2; - 第4行：
Cor1 += xb * xd4; ...Cor4 += xb * xd3;
通过这种旋转，xd1-xd4在每个样本位置都被充分利用，且每行只需要加载两个新操作数（一个xb，一个xd?），完美匹配硬件限制。基础内核循环WindowSize/4次，即可完成一个滞后量n的四个相关值的计算。
5.3 利用四操作数加载进一步优化
对于像StarCore SC140这样支持单指令加载四个连续内存字到四个寄存器的DSP，我们可以进行更激进的优化。我们可以一次性加载四个基础数据xb1:xb2:xb3:xb4和四个偏移数据xd1:xd2:xd3:xd4。
但这带来了新的挑战：数据的生命周期。xd4在被使用后，其“寿命”可能延续到后续多次计算中，如果立即加载下一组四个偏移数据，会覆盖掉还需要用的寄存器。解决方案是使用两组寄存器交替缓冲。
我们引入两组偏移数据寄存器：xd1-xd4和xd5-xd8。基础内核扩大为8行。在第一个4行循环中，使用xb1-xb4和xd1-xd4进行计算，并加载xd5-xd8。在第二个4行循环中，使用xb1-xb4和xd5-xd8进行计算，并加载下一组xd1-xd4。如此交替，确保了在需要使用某个数据时，它不会被覆盖。
这种优化将内存访问效率提升到了极致：每次循环迭代执行16次MAC操作（4个xb * 4个xd），但只进行了2次四操作数加载（共8个数据字）。平均下来，每次MAC操作仅需0.5次内存加载，极大地释放了数据总线带宽，让计算单元持续饱和工作。
实操要点：相关函数的优化，核心在于最大化数据复用，最小化数据移动。在实现时，要仔细绘制数据依赖图，明确每个数据的“生存期”。使用双缓冲甚至多缓冲寄存器组是解决数据生存期重叠问题的标准手法。在SIMD（单指令多数据）架构日益普及的今天，这种将多个标量计算重组为向量化、数据流清晰的内核的设计思想，具有极高的普适价值。
6. 从C仿真到汇编：StarCore SC140的实现精要
将优化后的算法映射到具体DSP硬件（如StarCore SC140）上，是最后也是最关键的一步。C代码描述了算法逻辑，而汇编代码则决定了硬件是否真的能以最高效率执行。
6.1 寄存器分配与指令并行
SC140拥有四个ALU（D0-D3常用于计算）和丰富的地址寄存器（R0-R3等）。我们的目标是在单条VLES（可变长指令集）指令中，让四个ALU都执行MAC操作，同时进行必要的数据加载。
以IIR滤波器为例，汇编内核的精妙之处在于将加载操作与计算操作紧密交织在并行指令中。例如：
[ mac d7,d8,d0 mac d6,d8,d1 mac d5,d8,d2 mac d4,d8,d3 move.f (r0)-,d7 move.f (r1)-,d8 ]这条指令在一个周期内完成了四件事：1）四个ALU分别用寄存器d7,d6,d5,d4与d8相乘，并累加到d0-d3；2）从r0指向的地址（系数数组）加载一个新系数到d7；3）从r1指向的地址（延迟线）加载一个新延迟值到d8；4）自动后递减指针r0和r1。这就是“加载与计算重叠”，完美隐藏了内存访问延迟。
6.2 模寻址与环形缓冲区管理
DSP处理流式数据时，延迟线通常用环形缓冲区实现。SC140的地址寄存器支持模寻址模式，通过设置模寄存器（M0），可以让指针在到达缓冲区边界时自动回绕，省去了边界检查的开销。在初始化时，设置M0 = IirSize*2 - 1（因为每个数据占2字节），并将指针指向缓冲区末尾，通过递减操作即可实现环形遍历。
move #Delay, r1 ; R1指向延迟线 move #(2*IirSize)-1, m0 ; 设置模寄存器，缓冲区大小为IirSize*2字节 ; 在循环中，使用 (r1)- 寻址，指针会自动在缓冲区循环6.3 循环控制与软件流水
SC140支持零开销硬件循环（DO循环）。我们需要设置两层循环：外层循环遍历数据块（DataBlockSize/4次），内层循环是IIR的基础内核（IirSize/4 - 1次）。使用dosetup和loopstart/loopend指令可以高效实现。
更高级的优化是软件流水，即手动安排指令，使得内层循环的加载、计算、存储操作在不同迭代间重叠，进一步填充指令流水线。这需要极其精细的指令调度，通常由编译器或经验丰富的工程师完成。在给出的示例代码中，交织的加载和计算指令已经体现了软件流水的思想。
6.4 定點數運算考量
在實際DSP中，大量使用定點數（Q格式）以提升速度、降低功耗。C仿真代码中的double类型在實作中需轉換為Word16（16位定點）或Word32（32位定點）。乘累加操作需使用飽和與舍入指令（如macr），並注意動態範圍管理，防止溢出。例如，在SC140彙編中，mac指令進行40位累加，最後用rnd指令進行舍入並存回32位寄存器。
踩坑實錄：從浮點C模型過渡到定點彙編時，最易忽略的是精度與溢出。務必先用定點C代碼驗證算法，仔細分析每個累加環節的位增長。在彙編中，要明確區分保護位（guard bits）的使用。另一個常見陷阱是指針對齊，SC140的四操作數加載要求數據地址必須是8字節對齊，否則會引發總線錯誤。在定義數據數組時，必須使用編譯器指令（如
align 8）確保對齊。
7. 性能評估與優化邊界
經過上述優化，我們可以量化性能收益：
- IIR濾波器：指令週期從每樣本 ~N 降至 ~N/4。記憶體訪問從每樣本 ~2N 次（讀係數、讀延遲）降至 ~N/2 次。
- 相關函數（使用四操作數加載）：指令週期降至每樣本 ~N/4，記憶體訪問降至每樣本 ~N/2，且由於加載合併，總線利用率更高。
然而，優化並非沒有代價：
- 代碼複雜度激增：並行化代碼比串行代碼難讀、難調試、難維護。
- 寄存器壓力巨大：為了保持數據在寄存器中復用，需要佔用大量寄存器，可能影響函數調用約定或需要額外的寄存器溢出/重載開銷。
- 對濾波器階數的約束：我們的優化假設濾波器階數
N是4的倍數。如果不是，需要進行邊界處理，會引入額外開銷。在設計系統時，應盡量將濾波器階數設為並行度的整數倍。 - 啟動與收尾開銷：並行處理需要“預熱”和“清理”階段（如我們代碼中循環前後的額外計算），對於處理極短數據塊，這些開銷可能抵消並行收益。
何時使用多採樣編程？
- 高階濾波器：N越大，並行化帶來的收益越能攤薄啟動開銷。
- 實時性要求極高的場景：如音頻編解碼、雷達信號處理。
- 功耗敏感場景：更快的計算意味著更早進入休眠狀態。
- 處理器具有顯著SIMD/多ALU能力：如現代DSP、ARM Neon、Intel AVX等。
進一步的優化維度：
- 記憶體佈局優化：將係數和延遲線數據組織為便於SIMD加載的結構（如交錯存儲）。
- 指令級並行（ILP）：在VLIW架構上，進一步挖掘指令間的獨立性，在一個超長指令字中安排更多無關聯的操作。
- 數據級並行（DLP）：如果硬件支持更寬的向量（如128位、256位），可以將4樣本並行擴展到8樣本甚至16樣本並行。
- 循環展開與軟件流水：對於小的內核循環，可以手動展開多次迭代，並重排指令以填充所有功能單元的流水線。
多採樣編程是將DSP算法性能推向極限的藝術。它要求開發者同時是算法專家、編譯器和硬件架構師。理解數據流、依賴關係和硬件限制是成功的基礎。從一個清晰正確的串行算法開始，逐步進行數學展開、依賴分析、內核設計、寄存器分配和指令調度，這個過程本身，就是對信號處理和計算機體系結構深刻理解的體現。當你看到優化後的代碼在示波器上跑出預期的實時性能時，那種成就感，是無與倫比的。