note: look inside the package for the author !
(its NOT done by me, joerg arndt)
Here is an n-dimensional FFT routine I wrote recently. It's based
on the Fortrash and Pascal routines in the book "Numerical Recipes"
(I forgot the author, and the book is at work.), but much of it has
been changed (to take advantage of C's logical operations for the bit
reversal, to use C's array index ordering and starting index of 0,
etc.)
If the preprocessor variable SPEED is defined, data[] must be an
array of fcomplex, not complex, which uses float's instead of double's.
I used the float's in the application I wrote this for to reduce the size
of the array, and to reduce the time taken by memory accesses, since high
accuracy wasn't needed.
The routine has been optimized for speed for MSC 4.0, compiling
for an 80286, using inline 8087/287 code (requires the 87/287), and it is
probably not optimal for other systems.
I have included 2 options for the complex multiplication, a complex
multiply using 4 real multiplies, and one using 3 real multiplies.
The 3 multiply method is used unless the preprocessor variable MULT4 is
defined (which is automatically done for MSC). Usually, the 3 multiply
method should be faster, but unfortunately, MSC 4.0 does almost no
optimization for the 80287, so it was storing t1 and t2 back to memory
and reading them back, which all took more time than doing the extra multiply!
Now for the usage (finally):
data[] is the array of complex numbers to be transformed,
nn[] is the array giving the dimensions (I mean size) of the array,
ndim is the number of dimensions of the array, and
isign is +1 for a forward transform, and -1 for an inverse transform.
data[] and nn[] are stored in the "natural" order for C:
nn[0] gives the number of elements along the leftmost index,
nn[ndim - 1] gives the number of elements along the rightmost index, and
data should be declared along the lines of
struct (f)complex data[nn[0], nn[1], ..., nn[ndim - 1]]
Additional notes: The routine does NO NORMALIZATION, so if you do a
forward, and then an inverse transform on an array, the result will
be identical to the original array MULTIPLIED BY THE NUMBER OF
ELEMENTS IN THE ARRAY. Also, of course, the dimensions of data[]
must all be powers of 2. I have also enclosed a version of the demo
program from the examples book for "Numerical Recipes".