Weird results from FFTW -- banding and large roundoff

I have been struggling to figure out some strange results from FFTW. I get
a angled banding in the DFT and the inverse gives an image that looks
generally OK except that the signal varies periodically, consistent with
those bands, by up to 1% above and below the original image signal.

I have put some picture and description here: 
http://www.mirametrics.com/pub/fftw_strange_results.htm

I have been yanking out my hair over this. Anyone have any ideas what is
wrong? THANKS!

Micahel

Michael48 wrote:
> I have been struggling to figure out some strange results from FFTW. I get
> a angled banding in the DFT and the inverse gives an image that looks
> generally OK except that the signal varies periodically, consistent with
> those bands, by up to 1% above and below the original image signal.
> 
> I have put some picture and description here: 
> http://www.mirametrics.com/pub/fftw_strange_results.htm

Looks as if you got your memory management/addressing wrong.
My first idea is that you gave FFTW the dimensions in the wrong
order. FFTW wants the dimension with consecutive memory as
the last dimension. So if your image is 1042x1032 and the
1042 pixels in a line are in consecutive memory locations,
you've got to plan a 1032x1042 FFT.

Maybe you could show some source?


Hendrik vdH

"Hendrik van der Heijden" <hvdh@gmx.de> wrote in message 
news:45da252d$0$6392$9b4e6d93@newsspool2.arcor-online.net...
> Michael48 wrote:
>> I have been struggling to figure out some strange results from FFTW. I 
>> get
>> a angled banding in the DFT and the inverse gives an image that looks
>> generally OK except that the signal varies periodically, consistent with
>> those bands, by up to 1% above and below the original image signal.
>>
>> I have put some picture and description here: 
>> http://www.mirametrics.com/pub/fftw_strange_results.htm
>
> Looks as if you got your memory management/addressing wrong.
> My first idea is that you gave FFTW the dimensions in the wrong
> order. FFTW wants the dimension with consecutive memory as
> the last dimension. So if your image is 1042x1032 and the
> 1042 pixels in a line are in consecutive memory locations,
> you've got to plan a 1032x1042 FFT.
>
> Maybe you could show some source?
>
>
> Hendrik vdH
>

It might also help to see the original image.  Starting from the DFTs is 
always questionable.

Fred 

That fixed the major problem but there is one remaining wrinkle to work
out. Thank you!!!

There's clearly omething I am missing. The docs say, for every example,
that the arguments are passed as nx, ny, nz. So, I guess that nx does not
mean "number of columns" and ny does not mean "number of rows" ? I am
confused: Should pass the args in the order nx, ny, nz but transpose the
data buffer before r2c and after c2r? Or should I reverse the order of
parameter passing and leave the data array intact? The 2D image is stored
as a single buffer with the column number (nx) varying most rapidly. And I
am calling the C libraries from C code. Something is wrong here!!

For the in-place case, I get great precision:

   (forward >> inverse DFT) - orignal < ~10^-18.

Here's the remaining wrinkle: If I use an auxiliary buffer rather than
in-place, then there's problem in how the data are packed and unpacked,
from the point of view of the column count. Maybe this is a transpose/arg
ordering issue like I mentioned above. I need to investigate this some
more.

Thank you again, for getting me over this first hurdle.


Michael

>Michael48 wrote:
>> I have been struggling to figure out some strange results from FFTW. I
get
>> a angled banding in the DFT and the inverse gives an image that looks
>> generally OK except that the signal varies periodically, consistent
with
>> those bands, by up to 1% above and below the original image signal.
>> 
>> I have put some picture and description here: 
>> http://www.mirametrics.com/pub/fftw_strange_results.htm
>
>Looks as if you got your memory management/addressing wrong.
>My first idea is that you gave FFTW the dimensions in the wrong
>order. FFTW wants the dimension with consecutive memory as
>the last dimension. So if your image is 1042x1032 and the
>1042 pixels in a line are in consecutive memory locations,
>you've got to plan a 1032x1042 FFT.
>
>Maybe you could show some source?
>
>
>Hendrik vdH
>
>

Michael48 wrote:
> There's clearly omething I am missing. The docs say, for every example,
> that the arguments are passed as nx, ny, nz. So, I guess that nx does not
> mean "number of columns" and ny does not mean "number of rows" ? I am
> confused: Should pass the args in the order nx, ny, nz but transpose the
> data buffer before r2c and after c2r?  Or should I reverse the order of
> parameter passing and leave the data array intact?

