c2dfftwcoefficients

2D R2C FFTW: Order of Fourier coefficients


I've used the 1D c2c transform several times, without any problems. The order of the Fourier coefficients for a transform with N grid points is: f_0, f_1, f_2, ..., f_N/2, f_-N/2+1, ...., f_-1.

I just can't figure out the order of the coefficients for the 2D R2C FFTW. I am using the following code. Using 2D_r2c, normalizing and then using 2D_c2r yields the original input so there should be no error.

void mFFTW2D(double INPUT[][STEPS], fftw_complex OUTPUT[][STEPS]){
    fftw_plan my_PLAN = fftw_plan_dft_r2c_2d(STEPS, 
                        STEPS,
                        *INPUT,
                        *OUTPUT,
                        FFTW_ESTIMATE);
    fftw_execute(my_PLAN);
    fftw_destroy_plan(my_PLAN);
}

void mIFFTW2D(fftw_complex INPUT[][STEPS], double OUTPUT[][STEPS]){
    fftw_plan my_PLAN = fftw_plan_dft_c2r_2d(STEPS,
                        STEPS,
                        *INPUT,
                        *OUTPUT,
                        FFTW_ESTIMATE);
    fftw_execute(my_PLAN);
    fftw_destroy_plan(my_PLAN);
    D2Norm(OUTPUT); //properly normalized: STEPS^-2
}

double INN[STEPS][STEPS];
fftw_complex OUTT[STEPS][STEPS];

// read in signal in INN
mFFTW2D(INN, OUTT);
// what is the order of the fourier coefficients in OUTT?
mIFFTW2D(OUTT, INN);

I used f(x,y)=sin(ax)*sin(ay) as a test input signal. 'a' was chosen in a manner, that the signal would be an integral multiple of one period of the sine(no leakage-effect). I was especially surprised that there was no symmetry of the Fourier coefficients for x and y.


Solution

  • The output of fftw_plan_dft_r2c_2d is not a STEP by STEP array of double. As the input is real, in the Fourier space, opposite frequencies are conjugate $X_{N-k} = X_k^*$. (http://en.wikipedia.org/wiki/Fast_Fourier_transform)

    fftw_c2r and fftw_r2c take account of that. Only half of the frequencies are stored and computations are reduced.

    http://www.fftw.org/fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format

    Therefore, you may rather use something like this (if it works) :

    mFFTW2D(double INPUT[][STEPS], fftw_complex OUTPUT[][STEPS/2+1])
    

    Watch the size of the array OUTT: you may reduce it as well and gain memory !

    fftw_complex OUTT[STEPS][STEPS/2+1];