FFTW Reference
**************

This chapter provides a complete reference for all sequential (i.e., one-processor) FFTW functions.  Parallel transforms are described in later chapters.

Data Types and Files
====================

All programs using FFTW should include its header file:

     #include <fftw3.h>

You must also link to the FFTW library.  On Unix, this means adding -lfftw3 -lm at the end of the link command.

Complex numbers
---------------

The default FFTW interface uses double precision for all floating-point numbers, and defines a fftw_complex type to hold complex numbers as:

     typedef double fftw_complex[2];

Here, the [0] element holds the real part and the [1] element holds the imaginary part.

Alternatively, if you have a C compiler (such as gcc) that supports the C99 revision of the ANSI C standard, you can use C's new native complex type (which is binary-compatible with the typedef above). In particular, if you #include <complex.h> before <fftw3.h>, then fftw_complex is defined to be the native complex type and you can manipulate it with ordinary arithmetic (e.g. x = y * (3+4*I), where x and y are fftw_complex and I is the standard symbol for the imaginary unit);

C++ has its own complex<T> template class, defined in the standard <complex> header file.  Reportedly, the C++ standards committee has recently agreed to mandate that the storage format used for this type be binary-compatible with the C99 type, i.e. an array T[2] with consecutive real [0] and imaginary [1] parts.  (See report http://www.open-std.org/jtc1/sc22/WG21/docs/papers/2002/n1388.pdf WG21/N1388.)  Although not part of the official standard as of this writing, the proposal stated that: ``This solution has been tested with all current major implementations of the standard library and shown to be working.''  To the extent that this is true, if you have a variable complex<double> *x, you can pass it directly to FFTW via reinterpret_cast<fftw_complex*>(x).

Precision
---------

You can install single and long-double precision versions of FFTW, which replace double with float and long double, respectively (see Installation and Customization).  To use these interfaces, you:

  * Link to the single/long-double libraries; on Unix, -lfftw3f or -lfftw3l instead of (or in addition to) -lfftw3.  (You can link to the different-precision libraries simultaneously.)

  * Include the same <fftw3.h> header file.

  * Replace all lowercase instances of fftw_ with fftwf_ or fftwl_ for single or long-double precision, respectively. (fftw_complex becomes fftwf_complex, fftw_execute becomes fftwf_execute, etcetera.)

  * Uppercase names, i.e. names beginning with FFTW_, remain the same.

  * Replace double with float or long double for subroutine parameters.

Depending upon your compiler and/or hardware, long double may not be any more precise than double (or may not be supported at all, although it is standard in C99).

We also support using the nonstandard __float128 quadruple-precision type provided by recent versions of gcc on 32- and 64-bit x86 hardware (see Installation and Customization). To use this type, link with -lfftw3q -lquadmath -lm (the libquadmath library provided by gcc is needed for quadruple-precision trigonometric functions) and use fftwq_ identifiers.

Memory Allocation
-----------------

     void *fftw_malloc(size_t n);
     void fftw_free(void *p);

These are functions that behave identically to malloc and free, except that they guarantee that the returned pointer obeys any special alignment restrictions imposed by any algorithm in FFTW (e.g. for SIMD acceleration).  See SIMD alignment and fftw_malloc.

Data allocated by fftw_malloc must be deallocated by fftw_free and not by the ordinary free.

These routines simply call through to your operating system's malloc or, if necessary, its aligned equivalent (e.g. memalign), so you normally need not worry about any significant time or space overhead.  You are not required to use them to allocate your data, but we strongly recommend it.

Note: in C++, just as with ordinary malloc, you must typecast the output of fftw_malloc to whatever pointer type you are allocating.

We also provide the following two convenience functions to allocate real and complex arrays with n elements, which are equivalent to (double *) fftw_malloc(sizeof(double) * n) and (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * n), respectively:

     double *fftw_alloc_real(size_t n);
     fftw_complex *fftw_alloc_complex(size_t n);

The equivalent functions in other precisions allocate arrays of n elements in that precision.  e.g. fftwf_alloc_real(n) is equivalent to (float *) fftwf_malloc(sizeof(float) * n).

Using Plans
===========

Plans for all transform types in FFTW are stored as type fftw_plan (an opaque pointer type), and are created by one of the various planning routines described in the following sections. An fftw_plan contains all information necessary to compute the transform, including the pointers to the input and output arrays.

     void fftw_execute(const fftw_plan plan);

This executes the plan, to compute the corresponding transform on the arrays for which it was planned (which must still exist).  The plan is not modified, and fftw_execute can be called as many times as desired.

To apply a given plan to a different array, you can use the new-array execute interface.  See New-array Execute Functions.

fftw_execute (and equivalents) is the only function in FFTW guaranteed to be thread-safe; see Thread safety.

This function:

     void fftw_destroy_plan(fftw_plan plan);

deallocates the plan and all its associated data.

FFTW's planner saves some other persistent data, such as the accumulated wisdom and a list of algorithms available in the current configuration.  If you want to deallocate all of that and reset FFTW to the pristine state it was in when you started your program, you can call:

     void fftw_cleanup(void);

After calling fftw_cleanup, all existing plans become undefined, and you should not attempt to execute them nor to destroy them.  You can however create and execute/destroy new plans, in which case FFTW starts accumulating wisdom information again.

fftw_cleanup does not deallocate your plans, however.  To prevent memory leaks, you must still call fftw_destroy_plan before executing fftw_cleanup.

Occasionally, it may useful to know FFTW's internal ``cost'' metric that it uses to compare plans to one another; this cost is proportional to an execution time of the plan, in undocumented units, if the plan was created with the FFTW_MEASURE or other timing-based options, or alternatively is a heuristic cost function for FFTW_ESTIMATE plans.  (The cost values of measured and estimated plans are not comparable, being in different units.  Also, costs from different FFTW versions or the same version compiled differently may not be in the same units.  Plans created from wisdom have a cost of 0 since no timing measurement is performed for them. Finally, certain problems for which only one top-level algorithm was possible may have required no measurements of the cost of the whole plan, in which case fftw_cost will also return 0.)  The cost metric for a given plan is returned by:

     double fftw_cost(const fftw_plan plan);

The following two routines are provided purely for academic purposes (that is, for entertainment).

     void fftw_flops(const fftw_plan plan,
                     double *add, double *mul, double *fma);

Given a plan, set add, mul, and fma to an exact count of the number of floating-point additions, multiplications, and fused multiply-add operations involved in the plan's execution.  The total number of floating-point operations (flops) is add + mul + 2*fma, or add + mul + fma if the hardware supports fused multiply-add instructions (although the number of FMA operations is only approximate because of compiler voodoo).  (The number of operations should be an integer, but we use double to avoid overflowing int for large transforms; the arguments are of type double even for single and long-double precision versions of FFTW.)

     void fftw_fprint_plan(const fftw_plan plan, FILE *output_file);
     void fftw_print_plan(const fftw_plan plan);
     char *fftw_sprint_plan(const fftw_plan plan);

This outputs a ``nerd-readable'' representation of the plan to the given file, to stdout, or two a newly allocated NUL-terminated string (which the caller is responsible for deallocating with free), respectively.

Basic Interface
===============

Recall that the FFTW API is divided into three parts (Gallia est omnis divisa in partes tres (Julius Caesar).): the basic interface computes a single transform of contiguous data, the advanced interface computes transforms of multiple or strided arrays, and the guru interface supports the most general data layouts, multiplicities, and strides.  This section describes the the basic interface, which we expect to satisfy the needs of most users.

Complex DFTs
------------

     fftw_plan fftw_plan_dft_1d(int n0,
                                fftw_complex *in, fftw_complex *out,
                                int sign, unsigned flags);
     fftw_plan fftw_plan_dft_2d(int n0, int n1,
                                fftw_complex *in, fftw_complex *out,
                                int sign, unsigned flags);
     fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2,
                                fftw_complex *in, fftw_complex *out,
                                int sign, unsigned flags);
     fftw_plan fftw_plan_dft(int rank, const int *n,
                             fftw_complex *in, fftw_complex *out,
                             int sign, unsigned flags);