Transposing a 1MP image takes some time.
Passing the args in a different order doesn't cost anything.

> For the in-place case, I get great precision:

> Here's the remaining wrinkle: If I use an auxiliary buffer rather than
> in-place, then there's problem in how the data are packed and unpacked,
> from the point of view of the column count.

If you do a R2C FFT out-of-place, the memory strides (# of bytes
between two line starts) of input and output image are different,
see

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

Another thing is that in out-of-place transforms, the input array
data will generally not be preserved unless you don't specify the
FFTW_PRESERVE_INPUT planner flag.
Maybe you just compared to you overwritten input array.


Hendrik vdH

Michael48 wrote:
 > There's clearly omething I am missing. The docs say, for every example,
 > that the arguments are passed as nx, ny, nz. So, I guess that nx does not
 > mean "number of columns" and ny does not mean "number of rows" ? I am
 > confused: Should pass the args in the order nx, ny, nz but transpose the
 > data buffer before r2c and after c2r?  Or should I reverse the order of
 > parameter passing and leave the data array intact?

Transposing a 1MP image takes some time.
Passing the args in a different order doesn't cost anything.

 > For the in-place case, I get great precision:

 > Here's the remaining wrinkle: If I use an auxiliary buffer rather than
 > in-place, then there's problem in how the data are packed and unpacked,
 > from the point of view of the column count.

If you do a R2C FFT out-of-place, the memory strides (# of bytes
between two line starts) of input and output image are different,
see

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


Hendrik vdH

Hendrik,

Thanks for your help. Now I get the forward/inverse pair for the 1042x1032
array in about 0.6 seconds and the result is clean to about 1 part in 10^14
or better. It is a significant improvement over my old FFT code.

I do appreciate the difference between the order of parameter passing and
whether or not to transpose twice. I know you can't speak for the FFTW
authors, but why on earth would the documentation specify "nx,ny, nz"
*everywhere*, when it really should always be "nz, ny, nx" ?? That is just
so arcane! Which is why it cost me most of 2 days of time and frustration
before I realized there was something I was simply not getting about using
this library. I am grateful to you for pointing me the right way!

Michael

>Michael48 wrote:
> > There's clearly omething I am missing. The docs say, for every
example,
> > that the arguments are passed as nx, ny, nz. So, I guess that nx does
not
> > mean "number of columns" and ny does not mean "number of rows" ? I am
> > confused: Should pass the args in the order nx, ny, nz but transpose
the
> > data buffer before r2c and after c2r?  Or should I reverse the order
of
> > parameter passing and leave the data array intact?
>
>Transposing a 1MP image takes some time.
>Passing the args in a different order doesn't cost anything.
>
> > For the in-place case, I get great precision:
>
> > Here's the remaining wrinkle: If I use an auxiliary buffer rather
than
> > in-place, then there's problem in how the data are packed and
unpacked,
> > from the point of view of the column count.
>
>If you do a R2C FFT out-of-place, the memory strides (# of bytes
>between two line starts) of input and output image are different,
>see
>
>http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html
>
>
>Hendrik vdH
>

On Feb 26, 10:10 am, "Michael48" <m...@mirametrics.com> wrote:
> I do appreciate the difference between the order of parameter passing and
> whether or not to transpose twice. I know you can't speak for the FFTW
> authors, but why on earth would the documentation specify "nx,ny, nz"
> *everywhere*, when it really should always be "nz, ny, nx" ?? That is just
> so arcane! Which is why it cost me most of 2 days of time and frustration
> before I realized there was something I was simply not getting about using
> this library. I am grateful to you for pointing me the right way!

I am one of the FFTW authors, and I can point out that FFTW's
dimension ordering conforms to the standard for C programs and for
most other modern programming languages other than Fortran.

In a C program, if you have an array[nx][ny], it is stored with the ny
dimension contiguous, and this is specified as an "nx, ny" array in
FFTW.  This is called "row-major order" for an nx by ny array.

The other convention is "column-major" order, in which an nx by ny
array is stored with the nx dimension contiguous.  This is used in
Fortran for an array(nx,ny), and in various other programs derived
from Fortran (such as Matlab).  A column-major array can be specified
in FFTW merely by passing the dimensions in reverse order.

