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

Question

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

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是递归的：不能考虑“以前的值”来“打包”double数组，对吗？

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

欢迎提供任何提示或建议！

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

Answer 1

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

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

从上一个等式：

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寄存器。现在很容易将其转换为imintrin指令。请注意，现在您不能更改任何样本的输入值，但有时可以更改四个样本的倍数。

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

Answer 2

如果前面的n步骤有一个闭式公式，您可以使用它来回避串行依赖关系。如果可以使用与1步相同的操作使用不同的系数来计算它，那么您只需要广播。

就像在这个例子中，z1=cz1*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变量，而不是数组引用的伪代码，那么c c*b和b^2都是常量。（因此，除了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);
}

使用浮动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个步骤。

通常，人们对SIMD使用浮点运算，因为每个向量的元素数是原来的两倍，内存带宽/缓存占用率是原来的一半，因此通常比浮点运算快2倍。（每个时钟的向量数/字节数相同。）

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

（或如果使用一个FMA而不是mul编译-

避开串行依赖性不仅对SIMD有用，而且对于避免FP add/mul（或FMA）延迟瓶颈也很有用，可以进一步提高速度，达到FP add mul延迟与add mul吞吐量的比率。

只需展开更多，并使用多个向量，例如zv0、zv1、zv2、zv3。这也增加了您一次执行的步骤数。例如，浮动的16字节向量，具有4个向量，将是4x4=16步。