Plan a complex input/output discrete Fourier transform (DFT) in zero or more dimensions, returning an fftw_plan (see Using Plans).

Once you have created a plan for a certain transform type and parameters, then creating another plan of the same type and parameters, but for different arrays, is fast and shares constant data with the first plan (if it still exists).

The planner returns NULL if the plan cannot be created.  In the standard FFTW distribution, the basic interface is guaranteed to return a non-NULL plan.  A plan may be NULL, however, if you are using a customized FFTW configuration supporting a restricted set of transforms.

Arguments
.........

  * rank is the rank of the transform (it should be the size of the array *n), and can be any non-negative integer.  (see Complex Multi-Dimensional DFTs, for the definition of ``rank''.)  The _1d, _2d, and _3d planners correspond to a rank of 1, 2, and 3, respectively.  The rank may be zero, which is equivalent to a rank-1 transform of size 1, i.e. a copy of one number from input to output.

  * n0, n1, n2, or n[0..rank-1] (as appropriate for each routine) specify the size of the transform dimensions.  They can be any positive integer.

       - Multi-dimensional arrays are stored in row-major order with dimensions: n0 x n1; or n0 x n1 x n2; or n[0] x n[1] x ... x n[rank-1]. See Multi-dimensional Array Format.

       - FFTW is best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, where e+f is either 0 or 1, and the other exponents are arbitrary.  Other sizes are computed by means of a slow, general-purpose algorithm (which nevertheless retains Onlogn performance even for prime sizes).  It is possible to customize FFTW for different array sizes; see Installation and Customization. Transforms whose sizes are powers of 2 are especially fast.

  * in and out point to the input and output arrays of the transform, which may be the same (yielding an in-place transform). These arrays are overwritten during planning, unless FFTW_ESTIMATE is used in the flags.  (The arrays need not be initialized, but they must be allocated.)

     If in == out, the transform is in-place and the input array is overwritten. If in != out, the two arrays must not overlap (but FFTW does not check for this condition).

  * sign is the sign of the exponent in the formula that defines the Fourier transform.  It can be -1 (= FFTW_FORWARD) or +1 (= FFTW_BACKWARD).

  * flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

FFTW computes an unnormalized transform: computing a forward followed by a backward transform (or vice versa) will result in the original data multiplied by the size of the transform (the product of the dimensions). For more information, see What FFTW Really Computes.

Planner Flags
-------------

All of the planner routines in FFTW accept an integer flags argument, which is a bitwise OR (|) of zero or more of the flag constants defined below.  These flags control the rigor (and time) of the planning process, and can also impose (or lift) restrictions on the type of transform algorithm that is employed.

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan.  The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

In all  cases, if  wisdom is  available for the  given problem  that was created  with equal-or-greater  planning rigor,  then the  more rigorous wisdom is used.  For example, in FFTW_ESTIMATE mode any available wisdom is used, whereas  in FFTW_PATIENT mode only wisdom created in patient or exhaustive mode can be used.  See Words of Wisdom-Saving Plans.

Planning-rigor flags
....................

  * FFTW_ESTIMATE specifies that, instead of actual measurements of different algorithms, a simple heuristic is used to pick a (probably sub-optimal) plan quickly.  With this flag, the input/output arrays are not overwritten during planning.

  * FFTW_MEASURE tells FFTW to find an optimized plan by actually computing several FFTs and measuring their execution time. Depending on your machine, this can take some time (often a few seconds).  FFTW_MEASURE is the default planning option.

  * FFTW_PATIENT is like FFTW_MEASURE, but considers a wider range of algorithms and often produces a ``more optimal'' plan (especially for large transforms), but at the expense of several times longer planning time (especially for large transforms).

  * FFTW_EXHAUSTIVE is like FFTW_PATIENT, but considers an even wider range of algorithms, including many that we think are unlikely to be fast, to produce the most optimal plan but with a substantially increased planning time.

  * FFTW_WISDOM_ONLY is a special planning mode in which the plan is only created if wisdom is available for the given problem, and otherwise a NULL plan is returned.  This can be combined with other flags, e.g. FFTW_WISDOM_ONLY | FFTW_PATIENT creates a plan only if wisdom is available that was created in FFTW_PATIENT or FFTW_EXHAUSTIVE mode.  The FFTW_WISDOM_ONLY flag is intended for users who need to detect whether wisdom is available; for example, if wisdom is not available one may wish to allocate new arrays for planning so that user data is not overwritten.

Algorithm-restriction flags
...........................

  * FFTW_DESTROY_INPUT specifies that an out-of-place transform is allowed to overwrite its input array with arbitrary data; this can sometimes allow more efficient algorithms to be employed.

  * FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array.  This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUT is the default.  In the latter cases, passing FFTW_PRESERVE_INPUT will attempt to use algorithms that do not destroy the input, at the expense of worse performance; for multi-dimensional c2r transforms, however, no input-preserving algorithms are implemented and the planner will return NULL if one is requested.

  * FFTW_UNALIGNED specifies that the algorithm may not impose any unusual alignment requirements on the input/output arrays (i.e. no SIMD may be used).  This flag is normally not necessary, since the planner automatically detects misaligned arrays.  The only use for this flag is if you want to use the new-array execute interface to execute a given plan on a different array that may not be aligned like the original.  (Using fftw_malloc makes this flag unnecessary even then.  You can also use fftw_alignment_of to detect whether two arrays are equivalently aligned.)

Limiting planning time
......................

     extern void fftw_set_timelimit(double seconds);

This function instructs FFTW to spend at most seconds seconds (approximately) in the planner.  If seconds == FFTW_NO_TIMELIMIT (the default value, which is negative), then planning time is unbounded.  Otherwise, FFTW plans with a progressively wider range of algorithms until the the given time limit is reached or the given range of algorithms is explored, returning the best available plan.

For example, specifying FFTW_PATIENT first plans in FFTW_ESTIMATE mode, then in FFTW_MEASURE mode, then finally (time permitting) in FFTW_PATIENT.  If FFTW_EXHAUSTIVE is specified instead, the planner will further progress to FFTW_EXHAUSTIVE mode.

Note that the seconds argument specifies only a rough limit; in practice, the planner may use somewhat more time if the time limit is reached when the planner is in the middle of an operation that cannot be interrupted.  At the very least, the planner will complete planning in FFTW_ESTIMATE mode (which is thus equivalent to a time limit of 0).

