c++ - 如何用 SSE 实现卷积算法？

Question

const int INPUT_SIGNAL_ARRAY_SIZE = 256896;
const int IMPULSE_RESPONSE_ARRAY_SIZE = 318264;
const int OUTPUT_SIGNAL_ARRAY_SIZE = INPUT_SIGNAL_ARRAY_SIZE + IMPULSE_RESPONSE_ARRAY_SIZE;

__declspec(align(16)) float inputSignal_dArray[INPUT_SIGNAL_ARRAY_SIZE];
__declspec(align(16)) float impulseResponse_dArray[IMPULSE_RESPONSE_ARRAY_SIZE];
__declspec(align(16)) float outputSignal_dArray[OUTPUT_SIGNAL_ARRAY_SIZE];

我有正确工作并产生一些正确值的方法：

//#pragma optimize( "", off )
void computeConvolutionOutputCPU(float* inputSignal, float* impulseResponse, float* outputSignal) {
    float* pInputSignal = inputSignal;
    float* pImpulseResponse = impulseResponse;
    float* pOutputSignal = outputSignal;

    #pragma loop(no_vector)
    for (int i = 0; i < OUTPUT_SIGNAL_ARRAY_SIZE; i++)
    {
        *(pOutputSignal + i) = 0;

        #pragma loop(no_vector)
        for (int j = 0; j < IMPULSE_RESPONSE_ARRAY_SIZE; j++)
        {
            if (i - j >= 0 && i - j < INPUT_SIGNAL_ARRAY_SIZE)
            {
                *(pOutputSignal + i) = *(pOutputSignal + i)  + *(pImpulseResponse + j) *  (*(pInputSignal + i - j));
            }
        }
    }
}
//#pragma optimize( "", on )

另一方面，我应该将函数与 SSE 一起使用。我尝试了以下代码，但输出不一样：

void computeConvolutionOutputSSE(float* inputSignal, float* impulseResponse, float* outputSignal) {
    __m128* pInputSignal = (__m128*) inputSignal;
    __m128* pImpulseResponse = (__m128*) impulseResponse;
    __m128* pOutputSignal = (__m128*) outputSignal;


    int nOuterLoop = OUTPUT_SIGNAL_ARRAY_SIZE / 4;
    int nInnerLoop = IMPULSE_RESPONSE_ARRAY_SIZE / 4;
    int quarterOfInputSignal = INPUT_SIGNAL_ARRAY_SIZE / 4;

    __m128 m0 = _mm_set_ps1(0);

    for (int i = 0; i < nOuterLoop; i++)
    {
        *(pOutputSignal + i) = m0;
        for (int j = 0; j < nInnerLoop; j++)
        {
            if ((i - j) >= 0 && (i - j) < quarterOfInputSignal)
            {
                *(pOutputSignal + i) = _mm_add_ps(
                    *(pOutputSignal + i), 
                    _mm_mul_ps(*(pImpulseResponse + j), *(pInputSignal + i - j))
                );
            }
        }
    }
}

当我查看结果时，我看到我的 SSE 变体跳过了一些值等等，我认为我没有正确使用 SSE 实现。算法应该是这样的：

"*(pInputSignal + i - j) 在 SSE 的情况下是不正确的，因为它不是远离当前值的 ij 偏移量，它是 (ij) * 4 。我记得，事情是这样使用指针的想法除非内在函数从那时起发生了变化，否则这是不正确的——在我的时代，在这种情况下，必须将值“加载”到 __m128 的实例中，因为 H(J) 和 X(IJ) 处于未对齐的位置（并且序列中断）。” - @Swift-FridayPie。我认为这是必不可少的信息。

标签： c++simdsseconvolution