Abstract

Introduction

FFTW(Fastest Fourier Transform in the West) is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data (as well as of even/odd data, i.e. the discrete cosine/sine transforms or DCT/DST). We believe that FFTW, which is free software, should become the FFT library of choice for most applications.

Install

ここから落として、 解凍⇒「./configure」⇒「make」⇒「make install」で「/usr/local/」以下にインストール。 もしくは、

tar zxvf fftw-3.1.2.tar.gz
mv fftw-3.1.2 fftw-3.1.2_00
cd  fftw-3.1.2_00
mkdir  -p /home/usrname/opt/pkgs/fftw-3.1.2
./configure CFLAGS=-m64 --with-pic --enable-threads --prefix=/home/usrname/opt/pkgs/fftw-3.1.2 
make
make check
 --------------------------------------------------------------
          FFTW transforms passed basic tests!
 --------------------------------------------------------------
make install
ln -s /home/usrname/opt/pkgs/fftw-3.1.2 /home/usrname/opt/fftw

「おまけ(?)」もインストールしたいのであれば、以下の様にする:

./configure CFLAGS=-m64 --with-pic --enable-shared --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install
make clean
./configure CFLAGS=-m64 --with-pic --enable-shared --enable-threads --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install
make clean
./configure CFLAGS=-m64 --with-pic --enable-shared --enable-single --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install
make clean
./configure CFLAGS=-m64 --with-pic --enable-shared --enable-single --enable-threads --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install
make clean
./configure CFLAGS=-m64 --with-pic --enable-shared --enable-long-double --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install
make clean
./configure CFLAGS=-m64 --with-pic --enable-shared --enable-long-double --enable-threads --prefix=/Users/username/opt/pkgs/fftw-3.2.2
make
make check
make install

Features

FFTW 3.1.2 is the latest official version of FFTW (refer to the release notes to find out what is new). Subscribe to the fftw-announce mailing list to receive announcements of future updates. Here is a list of some of FFTW's more interesting features:

  • Speed. (Supports SSE/SSE2/3dNow!/Altivec, since version 3.0.)
  • Both one-dimensional and multi-dimensional transforms.
  • Arbitrary-size transforms. (Sizes with small prime factors are best, but FFTW uses O(N log N) algorithms even for prime sizes.)
  • Fast transforms of purely real input or output data.
  • Transforms of real even/odd data: the discrete cosine transform (DCT) and the discrete sine transform (DST), types I-IV. (Version 3.0 or later.)
  • Efficient handling of multiple, strided transforms. (This lets you do things like transform multiple arrays at once, transform one dimension of a multi-dimensional array, or transform one field of a multi-component array.)
  • Parallel transforms: parallelized code for platforms with Cilk or for SMP machines with some flavor of threads (e.g. POSIX). An MPI version for distributed-memory transforms is also available, currently only as part of FFTW 2.1.5.
  • Portable to any platform with a C compiler. Documentation in HTML and other formats.
  • Both C and Fortran interfaces.
  • Free software, released under the GNU General Public License (GPL, see FFTW license). (Non-free licenses may also be purchased from MIT, for users who do not want their programs protected by the GPL. Contact us for details.) (Also see the FAQ.)

注意

from FAQ

Question 1.5. In the West? I thought MIT was in the East?
Not to an Italian. You could say that we're a Spaghetti Western (with apologies to Sergio Leone).
Question 2.9. Can I call FFTW from C++?
Most definitely. FFTW should compile and/or link under any C++ compiler. Moreover, it is likely that the C++ <complex> template class is bit-compatible with FFTW's complex-number format (see the FFTW manual for more details).
Question 3.10. Why does your inverse transform return a scaled result?
Computing the forward transform followed by the backward transform (or vice versa) yields the original array scaled by the size of the array. (For multi-dimensional transforms, the size of the array is the product of the dimensions.) We could, instead, have chosen a normalization that would have returned the unscaled array. Or, to accomodate the many conventions in this matter, the transform routines could have accepted a "scale factor" parameter. We did not do this, however, for two reasons. First, we didn't want to sacrifice performance in the common case where the scale factor is 1. Second, in real applications the FFT is followed or preceded by some computation on the data, into which the scale factor can typically be absorbed at little or no cost.
Question 3.11. How can I make FFTW put the origin (zero frequency) at the center of its output?
For human viewing of a spectrum, it is often convenient to put the origin in frequency space at the center of the output array, rather than in the zero-th element (the default in FFTW). If all of the dimensions of your array are even, you can accomplish this by simply multiplying each element of the input array by (-1)^(i + j + ...), where i, j, etcetera are the indices of the element. (This trick is a general property of the DFT, and is not specific to FFTW.)
Question 3.17. The output of FFTW's transform is all zeros.
You should initialize your input array after creating the plan, unless you use FFTW_ESTIMATE: planning with FFTW_MEASURE or FFTW_PATIENT overwrites the input/output arrays, as described in the manual.

