FIRSR多种方法实现，FFTW与KISSFFT性能对比

/*
* FIRSR
*/
void FIRSR(float* pSrc, float* pTaps, float* pDst, int numIters, int taps_len)
{
    int index;
    double y;
    int i, j, k;

    float* pWorkBuf = (float*)malloc(sizeof(float) * (taps_len - 1));

    index = taps_len - 2;
    for (k = 0; k < taps_len - 1; k++)
    {
        pWorkBuf[k] = 0;
    }

    for (i = 0; i < numIters; i++)
    {
        y = pTaps[0] * pSrc[i];
        for (j = index + 1; j < taps_len - 1; j++)
        {
            y += pTaps[j - index] * pWorkBuf[j];
        }

        for (k = 0; k <= index; k++)
        {
            y += pTaps[taps_len - 1 + k - index] * pWorkBuf[k];
        }

        pDst[i] = y;
        pWorkBuf[index] = pSrc[i];

        index = index == 0 ? taps_len - 2 : index - 1;
    }
    free(pWorkBuf);
}

#include "kiss_fft.h"
/*
* 使用kiss-FFT实现FIRSR
*/
void FIRSR_KISSFFT(float* pSrc, float* pTaps, float* pDst, int numIters, int taps_len)
{
    int i, j;

    // 对输入信号进行零填充
    kiss_fft_cpx* x_padded = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));
    memset(x_padded, 0, sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));

    // 对滤波器系数进行零填充
    kiss_fft_cpx* h_padded = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));
    memset(h_padded, 0, sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));

    // 转为复数类型
    for (int i = 0; i < numIters; i++)
    {
        x_padded[i].r = pSrc[i];
    }
    for (int i = 0; i < taps_len; i++)
    {
        h_padded[i].r = pTaps[i];
    }

    kiss_fft_cpx* X = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));
    kiss_fft_cpx* H = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));
    kiss_fft_cpx* Y = (kiss_fft_cpx*)malloc(sizeof(kiss_fft_cpx) * (numIters + taps_len - 1));
    // 正向FFT
    kiss_fft_cfg cfg = kiss_fft_alloc(numIters + taps_len - 1, 0, NULL, NULL);
    // 逆向FFT
    kiss_fft_cfg cfg_ifft = kiss_fft_alloc(numIters + taps_len - 1, 1, NULL, NULL);

    // 转换为频域
    kiss_fft(cfg, x_padded, X);
    kiss_fft(cfg, h_padded, H);

    // 对频域信号进行乘法运算
    for (i = 0; i < numIters + taps_len - 1; i++) {
        Y[i].r = X[i].r * H[i].r - X[i].i * H[i].i;
        Y[i].i = X[i].r * H[i].i + X[i].i * H[i].r;
    }

    for (i = 0; i < numIters + taps_len - 1; i++)
    {
        Y[i].r /= (numIters + taps_len - 1);
        Y[i].i /= (numIters + taps_len - 1);
    }

    // 使用IFFT得到滤波后的信号
    kiss_fft(cfg_ifft, Y, x_padded);

    for (int i = 0; i < numIters; i++)
    {
        pDst[i] = x_padded[i].r;
    }

    // 释放内存
    free(x_padded);
    free(h_padded);
    free(X);
    free(H);
    free(Y);
    kiss_fft_free(cfg);
    kiss_fft_free(cfg_ifft);
}

