如何使用SSE指令集实现卷积算法并修复SSE版本输出结果不一致问题

现有基础定义 #

你已经定义的16字节对齐数组是正确的，无需修改：

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版本的问题分析 #

你写的SSE版本存在三个核心问题：

1. 计算逻辑完全错位：把i、j都转换为4元素块的索引后，你的计算逻辑等价于同时计算4条独立的卷积链路，和标量版本单通道卷积的计算逻辑完全不匹配
2. 寻址逻辑错误：__m128指针偏移1对应4个float的偏移，你直接套用标量的i-j判断逻辑，实际访问的输入输出位置和标量版本完全不一致
3. 边界判断逻辑错误：块级的边界判断无法覆盖单元素的边界场景，会导致大量有效计算被跳过或者越界访问

正确SSE卷积实现 #

我们采用「单j遍历，向量化4个连续输出」的逻辑实现，既保证和标量计算逻辑一致，又能充分利用SSE的并行计算能力：

#include <smmintrin.h>

void computeConvolutionOutputSSE(float* inputSignal, float* impulseResponse, float* outputSignal) {
    // 先把输出全部初始化为0
    __m128 zero = _mm_setzero_ps();
    int out_size_aligned = OUTPUT_SIGNAL_ARRAY_SIZE / 4 * 4;
    for (int i = 0; i < out_size_aligned; i += 4) {
        _mm_store_ps(outputSignal + i, zero);
    }
    // 处理末尾不足4个的输出元素
    for (int i = out_size_aligned; i < OUTPUT_SIGNAL_ARRAY_SIZE; i++) {
        outputSignal[i] = 0.0f;
    }

    // 遍历卷积核的每个元素j
    for (int j = 0; j < IMPULSE_RESPONSE_ARRAY_SIZE; j++) {
        // 把当前卷积核值广播为4个相同的元素，用于并行计算
        __m128 h_j = _mm_set1_ps(impulseResponse[j]);
        // 计算当前j对应的输出有效范围
        int out_start = j;
        int out_end = j + INPUT_SIGNAL_ARRAY_SIZE;
        // 对齐到4的边界，保证向量化访问的效率
        int out_end_aligned = out_start + (out_end - out_start) / 4 * 4;

        // 向量化处理4个连续的输出
        for (int i = out_start; i < out_end_aligned; i += 4) {
            // 加载对应位置的4个输入元素，非对齐访问用loadu
            __m128 x_vals = _mm_loadu_ps(inputSignal + i - j);
            // 加载当前的4个输出值
            __m128 out_vals = _mm_load_ps(outputSignal + i);
            // 计算 h[j] * x[i-j] + out[i]，并行处理4个元素
            out_vals = _mm_add_ps(out_vals, _mm_mul_ps(h_j, x_vals));
            // 存回输出
            _mm_store_ps(outputSignal + i, out_vals);
        }

        // 处理末尾不足4个的输出元素，用标量计算
        for (int i = out_end_aligned; i < out_end; i++) {
            outputSignal[i] += impulseResponse[j] * inputSignal[i - j];
        }
    }
}

实现说明 #

• 循环顺序调整为外层遍历卷积核j，内层并行处理4个连续输出，完全对齐标量版本的计算逻辑，结果完全一致
• 输出数组本身是16字节对齐的，所以用_mm_load_ps/_mm_store_ps做对齐访问，性能更高
• 输入元素因为i-j的偏移可能非对齐，所以用_mm_loadu_ps做非对齐访问，兼容性更好
• 增加了末尾不足4个元素的标量兜底处理，保证所有边界场景计算正确
• 去掉了内层的if判断，直接通过计算有效输出范围跳过无效计算，比每次判断性能更高

内容的提问来源于stack exchange，提问作者mkip