fftw3: Complex One-Dimensional DFTs

 
 2.1 Complex One-Dimensional DFTs
 ================================
 
      Plan: To bother about the best method of accomplishing an
      accidental result.  [Ambrose Bierce, 'The Enlarged Devil's
      Dictionary'.]
 
    The basic usage of FFTW to compute a one-dimensional DFT of size 'N'
 is simple, and it typically looks something like this code:
 
      #include <fftw3.h>
      ...
      {
          fftw_complex *in, *out;
          fftw_plan p;
          ...
          in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
          out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
          p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
          ...
          fftw_execute(p); /* repeat as needed */
          ...
          fftw_destroy_plan(p);
          fftw_free(in); fftw_free(out);
      }
 
    You must link this code with the 'fftw3' library.  On Unix systems,
 link with '-lfftw3 -lm'.
 
    The example code first allocates the input and output arrays.  You
 can allocate them in any way that you like, but we recommend using
 'fftw_malloc', which behaves like 'malloc' except that it properly
 aligns the array when SIMD instructions (such as SSE and Altivec) are
 available (SeeSIMD alignment and fftw_malloc).  [Alternatively, we
 provide a convenient wrapper function 'fftw_alloc_complex(N)' which has
 the same effect.]
 
    The data is an array of type 'fftw_complex', which is by default a
 'double[2]' composed of the real ('in[i][0]') and imaginary ('in[i][1]')
 parts of a complex number.
 
    The next step is to create a "plan", which is an object that contains
 all the data that FFTW needs to compute the FFT. This function creates
 the plan:
 
      fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out,
                                 int sign, unsigned 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 (although
 prime sizes still use an O(n log n) algorithm).
 
    The next two arguments are pointers to the input and output arrays of
 the transform.  These pointers can be equal, indicating an "in-place"
 transform.
 
    The fourth argument, 'sign', can be either 'FFTW_FORWARD' ('-1') or
 'FFTW_BACKWARD' ('+1'), and indicates the direction of the transform you
 are interested in; technically, it is the sign of the exponent in the
 transform.
 
    The 'flags' argument is usually either 'FFTW_MEASURE' or
 'FFTW_ESTIMATE'.  'FFTW_MEASURE' instructs FFTW to run and measure the
 execution time of several FFTs in order to find the best way to compute
 the transform of size 'n'.  This process takes some time (usually a few
 seconds), depending on your machine and on the size of the transform.
 'FFTW_ESTIMATE', on the contrary, does not run any computation and just
 builds a reasonable plan that is probably sub-optimal.  In short, if
 your program performs many transforms of the same size and
 initialization time is not important, use 'FFTW_MEASURE'; otherwise use
 the estimate.
 
    _You must create the plan before initializing the input_, because
 'FFTW_MEASURE' overwrites the 'in'/'out' arrays.  (Technically,
 'FFTW_ESTIMATE' does not touch your arrays, but you should always create
 plans first just to be sure.)
 
    Once the plan has been created, you can use it as many times as you
 like for transforms on the specified 'in'/'out' arrays, computing the
 actual transforms via 'fftw_execute(plan)':
      void fftw_execute(const fftw_plan plan);
 
    The DFT results are stored in-order in the array 'out', with the
 zero-frequency (DC) component in 'out[0]'.  If 'in != out', the
 transform is "out-of-place" and the input array 'in' is not modified.
 Otherwise, the input array is overwritten with the transform.
 
    If you want to transform a _different_ array of the same size, you
 can create a new plan with 'fftw_plan_dft_1d' and FFTW automatically
 reuses the information from the previous plan, if possible.
 Alternatively, with the "guru" interface you can apply a given plan to a
 different array, if you are careful.  SeeFFTW Reference.
 
    When you are done with the plan, you deallocate it by calling
 'fftw_destroy_plan(plan)':
      void fftw_destroy_plan(fftw_plan plan);
    If you allocate an array with 'fftw_malloc()' you must deallocate it
 with 'fftw_free()'.  Do not use 'free()' or, heaven forbid, 'delete'.
 
    FFTW computes an _unnormalized_ DFT. Thus, computing a forward
 followed by a backward transform (or vice versa) results in the original
 array scaled by 'n'.  For the definition of the DFT, see SeeWhat FFTW
 Really Computes.
 
    If you have a C compiler, such as 'gcc', that supports the C99
 standard, and you '#include <complex.h>' _before_ '<fftw3.h>', then
 'fftw_complex' is the native double-precision complex type and you can
 manipulate it with ordinary arithmetic.  Otherwise, FFTW defines its own
 complex type, which is bit-compatible with the C99 complex type.  See
 Complex numbers.  (The C++ '<complex>' template class may also be
 usable via a typecast.)
 
    To use single or long-double precision versions of FFTW, replace the
 'fftw_' prefix by 'fftwf_' or 'fftwl_' and link with '-lfftw3f' or
 '-lfftw3l', but use the _same_ '<fftw3.h>' header file.
 
    Many more flags exist besides 'FFTW_MEASURE' and 'FFTW_ESTIMATE'.
 For example, use 'FFTW_PATIENT' if you're willing to wait even longer
 for a possibly even faster plan (SeeFFTW Reference).  You can also
 save plans for future use, as described by SeeWords of Wisdom-Saving
 Plans.