fftw_2.html

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

<P>
(Note that <CODE>fftw_real</CODE> 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:
<A NAME="IDX60"></A>


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

<P>
Here,
r<sub>k</sub>
is the real part of the kth output, and
i<sub>k</sub>
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 <CODE>hc[]</CODE>, the kth component has its real part in
<CODE>hc[k]</CODE> and its imaginary part in <CODE>hc[n-k]</CODE>, with the
exception of <CODE>k</CODE> <CODE>==</CODE> <CODE>0</CODE> or <CODE>n/2</CODE> (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</CODE> real values is a halfcomplex array of length
<CODE>n</CODE>, and vice versa.  <A NAME="DOCF1" HREF="fftw_foot.html#FOOT1">(1)</A>  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.


<P>
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</CODE>.
<A NAME="IDX61"></A>


<P>
Let us reiterate here our warning that an <CODE>FFTW_COMPLEX_TO_REAL</CODE>
transform has the side-effect of destroying its (halfcomplex) input.
The <CODE>FFTW_REAL_TO_COMPLEX</CODE> transform, however, leaves its (real)
input untouched, just as you would hope.


<P>
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):
<A NAME="IDX62"></A>



<PRE>
#include &#60;rfftw.h&#62;
...
{
     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 &#60; (N+1)/2; ++k)  /* (k &#60; 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);
}
</PRE>

<P>
Programs using RFFTW should link with <CODE>-lrfftw -lfftw -lm</CODE> on Unix,
or with the FFTW, RFFTW, and math libraries in general.
<A NAME="IDX63"></A>




<H2><A NAME="SEC6">Real Multi-dimensional Transforms Tutorial</A></H2>
<P>
<A NAME="IDX64"></A>


<P>
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.  


<P>
To use the real multi-dimensional transforms, you first <CODE>#include
&#60;rfftw.h&#62;</CODE> and then create a plan for the size and direction of
transform that you are interested in:



<PRE>
rfftwnd_plan rfftwnd_create_plan(int rank, const int *n,
                                 fftw_direction dir, int flags);
</PRE>

<P>
<A NAME="IDX65"></A>
<A NAME="IDX66"></A>


<P>
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</CODE>.  Just as
for fftwnd, there are two alternate versions of this routine,
<CODE>rfftw2d_create_plan</CODE> and <CODE>rfftw3d_create_plan</CODE>, that are
sometimes more convenient for two- and three-dimensional transforms.
<A NAME="IDX67"></A>
<A NAME="IDX68"></A>
Also as in fftwnd, rfftwnd supports true in-place transforms, specified
by including <CODE>FFTW_IN_PLACE</CODE> in the flags.


<P>
Once created, a plan can be used for any number of transforms, and is
deallocated by calling <CODE>rfftwnd_destroy_plan(plan)</CODE>.


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



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

<P>
<A NAME="IDX69"></A>
<A NAME="IDX70"></A>


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


<P>
The format of the complex arrays deserves careful attention.
<A NAME="IDX71"></A>
Suppose that the real data has dimensions
n<sub>1</sub> x n<sub>2</sub> x ... x n<sub>d</sub>
(in row-major order).  Then, after a real-to-complex transform, the
output is an
n<sub>1</sub> x n<sub>2</sub> x ... x (n<sub>d</sub>/2+1)
array of <CODE>fftw_complex</CODE> 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.


<P>
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.
<A NAME="IDX72"></A>
That is, the last dimension of the real data must physically contain
2 * (n<sub>d</sub>/2+1)
<CODE>fftw_real</CODE> values (exactly enough to hold the complex data).
This physical array size does not, however, change the <EM>logical</EM>
array size--only
n<sub>d</sub>
values are actually stored in the last dimension, and
n<sub>d</sub>
is the last dimension passed to <CODE>rfftwnd_create_plan</CODE>.


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


<P>
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.
<A NAME="IDX73"></A>


<P>
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</CODE> and <CODE>b</CODE> (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.
<A NAME="IDX74"></A>
<A NAME="IDX75"></A>



<PRE>
#include &#60;rfftw.h&#62;
...
{
     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*) &#38;a[0][0];
     B = (fftw_complex*) &#38;b[0][0];
     ...
     for (i = 0; i &#60; M; ++i)
          for (j = 0; j &#60; N; ++j) {
               a[i][j] = ... ;
               b[i][j] = ... ;
          }
     ...
     rfftwnd_one_real_to_complex(p, &#38;a[0][0], NULL);
     rfftwnd_one_real_to_complex(p, &#38;b[0][0], NULL);

     for (i = 0; i &#60; M; ++i)
          for (j = 0; j &#60; 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, &#38;C[0][0], &#38;c[0][0]);
     ...
     rfftwnd_destroy_plan(p);
     rfftwnd_destroy_plan(pinv);
}
</PRE>

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


<P>
As with the one-dimensional RFFTW, an out-of-place
<CODE>FFTW_COMPLEX_TO_REAL</CODE> 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.




<H2><A NAME="SEC7">Multi-dimensional Array Format</A></H2>

<P>
This section describes the format in which multi-dimensional arrays are
stored.  We felt that a detailed discussion of this topic was necessary,
since it is often a source of confusion among users and several
different formats are common.  Although the comments below refer to
<CODE>fftwnd</CODE>, they are also applicable to the <CODE>rfftwnd</CODE> routines.




<H3><A NAME="SEC8">Row-major Format</A></H3>
<P>
<A NAME="IDX78"></A>


<P>
The multi-dimensional arrays passed to <CODE>fftwnd</CODE> are expected to be
stored as a single contiguous block in <EM>row-major</EM> order (sometimes
called "C order").  Basically, this means that as you step through
adjacent memory locations, the first dimension's index varies most
slowly and the last dimension's index varies most quickly.


<P>
To be more explicit, let us consider an array of rank d whose
dimensions are
n<sub>1</sub> x n<sub>2</sub> x n<sub>3</sub> x ... x n<sub>d</sub>.
Now, we specify a location in the array by a sequence of (zero-based) indices,
one for each dimension:
(i<sub>1</sub>, i<sub>2</sub>, i<sub>3</sub>,..., i<sub>d</sub>).
If the array is stored in row-major
order, then this element is located at the position
i<sub>d</sub> + n<sub>d</sub> * (i<sub>d-1</sub> + n<sub>d-1</sub> * (... + n<sub>2</sub> * i<sub>1</sub>)).


<P>
Note that each element of the array must be of type <CODE>fftw_complex</CODE>;
i.e. a (real, imaginary) pair of (double-precision) numbers. Note also
that, in <CODE>fftwnd</CODE>, the expression above is multiplied by the stride
to get the actual array index--this is useful in situations where each
element of the multi-dimensional array is actually a data structure or
another array, and you just want to transform a single field. In most
cases, however, you use a stride of 1.
<A NAME="IDX79"></A>




<H3><A NAME="SEC9">Column-major Format</A></H3>
<P>
<A NAME="IDX80"></A>


<P>
Readers from the Fortran world are used to arrays stored in
<EM>column-major</EM> order (sometimes called "Fortran order").  This is
essentially the exact opposite of row-major order in that, here, the
<EM>first</EM> dimension's index varies most quickly.


<P>
If you have an array stored in column-major order and wish to transform
it using <CODE>fftwnd</CODE>, it is quite easy to do.  When creating the plan,
simply pass the dimensions of the array to <CODE>fftwnd_create_plan</CODE> in
<EM>reverse order</EM>.  For example, if your array is a rank three
<CODE>N x M x L</CODE> matrix in column-major order, you should pass the
dimensions of the array as if it were an <CODE>L x M x N</CODE> matrix (which
it is, from the perspective of <CODE>fftwnd</CODE>).  This is done for you
automatically by the FFTW Fortran wrapper routines (see Section <A HREF="fftw_5.html#SEC62">Calling FFTW from Fortran</A>).
<A NAME="IDX81"></A>