FFTW: successive real2complex 1-D FFTs fail

Hi everybody,

I've been experimenting lately with the great FFTW library in C, and I've
just encountered some very strange behavior.

I have a 1024x1024 image (w=1024, h=1024), on which I want to do the
following processing:

a) Perform 1-D FFTs on every column of the image array
b) Perform 1-D FFTs on every row of the same image array
c) more computations using the outputs of the two previous steps ...

The image is stored in RAM following a row-major order.

I've read the FFTW documentation and I came up with the following two
plans that make use of the advanced API:

fftw_p = fftw_plan_many_dft_r2c (1, &h, w, fftw_in, NULL, w, 1, fftw_out,
NULL, w, 1, FFTW_MEASURE);

fftw_p_T = fftw_plan_many_dft_r2c (1, &w, h, fftw_in, NULL, 1, fftw_out_T,
NULL, 1, h, FFTW_MEASURE);

So fftw_in is real (double), of size w*h and there are two output arrays,
fftw_out and fftw_out_T, both of size w*h/2+1.

I first create the plans, then I load the image into fftw_in, and then I
run the two plans, one after the other:

fftw_execute(fftw_p);
fftw_execute(fftw_p_T);

The strange thing is that the second transform produces wrong results. If
I comment out the first transform, then the second one works fine.

My original version used a complex2complex transform. That one did work
like a charm. So it must be something related to the real2complex version
of the plans.

I have tried the following:

1. switch to standard malloc instead of fftw_malloc for the input array
2. Use FFTW_ESTIMATE instead of FFTW_MEASURE
3. As ub 2m but also move the creation of the 2nd plan after the execution
of the 1st transform
4. As in 3, but also destroy the 1st plan before executing the 2nd one
5. As in 4, but also call fftw_cleanup() in between
6. I checked out the values of w and h before and after executing the 1st
transform, and they are always equal to 1024.

None of these steps helped me in any way. If I execute the 1st transform,
the 2nd one goes always wrong. If I dont, the 2nd one goes fine.

Furthermore, the whole idea of planning first, is that I want to reuse the
plans for real-time processing of a stream of images, so destroying and
cleaning up is not a realistic option anyway.

Does anybody have an idea what am I doing wrong? I'm very tempted to start
shouting that I've discovered a bug, but simple reasoning dictates that I'm
just a newbie in FFTW and I'm doing something incredibly stupid.

Just for the record, I'm working on an Intel Core 2 Duo T7200, running
32-bit Debian GNU/Linux, GCC version 4.3.3, and FFTW 3.1.2.

Thank you,

Dimitris


.

Relevant Pages

  • Re: FFTW: successive real2complex 1-D FFTs fail
    ... Perform 1-D FFTs on every row of the same image array ... I've read the FFTW documentation and I came up with the following two ... I first create the plans, then I load the image into fftw_in, and then I ...
    (comp.dsp)
  • fftw - loosing accuracy in cosine transform
    ... Sorry if this is a wrong place for posting fftw related question. ... transform as opposed to the fouerier transform? ... in cosine transform it is O. ... sample array contains input to the FFTW program ...
    (comp.lang.fortran)
  • Re: FFTW versus NRC FFT
    ... Before exchanging the FFT routine known from the NRC ... Applying both codes for that function the FFTW algorithm ... what the FFT NRC code does. ... In particular, if I remember correctly, the "sinft" transform that NRC ...
    (sci.math.num-analysis)
  • Re: FFTW Fortran Real2Complex R2C output array structure? fftw_plan_dft_r2c_3d?
    ... So, in the C FFTW ... for an r2c transform the last dimension is cut ...
    (comp.dsp)
  • Re: Confused by FFTW output
    ... > I am trying to get acquainted with the latest FFTW and wrote a small ... the FFT computes the discrete Fourier ... transform, which is *not the same* as the continuous Fourier ... (This periodicity is also why you ...
    (sci.math.num-analysis)