fftw-faq.ascii

FFTW, a collection of fast C routines to compute the Discrete Fourier Transform in one or more dime

If you don't need to compute many transforms and the time for the planner
is significant, you have two options.  First, you can use the
FFTW_ESTIMATE option in the planner, which uses heuristics instead of
runtime measurements and produces a good plan in a short time.  Second,
you can use the wisdom feature to precompute the plan; see Q3.7 `Can I
save FFTW's plans?'

-------------------------------------------------------------------------------

Question 3.2.  FFTW slows down after repeated calls.

Probably, NaNs or similar are creeping into your data, and the slowdown is
due to the resulting floating-point exceptions.  For example, be aware
that repeatedly FFTing the same array is a diverging process (because FFTW
computes the unnormalized transform).

-------------------------------------------------------------------------------

Question 3.3.  An FFTW routine is crashing when I call it.

You almost certainly have a bug in your code.  For example, you could be
passing invalid arguments (such as wrongly-sized arrays) to FFTW, or you
could simply have memory corruption elsewhere in your program that causes
random crashes later on.  Learn to debug, and don't complain to us unless
you can come up with a minimal program (preferably under 30 lines) that
illustrates the problem.

-------------------------------------------------------------------------------

Question 3.4.  My Fortran program crashes when calling FFTW.

As described in the manual, on 64-bit machines you must store the plans in
variables large enough to hold a pointer, for example integer*8.

-------------------------------------------------------------------------------

Question 3.5.  FFTW gives results different from my old FFT.

People follow many different conventions for the DFT, and you should be
sure to know the ones that we use (described in the FFTW manual).  In
particular, you should be aware that the FFTW_FORWARD/FFTW_BACKWARD
directions correspond to signs of -1/+1 in the exponent of the DFT
definition.  (*Numerical Recipes* uses the opposite convention.)

You should also know that we compute an unnormalized transform.  In
contrast, Matlab is an example of program that computes a normalized
transform.  See Q3.8 `Why does your inverse transform return a scaled
result?'.

Finally, note that floating-point arithmetic is not exact, so different
FFT algorithms will give slightly different results (on the order of the
numerical accuracy; typically a fractional difference of 1e-15 or so).

-------------------------------------------------------------------------------

Question 3.6.  Your in-place transform gives incorrect results.

As described in the FFTW manual, the output array argument has a special
meaning for FFTW_INPLACE transforms; you should not pass the input array
for this argument.

-------------------------------------------------------------------------------

Question 3.7.  Can I save FFTW's plans?