We can't help it that there are two common conventions for laying out
multidimensional arrays in memory; we have to pick one ordering, and
it is a C library so we pick the C convention.  I will note, however,
that there is an *entire section* of the FFTW manual discussing
multidimensional array formats in great detail.  I'm not sure what
more you can expect.

Regards,
Steven G. Johnson

Steve,

I really do appreciate your work of FFTW, and I understand how array
ordering works in Fortran vs C. But that wasn't my point. My point was
that your doc explicltly states the parameter passing order as "foo( nx,
ny, nz, ..."  So my mistake was to associate (nx) with the number of
columns and (ny) with number of rows. If you are looking at a 2-D image on
a screen, then don't you think that is somewhat hueristic? If I have a C
program, with its multi-D arrays in a native row-major order, and my code
calls your C-library and the docs show the parameter order as nx, ny, nz,
then can you see that there's room for confusion at the implementation
end? 

Michael

>On Feb 26, 10:10 am, "Michael48" <m...@mirametrics.com> wrote:
>> I do appreciate the difference between the order of parameter passing
and
>> whether or not to transpose twice. I know you can't speak for the FFTW
>> authors, but why on earth would the documentation specify "nx,ny, nz"
>> *everywhere*, when it really should always be "nz, ny, nx" ?? That is
just
>> so arcane! Which is why it cost me most of 2 days of time and
frustration
>> before I realized there was something I was simply not getting about
using
>> this library. I am grateful to you for pointing me the right way!
>
>I am one of the FFTW authors, and I can point out that FFTW's
>dimension ordering conforms to the standard for C programs and for
>most other modern programming languages other than Fortran.
>
>In a C program, if you have an array[nx][ny], it is stored with the ny
>dimension contiguous, and this is specified as an "nx, ny" array in
>FFTW.  This is called "row-major order" for an nx by ny array.
>
>The other convention is "column-major" order, in which an nx by ny
>array is stored with the nx dimension contiguous.  This is used in
>Fortran for an array(nx,ny), and in various other programs derived
>from Fortran (such as Matlab).  A column-major array can be specified
>in FFTW merely by passing the dimensions in reverse order.
>
>We can't help it that there are two common conventions for laying out
>multidimensional arrays in memory; we have to pick one ordering, and
>it is a C library so we pick the C convention.  I will note, however,
>that there is an *entire section* of the FFTW manual discussing
>multidimensional array formats in great detail.  I'm not sure what
>more you can expect.
>
>Regards,
>Steven G. Johnson
>
>

On Feb 26, 12:37 pm, "Michael48" <m...@mirametrics.com> wrote:
> I really do appreciate your work of FFTW, and I understand how array
> ordering works in Fortran vs C. But that wasn't my point. My point was
> that your doc explicltly states the parameter passing order as "foo( nx,
> ny, nz, ..."  So my mistake was to associate (nx) with the number of
> columns and (ny) with number of rows. If you are looking at a 2-D image on
> a screen, then don't you think that is somewhat hueristic? If I have a C
> program, with its multi-D arrays in a native row-major order, and my code
> calls your C-library and the docs show the parameter order as nx, ny, nz,
> then can you see that there's room for confusion at the implementation
> end?

There is room for confusion whether we list them as (nx,ny,nz) or
(nz,ny,nx): that's why you have to read the manual instead of
guessing.

In English, it is conventional to name the dimensions "x", "y", "z"
listed from left to right.  It is also conventional in discussing row-
major order to read the dimensions from left to right, with the right-
most being the contiguous dimension.

Conversely, it's not reasonable to expect FFTW's data formats to be
determined by how you happen to display your 2d data on a screen.

The fact is, when you use any library that processes multidimensional
data, you have to find out how the data must be stored---it is foolish
to guess.   And we say *over and over* in the manual that the data is
in row-major order.  In the tutorial.  In the reference section.  In
the section on multi-dimensional array formats.  As far as I can tell,
in every place the dimensions are discussed.

I'm not sure what you expect us to do.  If we label the dimensions
(nz,ny,nx), then at *least* as many people will be confused in the
opposite direction.  And for people who just read the function
prototype and don't read the accompanying explanations in the manual,
nothing we write will make any difference...

Regards,
Steven G. Johnson