Real-data DFTs
--------------

     fftw_plan fftw_plan_dft_r2c_1d(int n0,
                                    double *in, fftw_complex *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1,
                                    double *in, fftw_complex *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2,
                                    double *in, fftw_complex *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_r2c(int rank, const int *n,
                                 double *in, fftw_complex *out,
                                 unsigned flags);

Plan a real-input/complex-output discrete Fourier transform (DFT) in zero or more dimensions, returning an fftw_plan (see Using Plans).

Once you have created a plan for a certain transform type and parameters, then creating another plan of the same type and parameters, but for different arrays, is fast and shares constant data with the first plan (if it still exists).

The planner returns NULL if the plan cannot be created.  A non-NULL plan is always returned by the basic interface unless you are using a customized FFTW configuration supporting a restricted set of transforms, or if you use the FFTW_PRESERVE_INPUT flag with a multi-dimensional out-of-place c2r transform (see below).

Arguments
.........

  * rank is the rank of the transform (it should be the size of the array *n), and can be any non-negative integer.  (see Complex Multi-Dimensional DFTs, for the definition of ``rank''.)  The _1d, _2d, and _3d planners correspond to a rank of 1, 2, and 3, respectively.  The rank may be zero, which is equivalent to a rank-1 transform of size 1, i.e. a copy of one real number (with zero imaginary part) from input to output.

  * n0, n1, n2, or n[0..rank-1], (as appropriate for each routine) specify the size of the transform dimensions.  They can be any positive integer.  This is different in general from the physical array dimensions, which are described in Real-data DFT Array Format.

       - FFTW is best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, where e+f is either 0 or 1, and the other exponents are arbitrary.  Other sizes are computed by means of a slow, general-purpose algorithm (which nevertheless retains Onlogn performance even for prime sizes).  (It is possible to customize FFTW for different array sizes; see Installation and Customization.) Transforms whose sizes are powers of 2 are especially fast, and it is generally beneficial for the last dimension of an r2c/c2r transform to be even.

  * in and out point to the input and output arrays of the transform, which may be the same (yielding an in-place transform). These arrays are overwritten during planning, unless FFTW_ESTIMATE is used in the flags.  (The arrays need not be initialized, but they must be allocated.)  For an in-place transform, it is important to remember that the real array will require padding, described in Real-data DFT Array Format.

  * flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

The inverse transforms, taking complex input (storing the non-redundant half of a logically Hermitian array) to real output, are given by:

     fftw_plan fftw_plan_dft_c2r_1d(int n0,
                                    fftw_complex *in, double *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_c2r_2d(int n0, int n1,
                                    fftw_complex *in, double *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_c2r_3d(int n0, int n1, int n2,
                                    fftw_complex *in, double *out,
                                    unsigned flags);
     fftw_plan fftw_plan_dft_c2r(int rank, const int *n,
                                 fftw_complex *in, double *out,
                                 unsigned flags);

The arguments are the same as for the r2c transforms, except that the input and output data formats are reversed.

FFTW computes an unnormalized transform: computing an r2c followed by a c2r transform (or vice versa) will result in the original data multiplied by the size of the transform (the product of the logical dimensions). An r2c transform produces the same output as a FFTW_FORWARD complex DFT of the same input, and a c2r transform is correspondingly equivalent to FFTW_BACKWARD.  For more information, see What FFTW Really Computes.

Real-data DFT Array Format
--------------------------

The output of a DFT of real data (r2c) contains symmetries that, in principle, make half of the outputs redundant (see What FFTW Really Computes).  (Similarly for the input of an inverse c2r transform.)  In practice, it is not possible to entirely realize these savings in an efficient and understandable format that generalizes to multi-dimensional transforms.  Instead, the output of the r2c transforms is slightly over half of the output of the corresponding complex transform.  We do not ``pack'' the data in any way, but store it as an ordinary array of fftw_complex values. In fact, this data is simply a subsection of what would be the array in the corresponding complex transform.

Specifically, for a real transform of d (= rank) dimensions ndims, the complex data is an array of fftw_complex values in row-major order (with the division rounded down).  That is, we only store the lower half (non-negative frequencies), plus one element, of the last dimension of the data from the ordinary complex transform.  (We could have instead taken half of any other dimension, but implementation turns out to be simpler if the last, contiguous, dimension is used.)

For an out-of-place transform, the real data is simply an array with physical dimensions in row-major order.

For an in-place transform, some complications arise since the complex data is slightly larger than the real data.  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 * (n[d-1]/2+1) double values (exactly enough to hold the complex data).  This physical array size does not, however, change the logical array size---only n[d-1] values are actually stored in the last dimension, and n[d-1] is the last dimension passed to the planner.

Real-to-Real Transforms
-----------------------

     fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out,
                                fftw_r2r_kind kind, unsigned flags);
     fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out,
                                fftw_r2r_kind kind0, fftw_r2r_kind kind1,
                                unsigned flags);
     fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2,
                                double *in, double *out,
                                fftw_r2r_kind kind0,
                                fftw_r2r_kind kind1,
                                fftw_r2r_kind kind2,
                                unsigned flags);
     fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out,
                             const fftw_r2r_kind *kind, unsigned flags);

Plan a real input/output (r2r) transform of various kinds in zero or more dimensions, returning an fftw_plan (see Using Plans).

Once you have created a plan for a certain transform type and parameters, then creating another plan of the same type and parameters, but for different arrays, is fast and shares constant data with the first plan (if it still exists).

The planner returns NULL if the plan cannot be created.  A non-NULL plan is always returned by the basic interface unless you are using a customized FFTW configuration supporting a restricted set of transforms, or for size-1 FFTW_REDFT00 kinds (which are not defined).

Arguments
.........

  * rank is the dimensionality of the transform (it should be the size of the arrays *n and *kind), and can be any non-negative integer.  The _1d, _2d, and _3d planners correspond to a rank of 1, 2, and 3, respectively.  A rank of zero is equivalent to a copy of one number from input to output.

  * n, or n0/n1/n2, or n[rank], respectively, gives the (physical) size of the transform dimensions. They can be any positive integer.

       - Multi-dimensional arrays are stored in row-major order with dimensions: n0 x n1; or n0 x n1 x n2; or n[0] x n[1] x ... x n[rank-1]. See Multi-dimensional Array Format.

       - FFTW is generally best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, where e+f is either 0 or 1, and the other exponents are arbitrary.  Other sizes are computed by means of a slow, general-purpose algorithm (which nevertheless retains Onlogn performance even for prime sizes).  (It is possible to customize FFTW for different array sizes; see Installation and Customization.) Transforms whose sizes are powers of 2 are especially fast.

       - For a REDFT00 or RODFT00 transform kind in a dimension of size n, it is n-1 or n+1, respectively, that should be factorizable in the above form.

  * in and out point to the input and output arrays of the transform, which may be the same (yielding an in-place transform). These arrays are overwritten during planning, unless FFTW_ESTIMATE is used in the flags.  (The arrays need not be initialized, but they must be allocated.)

  * kind, or kind0/kind1/kind2, or kind[rank], is the kind of r2r transform used for the corresponding dimension.  The valid kind constants are described in Real-to-Real Transform Kinds.  In a multi-dimensional transform, what is computed is the separable product formed by taking each transform kind along the corresponding dimension, one dimension after another.

  * flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

Real-to-Real Transform Kinds
----------------------------

FFTW currently supports 11 different r2r transform kinds, specified by one of the constants below.  For the precise definitions of these transforms, see What FFTW Really Computes.  For a more colloquial introduction to these transform kinds, see More DFTs of Real Data.