Yes. Starting with version 1.2, FFTW provides the  wisdom mechanism for
saving plans.  See Q4.3 `What is this wisdom thing?' and the FFTW manual.

-------------------------------------------------------------------------------

Question 3.8.  Why does your inverse transform return a scaled result?

Computing the forward transform followed by the backward transform (or
vice versa) yields the original array scaled by the size of the array.
(For multi-dimensional transforms, the size of the array is the product of
the dimensions.)  We could, instead, have chosen a normalization that
would have returned the unscaled array. Or, to accomodate the many
conventions in this matter, the transform routines could have accepted a
"scale factor" parameter. We did not do this, however, for two reasons.
First, we didn't want to sacrifice performance in the common case where
the scale factor is 1. Second, in real applications the FFT is followed or
preceded by some computation on the data, into which the scale factor can
typically be absorbed at little or no cost.

-------------------------------------------------------------------------------

Question 3.9.  How can I make FFTW put the origin (zero frequency) at the center of its output?

For human viewing of a spectrum, it is often convenient to put the origin
in frequency space at the center of the output array, rather than in the
zero-th element (the default in FFTW).  If all of the dimensions of your
array are even, you can accomplish this by simply multiplying each element
of the input array by (-1)^(i + j + ...), where i, j, etcetera are the
indices of the element.  (This trick is a general property of the DFT, and
is not specific to FFTW.)

-------------------------------------------------------------------------------

Question 3.10.  How do I FFT an image/audio file in *foobar* format?

FFTW performs an FFT on an array of floating-point values.  You can
certainly use it to compute the transform of an image or audio stream, but
you are responsible for figuring out your data format and converting it to
the form FFTW requires.

-------------------------------------------------------------------------------

Question 3.11.  My program does not link (on Unix).

Please use the exact order in which libraries are specified by the FFTW
manual (e.g. -lrfftw -lfftw -lm).  Also, note that the libraries must be
listed after your program sources/objects.  (The general rule is that if
*A* uses *B*, then *A* must be listed before *B* in the link command.).
For example, switching the order to -lfftw -lrfftw -lm will fail.

-------------------------------------------------------------------------------

Question 3.12.  My program crashes, complaining about stack space.

You cannot declare large arrays statically; you should use malloc (or
equivalent) to allocate the arrays you want to transform if they are
larger than a few hundred elements.

===============================================================================

Section 4.  Internals of FFTW

 Q4.1        How does FFTW work?
 Q4.2        Why is FFTW so fast?
 Q4.3        What is this wisdom thing?
 Q4.4        Why do you use wisdom? I just wanted to save a plan.

-------------------------------------------------------------------------------

Question 4.1.  How does FFTW work?

The innovation (if it can be so called) in FFTW consists in having an
interpreter execute the transform.  The program for the interpreter (the
*plan*) is computed at runtime according to the characteristics of your
machine/compiler.  This peculiar software architecture allows FFTW to
adapt itself to almost any machine.

For more details, see the paper "The Fastest Fourier Transform in the
West", by M. Frigo and S. G. Johnson, available at the FFTW web page.  See
also "FFTW: An Adaptive Software Architecture for the FFT", in ICASSP '98.

-------------------------------------------------------------------------------

Question 4.2.  Why is FFTW so fast?

This is a complex question, and there is no simple answer.  In fact, the
authors do not fully know the answer, either.  In addition to many small
performance hacks throughout FFTW, there are three general reasons for
FFTW's speed.

* 	FFTW uses an internal interpreter to adapt itself to a machine.  See
  Q4.1 `How does FFTW work?'.
* 	FFTW uses a code generator to produce highly-optimized routines for
  computing small transforms.
* 	FFTW uses explicit divide-and-conquer to take advantage of the memory
  hierarchy.

For more details on these three topics, see the paper "The Fastest Fourier
Transform in the West", by M. Frigo and S. G. Johnson, available at the
FFTW web page.

-------------------------------------------------------------------------------

Question 4.3.  What is this wisdom thing?

wisdom is the name of the mechanism that FFTW uses to save and restore
plans.  Rather than just saving plans, FFTW remembers what it learns about
your machine, and becomes wiser and wiser as time passes by.  You can save
wisdom for later use.

-------------------------------------------------------------------------------

Question 4.4.  Why do you use wisdom? I just wanted to save a plan.

wisdom could be implemented with less effort than a general plan-saving
mechanism would have required.  In addition, wisdom provides additional
benefits.  For example, if you are planning transforms of size 1024, and
later you want a transform of size 2048, most of the calculations of the
1024 case can be reused.

In short, wisdom does more things with less effort, and seemed like The
Right Thing to do.

===============================================================================

Section 5.  Known bugs

 Q5.1        FFTW 1.1 crashes in rfftwnd on Linux.
 Q5.2        The MPI transforms in FFTW 1.2 give incorrect results/leak memory.
 Q5.3        The test programs in FFTW 1.2.1 fail when I change FFTW to use sin
 Q5.4        The test program in FFTW 1.2.1 fails for n > 46340.
 Q5.5        The threaded code fails on Linux Redhat 5.0
 Q5.6        FFTW 2.0's rfftwnd fails for rank > 1 transforms with a final dime
 Q5.7        FFTW 2.0's complex transforms give the wrong results with prime fa
 Q5.8        FFTW 2.1.1's MPI test programs crash with MPICH.
 Q5.9        FFTW 2.1.2's multi-threaded transforms don't work on AIX.
 Q5.10       FFTW 2.1.2's complex transforms give incorrect results for large p
 Q5.11       FFTW 2.1.3 crashes on AIX

