GNU Info

(fftw.info)Real Multi-dimensional Transforms Tutorial

---

Next: Multi-dimensional Array Format Prev: Real One-dimensional Transforms Tutorial Up: Tutorial

---

Enter node , (file) or (file)node

Real Multi-dimensional Transforms Tutorial
==========================================

   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 `#include
<rfftw.h>' and then create a plan for the size and direction of
transform that you are interested in:

     rfftwnd_plan rfftwnd_create_plan(int rank, const int *n,
                                      fftw_direction dir, int flags);

   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 `rfftw_create_plan'.  Just as for
fftwnd, there are two alternate versions of this routine,
`rfftw2d_create_plan' and `rfftw3d_create_plan', that are sometimes
more convenient for two- and three-dimensional transforms.  Also as in
fftwnd, rfftwnd supports true in-place transforms, specified by
including `FFTW_IN_PLACE' in the flags.

   Once created, a plan can be used for any number of transforms, and is
deallocated by calling `rfftwnd_destroy_plan(plan)'.

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

     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);

   As is clear from their names and parameter types, the former
function is for `FFTW_REAL_TO_COMPLEX' transforms and the latter is for
`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 `out' parameter is ignored for in-place transforms.

   The format of the complex arrays deserves careful attention.
Suppose that the real data has dimensions n1 x n2 x ... x nd (in
row-major order).  Then, after a real-to-complex transform, the output
is an n1 x n2 x ... x (nd/2+1) array of `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.  That is, the last dimension of
the real data must physically contain 2 * (nd/2+1) `fftw_real' values
(exactly enough to hold the complex data).  This physical array size
does not, however, change the *logical* array size--only nd values are
actually stored in the last dimension, and nd is the last dimension
passed to `rfftwnd_create_plan'.

   For example, consider the transform of a two-dimensional real array
of size `nx' by `ny'.  The output of the `rfftwnd' transform is a
two-dimensional real array of size `nx' by `ny/2+1', where the `y'
dimension has been cut nearly in half because of redundancies in the
output.  Because `fftw_complex' is twice the size of `fftw_real', the
output array is slightly bigger than the input array.  Thus, if we want
to compute the transform in place, we must *pad* the input array so
that it is of size `nx' by `2*(ny/2+1)'.  If `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).

   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.

   Below, we illustrate the use of RFFTWND by showing how you might use
it to compute the (cyclic) convolution of two-dimensional real arrays
`a' and `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.

     #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);
     }

   We access the complex outputs of the in-place transforms by casting
each real array to a `fftw_complex' pointer.  Because this is a "flat"
pointer, we have to compute the row-major index `ij' explicitly in the
convolution product loop.  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.  Notice the limits of the loops and
the dimensions of the various arrays.

   As with the one-dimensional RFFTW, an out-of-place
`FFTW_COMPLEX_TO_REAL' transform has the side-effect of overwriting its
input array.  (The real-to-complex transform, on the other hand, leaves
its input array untouched.)  If you use RFFTWND for a rank-one
transform, however, this side-effect does not occur.  Because of this
fact (and the simpler output format), users may find the RFFTWND
interface more convenient than RFFTW for one-dimensional transforms.
However, RFFTWND in one dimension is slightly slower than RFFTW because
RFFTWND uses an extra buffer array internally.