GNU Info

(fftw.info)Real One-dimensional Transforms Tutorial

---

Next: Real Multi-dimensional Transforms Tutorial Prev: Complex Multi-dimensional Transforms Tutorial Up: Tutorial

---

Enter node , (file) or (file)node

Real One-dimensional Transforms Tutorial
========================================

   If the input data are purely real, you can save roughly a factor of
two in both time and storage by using the "rfftw" transforms, which are
FFTs specialized for real data.  The output of a such a transform is a
"halfcomplex" array, which consists of only half of the complex DFT
amplitudes (since the negative-frequency amplitudes for real data are
the complex conjugate of the positive-frequency amplitudes).

   In exchange for these speed and space advantages, the user sacrifices
some of the simplicity of FFTW's complex transforms.  First of all, to
allow maximum performance, the output format of the one-dimensional real
transforms is different from that used by the multi-dimensional
transforms.  Second, the inverse transform (halfcomplex to real) has the
side-effect of destroying its input array.  Neither of these
inconveniences should pose a serious problem for users, but it is
important to be aware of them.  (Both the inconvenient output format
and the side-effect of the inverse transform can be ameliorated for
one-dimensional transforms, at the expense of some performance, by using
instead the multi-dimensional transform routines with a rank of one.)

   The computation of the plan is similar to that for the complex
transforms.  First, you `#include <rfftw.h>'.  Then, you create a plan
(of type `rfftw_plan') by calling:

     rfftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags);

   `n' is the length of the *real* array in the transform (even for
halfcomplex-to-real transforms), and can be any positive integer
(although sizes with small factors are transformed more efficiently).
`dir' is either `FFTW_REAL_TO_COMPLEX' or `FFTW_COMPLEX_TO_REAL'.  The
`flags' parameter is the same as in `fftw_create_plan'.

   Once created, a plan can be used for any number of transforms, and is
deallocated when you are done with it by calling
`rfftw_destroy_plan(plan)'.

   Given a plan, a real-to-complex or complex-to-real transform is
computed by calling:

     void rfftw_one(rfftw_plan plan, fftw_real *in, fftw_real *out);

   (Note that `fftw_real' is an alias for the floating-point type for
which FFTW was compiled.)  Depending upon the direction of the plan,
either the input or the output array is halfcomplex, and is stored in
the following format:

   r0, r1, r2, r(n/2), i((n+1)/2-1), ..., i2, i1

   Here, rk is the real part of the kth output, and ik is the imaginary
part.  (We follow here the C convention that integer division is
rounded down, e.g. 7 / 2 = 3.) For a halfcomplex array `hc[]', the kth
component has its real part in `hc[k]' and its imaginary part in
`hc[n-k]', with the exception of `k' `==' `0' or `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 `n' real values is a halfcomplex array of length `n', and
vice versa.  (1)  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, `n'.

   Let us reiterate here our warning that an `FFTW_COMPLEX_TO_REAL'
transform has the side-effect of destroying its (halfcomplex) input.
The `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):

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

   Programs using RFFTW should link with `-lrfftw -lfftw -lm' on Unix,
or with the FFTW, RFFTW, and math libraries in general.

   ---------- Footnotes ----------

   (1) The output for the multi-dimensional rfftw is a
more-conventional array of `fftw_complex' values, but the format here
permitted us greater efficiency in one dimension.