For dimension of size n, there is a corresponding ``logical'' dimension N that determines the normalization (and the optimal factorization); the formula for N is given for each kind below. Also, with each transform kind is listed its corrsponding inverse transform.  FFTW computes unnormalized transforms: a transform followed by its inverse will result in the original data multiplied by N (or the product of the N's for each dimension, in multi-dimensions).

  * FFTW_R2HC computes a real-input DFT with output in ``halfcomplex'' format, i.e. real and imaginary parts for a transform of size n stored as: r0, r1, r2, r(n/2), i((n+1)/2-1), ..., i2, i1 (Logical N=n, inverse is FFTW_HC2R.)

  * FFTW_HC2R computes the reverse of FFTW_R2HC, above. (Logical N=n, inverse is FFTW_R2HC.)

  * FFTW_DHT computes a discrete Hartley transform. (Logical N=n, inverse is FFTW_DHT.)

  * FFTW_REDFT00 computes an REDFT00 transform, i.e. a DCT-I. (Logical N=2*(n-1), inverse is FFTW_REDFT00.)

  * FFTW_REDFT10 computes an REDFT10 transform, i.e. a DCT-II (sometimes called ``the'' DCT). (Logical N=2*n, inverse is FFTW_REDFT01.)

  * FFTW_REDFT01 computes an REDFT01 transform, i.e. a DCT-III (sometimes called ``the'' IDCT, being the inverse of DCT-II). (Logical N=2*n, inverse is FFTW_REDFT=10.)

  * FFTW_REDFT11 computes an REDFT11 transform, i.e. a DCT-IV. (Logical N=2*n, inverse is FFTW_REDFT11.)

  * FFTW_RODFT00 computes an RODFT00 transform, i.e. a DST-I. (Logical N=2*(n+1), inverse is FFTW_RODFT00.)

  * FFTW_RODFT10 computes an RODFT10 transform, i.e. a DST-II. (Logical N=2*n, inverse is FFTW_RODFT01.)

  * FFTW_RODFT01 computes an RODFT01 transform, i.e. a DST-III. (Logical N=2*n, inverse is FFTW_RODFT=10.)

  * FFTW_RODFT11 computes an RODFT11 transform, i.e. a DST-IV. (Logical N=2*n, inverse is FFTW_RODFT11.)

Advanced Interface
==================

FFTW's ``advanced'' interface supplements the basic interface with four new planner routines, providing a new level of flexibility: you can plan a transform of multiple arrays simultaneously, operate on non-contiguous (strided) data, and transform a subset of a larger multi-dimensional array.  Other than these additional features, the planner operates in the same fashion as in the basic interface, and the resulting fftw_plan is used in the same way (see Using Plans).

Advanced Complex DFTs
---------------------

     fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
                                  fftw_complex *in, const int *inembed,
                                  int istride, int idist,
                                  fftw_complex *out, const int *onembed,
                                  int ostride, int odist,
                                  int sign, unsigned flags);

This routine plans multiple multidimensional complex DFTs, and it extends the fftw_plan_dft routine (see Complex DFTs) to compute howmany transforms, each having rank rank and size n.  In addition, the transform data need not be contiguous, but it may be laid out in memory with an arbitrary stride.  To account for these possibilities, fftw_plan_many_dft adds the new parameters howmany, {i,o}nembed, {i,o}stride, and {i,o}dist.  The FFTW basic interface (see Complex DFTs) provides routines specialized for ranks 1, 2, and 3, but the advanced interface handles only the general-rank case.

howmany is the (nonnegative) number of transforms to compute.  The resulting plan computes howmany transforms, where the input of the k-th transform is at location in+k*idist (in C pointer arithmetic), and its output is at location out+k*odist.  Plans obtained in this way can often be faster than calling FFTW multiple times for the individual transforms.  The basic fftw_plan_dft interface corresponds to howmany=1 (in which case the dist parameters are ignored).

Each of the howmany transforms has rank rank and size n, as in the basic interface.  In addition, the advanced interface allows the input and output arrays of each transform to be row-major subarrays of larger rank-rank arrays, described by inembed and onembed parameters, respectively. {i,o}nembed must be arrays of length rank, and n should be elementwise less than or equal to {i,o}nembed.  Passing NULL for an nembed parameter is equivalent to passing n (i.e. same physical and logical dimensions, as in the basic interface.)

The stride parameters indicate that the j-th element of the input or output arrays is located at j*istride or j*ostride, respectively.  (For a multi-dimensional array, j is the ordinary row-major index.)  When combined with the k-th transform in a howmany loop, from above, this means that the (j,k)-th element is at j*stride+k*dist. (The basic fftw_plan_dft interface corresponds to a stride of 1.)

For in-place transforms, the input and output stride and dist parameters should be the same; otherwise, the planner may return NULL.

Arrays n, inembed, and onembed are not used after this function returns.  You can safely free or reuse them.

Examples: One transform of one 5 by 6 array contiguous in memory:

        int rank = 2;
        int n[] = {5, 6};
        int howmany = 1;
        int idist = odist = 0; /* unused because howmany = 1 */
        int istride = ostride = 1; /* array is contiguous in memory */
        int *inembed = n, *onembed = n;

Transform of three 5 by 6 arrays, each contiguous in memory, stored in memory one after another:

        int rank = 2;
        int n[] = {5, 6};
        int howmany = 3;
        int idist = odist = n[0]*n[1]; /* = 30, the distance in memory
                                          between the first element
                                          of the first array and the
                                          first element of the second array */
        int istride = ostride = 1; /* array is contiguous in memory */
        int *inembed = n, *onembed = n;

Transform each column of a 2d array with 10 rows and 3 columns:

        int rank = 1; /* not 2: we are computing 1d transforms */
        int n[] = {10}; /* 1d transforms of length 10 */
        int howmany = 3;
        int idist = odist = 1;
        int istride = ostride = 3; /* distance between two elements in
                                      the same column */
        int *inembed = n, *onembed = n;

Advanced Real-data DFTs
-----------------------

     fftw_plan fftw_plan_many_dft_r2c(int rank, const int *n, int howmany,
                                      double *in, const int *inembed,
                                      int istride, int idist,
                                      fftw_complex *out, const int *onembed,
                                      int ostride, int odist,
                                      unsigned flags);
     fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany,
                                      fftw_complex *in, const int *inembed,
                                      int istride, int idist,
                                      double *out, const int *onembed,
                                      int ostride, int odist,
                                      unsigned flags);

Like fftw_plan_many_dft, these two functions add howmany, nembed, stride, and dist parameters to the fftw_plan_dft_r2c and fftw_plan_dft_c2r functions, but otherwise behave the same as the basic interface.

The interpretation of howmany, stride, and dist are the same as for fftw_plan_many_dft, above.  Note that the stride and dist for the real array are in units of double, and for the complex array are in units of fftw_complex.

If an nembed parameter is NULL, it is interpreted as what it would be in the basic interface, as described in Real-data DFT Array Format.  That is, for the complex array the size is assumed to be the same as n, but with the last dimension cut roughly in half. For the real array, the size is assumed to be n if the transform is out-of-place, or n with the last dimension ``padded'' if the transform is in-place.

If an nembed parameter is non-NULL, it is interpreted as the physical size of the corresponding array, in row-major order, just as for fftw_plan_many_dft.  In this case, each dimension of nembed should be >= what it would be in the basic interface (e.g. the halved or padded n).

Arrays n, inembed, and onembed are not used after this function returns.  You can safely free or reuse them.

