fftw.texi

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

@ifhtml
<p align=center>
r<sub>0</sub>, r<sub>1</sub>, r<sub>2</sub>, ..., r<sub>n/2</sub>, i<sub>(n+1)/2-1</sub>, ..., i<sub>2</sub>, i<sub>1</sub>
</p>
@end ifhtml

Here,
@ifinfo
rk
@end ifinfo
@tex
$r_k$
@end tex
@ifhtml
r<sub>k</sub>
@end ifhtml
is the real part of the @math{k}th output, and
@ifinfo
ik
@end ifinfo
@tex
$i_k$
@end tex
@ifhtml
i<sub>k</sub>
@end ifhtml
is the imaginary part.  (We follow here the C convention that integer
division is rounded down, e.g. @math{7 / 2 = 3}.) For a halfcomplex
array @code{hc[]}, the @math{k}th component has its real part in
@code{hc[k]} and its imaginary part in @code{hc[n-k]}, with the
exception of @code{k} @code{==} @code{0} or @code{n/2} (the latter only
if n is even)---in these two cases, the imaginary part is zero due to
symmetries of the real-complex transform, and is not stored.  Thus, the
transform of @code{n} real values is a halfcomplex array of length
@code{n}, and vice versa.  @footnote{The output for the
multi-dimensional rfftw is a more-conventional array of
@code{fftw_complex} values, but the format here permitted us greater
efficiency in one dimension.}  This is actually only half of the DFT
spectrum of the data.  Although the other half can be obtained by
complex conjugation, it is not required by many applications such as
convolution and filtering.

Like the complex transforms, the RFFTW transforms are unnormalized, so a
forward followed by a backward transform (or vice-versa) yields the
original data scaled by the length of the array, @code{n}.
@cindex normalization

Let us reiterate here our warning that an @code{FFTW_COMPLEX_TO_REAL}
transform has the side-effect of destroying its (halfcomplex) input.
The @code{FFTW_REAL_TO_COMPLEX} transform, however, leaves its (real)
input untouched, just as you would hope.

As an example, here is an outline of how you might use RFFTW to compute
the power spectrum of a real array (i.e. the squares of the absolute
values of the DFT amplitudes):
@cindex power spectrum

@example
#include <rfftw.h>
...
@{
     fftw_real in[N], out[N], power_spectrum[N/2+1];
     rfftw_plan p;
     int k;
     ...
     p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
     ...
     rfftw_one(p, in, out);
     power_spectrum[0] = out[0]*out[0];  /* DC component */
     for (k = 1; k < (N+1)/2; ++k)  /* (k < N/2 rounded up) */
          power_spectrum[k] = out[k]*out[k] + out[N-k]*out[N-k];
     if (N % 2 == 0) /* N is even */
          power_spectrum[N/2] = out[N/2]*out[N/2];  /* Nyquist freq. */
     ...
     rfftw_destroy_plan(p);
@}
@end example

Programs using RFFTW should link with @code{-lrfftw -lfftw -lm} on Unix,
or with the FFTW, RFFTW, and math libraries in general.
@cindex linking on Unix

@node Real Multi-dimensional Transforms Tutorial, Multi-dimensional Array Format, Real One-dimensional Transforms Tutorial, Tutorial
@section Real Multi-dimensional Transforms Tutorial
@cindex real multi-dimensional transform

FFTW includes multi-dimensional transforms for real data of any rank.
As with the one-dimensional real transforms, they save roughly a factor
of two in time and storage over complex transforms of the same size.
Also as in one dimension, these gains come at the expense of some
increase in complexity---the output format is different from the
one-dimensional RFFTW (and is more similar to that of the complex FFTW)
and the inverse (complex to real) transforms have the side-effect of
overwriting their input data.  

To use the real multi-dimensional transforms, you first @code{#include
<rfftw.h>} and then create a plan for the size and direction of
transform that you are interested in:

@example
rfftwnd_plan rfftwnd_create_plan(int rank, const int *n,
                                 fftw_direction dir, int flags);
@end example
@tindex rfftwnd_plan
@findex rfftwnd_create_plan

