使用 NEON 指令加速级联双二阶 - 它是如何工作的？

Posted 2023-02-16

【中文标题】使用 NEON 指令加速级联双二阶 - 它是如何工作的？

【英文标题】：Using NEON instructions to speed up cascaded biquads - how it works?

【发布时间】：2020-07-16 17:02:07

【问题描述】：

我试图了解如何使用 Neon 扩展为 CMSIS 中的 Arm 处理器优化级联双二阶滤波。 代码在#if defined(ARM_MATH_NEON)here下ifdefed，文档为here。

NEON intrinsics 在级联超过 4 个双二阶时使用。如果一个biduaq的输出作为输入馈送到下一个biduaq，我很困惑如何执行任何类型的并行指令？谁能解释一下在代码的和平中并行完成了什么？

【问题讨论】：

SO 问题需要独立/抗链接腐烂；除了链接到源代码之外，还要引用您要询问的部分。

我建议使用编译器的向量扩展。这将在必要时自动使用霓虹灯并轻松移植到其他架构。

【参考方案1】：

双二阶级联可以通过及时偏移来实现并行化。

如果您一次计算 4 个 biquad，则最后一个级联 biquad 不会对同一批次 4 中前一个 biquad 的结果进行运算，而是对前一批 4 中保存的结果进行运算。这消除了内部的依赖关系每批。因此，从第一个双二阶对角线传播数据到最后一个双二阶需要 4 步延迟，但吞吐量在每个时间步完成 4 个双二阶，或者说吞吐量比一次计算一个双二阶高 4 倍。

【参考方案2】：

这是文档中的公式：

y[ n ] = b0 * x[ n ] + d1;
d1 = b1 * x[ n ] + a1 * y[ n ] + d2;
d2 = b2 * x[ n ] + a2 * y[ n ];

让我们通过重命名变量来摆脱可变状态，循环进行 2 次迭代：

// Iteration 1
y[ n ] = b0 * x[ n ] + d1_0;
const float d1_1 = b1 * x[ n ] + a1 * y[ n ] + d2_0;
const float d2_1 = b2 * x[ n ] + a2 * y[ n ];

// Iteration 2
y[ n + 1 ] = b0 * x[ n + 1 ] + d1_1;
const float d1_2 = b1 * x[ n + 1 ] + a1 * y[ n + 1 ] + d2_1;
const float d2_2 = b2 * x[ n + 1 ] + a2 * y[ n + 1 ];

当它这样写时，很明显你可以替换变量，并并行计算 2 次迭代，方法如下：

// Rewriting iterations to only use data available before the #1
y[ n ] = b0 * x[ n ] + d1_0;
y[ n + 1 ] = b0 * x[ n + 1 ] + b1 * x[ n ] + a1 * b0 * x[ n ] + d1_0 + d2_0;
const float d1_2 = b1 * x[ n + 1 ] + a1 * y[ n + 1 ] + b2 * x[ n ] + a2 * y[ n ];
const float d2_2 = b2 * x[ n + 1 ] + a2 * y[ n + 1 ];

很确定我把上面的代数搞砸了，但我希望你明白了。该方法以更多计算为代价消除了数据依赖性。

该特定实现通过移动向量和进行大量额外计算来完成 4 次迭代而不是 2 次。这是带有 HLSL 样式的 cmets 的主要 NEON 循环，关于 YnV SIMD 向量的通道发生了什么。

float32x4_t YnV = s;
// YnV.w += t1.w * dV.val[ 0 ].x;
s = vextq_f32( zeroV, dV.val[ 0 ], 3 );
YnV = vmlaq_f32( YnV, t1, s );

// YnV.zw += t2.zw * dV.val[ 0 ].xy;
s = vextq_f32( zeroV, dV.val[ 0 ], 2 );
YnV = vmlaq_f32( YnV, t2, s );

// YnV.yzw += t3.yzw * dV.val[ 0 ].xyz
s = vextq_f32( zeroV, dV.val[ 0 ], 1 );
YnV = vmlaq_f32( YnV, t3, s );

// And finally the all-lanes version without shifts:
// YnV.xyzw += t4.xyzw * XnV.xyzw
YnV = vmlaq_f32( YnV, t4, XnV );

【讨论】：

好的，我明白了，蛮力代数。但是...您在编写时考虑了单个双二阶（是吗？），并且我提到的代码不使用任何 NEON，除非有 4 个双二阶级联。你认为为什么会这样？如果只有单个双二阶，CMSIS 库仅使用循环展开（检查上面的链接），仅此而已，没有 NEON。