基于x86 Intrinsic的FIR滤波器内外循环优化问询

原代码 #

void Fir(float* pIn, float* pOut, float* pCoeff, float* pStage, uint32_t N, uint32_t FilterLength)
{
    int n, k;
    float* pSrc;
    float* pCoeffSrc = pCoeff;
    float* pDst = pOut;
    float s0, s1, s2, s3;
    __m128 Vec, Mul;
    __m128 Sum0,Sum1,Sum2,Sum3;

    __m128 Zero = _mm_set_ps1(0);
    memcpy(&pStage[FilterLength - 1], pIn, N * sizeof(float));

    for (n = 0; n < N; n+=4)
    {
        //Sum0
        pSrc = &pStage[n];
        Sum0 = _mm_set_ps1(0);
        pCoeffSrc = pCoeff;

        for (k = 0; k < FilterLength >> 2; k++)
        {
            __m128 Coeff = _mm_load_ps(pCoeffSrc);
            Vec = _mm_load_ps(pSrc); 
            Sum0  = _mm_fmadd_ps(Coeff, Vec, Sum0);
            pCoeffSrc += 4;
            pSrc += 4;
        }

        Sum0 = _mm_hadd_ps(Sum0, Zero);
        Sum0 = _mm_hadd_ps(Sum0, Zero);

        //Sum1
        pSrc = &pStage[n+1];
        Sum1 = _mm_set_ps1(0);
        pCoeffSrc = pCoeff;

        for (k = 0; k < FilterLength >> 2; k++)
        {
            __m128 Coeff = _mm_load_ps(pCoeffSrc);
            Vec = _mm_load_ps(pSrc);
            Sum1 = _mm_fmadd_ps(Coeff, Vec, Sum1);
            pCoeffSrc += 4;
            pSrc += 4;
        }

        Sum1 = _mm_hadd_ps(Sum1, Zero);
        Sum1 = _mm_hadd_ps(Sum1, Zero);

        //Sum2
        pSrc = &pStage[n+2];
        Sum2 = _mm_set_ps1(0);
        pCoeffSrc = pCoeff;

        for (k = 0; k < FilterLength >> 2; k++)
        {
            __m128 Coeff = _mm_load_ps(pCoeffSrc);
            Vec = _mm_load_ps(pSrc);
            Sum2 = _mm_fmadd_ps(Coeff, Vec, Sum2);
            pCoeffSrc += 4;
            pSrc += 4;
        }

        Sum2 = _mm_hadd_ps(Sum2, Zero);
        Sum2 = _mm_hadd_ps(Sum2, Zero);

        //Sum3
        pSrc = &pStage[n+3];
        Sum3 = _mm_set_ps1(0);
        pCoeffSrc = pCoeff;

        for (k = 0; k < FilterLength >> 2; k++)
        {
            __m128 Coeff = _mm_load_ps(pCoeffSrc);
            Vec = _mm_load_ps(pSrc);
            Sum3 = _mm_fmadd_ps(Coeff, Vec, Sum3);
            pCoeffSrc += 4;
            pSrc += 4;
        }

        Sum3 = _mm_hadd_ps(Sum3, Zero);
        Sum3 = _mm_hadd_ps(Sum3, Zero);

        Vec = _mm_set_ps(Sum3.m128_f32[0], Sum2.m128_f32[0], Sum1.m128_f32[0], Sum0.m128_f32[0]);
        _mm_store_ps(pDst, Vec);
        pDst+=4;
    }
}

原代码存在重复内循环、冗余系数加载、低效标量重组等问题，完全可以针对性优化，具体方向如下：

1. 合并内循环，复用系数加载 #

原代码对4个输出样本分别执行独立内循环，每次都重复加载滤波器系数，浪费大量内存带宽和指令周期。可以将4个Sum的计算合并到同一个内循环中，只加载一次系数，同时处理4组输入数据：