The first two parameters describe the size of the real data (not the
halfcomplex data, which will have different dimensions).  The last two
parameters are the same as those for @code{rfftw_create_plan}.  Just as
for fftwnd, there are two alternate versions of this routine,
@code{rfftw2d_create_plan} and @code{rfftw3d_create_plan}, that are
sometimes more convenient for two- and three-dimensional transforms.
@findex rfftw2d_create_plan
@findex rfftw3d_create_plan
Also as in fftwnd, rfftwnd supports true in-place transforms, specified
by including @code{FFTW_IN_PLACE} in the flags.

Once created, a plan can be used for any number of transforms, and is
deallocated by calling @code{rfftwnd_destroy_plan(plan)}.

Given a plan, the transform is computed by calling one of the following
two routines:

@example
void rfftwnd_one_real_to_complex(rfftwnd_plan plan,
                                 fftw_real *in, fftw_complex *out);
void rfftwnd_one_complex_to_real(rfftwnd_plan plan,
                                 fftw_complex *in, fftw_real *out);
@end example
@findex rfftwnd_one_real_to_complex
@findex rfftwnd_one_complex_to_real

As is clear from their names and parameter types, the former function is
for @code{FFTW_REAL_TO_COMPLEX} transforms and the latter is for
@code{FFTW_COMPLEX_TO_REAL} transforms.  (We could have used only a
single routine, since the direction of the transform is encoded in the
plan, but we wanted to correctly express the datatypes of the
parameters.)  The latter routine, as we discuss elsewhere, has the
side-effect of overwriting its input (except when the rank of the array
is one).  In both cases, the @code{out} parameter is ignored for
in-place transforms.

The format of the complex arrays deserves careful attention.
@cindex rfftwnd array format
Suppose that the real data has dimensions
@tex
$n_1 \times n_2 \times \cdots \times n_d$
@end tex
@ifinfo
n1 x n2 x ... x nd
@end ifinfo
@ifhtml
n<sub>1</sub> x n<sub>2</sub> x ... x n<sub>d</sub>
@end ifhtml
(in row-major order).  Then, after a real-to-complex transform, the
output is an
@tex
$n_1 \times n_2 \times \cdots \times (n_d/2+1)$
@end tex
@ifinfo
n1 x n2 x ... x (nd/2+1)
@end ifinfo
@ifhtml
n<sub>1</sub> x n<sub>2</sub> x ... x (n<sub>d</sub>/2+1)
@end ifhtml
array of @code{fftw_complex} values in row-major order, corresponding to
slightly over half of the output of the corresponding complex transform.
(Note that the division is rounded down.)  The ordering of the data is
otherwise exactly the same as in the complex case.  (In principle, the
output could be exactly half the size of the complex transform output,
but in more than one dimension this requires too complicated a format to
be practical.)  Note that, unlike the one-dimensional RFFTW, the real
and imaginary parts of the DFT amplitudes are here stored together in
the natural way.

Since the complex data is slightly larger than the real data, some
complications arise for in-place transforms.  In this case, the final
dimension of the real data must be padded with extra values to
accommodate the size of the complex data---two extra if the last
dimension is even and one if it is odd.
@cindex padding
That is, the last dimension of the real data must physically contain
@tex
$2 (n_d/2+1)$
@end tex
@ifinfo
2 * (nd/2+1)
@end ifinfo
@ifhtml
2 * (n<sub>d</sub>/2+1)
@end ifhtml
@code{fftw_real} values (exactly enough to hold the complex data).
This physical array size does not, however, change the @emph{logical}
array size---only
@tex
$n_d$
@end tex
@ifinfo
nd
@end ifinfo
@ifhtml
n<sub>d</sub>
@end ifhtml
values are actually stored in the last dimension, and
@tex
$n_d$
@end tex
@ifinfo
nd
@end ifinfo
@ifhtml
n<sub>d</sub>
@end ifhtml
is the last dimension passed to @code{rfftwnd_create_plan}.

