- publishing free software manuals
 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. RossiPaperback (6"x9"), 592 pages, 60 figuresISBN 0954612078RRP £24.95 (\$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;
}
```
 ISBN 0954612078 GNU Scientific Library Reference Manual - Third Edition (v1.12) See the print edition