是否可以在计算中对串行依赖性使用SIMD，如指数移动平均滤波器

Question

我正在处理音频应用程序中不同参数的多个(独立的(指数移动平均单极滤波器，目的是以音频速率平滑每个参数值：

for (int i = 0; i < mParams.GetSize(); i++) {
mParams.Get(i)->SmoothBlock(blockSize);
}
...
inline void SmoothBlock(int blockSize) {
double inputA0 = mValue * a0;
for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
mSmoothedValues[sampleIndex] = z1 = inputA0 + z1 * b1;
}
}   

我想利用CPUSIMD指令，并行处理它们，但我真的不确定如何实现这一点。

事实上，z1是递归的：不能考虑"以前的值"来"打包"二重数组，对吧？

也许有没有一种方法可以正确地组织不同过滤器的数据并并行处理它们？

欢迎任何提示或建议！

请注意：我没有几个信号路径。任何参数表示对(唯一的(处理信号的不同控制。假设我有一个sin信号：参数1将影响增益，参数2间距，参数3滤波器截止，参数4平移，等等

Answer 1

如果前面有n步骤的闭合形式公式，则可以使用它来避开串行依赖关系。如果可以使用与1步相同的操作，仅使用不同的系数来计算，那么广播就是您所需要的。

就像在这种情况下，z1 = c + z1 * b，所以应用两次，我们得到

# I'm using Z0..n as the elements in the series your loop calculates
Z2 = c + (c+Z0*b)*b
= c + c*b + Z0*b^2

c + c*b和b^2都是常量，如果我正确理解您的代码，所有的C变量实际上只是C变量，而不是数组引用的伪代码。(所以除了z1之外的所有内容都是循环不变的(。

---

因此，如果我们有一个由2个元素组成的SIMD向量，从Z0和Z1开始，我们可以将它们中的每个向前推进2，得到Z2和Z3。

void SmoothBlock(int blockSize, double b, double c, double z_init) {
// z1 = inputA0 + z1 * b1;
__m128d zv = _mm_setr_pd(z_init, z_init*b + c);
__m128d step2_mul = _mm_set1_pd(b*b);
__m128d step2_add = _mm_set1_pd(c + c*b);
for (int i = 0; i < blockSize-1; i+=2) {
_mm_storeu_pd(mSmoothedValues + i, zv);
zv = _mm_mul_pd(zv, step2_mul);
zv = _mm_add_pd(zv, step2_add);
// compile with FMA + fast-math for the compiler to fold the mul/add into one FMA
}
// handle final odd element if necessary
if(blockSize%2 != 0)
_mm_store_sd(mSmoothedValues+blockSize-1, zv);
}

有了float+AVX(每个矢量8个元素(，就有了

__m256 zv = _mm256_setr_ps(z_init, c + z_init*b,
c + c*b + z_init*b*b,
c + c*b + c*b*b + z_init*b*b*b, ...);
// Z2 = c + c*b + Z0*b^2
// Z3 = c + c*b + (c + Z0*b) * b^2
// Z3 = c + c*b + c*b^2 + Z0*b^3

并且加法/复数因子将用于8个步骤。

通常，人们将float用于SIMD，因为每个向量的元素数量是原来的两倍，内存带宽/缓存占用面积是原来的一半，所以通常比float加速2倍。(每个时钟的矢量/字节数相同。(

---

例如，Haswell或Sandybridge CPU上的上述循环将以每8个周期一个矢量的速度运行，受到mulpd(5个周期(+addps(3个周期(延迟的限制我们为每个向量生成2个double结果，但与每个时钟1个mul和1个add的吞吐量相比，这仍然是一个巨大的瓶颈。我们错过了8倍的吞吐量。

(或者，如果用一个FMA而不是mul->add编译，那么我们有5个周期的延迟(。

并行串行依赖性不仅对SIMD有用，而且对于避免FP add/mul(或FMA(延迟上的瓶颈也很有用，这将带来进一步的加速，最高可达FP add/mol延迟与add+mul吞吐量的比率

只需展开更多，并使用多个矢量，如zv0、zv1、zv2、zv3。这也增加了您同时执行的步骤数。因此，例如，float的16字节矢量，具有4个矢量，将是4x4=16步。

Answer 2

您有一个特殊情况，其中输入信号是Heaviside阶跃函数。您希望获得对此函数的筛选器响应，称为Step响应。在这种情况下可以消除递归。首先，让我们展开递归几个步骤。

z[1] = in + z[0]*b
z[2] = in + z[1]*b = in + (in + z[0]*b)*b  = in*(1 + b) + z[0]*b^2
z[3] = in + z[2]*b = in*(1 + b + b^2) + z[0]*b^3
z[4] = in + z[3]*b = in*(1 + b + b^2 + b^3) + z[0]*b^4

来自上一个eq:

z[1] = in*(1 + b + b^2 + b^3) + z[-3]*b^4
z[2] = in*(1 + b + b^2 + b^3) + z[-2]*b^4
z[3] = in*(1 + b + b^2 + b^3) + z[-1]*b^4
z[4] = in*(1 + b + b^2 + b^3) + z[0]*b^4

现在很容易以矢量化的形式重写它

in' = {in, in, in, in};
z' = in' * (1 + b + b^2 + b^3) + z'*b^4

其中"表示矢量或单个SIMD寄存器。现在，将其翻译成immontrin指令将很容易。请注意，现在您不能更改任何样本的输入值，但有时可以更改四个样本的倍数。

此外，您可以将两个或多个SIMD寄存器表示为一个向量，并进一步扩展递归。这将提高性能，因为管道利用率更高，但不要过度使用，否则您将没有足够的寄存器。