Advanced Real-to-real Transforms
--------------------------------

     fftw_plan fftw_plan_many_r2r(int rank, const int *n, int howmany,
                                  double *in, const int *inembed,
                                  int istride, int idist,
                                  double *out, const int *onembed,
                                  int ostride, int odist,
                                  const fftw_r2r_kind *kind, unsigned flags);

Like fftw_plan_many_dft, this functions adds howmany, nembed, stride, and dist parameters to the fftw_plan_r2r function, but otherwise behave the same as the basic interface.  The interpretation of those additional parameters are the same as for fftw_plan_many_dft.  (Of course, the stride and dist parameters are now in units of double, not fftw_complex.)

Arrays n, inembed, onembed, and kind are not used after this function returns.  You can safely free or reuse them.

Guru Interface
==============

The ``guru'' interface to FFTW is intended to expose as much as possible of the flexibility in the underlying FFTW architecture.  It allows one to compute multi-dimensional ``vectors'' (loops) of multi-dimensional transforms, where each vector/transform dimension has an independent size and stride. One can also use more general complex-number formats, e.g. separate real and imaginary arrays.

For those users who require the flexibility of the guru interface, it is important that they pay special attention to the documentation lest they shoot themselves in the foot.

Interleaved and split arrays
----------------------------

The guru interface supports two representations of complex numbers, which we call the interleaved and the split format.

The interleaved format is the same one used by the basic and advanced interfaces, and it is documented in Complex numbers. In the interleaved format, you provide pointers to the real part of a complex number, and the imaginary part understood to be stored in the next memory location.

The split format allows separate pointers to the real and imaginary parts of a complex array.

Technically, the interleaved format is redundant, because you can always express an interleaved array in terms of a split array with appropriate pointers and strides.  On the other hand, the interleaved format is simpler to use, and it is common in practice.  Hence, FFTW supports it as a special case.

Guru vector and transform sizes
-------------------------------

The guru interface introduces one basic new data structure, fftw_iodim, that is used to specify sizes and strides for multi-dimensional transforms and vectors:

     typedef struct {
          int n;
          int is;
          int os;
     } fftw_iodim;

Here, n is the size of the dimension, and is and os are the strides of that dimension for the input and output arrays.  (The stride is the separation of consecutive elements along this dimension.)

The meaning of the stride parameter depends on the type of the array that the stride refers to.  If the array is interleaved complex, strides are expressed in units of complex numbers (fftw_complex).  If the array is split complex or real, strides are expressed in units of real numbers (double).  This convention is consistent with the usual pointer arithmetic in the C language.  An interleaved array is denoted by a pointer p to fftw_complex, so that p+1 points to the next complex number.  Split arrays are denoted by pointers to double, in which case pointer arithmetic operates in units of sizeof(double).

The guru planner interfaces all take a (rank, dims[rank]) pair describing the transform size, and a (howmany_rank, howmany_dims[howmany_rank]) pair describing the ``vector'' size (a multi-dimensional loop of transforms to perform), where dims and howmany_dims are arrays of fftw_iodim.  Each n field must be positive for dims and nonnegative for howmany_dims, while both rank and howmany_rank must be nonnegative.

For example, the howmany parameter in the advanced complex-DFT interface corresponds to howmany_rank = 1, howmany_dims[0].n = howmany, howmany_dims[0].is = idist, and howmany_dims[0].os = odist. (To compute a single transform, you can just use howmany_rank = 0.)

A row-major multidimensional array with dimensions n[rank] (see Row-major Format) corresponds to dims[i].n = n[i] and the recurrence dims[i].is = n[i+1] * dims[i+1].is (similarly for os).  The stride of the last (i=rank-1) dimension is the overall stride of the array. e.g. to be equivalent to the advanced complex-DFT interface, you would have dims[rank-1].is = istride and dims[rank-1].os = ostride.

In general, we only guarantee FFTW to return a non-NULL plan if the vector and transform dimensions correspond to a set of distinct indices, and for in-place transforms the input/output strides should be the same.

Guru Complex DFTs
-----------------

     fftw_plan fftw_plan_guru_dft(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          fftw_complex *in, fftw_complex *out,
          int sign, unsigned flags);

     fftw_plan fftw_plan_guru_split_dft(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          double *ri, double *ii, double *ro, double *io,
          unsigned flags);

These two functions plan a complex-data, multi-dimensional DFT for the interleaved and split format, respectively. Transform dimensions are given by (rank, dims) over a multi-dimensional vector (loop) of dimensions (howmany_rank, howmany_dims).  dims and howmany_dims should point to fftw_iodim arrays of length rank and howmany_rank, respectively.

flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

In the fftw_plan_guru_dft function, the pointers in and out point to the interleaved input and output arrays, respectively.  The sign can be either -1 (= FFTW_FORWARD) or +1 (= FFTW_BACKWARD).  If the pointers are equal, the transform is in-place.

In the fftw_plan_guru_split_dft function, ri and ii point to the real and imaginary input arrays, and ro and io point to the real and imaginary output arrays.  The input and output pointers may be the same, indicating an in-place transform.  For example, for fftw_complex pointers in and out, the corresponding parameters are:

     ri = (double *) in;
     ii = (double *) in + 1;
     ro = (double *) out;
     io = (double *) out + 1;

Because fftw_plan_guru_split_dft accepts split arrays, strides are expressed in units of double.  For a contiguous fftw_complex array, the overall stride of the transform should be 2, the distance between consecutive real parts or between consecutive imaginary parts; see Guru vector and transform sizes.  Note that the dimension strides are applied equally to the real and imaginary parts; real and imaginary arrays with different strides are not supported.

There is no sign parameter in fftw_plan_guru_split_dft. This function always plans for an FFTW_FORWARD transform.  To plan for an FFTW_BACKWARD transform, you can exploit the identity that the backwards DFT is equal to the forwards DFT with the real and imaginary parts swapped.  For example, in the case of the fftw_complex arrays above, the FFTW_BACKWARD transform is computed by the parameters:

     ri = (double *) in + 1;
     ii = (double *) in;
     ro = (double *) out + 1;
     io = (double *) out;

Guru Real-data DFTs
-------------------

     fftw_plan fftw_plan_guru_dft_r2c(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          double *in, fftw_complex *out,
          unsigned flags);

     fftw_plan fftw_plan_guru_split_dft_r2c(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          double *in, double *ro, double *io,
          unsigned flags);

     fftw_plan fftw_plan_guru_dft_c2r(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          fftw_complex *in, double *out,
          unsigned flags);

     fftw_plan fftw_plan_guru_split_dft_c2r(
          int rank, const fftw_iodim *dims,
          int howmany_rank, const fftw_iodim *howmany_dims,
          double *ri, double *ii, double *out,
          unsigned flags);

Plan a real-input (r2c) or real-output (c2r), multi-dimensional DFT with transform dimensions given by (rank, dims) over a multi-dimensional vector (loop) of dimensions (howmany_rank, howmany_dims).  dims and howmany_dims should point to fftw_iodim arrays of length rank and howmany_rank, respectively.  As for the basic and advanced interfaces, an r2c transform is FFTW_FORWARD and a c2r transform is FFTW_BACKWARD.

The last dimension of dims is interpreted specially: that dimension of the real array has size dims[rank-1].n, but that dimension of the complex array has size dims[rank-1].n/2+1 (division rounded down).  The strides, on the other hand, are taken to be exactly as specified.  It is up to the user to specify the strides appropriately for the peculiar dimensions of the data, and we do not guarantee that the planner will succeed (return non-NULL) for any dimensions other than those described in Real-data DFT Array Format and generalized in Advanced Real-data DFTs.  (That is, for an in-place transform, each individual dimension should be able to operate in place.)

