How to implement convolution algorithm with 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];

I have method that works correct and produces some correct values:

//#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 )

On the other hand I should use function with SSE. I tried the following code, but the output not the same:

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))
                );
            }
        }
    }
}

When I was viewing the results, I saw that my SSE variant skipped some values and so on, and I think, that I'm not using properly implementation with SSE. The algorithm should be this:

"*(pInputSignal + i - j) is incorrect in case of SSE, because it's not an i-j offset away from current value, it's (i-j) * 4 . THe thing is, as I remember it, the idea of using pointer that way is incorrect unless intrinsics had changed since then - in my time one had to "load" values into an instance of __m128 in this case, as H(J) and X(I-J) are in unaligned location (and sequence breaks)." - @Swift-FridayPie. I think that this is essential message .

How to implement convolution algorithm with SSE?

0 Answers