C - 使用 FFTW, fftw_plan_dft_1d 做卷积

Question

已解决：结果证明我的代码大部分是正确的，但我忘记了卷积会从搜索中翻转内核，所以我误用了该函数。查看我对当前代码的回答，并提供了一个类似的非 FFT 函数

似乎有一些人提出了一个问题，即如何实际进行一维卷积，特别是使用 FFTW。我正在尝试使用 FFTW3，但绝对没有成功。

其他人也提出了类似的问题，但将不相关的语言回复发布为“解决方案”。如果我解决了这个问题，我将发布一个实际的 C 解决方案！

请参阅：How to multiply 2 fftw_complex arrays 和 Calculating convolution of two functions using FFT (FFTW)

void Convolve( double * data, double * kernel, double * convout, int size )
{
    int i;

    size *= 2;  //Create zero-padded arrays.

    fftw_complex in_sequence[size], freq_sequence[size];
    fftw_complex in_data[size], freq_data[size];
    fftw_complex rev_data[size], time_data[size];
    fftw_plan p1 = fftw_plan_dft_1d(size, in_sequence, freq_sequence, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan p2 = fftw_plan_dft_1d(size, in_data, freq_data,         FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan rev = fftw_plan_dft_1d(size, rev_data, time_data,       FFTW_BACKWARD, FFTW_ESTIMATE);

    //Load our real data into the complex arrays.
    for( i = 0; i < size/2; i++ )
    {
        in_sequence[i][0] = kernel[i];
        in_sequence[i][1] = 0;

        in_data[i][0] = data[i];
        in_data[i][1] = 0;
    }
    for( ; i < size; i++ )
    {
        in_sequence[i][0] = in_sequence[i][1] = 0;
        in_data[i][0] = in_data[i][1] = 0;
    }

    fftw_execute(p1);
    fftw_execute(p2);

    for( i = 0; i < size; i++ )
    {
        double realD = freq_data[i][0];
        double imagD = freq_data[i][1];
        double realS = freq_sequence[i][0];
        double imagS = freq_sequence[i][1];
        rev_data[i][0] = (realD * realS - imagD * imagS)/size;
        rev_data[i][1] = (realD * imagS + imagD * realS)/size;
    }

    fftw_execute(rev);

    for( i = 0; i < size/2; i++ )
    {
        convout[i] = (time_data[i][0]*time_data[i][0]-time_data[i][1]*time_data[i][1]);
    }
}

此代码的输出有数字，但这些数字不对应于信号的卷积。

Answer 1

因此，如问题更新中所述，我的问题是您正在与之进行卷积的内核在执行卷积时被翻转。我已经更新了这里的代码和使用常规 C 的该函数的副本。我已经验证了这些函数会产生类似的输出。

void Convolve( double * data, double * kernel, double * convout, int size )
{
    int a, b;
    for( a = 0; a < size*2; a++ )
    {
        double cvo = 0;
        for( b = 0; b < size; b++ )
        {
            if( a-b<0 ) continue;
            cvo += data[b] * kernel[a-b];
        }
        convout[a] = cvo;
    }
}

//Convout's array size must be 2*size.
void FFTConvolve( double * data, double * kernel, double * convout, int size )
{
    struct timeval tvs, tve;
    int i;

    //Pad the 2nd half of the data with 0's.
    size *= 2;

    fftw_complex in_sequence[size], freq_sequence[size];
    fftw_complex in_data[size], freq_data[size];
    fftw_complex rev_data[size], time_data[size];
    fftw_plan p1 = fftw_plan_dft_1d(size, in_sequence, freq_sequence, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan p2 = fftw_plan_dft_1d(size, in_data, freq_data,         FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan rev = fftw_plan_dft_1d(size, rev_data, time_data,       FFTW_BACKWARD, FFTW_ESTIMATE);

    for( i = 0; i < size/2; i++ )
    {
        in_sequence[i][0] = kernel[i];
        in_sequence[i][1] = 0;

        in_data[i][0] = data[i];
        in_data[i][1] = 0;
    }
    for( ; i < size; i++ )
    {
        in_sequence[i][0] = in_sequence[i][1] = 0;
        in_data[i][0] = in_data[i][1] = 0;
    }

    fftw_execute(p1);
    fftw_execute(p2);

    double mux = 1.0/size;

    for( i = 0; i < size; i++ )
    {
        double realD = freq_data[i][0];
        double imagD = freq_data[i][1];
        double realS = freq_sequence[i][0];
        double imagS = freq_sequence[i][1];
        rev_data[i][0] = (realD * realS - imagD * imagS)*mux;
        rev_data[i][1] = (realD * imagS + imagD * realS)*mux;
    }

    fftw_execute(rev);

    for( i = 0; i < size; i++ )
    {
        convout[i] = time_data[i][0]; // The i term should be non-existent.
    }
}