This section describes mixed-radix FFT algorithms for complex data. The mixed-radix 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 Self-sorting Mixed-radix FFTs by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK.
The mixed-radix algorithm is based on sub-transform 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 @math{2*2} and @math{2*3}.
For factors which are not implemented as modules there is a fall-back to a general length-@math{n} module which uses Singleton's method for efficiently computing a DFT. This module is @math{O(n^2)}, and slower than a dedicated module would be but works for any length @math{n}. Of course, lengths which use the general length-@math{n} module will still be factorized as much as possible. For example, a length of 143 will be factorized into @math{11*13}. Large prime factors are the worst case scenario, e.g. as found in @math{n=2*3*99991}, and should be avoided because their @math{O(n^2)} scaling will dominate the run-time (consult the document GSL FFT Algorithms included in the GSL distribution if you encounter this problem).
The mixed-radix initialization function gsl_fft_complex_wavetable_alloc
returns the list of factors chosen by the library for a given length
@math{N}. It can be used to check how well the length has been
factorized, and estimate the run-time. To a first approximation the
run-time scales as @math{N \sum f_i}, where the @math{f_i} are the
factors of @math{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 small-prime modules consult
GSL FFT Algorithms for details on adding support for other
factors.
All these functions are declared in the header file `gsl_fft_complex.h'.
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 to sin
and cos
, 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 one-off 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.
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 run-time or
numerical error.
The wavetable structure is declared in the header file `gsl_fft_complex.h'.
size_t n
size_t nf
n
was decomposed into.
size_t factor[64]
nf
elements are
used.
gsl_complex * trig
n
complex elements.
gsl_complex * twiddle[64]
trig
, giving the twiddle
factors for each pass.
The mixed radix algorithms require an additional working space to hold the intermediate steps of the transform.
The following functions compute the transform,
These functions compute forward, backward and inverse FFTs of length n with stride stride, on the packed complex array data, using a mixed radix decimation-in-frequency 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, @math{O(n^2)}, general-@math{n} module. The caller must supply a wavetable containing the trigonometric lookup tables and a workspace work.
The functions return a value of 0
if no errors were detected. The
following gsl_errno
conditions are defined for these functions:
GSL_EDOM
GSL_EINVAL
Here is an example program which computes the FFT of a short pulse in a sample of length 630 (@math{=2*3*3*5*7}) using the mixed-radix 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,n-i) = 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.