#include "fftw3.h"
/*
* 使用FFTW实现FIRSR
*/
void FIRSR_FFTW(float* pSrc, float* pTaps, float* pDst, int numIters, int taps_len)
{
    int i;
    int size = numIters + taps_len - 1;
    float* pSrc_padded, * pTaps_padded;
    fftwf_plan p1, p2, p3;

    // 填充输入信号与滤波器系数
    pSrc_padded  = (float*)calloc(size, sizeof(float));
    pTaps_padded = (float*)calloc(size, sizeof(float));
    memcpy(pSrc_padded, pSrc, numIters * sizeof(float));
    memcpy(pTaps_padded, pTaps, taps_len * sizeof(float));

    fftwf_complex* X = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * size);
    fftwf_complex* H = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * size);
    fftwf_complex* Y = (fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex) * size);

    // FFT-时域转频域
    p1 = fftwf_plan_dft_r2c_1d(size, pSrc_padded, X, FFTW_ESTIMATE);
    p2 = fftwf_plan_dft_r2c_1d(size, pTaps_padded, H, FFTW_ESTIMATE);
    fftwf_execute(p1);
    fftwf_execute(p2);

    // FIRSR-频域相乘并归一化
    for (i = 0; i < size; i++)
    {
        Y[i][0] = (X[i][0] * H[i][0] - X[i][1] * H[i][1]) / size;
        Y[i][1] = (X[i][0] * H[i][1] + X[i][1] * H[i][0]) / size;
    }

    // 频域转时域
    p3 = fftwf_plan_dft_c2r_1d(size, Y, pDst, FFTW_ESTIMATE);
    fftwf_execute(p3);

    fftwf_destroy_plan(p1);
    fftwf_destroy_plan(p2);
    fftwf_destroy_plan(p3);
    fftwf_free(X);
    fftwf_free(H);
    fftwf_free(Y);
}

测试函数：

void FIRSR_TEST()
{
    int i;
    // 生成随机的输入信号和滤波器系数
    int M = 4096;  // 输入信号的长度
    int N = 1024;   // 滤波器系数的长度
    //int M = 32;  // 输入信号的长度
    //int N = 8;   // 滤波器系数的长度
    const int loop = 20;
    float* x = (float*)malloc(sizeof(float) * M);
    float* h = (float*)malloc(sizeof(float) * N);

    // 调用FIR_filter函数进行滤波
    float* y = (float*)malloc(sizeof(float) * (M + N - 1));
    float* y_ = (float*)malloc(sizeof(float) * (M + N - 1));
    clock_t start_time, end_time;
    double t1, t2, t3;
    t1 = 0;
    t2 = 0;
    t3 = 0;

    // accuracy
    //for (i = 0; i < M; i++)
    //{
    //    x[i] = rand() / (float)RAND_MAX;
    //}
    //for (i = 0; i < N; i++)
    //{
    //    h[i] = rand() / (float)RAND_MAX;
    //}
    //FIRSR(x, h, y, M, N);
    //FIRSR_FFTW(x, h, y_, M, N);
    FIRSR_KISSFFT(x, h, y_, M, N);
    //for (i = 0; i < M; i++)
    //{
    //    cout << y[i] << "     " << y_[i] << endl;
    //}

    //performance
    for (int k = 0; k < loop; k++)
    {
        for (i = 0; i < M; i++)
        {
            x[i] = rand() / (float)RAND_MAX;
        }
        for (i = 0; i < N; i++)
        {
            h[i] = rand() / (float)RAND_MAX;
        }

        start_time = clock();
        FIRSR_KISSFFT(x, h, y, M, N);
        end_time = clock();
        t1 += double(end_time - start_time);
    }

    for (int k = 0; k < loop; k++)
    {
        for (i = 0; i < M; i++)
        {
            x[i] = rand() / (float)RAND_MAX;
        }
        for (i = 0; i < N; i++)
        {
            h[i] = rand() / (float)RAND_MAX;
        }

        start_time = clock();
        FIRSR(x, h, y_, M, N);
        end_time = clock();
        t2 += double(end_time - start_time);
    }

    for (int k = 0; k < loop; k++)
    {
        for (i = 0; i < M; i++)
        {
            x[i] = rand() / (float)RAND_MAX;
        }
        for (i = 0; i < N; i++)
        {
            h[i] = rand() / (float)RAND_MAX;
        }

        start_time = clock();
        FIRSR_FFTW(x, h, y_, M, N);
        end_time = clock();
        t3 += double(end_time - start_time);
    }

    cout << "FIRSR_KISSFFT TIME = " << t1 / CLOCKS_PER_SEC * 1000 << "ms" << endl;
    cout << "FIRSR         TIME = " << t2 / CLOCKS_PER_SEC * 1000 << "ms" << endl;
    cout << "FIRSR_FFTW    TIME = " << t3 / CLOCKS_PER_SEC * 1000 << "ms" << endl;
}

运行结果：

FFTW >> KISSFFT