Next: FFT References and Further Reading, Previous: Radix-2 FFT routines for real data, Up: Fast Fourier Transforms [Index]

This section describes mixed-radix FFT algorithms for real data. The mixed-radix functions work for FFTs of any length. They are a reimplementation of the real-FFT routines in the Fortran FFTPACK library by Paul Swarztrauber. The theory behind the algorithm is explained in the article Fast Mixed-Radix 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 half-complex
sequences. In this convention the half-complex 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 zero-frequency 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_{n-k}^** implies that this is
purely real too.

The storage scheme is best shown by some examples. The table below
shows the output for an odd-length sequence, *n=5*. The two columns
give the correspondence between the 5 values in the half-complex
sequence returned by `gsl_fft_real_transform`

, `halfcomplex`[] and the
values `complex`[] that would be returned if the same real input
sequence were passed to `gsl_fft_complex_backward`

as a complex
sequence (with imaginary parts set to `0`

),

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, `complex[3]`

and
`complex[4]`

are filled in using the symmetry condition. The
imaginary part of the zero-frequency term `complex[0].imag`

is
known to be zero by the symmetry.

The next table shows the output for an even-length sequence, *n=6*.
In the even case 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, `complex[4]`

and
`complex[5]`

are filled in using the symmetry condition. Both
`complex[0].imag`

and `complex[3].imag`

are known to be zero.

All these functions are declared in the header files
`gsl_fft_real.h` and `gsl_fft_halfcomplex.h`.

- Function:
*gsl_fft_real_wavetable ****gsl_fft_real_wavetable_alloc***(size_t*`n`) - Function:
*gsl_fft_halfcomplex_wavetable ****gsl_fft_halfcomplex_wavetable_alloc***(size_t*`n`) -
These functions prepare trigonometric lookup tables for an FFT of size

*n*real elements. The functions return a pointer to the newly allocated struct 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 computing the wavetable 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 appropriate type of wavetable must be used for forward real or inverse half-complex transforms.

- Function:
*void***gsl_fft_real_wavetable_free***(gsl_fft_real_wavetable **`wavetable`) - Function:
*void***gsl_fft_halfcomplex_wavetable_free***(gsl_fft_halfcomplex_wavetable **`wavetable`) These functions free the memory associated with the wavetable

`wavetable`. The wavetable can be freed if no further FFTs of the same length will be needed.

The mixed radix algorithms require additional working space to hold the intermediate steps of the transform,

- Function:
*gsl_fft_real_workspace ****gsl_fft_real_workspace_alloc***(size_t*`n`) -
This function allocates a workspace for a real transform of length

`n`. The same workspace can be used for both forward real and inverse halfcomplex transforms.

- Function:
*void***gsl_fft_real_workspace_free***(gsl_fft_real_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 transforms of real and half-complex data,

- Function:
*int***gsl_fft_real_transform***(double*`data`[], size_t`stride`, size_t`n`, const gsl_fft_real_wavetable *`wavetable`, gsl_fft_real_workspace *`work`) - Function:
*int***gsl_fft_halfcomplex_transform***(double*`data`[], size_t`stride`, size_t`n`, const gsl_fft_halfcomplex_wavetable *`wavetable`, gsl_fft_real_workspace *`work`) These functions compute the FFT of

`data`, a real or half-complex array of length`n`, using a mixed radix decimation-in-frequency algorithm. For`gsl_fft_real_transform`

`data`is an array of time-ordered real data. For`gsl_fft_halfcomplex_transform`

`data`contains Fourier coefficients in the half-complex ordering described above. There is no restriction on the length`n`. Efficient modules are provided for subtransforms of length 2, 3, 4 and 5. Any remaining factors are computed with a slow,*O(n^2)*, general-n module. The caller must supply a`wavetable`containing trigonometric lookup tables and a workspace`work`.

- Function:
*int***gsl_fft_real_unpack***(const double*`real_coefficient`[], gsl_complex_packed_array`complex_coefficient`, size_t`stride`, size_t`n`) -
This function converts a single real array,

`real_coefficient`into an equivalent complex array,`complex_coefficient`, (with imaginary part set to zero), suitable for`gsl_fft_complex`

routines. The algorithm for the conversion is simply,for (i = 0; i < n; i++) { complex_coefficient[i*stride].real = real_coefficient[i*stride]; complex_coefficient[i*stride].imag = 0.0; }

- Function:
*int***gsl_fft_halfcomplex_unpack***(const double*`halfcomplex_coefficient`[], gsl_complex_packed_array`complex_coefficient`, size_t`stride`, size_t`n`) -
This function converts

`halfcomplex_coefficient`, an array of half-complex coefficients as returned by`gsl_fft_real_transform`

, into an ordinary complex array,`complex_coefficient`. It fills in the complex array using the symmetry*z_k = z_{n-k}^**to reconstruct the redundant elements. The algorithm for the conversion is,complex_coefficient[0].real = halfcomplex_coefficient[0]; complex_coefficient[0].imag = 0.0; for (i = 1; i < n - i; i++) { double hc_real = halfcomplex_coefficient[(2 * i - 1)*stride]; double hc_imag = halfcomplex_coefficient[(2 * i)*stride]; complex_coefficient[i*stride].real = hc_real; complex_coefficient[i*stride].imag = hc_imag; complex_coefficient[(n - i)*stride].real = hc_real; complex_coefficient[(n - i)*stride].imag = -hc_imag; } if (i == n - i) { complex_coefficient[i*stride].real = halfcomplex_coefficient[(n - 1)*stride]; complex_coefficient[i*stride].imag = 0.0; }

Here is an example program using `gsl_fft_real_transform`

and
`gsl_fft_halfcomplex_inverse`

. It generates a real signal in the
shape of a square pulse. The pulse is Fourier transformed to frequency
space, and all but the lowest ten frequency components are removed from
the array of Fourier coefficients returned by
`gsl_fft_real_transform`

.

The remaining Fourier coefficients are transformed back to the time-domain, to give a filtered version of the square pulse. Since Fourier coefficients are stored using the half-complex 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; }