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) Get a printed copy>>>

GNU Scientific Library Reference Manual - Third Edition (v1.12) Buy the book here >>> $39.95

15.6 Radix-2 FFT routines for real data

This section describes radix-2 FFT algorithms for real data. They use the Cooley-Tukey algorithm to compute in-place FFTs for lengths which are a power of 2.

The radix-2 FFT functions for real data are declared in the header files ‘gsl_fft_real.h’

Function: int gsl_fft_real_radix2_transform (double data[], size_t stride, size_t n)

This function computes an in-place radix-2 FFT of length n and stride stride on the real array data. The output is a half-complex sequence, which is stored in-place. The arrangement of the half-complex terms uses the following scheme: for k < N/2 the real part of the k-th term is stored in location k, and the corresponding imaginary part is stored in location N-k. Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and count as a special case. Their real parts are stored in locations 0 and N/2 respectively, while their imaginary parts which are zero are not stored.

The following table shows the correspondence between the output data and the equivalent results obtained by considering the input data as a complex sequence with zero imaginary part (assuming stride=1),

complex[0].real    =    data[0] 
complex[0].imag    =    0 
complex[1].real    =    data[1] 
complex[1].imag    =    data[N-1]
...............         ................
complex[k].real    =    data[k]
complex[k].imag    =    data[N-k] 
...............         ................
complex[N/2].real  =    data[N/2]
complex[N/2].imag  =    0
...............         ................
complex[k'].real   =    data[k]        k' = N - k
complex[k'].imag   =   -data[N-k] 
...............         ................
complex[N-1].real  =    data[1]
complex[N-1].imag  =   -data[N-1]

Note that the output data can be converted into the full complex sequence using the function gsl_fft_halfcomplex_radix2_unpack described below.

The radix-2 FFT functions for halfcomplex data are declared in the header file ‘gsl_fft_halfcomplex.h’.

Function: int gsl_fft_halfcomplex_radix2_inverse (double data[], size_t stride, size_t n)

Function: int gsl_fft_halfcomplex_radix2_backward (double data[], size_t stride, size_t n)

These functions compute the inverse or backwards in-place radix-2 FFT of length n and stride stride on the half-complex sequence data stored according the output scheme used by gsl_fft_real_radix2. The result is a real array stored in natural order.

Function: int gsl_fft_halfcomplex_radix2_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_radix2_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[i*stride];
    double hc_imag 
      = halfcomplex_coefficient[(n-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;
  }