For example, consider the transform of a two-dimensional real array of
size @code{nx} by @code{ny}.  The output of the @code{rfftwnd} transform
is a two-dimensional complex array of size @code{nx} by @code{ny/2+1},
where the @code{y} dimension has been cut nearly in half because of
redundancies in the output.  Because @code{fftw_complex} is twice the
size of @code{fftw_real}, the output array is slightly bigger than the
input array.  Thus, if we want to compute the transform in place, we
must @emph{pad} the input array so that it is of size @code{nx} by
@code{2*(ny/2+1)}.  If @code{ny} is even, then there are two padding
elements at the end of each row (which need not be initialized, as they
are only used for output).
@ifhtml
The following illustration depicts the input and output arrays just
described, for both the out-of-place and in-place transforms (with the
arrows indicating consecutive memory locations):

<p align=center><img src="rfftwnd.gif" width=389 height=583>
@end ifhtml

@tex
Figure 1 depicts the input and output arrays just
described, for both the out-of-place and in-place transforms (with the
arrows indicating consecutive memory locations).

{
\pageinsert
\vfill
\vskip405pt
\hskip40pt
\special{psfile="rfftwnd.eps" 
}
\vskip 24pt
Figure 1: Illustration of the data layout for real to complex
transforms.  
\vfill
\endinsert}
@end tex

The RFFTWND transforms are unnormalized, so a forward followed by a
backward transform will result in the original data scaled by the number
of real data elements---that is, the product of the (logical) dimensions
of the real data.
@cindex normalization

Below, we illustrate the use of RFFTWND by showing how you might use it
to compute the (cyclic) convolution of two-dimensional real arrays
@code{a} and @code{b} (using the identity that a convolution corresponds
to a pointwise product of the Fourier transforms).  For variety,
in-place transforms are used for the forward FFTs and an out-of-place
transform is used for the inverse transform.
@cindex convolution
@cindex cyclic convolution

@example
#include <rfftw.h>
...
@{
     fftw_real a[M][2*(N/2+1)], b[M][2*(N/2+1)], c[M][N];
     fftw_complex *A, *B, C[M][N/2+1];
     rfftwnd_plan p, pinv;
     fftw_real scale = 1.0 / (M * N);
     int i, j;
     ...
     p    = rfftw2d_create_plan(M, N, FFTW_REAL_TO_COMPLEX,
                                FFTW_ESTIMATE | FFTW_IN_PLACE);
     pinv = rfftw2d_create_plan(M, N, FFTW_COMPLEX_TO_REAL,
                                FFTW_ESTIMATE);

     /* aliases for accessing complex transform outputs: */
     A = (fftw_complex*) &a[0][0];
     B = (fftw_complex*) &b[0][0];
     ...
     for (i = 0; i < M; ++i)
          for (j = 0; j < N; ++j) @{
               a[i][j] = ... ;
               b[i][j] = ... ;
          @}
     ...
     rfftwnd_one_real_to_complex(p, &a[0][0], NULL);
     rfftwnd_one_real_to_complex(p, &b[0][0], NULL);

     for (i = 0; i < M; ++i)
          for (j = 0; j < N/2+1; ++j) @{
               int ij = i*(N/2+1) + j;
               C[i][j].re = (A[ij].re * B[ij].re
                             - A[ij].im * B[ij].im) * scale;
               C[i][j].im = (A[ij].re * B[ij].im
                             + A[ij].im * B[ij].re) * scale;
          @}

     /* inverse transform to get c, the convolution of a and b;
        this has the side effect of overwriting C */
     rfftwnd_one_complex_to_real(pinv, &C[0][0], &c[0][0]);
     ...
     rfftwnd_destroy_plan(p);
     rfftwnd_destroy_plan(pinv);
@}
@end example

We access the complex outputs of the in-place transforms by casting
each real array to a @code{fftw_complex} pointer.  Because this is a
``flat'' pointer, we have to compute the row-major index @code{ij}
explicitly in the convolution product loop.
@cindex row-major
In order to normalize the convolution, we must multiply by a scale
factor---we can do so either before or after the inverse transform, and
choose the former because it obviates the necessity of an additional
loop.
@cindex normalization
Notice the limits of the loops and the dimensions of the various arrays.

As with the one-dimensional RFFTW, an out-of-place
@code{FFTW_COMPLEX_TO_REAL} transform has the side-effect of overwriting