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.
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];