DSPRelated.com
Forums

FFTW: a problem with in-place 2D r2c transform

Started by young March 13, 2006
Dear FFTW experts,

I get very strange results when I try to use in-place 2D r2c transforms.
Somehow, results are different for in-place and out-of-place transforms of
the same array, although an inverse transform always takes me back to my
original array (after scaling). For 1D there is no difference. Please,
tell me, what I'm doing wrong?!
My script is as follows (I know I allocate more space than is needed):

int main()
{
   fftw_complex *out;
   double *in;
   fftw_plan p;
   int i, j;
   int n[]={5, 5};

   in = (double*)malloc(sizeof(double)*(n[0]*n[1]+10));
   out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (n[0]*n[1]));


   p = fftw_plan_dft_r2c(2, n, in, (fftw_complex*)in,FFTW_MEASURE );
   //p = fftw_plan_dft_r2c(2, n, in, out,FFTW_MEASURE );

   in[0] = 12;
   in[1] = 14;
   in[2] = 11; 
   in[3] = 15;
   in[4] = 15 ;
   in[5] = 17;
   in[6] = 32;
   in[7] = 28;
   in[8] = 17;
   in[9] = 14;
   in[10] = 17;
   in[11] = 23;
   in[12] = 14;
   in[13] = 29;
   in[14] = 21;
   in[15] = 10;
   in[16] = 21; 
   in[17] = 21;
   in[18] = 10;
   in[19] = 16; 
   in[20] = 11;
   in[21] = 16 ;
   in[22] = 12;
   in[23] = 20;
   in[24] = 16;

   fftw_execute(p);

   printf("in[0]=%f, out[0][0]=%f\n", in[0], out[0][0]);
   free(in); fftw_free(out);
   fftw_destroy_plan(p); 
}



The array format for r2c transforms is described in detail at:


http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html

In particular, note that

1) the output is n[0] * (n[1]/2 + 1) in size, not n[0]*n[1]

2) the input is n[0] * n[1] in size for out-of-place and n[0] * 2 *
(n[1]/2+1) in size for in-place, where the latter reflects padding to
make room for the complex output.

This is all described in the manual.

Cordially,
Steven G. Johnson

The array format for r2c transforms is described in detail at:


http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html

In particular, note that

1) the output is n[0] * (n[1]/2 + 1) in size, not n[0]*n[1]

2) the input is n[0] * n[1] in size for out-of-place and n[0] * 2 *
(n[1]/2+1) in size for in-place, where the latter reflects padding to
make room for the complex output.

This is all described in the manual.

Cordially,
Steven G. Johnson

stevenj@alum.mit.edu wrote:
> The array format for r2c transforms is described in detail at: > > > http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html > > In particular, note that > > 1) the output is n[0] * (n[1]/2 + 1) in size, not n[0]*n[1] > > 2) the input is n[0] * n[1] in size for out-of-place and n[0] * 2 * > (n[1]/2+1) in size for in-place, where the latter reflects padding to > make room for the complex output. > > This is all described in the manual. > > Cordially,
Steven, You are a gem of patience and good will! Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Thanks for the kind words, Jerry.

Thanks a lot for your kind reply, Steven. I feel very flattered to be
answered by the author himself, and even more confused to continue asking
my dumb questions, but:
1) I read this part of the manual. I didn't realize that allocating more
than is needed was a mistake.
2) I've changed the array sizes to the exact values in your reply. I still
get zero-frequency of 351 for in-place, and 432 for out-of-place.
What am I doing wrong?? 

Thanks again

>The array format for r2c transforms is described in detail at: > > >http://www.fftw.org/doc/Multi_002dDimensional-DFTs-of-Real-Data.html > >In particular, note that > >1) the output is n[0] * (n[1]/2 + 1) in size, not n[0]*n[1] > >2) the input is n[0] * n[1] in size for out-of-place and n[0] * 2 * >(n[1]/2+1) in size for in-place, where the latter reflects padding to >make room for the complex output. > >This is all described in the manual. > >Cordially, >Steven G. Johnson > >
Probably you are writing or reading the input and/or output
incorrectly.

Allocating more space than necessary in itself is harmless.  However,
how you *access* the array elements depends on the array size.  If you
have an NxM array in row-major order, then the (i,j) element is
accessed as M*i+j.

Steven

Thanks again for your patience, Steven, and for your great program :)
I understand now what I had wrong - the padding should be at the end of
each line, not at the end of the whole array! So the same input array can
not be used for both in-place and out-of-place transforms, even if it's
big enough.
Sorry for bothering you with something so stupid.

- 

>Probably you are writing or reading the input and/or output >incorrectly. > >Allocating more space than necessary in itself is harmless. However, >how you *access* the array elements depends on the array size. If you >have an NxM array in row-major order, then the (i,j) element is >accessed as M*i+j. > >Steven > >
young wrote:
> Thanks again for your patience, Steven, and for your great program :) > I understand now what I had wrong - the padding should be at the end of > each line, not at the end of the whole array! So the same input array can > not be used for both in-place and out-of-place transforms, even if it's > big enough.
Actually, the same input array can be used for both in-place and out-of-place, but you have to use the advanced API - the advanced API allows you to specify an arbitrary amount of "padding" in your array.