This chapter describes the basic usage of FFTW for computing complex multi-dimensional Fourier Transforms. See section FFTW Reference, for more information.
The syntax for performing multi-dimensional transforms in FFTW is very
similar to the syntax for one-dimensional transforms; the main
difference is that `fftw_' gets replaced by `fftwnd_'.
(fftwnd stands for "fftw in N dimensions.")
As before, include <fftw.h>
at the top of your program:
#include <fftw.h>
Then, just as in one dimension, you must create a plan before computing
a transform of any given size. The plan is created by
fftwnd_create_plan.
fftwnd_plan fftwnd_create_plan(int rank, const int *n,
fftw_direction dir, int flags);
The arguments of this function are analogous to those of
fftw_create_plan. rank is the dimensionality of the
array, and can be any non-negative integer. The next argument,
n, is a pointer to an array of length rank containing the
(positive) sizes of each dimension of the array. (Note that the array
will be stored in row-major order. See section Multi-dimensional Array Format, for more information.) dir, as in the one-dimensional
transform, specifies the direction of the transform, and can be
FFTW_FORWARD or FFTW_BACKWARD (or, equivalently, the sign,
or respectively, of the exponent of the transform).
The last argument, flags, is a bitwise OR (`|') of options.
Two possible options are FFTW_ESTIMATE (the default) and
FFTW_MEASURE; just as in fftw_create_plan, these specify
whether the code just estimates the optimal plan or actually performs
timing measurements. Also, if you intend to perform in-place
transforms, you must include the option FFTW_IN_PLACE in the
flags (the default is to perform out-of-place transforms). (An
in-place transform is one in which the output data overwrite the
input data. It thus requires half as much memory as--and is often
faster than--its opposite, an out-of-place transform.)
For two and three dimensional transforms, there are optional,
alternative routines that allow you to simply pass the sizes of each
dimension directly, rather than specifying a rank and an array of sizes.
These are otherwise identical to fftwnd_create_plan, and are
sometimes more convenient:
fftwnd_plan fftw2d_create_plan(int nx, int ny,
fftw_direction dir, int flags);
fftwnd_plan fftw3d_create_plan(int nx, int ny, int nz,
fftw_direction dir, int flags);
Once the plan has been created, you can use it any number of times for
transforms of the same size. The subroutine that performs the transform
is fftwnd:
void fftwnd(fftwnd_plan plan, int howmany,
FFTW_COMPLEX *in, int istride, int idist,
FFTW_COMPLEX *out, int ostride, int odist);
This means: compute howmany multi-dimensional FFTs according to
the plan plan (which contains the size and rank of the
transform). The first (0th) element of the input array is pointed to by
in, and consists of howmany subarrays, placed at a
distance idist from each other, and whose elements are
istride apart and are stored in row-major order.
(See section Multi-dimensional Array Format, for more detail.) In other
words, the element at the i-th location in the j-th
subarray will be in position in[i * istride + j * idist], where
i and j are both zero-based. (Note that i here
refers to the actual index into the row-major-ordered array, and not an
index in any particular dimension.) If the plan was created with the
FFTW_IN_PLACE option, the results are stored in the in
array with the same stride, etcetera, and the last three arguments are
ignored. Otherwise, the result will be stored in the out array,
which consists of howmany subarrays, placed at a distance
odist from each other, and whose elements are ostride
apart.
Finally, don't forget to destroy the plan once you are done using it:
void fftwnd_destroy_plan(fftwnd_plan plan);
This frees all the memory and data structures associated with the plan.
As a simple example, let us compute an out-of-place, forward transform
of a four-dimensional array. In this case, the strides are all 1
and, since we are only computing a single transform (howmany is
1), the idist and odist parameters are ignored.
/* declare static in/out arrays */
FFTW_COMPLEX in[4][5][6][7], out[4][5][6][7];
/* the dimensions of the array */
int n[4] = { 4, 5, 6, 7 };
fftwnd_plan plan;
/* Create the plan (just estimate the optimal strategy): */
plan = fftwnd_create_plan(4, n, FFTW_FORWARD,
FFTW_ESTIMATE);
/* ... initialize the in array ... */
/* Compute the transform, writing the result to out: */
fftwnd(plan, 1, (FFTW_COMPLEX *) in, 1, 0, (FFTW_COMPLEX *) out, 1, 0);
fftwnd_destroy_plan(plan);
Note that the type-casts were necessary in order to get an ordinary
pointer to the first element. Alternatively, we could have used
&in[0][0][0][0] and &out[0][0][0][0].
The previous example is probably typical of the needs of most users.
However, more sophisticated uses of fftwnd are possible. For a
more complicated example involving strides, consider the following:
Suppose that we have a three dimensional array, but the array elements,
rather than being single FFTW_COMPLEX values, are "vectors," or
arrays of three FFTW_COMPLEX values. (Thus, our data form a
complex vector field in three dimensions.) We will declare the
(10x11x12) array as follows:
typedef COMPLEX cvector[3]; cvector vector_field[10][11][12];
Now, we will create the plan for an in-place transform (using
fftw3d_create_plan instead of fftwnd_create_plan just for
variety):
fftwnd_plan plan;
plan = fftw3d_create_plan(10, 11, 12, FFTW_FORWARD,
FFTW_ESTIMATE | FFTW_IN_PLACE);
Finally, let us transform each component of the vector field
separately (i.e. three 3D transforms, and not one 4D transform),
using three separate calls to fftwnd (note that the last three
parameters to fftwnd are ignored since we are transforming
in-place). We use an istride of 3 since each vector has three
components that are stored in adjacent locations in memory.
/* transform 1st component */ fftwnd(plan, 1, &in[0][0][0][0], 3, 0, 0, 0, 0); /* transform 2nd component */ fftwnd(plan, 1, &in[0][0][0][1], 3, 0, 0, 0, 0); /* transform 3rd component */ fftwnd(plan, 1, &in[0][0][0][2], 3, 0, 0, 0, 0);
Alternatively, we could use the howmany parameter to perform all
three transforms with one subroutine call. Here, we use an idist
of 1 since the components of the first element of the array are stored
in adjacent locations. howmany is 3 since there are three arrays
to transform, and the stride is unchanged:
/* transform all components at once */ fftwnd(plan, 3, &in[0][0][0][0], 3, 1, 0, 0, 0);
Note that we used the same plan in all instances, since the dimensions of the data to transform were the same, even though other parameters were different.