in and out point to the input and output arrays for r2c and c2r transforms, respectively.  For split arrays, ri and ii point to the real and imaginary input arrays for a c2r transform, and ro and io point to the real and imaginary output arrays for an r2c transform.  in and ro or ri and out may be the same, indicating an in-place transform.   (In-place transforms where in and io or ii and out are the same are not currently supported.)

flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

In-place transforms of rank greater than 1 are currently only supported for interleaved arrays.  For split arrays, the planner will return NULL.

Guru Real-to-real Transforms
----------------------------

     fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims,
                                  int howmany_rank,
                                  const fftw_iodim *howmany_dims,
                                  double *in, double *out,
                                  const fftw_r2r_kind *kind,
                                  unsigned flags);

Plan a real-to-real (r2r) multi-dimensional FFTW_FORWARD transform with transform dimensions given by (rank, dims) over a multi-dimensional vector (loop) of dimensions (howmany_rank, howmany_dims).  dims and howmany_dims should point to fftw_iodim arrays of length rank and howmany_rank, respectively.

The transform kind of each dimension is given by the kind parameter, which should point to an array of length rank.  Valid fftw_r2r_kind constants are given in Real-to-Real Transform Kinds.

in and out point to the real input and output arrays; they may be the same, indicating an in-place transform.

flags is a bitwise OR (|) of zero or more planner flags, as defined in Planner Flags.

64-bit Guru Interface
---------------------

When compiled in 64-bit mode on a 64-bit architecture (where addresses are 64 bits wide), FFTW uses 64-bit quantities internally for all transform sizes, strides, and so on---you don't have to do anything special to exploit this.  However, in the ordinary FFTW interfaces, you specify the transform size by an int quantity, which is normally only 32 bits wide.  This means that, even though FFTW is using 64-bit sizes internally, you cannot specify a single transform dimension larger than 2^31-1 numbers.

We expect that few users will require transforms larger than this, but, for those who do, we provide a 64-bit version of the guru interface in which all sizes are specified as integers of type ptrdiff_t instead of int.  (ptrdiff_t is a signed integer type defined by the C standard to be wide enough to represent address differences, and thus must be at least 64 bits wide on a 64-bit machine.)  We stress that there is no performance advantage to using this interface---the same internal FFTW code is employed regardless---and it is only necessary if you want to specify very large transform sizes.

In particular, the 64-bit guru interface is a set of planner routines that are exactly the same as the guru planner routines, except that they are named with guru64 instead of guru and they take arguments of type fftw_iodim64 instead of fftw_iodim. For example, instead of fftw_plan_guru_dft, we have fftw_plan_guru64_dft.

     fftw_plan fftw_plan_guru64_dft(
          int rank, const fftw_iodim64 *dims,
          int howmany_rank, const fftw_iodim64 *howmany_dims,
          fftw_complex *in, fftw_complex *out,
          int sign, unsigned flags);

The fftw_iodim64 type is similar to fftw_iodim, with the same interpretation, except that it uses type ptrdiff_t instead of type int.

     typedef struct {
          ptrdiff_t n;
          ptrdiff_t is;
          ptrdiff_t os;
     } fftw_iodim64;

Every other fftw_plan_guru function also has a fftw_plan_guru64 equivalent, but we do not repeat their documentation here since they are identical to the 32-bit versions except as noted above.

New-array Execute Functions
===========================

Normally, one executes a plan for the arrays with which the plan was created, by calling fftw_execute(plan) as described in Using Plans. However, it is possible for sophisticated users to apply a given plan to a different array using the ``new-array execute'' functions detailed below, provided that the following conditions are met:

  * The array size, strides, etcetera are the same (since those are set by the plan).

  * The input and output arrays are the same (in-place) or different (out-of-place) if the plan was originally created to be in-place or out-of-place, respectively.

  * For split arrays, the separations between the real and imaginary parts, ii-ri and io-ro, are the same as they were for the input and output arrays when the plan was created.  (This condition is automatically satisfied for interleaved arrays.)

  * The alignment of the new input/output arrays is the same as that of the input/output arrays when the plan was created, unless the plan was created with the FFTW_UNALIGNED flag. Here, the alignment is a platform-dependent quantity (for example, it is the address modulo 16 if SSE SIMD instructions are used, but the address modulo 4 for non-SIMD single-precision FFTW on the same machine).  In general, only arrays allocated with fftw_malloc are guaranteed to be equally aligned (see SIMD alignment and fftw_malloc).

The alignment issue is especially critical, because if you don't use fftw_malloc then you may have little control over the alignment of arrays in memory.  For example, neither the C++ new function nor the Fortran allocate statement provide strong enough guarantees about data alignment.  If you don't use fftw_malloc, therefore, you probably have to use FFTW_UNALIGNED (which disables most SIMD support).  If possible, it is probably better for you to simply create multiple plans (creating a new plan is quick once one exists for a given size), or better yet re-use the same array for your transforms.

For rare circumstances in which you cannot control the alignment of allocated memory, but wish to determine where a given array is aligned like the original array for which a plan was created, you can use the fftw_alignment_of function:

     int fftw_alignment_of(double *p);

Two arrays have equivalent alignment (for the purposes of applying a plan) if and only if fftw_alignment_of returns the same value for the corresponding pointers to their data (typecast to double* if necessary).

If you are tempted to use the new-array execute interface because you want to transform a known bunch of arrays of the same size, you should probably go use the advanced interface instead (see Advanced Interface)).

The new-array execute functions are:

     void fftw_execute_dft(
          const fftw_plan p,
          fftw_complex *in, fftw_complex *out);

     void fftw_execute_split_dft(
          const fftw_plan p,
          double *ri, double *ii, double *ro, double *io);

     void fftw_execute_dft_r2c(
          const fftw_plan p,
          double *in, fftw_complex *out);

     void fftw_execute_split_dft_r2c(
          const fftw_plan p,
          double *in, double *ro, double *io);

     void fftw_execute_dft_c2r(
          const fftw_plan p,
          fftw_complex *in, double *out);

     void fftw_execute_split_dft_c2r(
          const fftw_plan p,
          double *ri, double *ii, double *out);

     void fftw_execute_r2r(
          const fftw_plan p,
          double *in, double *out);

These execute the plan to compute the corresponding transform on the input/output arrays specified by the subsequent arguments.  The input/output array arguments have the same meanings as the ones passed to the guru planner routines in the preceding sections.  The plan is not modified, and these routines can be called as many times as desired, or intermixed with calls to the ordinary fftw_execute.

The plan must have been created for the transform type corresponding to the execute function, e.g. it must be a complex-DFT plan for fftw_execute_dft.  Any of the planner routines for that transform type, from the basic to the guru interface, could have been used to create the plan, however.

Wisdom
======

This section documents the FFTW mechanism for saving and restoring plans from disk.  This mechanism is called wisdom.

Wisdom Export
-------------

     int fftw_export_wisdom_to_filename(const char *filename);
     void fftw_export_wisdom_to_file(FILE *output_file);
     char *fftw_export_wisdom_to_string(void);
     void fftw_export_wisdom(void (*write_char)(char c, void *), void *data);

These functions allow you to export all currently accumulated wisdom in a form from which it can be later imported and restored, even during a separate run of the program. (see Words of Wisdom-Saving Plans.)  The current store of wisdom is not affected by calling any of these routines.

