C-使用FFTW、fftw_plan_dft_1d进行卷积

Question

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

似乎有几个人问了一个问题，如何实际做一个一维卷积，特别是与FFTW。我试着用FFTW3来做这件事，但完全没有成功。

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

见：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.
    }
}