
Go to the first, previous, next, last section, table of contents.
Fast Fourier Transforms (FFTs)This chapter describes functions for performing Fast Fourier Transforms (FFTs). The library includes radix2 routines (for lengths which are a power of two) and mixedradix routines (which work for any length). For efficiency there are separate versions of the routines for real data and for complex data. The mixedradix routines are a reimplementation of the FFTPACK library by Paul Swarztrauber. Fortran code for FFTPACK is available on Netlib (FFTPACK also includes some routines for sine and cosine transforms but these are currently not available in GSL). For details and derivations of the underlying algorithms consult the document GSL FFT Algorithms (see section References and Further Reading) Mathematical DefinitionsFast Fourier Transforms are efficient algorithms for calculating the discrete fourier transform (DFT), x_j = \sum_{k=0}^{N1} z_k \exp(2\pi i j k / N) The DFT usually arises as an approximation to the continuous fourier transform when functions are sampled at discrete intervals in space or time. The naive evaluation of the discrete fourier transform is a matrixvector multiplication W\vec{z}. A general matrixvector multiplication takes O(N^2) operations for N datapoints. Fast fourier transform algorithms use a divideandconquer strategy to factorize the matrix W into smaller submatrices, corresponding to the integer factors of the length N. If N can be factorized into a product of integers f_1 f_2 ... f_n then the DFT can be computed in O(N \sum f_i) operations. For a radix2 FFT this gives an operation count of O(N \log_2 N). All the FFT functions offer three types of transform: forwards, inverse and backwards, based on the same mathematical definitions. The definition of the forward fourier transform, x = FFT(z), is, x_j = \sum_{k=0}^{N1} z_k \exp(2\pi i j k / N) and the definition of the inverse fourier transform, x = IFFT(z), is, z_j = {1 \over N} \sum_{k=0}^{N1} x_k \exp(2\pi i j k / N).
The factor of 1/N makes this a true inverse. For example, a call
to In general there are two possible choices for the sign of the exponential in the transform/ inversetransform pair. GSL follows the same convention as FFTPACK, using a negative exponential for the forward transform. The advantage of this convention is that the inverse transform recreates the original function with simple fourier synthesis. Numerical Recipes uses the opposite convention, a positive exponential in the forward transform. The backwards FFT is simply our terminology for an unscaled version of the inverse FFT, z^{backwards}_j = \sum_{k=0}^{N1} x_k \exp(2\pi i j k / N). When the overall scale of the result is unimportant it is often convenient to use the backwards FFT instead of the inverse to save unnecessary divisions. Overview of complex data FFTsThe inputs and outputs for the complex FFT routines are packed arrays of floating point numbers. In a packed array the real and imaginary parts of each complex number are placed in alternate neighboring elements. For example, the following definition of a packed array of length 6, gsl_complex_packed_array data[6];
can be used to hold an array of three complex numbers, data[0] = Re(z[0]) data[1] = Im(z[0]) data[2] = Re(z[1]) data[3] = Im(z[1]) data[4] = Re(z[2]) data[5] = Im(z[2])
A stride parameter allows the user to perform transforms on the
elements The array indices have the same ordering as those in the definition of the DFT  i.e. there are no index transformations or permutations of the data. For physical applications it is important to remember that the index appearing in the DFT does not correspond directly to a physical frequency. If the timestep of the DFT is \Delta then the frequencydomain includes both positive and negative frequencies, ranging from 1/(2\Delta) through 0 to +1/(2\Delta). The positive frequencies are stored from the beginning of the array up to the middle, and the negative frequencies are stored backwards from the end of the array. Here is a table which shows the layout of the array data, and the correspondence between the timedomain data z, and the frequencydomain data x. index z x = FFT(z) 0 z(t = 0) x(f = 0) 1 z(t = 1) x(f = 1/(N Delta)) 2 z(t = 2) x(f = 2/(N Delta)) . ........ .................. N/2 z(t = N/2) x(f = +1/(2 Delta), 1/(2 Delta)) . ........ .................. N3 z(t = N3) x(f = 3/(N Delta)) N2 z(t = N2) x(f = 2/(N Delta)) N1 z(t = N1) x(f = 1/(N Delta)) When N is even the location N/2 contains the most positive and negative frequencies +1/(2 \Delta), 1/(2 \Delta)) which are equivalent. If N is odd then general structure of the table above still applies, but N/2 does not appear. Radix2 FFT routines for complex dataThe radix2 algorithms described in this section are simple and compact, although not necessarily the most efficient. They use the CooleyTukey algorithm to compute inplace complex FFTs for lengths which are a power of 2  no additional storage is required. The corresponding selfsorting mixedradix routines offer better performance at the expense of requiring additional working space. All these functions are declared in the header file `gsl_fft_complex.h'.
Here is an example program which computes the FFT of a short pulse in a sample of length 128. To make the resulting fourier transform real the pulse is defined for equal positive and negative times (10 ... 10), where the negative times wrap around the end of the array. #include <stdio.h> #include <math.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_fft_complex.h> #define REAL(z,i) ((z)[2*(i)]) #define IMAG(z,i) ((z)[2*(i)+1]) int main (void) { int i; double data[2*128]; for (i = 0; i < 128; i++) { REAL(data,i) = 0.0; IMAG(data,i) = 0.0; } REAL(data,0) = 1.0; for (i = 1; i <= 10; i++) { REAL(data,i) = REAL(data,128i) = 1.0; } for (i = 0; i < 128; i++) { printf ("%d %e %e\n", i, REAL(data,i), IMAG(data,i)); } printf ("\n"); gsl_fft_complex_radix2_forward (data, 1, 128); for (i = 0; i < 128; i++) { printf ("%d %e %e\n", i, REAL(data,i)/sqrt(128), IMAG(data,i)/sqrt(128)); } return 0; }
Note that we have assumed that the program is using the default error
handler (which calls The transformed data is rescaled by 1/\sqrt N so that it fits on the same plot as the input. Only the real part is shown, by the choice of the input data the imaginary part is zero. Allowing for the wraparound of negative times at t=128, and working in units of k/N, the DFT approximates the continuum fourier transform, giving a modulated \sin function. Mixedradix FFT routines for complex dataThis section describes mixedradix FFT algorithms for complex data. The mixedradix functions work for FFTs of any length. They are a reimplementation of the Fortran FFTPACK library by Paul Swarztrauber. The theory is explained in the review article Selfsorting Mixedradix FFTs by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK. The mixedradix algorithm is based on subtransform modules  highly optimized small length FFTs which are combined to create larger FFTs. There are efficient modules for factors of 2, 3, 4, 5, 6 and 7. The modules for the composite factors of 4 and 6 are faster than combining the modules for 2*2 and 2*3. For factors which are not implemented as modules there is a fallback to a general lengthn module which uses Singleton's method for efficiently computing a DFT. This module is O(n^2), and slower than a dedicated module would be but works for any length n. Of course, lengths which use the general lengthn module will still be factorized as much as possible. For example, a length of 143 will be factorized into 11*13. Large prime factors are the worst case scenario, e.g. as found in n=2*3*99991, and should be avoided because their O(n^2) scaling will dominate the runtime (consult the document GSL FFT Algorithms included in the GSL distribution if you encounter this problem).
The mixedradix initialization function All these functions are declared in the header file `gsl_fft_complex.h'.
gsl_fft_complex_wavetable structure
which contains internal parameters for the FFT. It is not necessary to
set any of the components directly but it can sometimes be useful to
examine them. For example, the chosen factorization of the FFT length
is given and can be used to provide an estimate of the runtime or
numerical error.
The wavetable structure is declared in the header file `gsl_fft_complex.h'.
The mixed radix algorithms require an additional working space to hold the intermediate steps of the transform.
The following functions compute the transform,
Here is an example program which computes the FFT of a short pulse in a sample of length 630 (=2*3*3*5*7) using the mixedradix algorithm. #include <stdio.h> #include <math.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_fft_complex.h> #define REAL(z,i) ((z)[2*(i)]) #define IMAG(z,i) ((z)[2*(i)+1]) int main (void) { int i; const int n = 630; double data[2*n]; gsl_fft_complex_wavetable * wavetable; gsl_fft_complex_workspace * workspace; for (i = 0; i < n; i++) { REAL(data,i) = 0.0; IMAG(data,i) = 0.0; } data[0].real = 1.0; for (i = 1; i <= 10; i++) { REAL(data,i) = REAL(data,ni) = 1.0; } for (i = 0; i < n; i++) { printf ("%d: %e %e\n", i, REAL(data,i), IMAG(data,i)); } printf ("\n"); wavetable = gsl_fft_complex_wavetable_alloc (n); workspace = gsl_fft_complex_workspace_alloc (n); for (i = 0; i < wavetable>nf; i++) { printf("# factor %d: %d\n", i, wavetable>factor[i]); } gsl_fft_complex_forward (data, 1, n, wavetable, workspace); for (i = 0; i < n; i++) { printf ("%d: %e %e\n", i, REAL(data,i), IMAG(data,i)); } gsl_fft_complex_wavetable_free (wavetable); gsl_fft_complex_workspace_free (workspace); return 0; }
Note that we have assumed that the program is using the default
Overview of real data FFTsThe functions for real data are similar to those for complex data. However, there is an important difference between forward and inverse transforms. The fourier transform of a real sequence is not real. It is a complex sequence with a special symmetry: z_k = z_{Nk}^*
A sequence with this symmetry is called conjugatecomplex or
halfcomplex. This different structure requires different
storage layouts for the forward transform (from real to halfcomplex)
and inverse transform (from halfcomplex back to real). As a
consequence the routines are divided into two sets: functions in
Functions in c_k = \sum_{j=0}^{N1} x_k \exp(2 \pi i j k /N)
Functions in x_j = {1 \over N} \sum_{k=0}^{N1} c_k \exp(2 \pi i j k /N) The symmetry of the halfcomplex sequence implies that only half of the complex numbers in the output need to be stored. The remaining half can be reconstructed using the halfcomplex symmetry condition. (This works for all lengths, even and odd. When the length is even the middle value, where k=N/2, is also real). Thus only N real numbers are required to store the halfcomplex sequence, and the transform of a real sequence can be stored in the same size array as the original data. The precise storage arrangements depend on the algorithm, and are different for radix2 and mixedradix routines. The radix2 function operates inplace, which constrain the locations where each element can be stored. The restriction forces real and imaginary parts to be stored far apart. The mixedradix algorithm does not have this restriction, and it stores the real and imaginary parts of a given term in neighboring locations. This is desirable for better locality of memory accesses. Radix2 FFT routines for real dataThis section describes radix2 FFT algorithms for real data. They use the CooleyTukey algorithm to compute inplace FFTs for lengths which are a power of 2. The radix2 FFT functions for real data are declared in the header files `gsl_fft_real.h'
The radix2 FFT functions for halfcomplex data are declared in the header file `gsl_fft_halfcomplex.h'.
Mixedradix FFT routines for real dataThis section describes mixedradix FFT algorithms for real data. The mixedradix functions work for FFTs of any length. They are a reimplementation of the realFFT routines in the Fortran FFTPACK library by Paul Swarztrauber. The theory behind the algorithm is explained in the article Fast MixedRadix Real Fourier Transforms by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK. The functions use the FFTPACK storage convention for halfcomplex sequences. In this convention the halfcomplex transform of a real sequence is stored with frequencies in increasing order, starting at zero, with the real and imaginary parts of each frequency in neighboring locations. When a value is known to be real the imaginary part is not stored. The imaginary part of the zerofrequency component is never stored. It is known to be zero (since the zero frequency component is simply the sum of the input data (all real)). For a sequence of even length the imaginary part of the frequency n/2 is not stored either, since the symmetry z_k = z_{Nk}^* implies that this is purely real too.
The storage scheme is best shown by some examples. The table below
shows the output for an oddlength sequence, n=5. The two columns
give the correspondence between the 5 values in the halfcomplex
sequence returned by complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[3] complex[3].imag = halfcomplex[4] complex[4].real = halfcomplex[1] complex[4].imag = halfcomplex[2]
The upper elements of the complex array, The next table shows the output for an evenlength sequence, n=5 In the even case both the there are two values which are purely real, complex[0].real = halfcomplex[0] complex[0].imag = 0 complex[1].real = halfcomplex[1] complex[1].imag = halfcomplex[2] complex[2].real = halfcomplex[3] complex[2].imag = halfcomplex[4] complex[3].real = halfcomplex[5] complex[3].imag = 0 complex[4].real = halfcomplex[3] complex[4].imag = halfcomplex[4] complex[5].real = halfcomplex[1] complex[5].imag = halfcomplex[2]
The upper elements of the complex array, All these functions are declared in the header files `gsl_fft_real.h' and `gsl_fft_halfcomplex.h'.
Here is an example program using The remaining fourier coefficients are transformed back to the timedomain, to give a filtered version of the square pulse. Since fourier coefficients are stored using the halfcomplex symmetry both positive and negative frequencies are removed and the final filtered signal is also real. #include <stdio.h> #include <math.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_fft_real.h> #include <gsl/gsl_fft_halfcomplex.h> int main (void) { int i, n = 100; double data[n]; gsl_fft_real_wavetable * real; gsl_fft_halfcomplex_wavetable * hc; gsl_fft_real_workspace * work; for (i = 0; i < n; i++) { data[i] = 0.0; } for (i = n / 3; i < 2 * n / 3; i++) { data[i] = 1.0; } for (i = 0; i < n; i++) { printf ("%d: %e\n", i, data[i]); } printf ("\n"); work = gsl_fft_real_workspace_alloc (n); real = gsl_fft_real_wavetable_alloc (n); gsl_fft_real_transform (data, 1, n, real, work); gsl_fft_real_wavetable_free (real); for (i = 11; i < n; i++) { data[i] = 0; } hc = gsl_fft_halfcomplex_wavetable_alloc (n); gsl_fft_halfcomplex_inverse (data, 1, n, hc, work); gsl_fft_halfcomplex_wavetable_free (hc); for (i = 0; i < n; i++) { printf ("%d: %e\n", i, data[i]); } gsl_fft_real_workspace_free (work); return 0; } References and Further ReadingA good starting point for learning more about the FFT is the review article Fast Fourier Transforms: A Tutorial Review and A State of the Art by Duhamel and Vetterli,
To find out about the algorithms used in the GSL routines you may want to consult the latex document GSL FFT Algorithms (it is included in GSL, as `doc/fftalgorithms.tex'). This has general information on FFTs and explicit derivations of the implementation for each routine. There are also references to the relevant literature. For convenience some of the more important references are reproduced below. There are several introductory books on the FFT with example programs, such as The Fast Fourier Transform by Brigham and DFT/FFT and Convolution Algorithms by Burrus and Parks,
Both these introductory books cover the radix2 FFT in some detail. The mixedradix algorithm at the heart of the FFTPACK routines is reviewed in Clive Temperton's paper,
The derivation of FFTs for realvalued data is explained in the following two articles,
In 1979 the IEEE published a compendium of carefullyreviewed Fortran FFT programs in Programs for Digital Signal Processing. It is a useful reference for implementations of many different FFT algorithms,
For serious FFT work we recommend the use of the dedicated FFTW library by Frigo and Johnson. The FFTW library is selfoptimizing  it automatically tunes itself for each hardware platform in order to achieve maximum performance. It is available under the GNU GPL.
Go to the first, previous, next, last section, table of contents. 