fftw_export_wisdom exports the wisdom to any output medium, as specified by the callback function write_char. write_char is a putc-like function that writes the character c to some output; its second parameter is the data pointer passed to fftw_export_wisdom.  For convenience, the following three ``wrapper'' routines are provided:

fftw_export_wisdom_to_filename writes wisdom to a file named filename (which is created or overwritten), returning 1 on success and 0 on failure.  A lower-level function, which requires you to open and close the file yourself (e.g. if you want to write wisdom to a portion of a larger file) is fftw_export_wisdom_to_file.  This writes the wisdom to the current position in output_file, which should be open with write permission; upon exit, the file remains open and is positioned at the end of the wisdom data.

fftw_export_wisdom_to_string returns a pointer to a NULL-terminated string holding the wisdom data. This string is dynamically allocated, and it is the responsibility of the caller to deallocate it with free when it is no longer needed.

All of these routines export the wisdom in the same format, which we will not document here except to say that it is LISP-like ASCII text that is insensitive to white space.

Wisdom Import
-------------

     int fftw_import_system_wisdom(void);
     int fftw_import_wisdom_from_filename(const char *filename);
     int fftw_import_wisdom_from_string(const char *input_string);
     int fftw_import_wisdom(int (*read_char)(void *), void *data);

These functions import wisdom into a program from data stored by the fftw_export_wisdom functions above. (see Words of Wisdom-Saving Plans.)  The imported wisdom replaces any wisdom already accumulated by the running program.

fftw_import_wisdom imports wisdom from any input medium, as specified by the callback function read_char. read_char is a getc-like function that returns the next character in the input; its parameter is the data pointer passed to fftw_import_wisdom. If the end of the input data is reached (which should never happen for valid data), read_char should return EOF (as defined in <stdio.h>).  For convenience, the following three ``wrapper'' routines are provided:

fftw_import_wisdom_from_filename reads wisdom from a file named filename.  A lower-level function, which requires you to open and close the file yourself (e.g. if you want to read wisdom from a portion of a larger file) is fftw_import_wisdom_from_file. This reads wisdom from the current position in input_file (which should be open with read permission); upon exit, the file remains open, but the position of the read pointer is unspecified.

fftw_import_wisdom_from_string reads wisdom from the NULL-terminated string input_string.

fftw_import_system_wisdom reads wisdom from an implementation-defined standard file (/etc/fftw/wisdom on Unix and GNU systems).

The return value of these import routines is 1 if the wisdom was read successfully and 0 otherwise. Note that, in all of these functions, any data in the input stream past the end of the wisdom data is simply ignored.

Forgetting Wisdom
-----------------

     void fftw_forget_wisdom(void);

Calling fftw_forget_wisdom causes all accumulated wisdom to be discarded and its associated memory to be freed. (New wisdom can still be gathered subsequently, however.)

Wisdom Utilities
----------------

FFTW includes two standalone utility programs that deal with wisdom.  We merely summarize them here, since they come with their own man pages for Unix and GNU systems (with HTML versions on our web site).

The first program is fftw-wisdom (or fftwf-wisdom in single precision, etcetera), which can be used to create a wisdom file containing plans for any of the transform sizes and types supported by FFTW.  It is preferable to create wisdom directly from your executable (see Caveats in Using Wisdom), but this program is useful for creating global wisdom files for fftw_import_system_wisdom.

The second program is fftw-wisdom-to-conf, which takes a wisdom file as input and produces a configuration routine as output.  The latter is a C subroutine that you can compile and link into your program, replacing a routine of the same name in the FFTW library, that determines which parts of FFTW are callable by your program. fftw-wisdom-to-conf produces a configuration routine that links to only those parts of FFTW needed by the saved plans in the wisdom, greatly reducing the size of statically linked executables (which should only attempt to create plans corresponding to those in the wisdom, however).

What FFTW Really Computes
=========================

In this section, we provide precise mathematical definitions for the transforms that FFTW computes.  These transform definitions are fairly standard, but some authors follow slightly different conventions for the normalization of the transform (the constant factor in front) and the sign of the complex exponent.  We begin by presenting the one-dimensional (1d) transform definitions, and then give the straightforward extension to multi-dimensional transforms.

The 1d Discrete Fourier Transform (DFT)
---------------------------------------

The forward (FFTW_FORWARD) discrete Fourier transform (DFT) of a 1d complex array X of size n computes an array Y, where:
Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(-2 pi j k sqrt(-1)/n) .
The backward (FFTW_BACKWARD) DFT computes:
Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(2 pi j k sqrt(-1)/n) .

FFTW computes an unnormalized transform, in that there is no coefficient in front of the summation in the DFT.  In other words, applying the forward and then the backward transform will multiply the input by n.

From above, an FFTW_FORWARD transform corresponds to a sign of -1 in the exponent of the DFT.  Note also that we use the standard ``in-order'' output ordering---the k-th output corresponds to the frequency k/n (or k/T, where T is your total sampling period).  For those who like to think in terms of positive and negative frequencies, this means that the positive frequencies are stored in the first half of the output and the negative frequencies are stored in backwards order in the second half of the output.  (The frequency -k/n is the same as the frequency (n-k)/n.)

The 1d Real-data DFT
--------------------

The real-input (r2c) DFT in FFTW computes the forward transform Y of the size n real array X, exactly as defined above, i.e.
Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(-2 pi j k sqrt(-1)/n) .
This output array Y can easily be shown to possess the ``Hermitian'' symmetry Y[k] = Y[n-k]*, where we take Y to be periodic so that Y[n] = Y[0].

As a result of this symmetry, half of the output Y is redundant (being the complex conjugate of the other half), and so the 1d r2c transforms only output elements 0...n/2 of Y (n/2+1 complex numbers), where the division by 2 is rounded down.

Moreover, the Hermitian symmetry implies that Y[0] and, if n is even, the Y[n/2] element, are purely real.  So, for the R2HC r2r transform, the halfcomplex format does not store the imaginary parts of these elements.

The c2r and H2RC r2r transforms compute the backward DFT of the complex array X with Hermitian symmetry, stored in the r2c/R2HC output formats, respectively, where the backward transform is defined exactly as for the complex case:
Y[k] = sum for j = 0 to (n - 1) of X[j] * exp(2 pi j k sqrt(-1)/n) .
The outputs Y of this transform can easily be seen to be purely real, and are stored as an array of real numbers.

Like FFTW's complex DFT, these transforms are unnormalized.  In other words, applying the real-to-complex (forward) and then the complex-to-real (backward) transform will multiply the input by n.

1d Real-even DFTs (DCTs)
------------------------

The Real-even symmetry DFTs in FFTW are exactly equivalent to the unnormalized forward (and backward) DFTs as defined above, where the input array X of length N is purely real and is also even symmetry.  In this case, the output array is likewise real and even symmetry.

For the case of REDFT00, this even symmetry means that X[j] = X[N-j], where we take X to be periodic so that X[N] = X[0]. Because of this redundancy, only the first n real numbers are actually stored, where N = 2(n-1).

The proper definition of even symmetry for REDFT10, REDFT01, and REDFT11 transforms is somewhat more intricate because of the shifts by 1/2 of the input and/or output, although the corresponding boundary conditions are given in Real even/odd DFTs (cosine/sine transforms).  Because of the even symmetry, however, the sine terms in the DFT all cancel and the remaining cosine terms are written explicitly below.  This formulation often leads people to call such a transform a discrete cosine transform (DCT), although it is really just a special case of the DFT.