from User's manual

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);
    }
(When you compile, you must also link with the fftw3 library, e.g. -lfftw3 -lm on Unix systems.)

First you allocate 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 (see SIMD alignment and fftw_malloc).

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. The data in the in/out arrays is overwritten during FFTW_MEASURE planning, so such planning should be done before the input is initialized by the user.

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);
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. See FFTW Reference.)

When you are done with the plan, you deallocate it by calling fftw_destroy_plan(plan):

    void fftw_destroy_plan(fftw_plan plan);

Arrays allocated with fftw_malloc should be deallocated by fftw_free rather than the ordinary free (or, heaven forbid, delete). 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.

Users should note that 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 What FFTW Really Computes. If you have a C compiler, such as gcc, that supports the recent 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.) Single and long-double precision versions of FFTW may be installed; to use them, 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 (see FFTW Reference). You can also save plans for future use, as described by Words of Wisdom-Saving Plans.

4.1.1 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 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).

使い方

Complex DFTs

    fftw_plan fftw_plan_dft_1d(int n,
                               fftw_complex *in, fftw_complex *out,
                               int sign, unsigned flags);
    fftw_plan fftw_plan_dft_2d(int nx, int ny,
                               fftw_complex *in, fftw_complex *out,
                               int sign, unsigned flags);
    fftw_plan fftw_plan_dft_3d(int nx, int ny, int nz,
                               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. 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.

Arguments

  • rank is the dimensionality of the transform (it should be the size of the array *n), 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 transform of size 1, i.e. a copy of one number from input to output.
  • n, or nx/ny/nz, or n[rank], respectively, gives the size of the transform dimensions. They can be any positive integer.
    • Multi-dimensional arrays are stored in row-major order with dimensions: nx x ny; or nx x ny x nz; 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 O(n log n)
      n12345678910
      2^n2481632641282565121024
      3^n392781243729218765611968359049
      5^n5251256253125156257812539062519531259765625
      7^n749343240116807117649823543576480140353607282475249

      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.

Real-data DFTs

    fftw_plan fftw_plan_dft_r2c_1d(int n,
                                   double *in, fftw_complex *out,
                                   unsigned flags);
    fftw_plan fftw_plan_dft_r2c_2d(int nx, int ny,
                                   double *in, fftw_complex *out,
                                   unsigned flags);
    fftw_plan fftw_plan_dft_r2c_3d(int nx, int ny, int nz,
                                   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 dimensionality of the transform (it should be the size of the array *n), 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 transform of size 1, i.e. a copy of one number (with zero imaginary part) from input to output.
  • n, or nx/ny/nz, or n[rank], respectively, gives the size of the logical 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 O(n log n)
      n12345678910
      2^n2481632641282565121024
      3^n392781243729218765611968359049
      5^n5251256253125156257812539062519531259765625
      7^n749343240116807117649823543576480140353607282475249
    • 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 n,
                                   fftw_complex *in, double *out,
                                   unsigned flags);
    fftw_plan fftw_plan_dft_c2r_2d(int nx, int ny,
                                   fftw_complex *in, double *out,
                                   unsigned flags);
    fftw_plan fftw_plan_dft_c2r_3d(int nx, int ny, int nz,
                                   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.