1
votes

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.

2
Sorry, I am facing a similar problem, but I don't understand your code at all. Where are you creating the Fourier coefficients?Tropilio

2 Answers

0
votes

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

Bye,

0
votes

The OUT array gets modified in when it is used in the inverse transform routine, so you have to check it before doing that. The wavenumbers are as the documentation then.