In each of the definitions below, we transform a real array X of length n to a real array Y of length n:

REDFT00 (DCT-I)
...............

An REDFT00 transform (type-I DCT) in FFTW is defined by: Y[k] = X[0] + (-1)^k X[n-1] + 2 (sum for j = 1 to n-2 of X[j] cos(pi jk /(n-1))). Note that this transform is not defined for n=1.  For n=2, the summation term above is dropped as you might expect.

REDFT10 (DCT-II)
................

An REDFT10 transform (type-II DCT, sometimes called ``the'' DCT) in FFTW is defined by: Y[k] = 2 (sum for j = 0 to n-1 of X[j] cos(pi (j+1/2) k / n)).

REDFT01 (DCT-III)
.................

An REDFT01 transform (type-III DCT) in FFTW is defined by: Y[k] = X[0] + 2 (sum for j = 1 to n-1 of X[j] cos(pi j (k+1/2) / n)). In the case of n=1, this reduces to Y[0] = X[0]. Up to a scale factor (see below), this is the inverse of REDFT10 (``the'' DCT), and so the REDFT01 (DCT-III) is sometimes called the ``IDCT''.

REDFT11 (DCT-IV)
................

An REDFT11 transform (type-IV DCT) in FFTW is defined by: Y[k] = 2 (sum for j = 0 to n-1 of X[j] cos(pi (j+1/2) (k+1/2) / n)).

Inverses and Normalization
..........................

These definitions correspond directly to the unnormalized DFTs used elsewhere in FFTW (hence the factors of 2 in front of the summations).  The unnormalized inverse of REDFT00 is REDFT00, of REDFT10 is REDFT01 and vice versa, and of REDFT11 is REDFT11.  Each unnormalized inverse results in the original array multiplied by N, where N is the logical DFT size.  For REDFT00, N=2(n-1) (note that n=1 is not defined); otherwise, N=2n.

In defining the discrete cosine transform, some authors also include additional factors of sqrt(2) (or its inverse) multiplying selected inputs and/or outputs.  This is a mostly cosmetic change that makes the transform orthogonal, but sacrifices the direct equivalence to a symmetric DFT.

1d Real-odd DFTs (DSTs)
-----------------------

The Real-odd symmetry DFTs in FFTW are exactly equivalent to the unnormalized forward (and backward) DFTs as defined above, where the input array X of length N is purely real and is also odd symmetry.  In this case, the output is odd symmetry and purely imaginary.

For the case of RODFT00, this odd symmetry means that X[j] = -X[N-j], where we take X to be periodic so that X[N] = X[0]. Because of this redundancy, only the first n real numbers starting at j=1 are actually stored (the j=0 element is zero), where N = 2(n+1).

The proper definition of odd symmetry for RODFT10, RODFT01, and RODFT11 transforms is somewhat more intricate because of the shifts by 1/2 of the input and/or output, although the corresponding boundary conditions are given in Real even/odd DFTs (cosine/sine transforms).  Because of the odd symmetry, however, the cosine terms in the DFT all cancel and the remaining sine terms are written explicitly below.  This formulation often leads people to call such a transform a discrete sine transform (DST), although it is really just a special case of the DFT.

In each of the definitions below, we transform a real array X of length n to a real array Y of length n:

RODFT00 (DST-I)
...............

An RODFT00 transform (type-I DST) in FFTW is defined by: Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1)(k+1) / (n+1))).

RODFT10 (DST-II)
................

An RODFT10 transform (type-II DST) in FFTW is defined by: Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1/2) (k+1) / n)).

RODFT01 (DST-III)
.................

An RODFT01 transform (type-III DST) in FFTW is defined by: Y[k] = (-1)^k X[n-1] + 2 (sum for j = 0 to n-2 of X[j] sin(pi (j+1) (k+1/2) / n)). In the case of n=1, this reduces to Y[0] = X[0].

RODFT11 (DST-IV)
................

An RODFT11 transform (type-IV DST) in FFTW is defined by: Y[k] = 2 (sum for j = 0 to n-1 of X[j] sin(pi (j+1/2) (k+1/2) / n)).

Inverses and Normalization
..........................

These definitions correspond directly to the unnormalized DFTs used elsewhere in FFTW (hence the factors of 2 in front of the summations).  The unnormalized inverse of RODFT00 is RODFT00, of RODFT10 is RODFT01 and vice versa, and of RODFT11 is RODFT11.  Each unnormalized inverse results in the original array multiplied by N, where N is the logical DFT size.  For RODFT00, N=2(n+1); otherwise, N=2n.

In defining the discrete sine transform, some authors also include additional factors of sqrt(2) (or its inverse) multiplying selected inputs and/or outputs.  This is a mostly cosmetic change that makes the transform orthogonal, but sacrifices the direct equivalence to an antisymmetric DFT.

1d Discrete Hartley Transforms (DHTs)
-------------------------------------

The discrete Hartley transform (DHT) of a 1d real array X of size n computes a real array Y of the same size, where:
Y[k] = sum for j = 0 to (n - 1) of X[j] * [cos(2 pi j k / n) + sin(2 pi j k / n)].

FFTW computes an unnormalized transform, in that there is no coefficient in front of the summation in the DHT.  In other words, applying the transform twice (the DHT is its own inverse) will multiply the input by n.

Multi-dimensional Transforms
----------------------------

The multi-dimensional transforms of FFTW, in general, compute simply the separable product of the given 1d transform along each dimension of the array.  Since each of these transforms is unnormalized, computing the forward followed by the backward/inverse multi-dimensional transform will result in the original array scaled by the product of the normalization factors for each dimension (e.g. the product of the dimension sizes, for a multi-dimensional DFT).

The definition of FFTW's multi-dimensional DFT of real data (r2c) deserves special attention.  In this case, we logically compute the full multi-dimensional DFT of the input data; since the input data are purely real, the output data have the Hermitian symmetry and therefore only one non-redundant half need be stored.  More specifically, for an multi-dimensional real-input DFT, the full (logical) complex output array Y[k[0], k[1], ..., k[d-1]] has the symmetry: Y[k[0], k[1], ..., k[d-1]] = Y[n[0] - k[0], n[1] - k[1], ..., n[d-1] - k[d-1]]* (where each dimension is periodic).  Because of this symmetry, we only store the k[d-1] = 0...n[d-1]/2 elements of the last dimension (division by 2 is rounded down).  (We could instead have cut any other dimension in half, but the last dimension proved computationally convenient.)  This results in the peculiar array format described in more detail by Real-data DFT Array Format.

The multi-dimensional c2r transform is simply the unnormalized inverse of the r2c transform.  i.e. it is the same as FFTW's complex backward multi-dimensional DFT, operating on a Hermitian input array in the peculiar format mentioned above and outputting a real array (since the DFT output is purely real).

We should remind the user that the separable product of 1d transforms along each dimension, as computed by FFTW, is not always the same thing as the usual multi-dimensional transform.  A multi-dimensional R2HC (or HC2R) transform is not identical to the multi-dimensional DFT, requiring some post-processing to combine the requisite real and imaginary parts, as was described in The Halfcomplex-format DFT.  Likewise, FFTW's multidimensional FFTW_DHT r2r transform is not the same thing as the logical multi-dimensional discrete Hartley transform defined in the literature, as discussed in The Discrete Hartley Transform.