C - 使用 FFTW,fftw_plan_dft_1d 做卷积

C - Using FFTW, fftw_plan_dft_1d to do convolution

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

似乎已经有一些人问过如何实际进行一维卷积的问题,特别是 FFTW。我正在尝试用 FFTW3 来做,但完全没有成功。

其他人也问过类似的问题,但 post 不相关的语言回答为 "solutions." 如果我解决了这个问题,我将 post 一个实际的 C 解决方案!

参见:How to multiply 2 fftw_complex arrays and 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]);
    }
}

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

因此,如问题更新中所述,我的问题是您正在使用的内核在执行卷积时被翻转了。我在这里更新了代码和使用常规 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.
    }
}