-------------------------------------------------------------------------------

Question 5.1.  FFTW 1.1 crashes in rfftwnd on Linux.

This bug was fixed in FFTW 1.2.  There was a bug in rfftwnd causing an
incorrect amount of memory to be allocated.  The bug showed up in Linux
with libc-5.3.12 (and nowhere else that we know of).

-------------------------------------------------------------------------------

Question 5.2.  The MPI transforms in FFTW 1.2 give incorrect results/leak memory.

These bugs were corrected in FFTW 1.2.1.  The MPI transforms (really, just
the transpose routines) in FFTW 1.2 had bugs that could cause errors in
some situations.

-------------------------------------------------------------------------------

Question 5.3.  The test programs in FFTW 1.2.1 fail when I change FFTW to use single precision.

This bug was fixed in FFTW 1.3.  (Older versions of FFTW did work in
single precision, but the test programs didn't--the error tolerances in
the tests were set for double precision.)

-------------------------------------------------------------------------------

Question 5.4.  The test program in FFTW 1.2.1 fails for n > 46340.

This bug was fixed in FFTW 1.3.  FFTW 1.2.1 produced the right answer, but
the test program was wrong.  For large n, n*n in the naive transform that
we used for comparison overflows 32 bit integer precision, breaking the
test.

-------------------------------------------------------------------------------

Question 5.5.  The threaded code fails on Linux Redhat 5.0

We had problems with glibc-2.0.5.  The code should work with glibc-2.0.7.

-------------------------------------------------------------------------------

Question 5.6.  FFTW 2.0's rfftwnd fails for rank > 1 transforms with a final dimension >= 65536.

This bug was fixed in FFTW 2.0.1.  (There was a 32-bit integer overflow
due to a poorly-parenthesized expression.)

-------------------------------------------------------------------------------

Question 5.7.  FFTW 2.0's complex transforms give the wrong results with prime factors 17 to 97.

There was a bug in the complex transforms that could cause incorrect
results under (hopefully rare) circumstances for lengths with
intermediate-size prime factors (17-97).  This bug was fixed in FFTW
2.1.1.

-------------------------------------------------------------------------------

Question 5.8.  FFTW 2.1.1's MPI test programs crash with MPICH.

This bug was fixed in FFTW 2.1.2.  The 2.1/2.1.1 MPI test programs crashed
when using the MPICH implementation of MPI with the ch_p4 device (TCP/IP);
the transforms themselves worked fine.

-------------------------------------------------------------------------------

Question 5.9.  FFTW 2.1.2's multi-threaded transforms don't work on AIX.

This bug was fixed in FFTW 2.1.3.  The multi-threaded transforms in
previous versions didn't work with AIX's pthreads implementation, which
idiosyncratically creates threads in detached (non-joinable) mode by
default.

-------------------------------------------------------------------------------

Question 5.10.  FFTW 2.1.2's complex transforms give incorrect results for large prime sizes.

This bug was fixed in FFTW 2.1.3.  FFTW's complex-transform algorithm for
prime sizes (in versions 2.0 to 2.1.2) had an integer overflow problem
that caused incorrect results for many primes greater than 32768 (on
32-bit machines).  (Sizes without large prime factors are not affected.)

-------------------------------------------------------------------------------

Question 5.11.  FFTW 2.1.3 crashes on AIX

The FFTW 2.1.3 configure script picked incorrect compiler flags for the
xlc compiler on newer IBM processors.  This is fixed in FFTW 2.1.4.