GNU Scientific Library Reference Manual  Third Edition (v1.12) by M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, F. Rossi Paperback (6"x9"), 592 pages, 60 figures ISBN 0954612078 RRP £24.95 ($39.95) 
15.4 Mixedradix FFT routines for complex data
This section describes mixedradix FFT algorithms for complex data. The mixedradix functions work for FFTs of any length. They are a reimplementation of Paul Swarztrauber's Fortran fftpack library. 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 moduleshighly 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 gsl_fft_complex_wavetable_alloc
returns the list of factors chosen by the library for a given length
N. It can be used to check how well the length has been
factorized, and estimate the runtime. To a first approximation the
runtime scales as N \sum f_i, where the f_i are the
factors of N. For programs under user control you may wish to
issue a warning that the transform will be slow when the length is
poorly factorized. If you frequently encounter data lengths which
cannot be factorized using the existing smallprime modules consult
GSL FFT Algorithms for details on adding support for other
factors.
All the functions described in this section are declared in the header file ‘gsl_fft_complex.h’.
 Function: gsl_fft_complex_wavetable * gsl_fft_complex_wavetable_alloc (size_t n)
 This function prepares a trigonometric lookup table for a complex FFT of
length n. The function returns a pointer to the newly allocated
gsl_fft_complex_wavetable
if no errors were detected, and a null pointer in the case of error. The length n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls tosin
andcos
, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then this computation is a oneoff overhead which does not affect the final throughput.The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The same wavetable can be used for both forward and backward (or inverse) transforms of a given length.
 Function: void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * wavetable)
 This function frees the memory associated with the wavetable wavetable. The wavetable can be freed if no further FFTs of the same length will be needed.
These functions operate on a 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’.
 Data Type: gsl_fft_complex_wavetable
 This is a structure that holds the factorization and trigonometric
lookup tables for the mixed radix fft algorithm. It has the following
components:
size_t n
 This is the number of complex data points
size_t nf

This is the number of factors that the length
n
was decomposed into. size_t factor[64]

This is the array of factors. Only the first
nf
elements are used. gsl_complex * trig

This is a pointer to a preallocated trigonometric lookup table of
n
complex elements. gsl_complex * twiddle[64]

This is an array of pointers into
trig
, giving the twiddle factors for each pass.
The mixed radix algorithms require additional working space to hold the intermediate steps of the transform.
 Function: gsl_fft_complex_workspace * gsl_fft_complex_workspace_alloc (size_t n)
 This function allocates a workspace for a complex transform of length n.
 Function: void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * workspace)
 This function frees the memory associated with the workspace workspace. The workspace can be freed if no further FFTs of the same length will be needed.
The following functions compute the transform,
 Function: int gsl_fft_complex_forward (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
 Function: int gsl_fft_complex_transform (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work, gsl_fft_direction sign)
 Function: int gsl_fft_complex_backward (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
 Function: int gsl_fft_complex_inverse (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
 These functions compute forward, backward and inverse FFTs of length
n with stride stride, on the packed complex array
data, using a mixed radix decimationinfrequency algorithm.
There is no restriction on the length n. Efficient modules are
provided for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining
factors are computed with a slow, O(n^2), generaln
module. The caller must supply a wavetable containing the
trigonometric lookup tables and a workspace work. For the
transform
version of the function the sign argument can be eitherforward
(1) orbackward
(+1).The functions return a value of
0
if no errors were detected. The followinggsl_errno
conditions are defined for these functions:GSL_EDOM
 The length of the data n is not a positive integer (i.e. n is zero).
GSL_EINVAL
 The length of the data n and the length used to compute the given wavetable do not match.
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] = 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
gsl
error handler (which calls abort
for any errors). If
you are not using a safe error handler you would need to check the
return status of all the gsl
routines.
ISBN 0954612078  GNU Scientific Library Reference Manual  Third Edition (v1.12)  See the print edition 