// 优化后的核心内循环逻辑
for (n = 0; n < N; n += 4)
{
    Sum0 = _mm_set_ps1(0.0f);
    Sum1 = _mm_set_ps1(0.0f);
    Sum2 = _mm_set_ps1(0.0f);
    Sum3 = _mm_set_ps1(0.0f);
    pCoeffSrc = pCoeff;

    // 合并四个内循环，仅加载一次系数
    for (k = 0; k < FilterLength >> 2; k++)
    {
        __m128 Coeff = _mm_load_ps(pCoeffSrc);
        uint32_t stageOffset = 4 * k;
        
        // 一次性加载4组相邻输入向量，对应n、n+1、n+2、n+3的偏移
        __m128 Vec0 = _mm_load_ps(&pStage[n + stageOffset]);
        __m128 Vec1 = _mm_load_ps(&pStage[n + 1 + stageOffset]);
        __m128 Vec2 = _mm_load_ps(&pStage[n + 2 + stageOffset]);
        __m128 Vec3 = _mm_load_ps(&pStage[n + 3 + stageOffset]);
        
        Sum0 = _mm_fmadd_ps(Coeff, Vec0, Sum0);
        Sum1 = _mm_fmadd_ps(Coeff, Vec1, Sum1);
        Sum2 = _mm_fmadd_ps(Coeff, Vec2, Sum2);
        Sum3 = _mm_fmadd_ps(Coeff, Vec3, Sum3);
        
        pCoeffSrc += 4;
    }

    // 后续求和与重组操作...
}

2. 用向量指令替代标量提取与重组 #

原代码通过SumX.m128_f32[0]提取标量再用_mm_set_ps重组，会触发不必要的标量-向量转换，效率极低。可以用_mm_hadd_ps的组合直接在向量域完成结果合并：

// 优化后的求和与合并逻辑
// 第一次水平加法，得到两个半和
Sum0 = _mm_hadd_ps(Sum0, Sum1);
Sum2 = _mm_hadd_ps(Sum2, Sum3);
// 第二次水平加法，得到四个结果的标量和，分别存在向量的四个分量中
Sum0 = _mm_hadd_ps(Sum0, Sum2);
// 直接存储合并后的向量结果
_mm_store_ps(pDst, Sum0);
pDst += 4;

3. 额外优化点 #

• 内存对齐：确保pIn、pOut、pCoeff、pStage都是16字节对齐（可使用_mm_malloc或编译器对齐属性），避免非对齐内存访问的性能损失。
• 循环展开：若FilterLength为固定值，可手动展开内循环，进一步减少循环控制开销。
• memcpy替换：若N较大，可用批量向量加载/存储替代memcpy，提升数据搬移效率。

完整优化代码示例 #

#include <emmintrin.h>
#include <smmintrin.h> // 针对_mm_fmadd_ps

void OptimizedFir(float* pIn, float* pOut, float* pCoeff, float* pStage, uint32_t N, uint32_t FilterLength)
{
    uint32_t n, k;
    float* pDst = pOut;
    __m128 Sum0, Sum1, Sum2, Sum3;
    __m128 Zero = _mm_set_ps1(0.0f);

    // 用向量指令替代memcpy（可选，视内存对齐情况调整）
    for (uint32_t i = 0; i < N; i += 4)
    {
        __m128 inVec = _mm_load_ps(&pIn[i]);
        _mm_store_ps(&pStage[FilterLength - 1 + i], inVec);
    }

    for (n = 0; n < N; n += 4)
    {
        Sum0 = Zero;
        Sum1 = Zero;
        Sum2 = Zero;
        Sum3 = Zero;
        float* pCoeffSrc = pCoeff;

        for (k = 0; k < FilterLength >> 2; k++)
        {
            __m128 Coeff = _mm_load_ps(pCoeffSrc);
            uint32_t stageOffset = 4 * k;
            
            __m128 Vec0 = _mm_load_ps(&pStage[n + stageOffset]);
            __m128 Vec1 = _mm_load_ps(&pStage[n + 1 + stageOffset]);
            __m128 Vec2 = _mm_load_ps(&pStage[n + 2 + stageOffset]);
            __m128 Vec3 = _mm_load_ps(&pStage[n + 3 + stageOffset]);
            
            Sum0 = _mm_fmadd_ps(Coeff, Vec0, Sum0);
            Sum1 = _mm_fmadd_ps(Coeff, Vec1, Sum1);
            Sum2 = _mm_fmadd_ps(Coeff, Vec2, Sum2);
            Sum3 = _mm_fmadd_ps(Coeff, Vec3, Sum3);
            
            pCoeffSrc += 4;
        }

        // 合并求和结果
        Sum0 = _mm_hadd_ps(Sum0, Sum1);
        Sum2 = _mm_hadd_ps(Sum2, Sum3);
        Sum0 = _mm_hadd_ps(Sum0, Sum2);
        
        _mm_store_ps(pDst, Sum0);
        pDst += 4;
    }
}

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