Go to the first, previous, next, last section, table of contents.


Tutorial

This chapter describes the basic usage of FFTW, i.e., how to compute the Fourier transform of a single array. This chapter tells the truth, but not the whole truth. (1) Specifically, FFTW implements additional routines and flags, providing extra functionality, that are not documented here. See Section FFTW Reference, for more complete information. (Note that you need to compile and install FFTW before you can use it in a program. See Section Installation and Customization, for the details of the installation.)

This chapter is divided into four sections. Section Complex One-dimensional Transforms Tutorial describes the basic usage of the one-dimensional transform of complex data. Section Complex Multi-dimensional Transforms Tutorial describes the basic usage of the multi-dimensional transform of complex data. Section Real One-dimensional Transforms Tutorial describes the one-dimensional transform of real data and its inverse. Finally, Section Real Multi-dimensional Transforms Tutorial describes the multi-dimensional transform of real data and its inverse. We recommend that you read these sections in the order that they are presented.

Complex One-dimensional Transforms Tutorial

The basic usage of FFTW is simple. A typical call to FFTW looks like:

#include <fftw.h>
...
{
     fftw_complex in[N], out[N];
     fftw_plan p;
     ...
     p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
     ...
     fftw_one(p, in, out);
     ...
     fftw_destroy_plan(p);  
}

The first thing we do is to create a plan, which is an object that contains all the data that FFTW needs to compute the FFT, using the following function:

fftw_plan fftw_create_plan(int n, fftw_direction dir, int flags);

The first argument, n, is the size of the transform you are trying to compute. The size n can be any positive integer, but sizes that are products of small factors are transformed most efficiently. The second argument, dir, can be either FFTW_FORWARD or FFTW_BACKWARD, and indicates the direction of the transform you are interested in. Alternatively, you can use the sign of the exponent in the transform, -1 or +1, which corresponds to FFTW_FORWARD or FFTW_BACKWARD respectively. The flags argument is either FFTW_MEASURE or FFTW_ESTIMATE. FFTW_MEASURE means that FFTW actually runs and measures the execution time of several FFTs in order to find the best way to compute the transform of size n. This may take some time, depending on your installation and on the precision of the timer in your machine. FFTW_ESTIMATE, on the contrary, does not run any computation, and just builds a reasonable plan, which may be sub-optimal. In other words, if your program performs many transforms of the same size and initialization time is not important, use FFTW_MEASURE; otherwise use the estimate. (A compromise between these two extremes exists. See Section Words of Wisdom.)

Once the plan has been created, you can use it as many times as you like for transforms on arrays of the same size. When you are done with the plan, you deallocate it by calling fftw_destroy_plan(plan).

The transform itself is computed by passing the plan along with the input and output arrays to fftw_one:

void fftw_one(fftw_plan plan, fftw_complex *in, fftw_complex *out);

Note that the transform is out of place, and operates on arrays of type fftw_complex, a data structure with real (in[i].re) and imaginary (in[i].im) floating-point components. The in and out arrays should have the length specified when the plan was created. An alternative function, fftw, allows you to efficiently perform multiple and/or strided transforms (see Section FFTW Reference).

The DFT results are stored in-order in the array out, with the zero-frequency (DC) component in out[0]. The array in is not modified. Users should note that FFTW computes an unnormalized DFT, the sign of whose exponent is given by the dir parameter of fftw_create_plan. Thus, computing a forward followed by a backward transform (or vice versa) results in the original array scaled by n. See Section What FFTW Really Computes, for the definition of DFT.

A program using FFTW should be linked with -lfftw -lm on Unix systems, or with the FFTW and standard math libraries in general.

Complex Multi-dimensional Transforms Tutorial

FFTW can also compute transforms of any number of dimensions (rank). The syntax is similar to that for the one-dimensional transforms, with `fftw_' replaced by `fftwnd_' (which stands for "fftw in N dimensions").

As before, we #include <fftw.h> and create a plan for the transforms, this time of type fftwnd_plan:

fftwnd_plan fftwnd_create_plan(int rank, const int *n,
                               fftw_direction dir, int flags);

rank is the dimensionality of the array, and can be any non-negative integer. The next argument, n, is a pointer to an integer 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 information on row-major order.) The last two parameters are the same as in fftw_create_plan. We now, however, have an additional possible flag, FFTW_IN_PLACE, since fftwnd supports true in-place transforms. Multiple flags are combined using a bitwise or (`|'). (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, FFTWND provides alternative routines that accept the sizes of each dimension directly, rather than indirectly through 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. When you do not need a plan anymore, you can deallocate the plan by calling fftwnd_destroy_plan(plan).

Given a plan, you can compute the transform of an array of data by calling:

void fftwnd_one(fftwnd_plan plan, fftw_complex *in, fftw_complex *out);

Here, in and out point to multi-dimensional arrays in row-major order, of the size specified when the plan was created. In the case of an in-place transform, the out parameter is ignored and the output data are stored in the in array. The results are stored in-order, unnormalized, with the zero-frequency component in out[0]. A forward followed by a backward transform (or vice-versa) yields the original data multiplied by the size of the array (i.e. the product of the dimensions). See Section What FFTWND Really Computes, for a discussion of what FFTWND computes.

For example, code to perform an in-place FFT of a three-dimensional array might look like:

#include <fftw.h>
...
{
     fftw_complex in[L][M][N];
     fftwnd_plan p;
     ...
     p = fftw3d_create_plan(L, M, N, FFTW_FORWARD,
                            FFTW_MEASURE | FFTW_IN_PLACE);
     ...
     fftwnd_one(p, &in[0][0][0], NULL);
     ...
     fftwnd_destroy_plan(p);  
}

Note that in is a statically-declared array, which is automatically in row-major order, but we must take the address of the first element in order to fit the type expected by fftwnd_one. (See Section Multi-dimensional Array Format.)

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, ..., rn/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. (2) 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.

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.

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.


Go to the first, previous, next, last section, table of contents.