Added MACS source

[htsworkflow.git] / htswanalysis / MACS / lib / gsl / gsl-1.11 / doc / fft.texi

@cindex FFT
@cindex Fast Fourier Transforms, see FFT
@cindex Fourier Transforms, see FFT
@cindex Discrete Fourier Transforms, see FFT
@cindex DFTs, see FFT

This chapter describes functions for performing Fast Fourier Transforms
(FFTs).  The library includes radix-2 routines (for lengths which are a
power of two) and mixed-radix routines (which work for any length).  For
efficiency there are separate versions of the routines for real data and
for complex data.  The mixed-radix routines are a reimplementation of the
@sc{fftpack} library of Paul Swarztrauber.  Fortran code for @sc{fftpack} is
available on Netlib (@sc{fftpack} also includes some routines for sine and
cosine transforms but these are currently not available in GSL).  For
details and derivations of the underlying algorithms consult the
document @cite{GSL FFT Algorithms} (@pxref{FFT References and Further Reading})

@menu
* Mathematical Definitions::    
* Overview of complex data FFTs::  
* Radix-2 FFT routines for complex data::  
* Mixed-radix FFT routines for complex data::  
* Overview of real data FFTs::  
* Radix-2 FFT routines for real data::  
* Mixed-radix FFT routines for real data::  
* FFT References and Further Reading::  
@end menu

@node Mathematical Definitions
@section Mathematical Definitions 
@cindex FFT mathematical definition

Fast Fourier Transforms are efficient algorithms for
calculating the discrete fourier transform (DFT),
@tex
\beforedisplay
$$
x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N) 
$$
\afterdisplay
@end tex
@ifinfo

@example
x_j = \sum_@{k=0@}^@{N-1@} z_k \exp(-2\pi i j k / N) 
@end example
@end ifinfo

The DFT usually arises as an approximation to the continuous fourier
transform when functions are sampled at discrete intervals in space or
time.  The naive evaluation of the discrete fourier transform is a
matrix-vector multiplication 
@c{$W\vec{z}$}
@math{W\vec@{z@}}. A general matrix-vector multiplication takes
@math{O(N^2)} operations for @math{N} data-points.  Fast fourier
transform algorithms use a divide-and-conquer strategy to factorize the
matrix @math{W} into smaller sub-matrices, corresponding to the integer
factors of the length @math{N}.  If @math{N} can be factorized into a
product of integers
@c{$f_1 f_2 \ldots f_n$} 
@math{f_1 f_2 ... f_n} then the DFT can be computed in @math{O(N \sum
f_i)} operations.  For a radix-2 FFT this gives an operation count of
@math{O(N \log_2 N)}.

All the FFT functions offer three types of transform: forwards, inverse
and backwards, based on the same mathematical definitions.  The
definition of the @dfn{forward fourier transform},
@c{$x = \hbox{FFT}(z)$}
@math{x = FFT(z)}, is,
@tex
\beforedisplay
$$
x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N) 
$$
\afterdisplay
@end tex
@ifinfo

@example
x_j = \sum_@{k=0@}^@{N-1@} z_k \exp(-2\pi i j k / N) 
@end example

@end ifinfo
@noindent
and the definition of the @dfn{inverse fourier transform},
@c{$x = \hbox{IFFT}(z)$}
@math{x = IFFT(z)}, is,
@tex
\beforedisplay
$$
z_j = {1 \over N} \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N).
$$
\afterdisplay
@end tex
@ifinfo

@example
z_j = @{1 \over N@} \sum_@{k=0@}^@{N-1@} x_k \exp(2\pi i j k / N).
@end example

@end ifinfo
@noindent
The factor of @math{1/N} makes this a true inverse.  For example, a call
to @code{gsl_fft_complex_forward} followed by a call to
@code{gsl_fft_complex_inverse} should return the original data (within
numerical errors).

In general there are two possible choices for the sign of the
exponential in the transform/ inverse-transform pair. GSL follows the
same convention as @sc{fftpack}, using a negative exponential for the forward
transform.  The advantage of this convention is that the inverse
transform recreates the original function with simple fourier
synthesis.  Numerical Recipes uses the opposite convention, a positive
exponential in the forward transform.

The @dfn{backwards FFT} is simply our terminology for an unscaled
version of the inverse FFT,
@tex
\beforedisplay
$$
z^{backwards}_j = \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N).
$$
\afterdisplay
@end tex
@ifinfo

@example
z^@{backwards@}_j = \sum_@{k=0@}^@{N-1@} x_k \exp(2\pi i j k / N).
@end example

@end ifinfo
@noindent
When the overall scale of the result is unimportant it is often
convenient to use the backwards FFT instead of the inverse to save
unnecessary divisions.

@node Overview of complex data FFTs
@section Overview of complex data FFTs
@cindex FFT, complex data

The inputs and outputs for the complex FFT routines are @dfn{packed
arrays} of floating point numbers.  In a packed array the real and
imaginary parts of each complex number are placed in alternate
neighboring elements.  For example, the following definition of a packed
array of length 6,

@example
double x[3*2];
gsl_complex_packed_array data = x;
@end example

@noindent
can be used to hold an array of three complex numbers, @code{z[3]}, in
the following way,

@example
data[0] = Re(z[0])
data[1] = Im(z[0])
data[2] = Re(z[1])
data[3] = Im(z[1])
data[4] = Re(z[2])
data[5] = Im(z[2])
@end example

@noindent
The array indices for the data have the same ordering as those
in the definition of the DFT---i.e. there are no index transformations
or permutations of the data.

A @dfn{stride} parameter allows the user to perform transforms on the
elements @code{z[stride*i]} instead of @code{z[i]}.  A stride greater
than 1 can be used to take an in-place FFT of the column of a matrix. A
stride of 1 accesses the array without any additional spacing between
elements.  

To perform an FFT on a vector argument, such as @code{gsl_vector_complex
* v}, use the following definitions (or their equivalents) when calling
the functions described in this chapter:

@example
gsl_complex_packed_array data = v->data;
size_t stride = v->stride;
size_t n = v->size;
@end example

For physical applications it is important to remember that the index
appearing in the DFT does not correspond directly to a physical
frequency.  If the time-step of the DFT is @math{\Delta} then the
frequency-domain includes both positive and negative frequencies,
ranging from @math{-1/(2\Delta)} through 0 to @math{+1/(2\Delta)}.  The
positive frequencies are stored from the beginning of the array up to
the middle, and the negative frequencies are stored backwards from the
end of the array.

Here is a table which shows the layout of the array @var{data}, and the
correspondence between the time-domain data @math{z}, and the
frequency-domain data @math{x}.

@example
index    z               x = FFT(z)

0        z(t = 0)        x(f = 0)
1        z(t = 1)        x(f = 1/(N Delta))
2        z(t = 2)        x(f = 2/(N Delta))
.        ........        ..................
N/2      z(t = N/2)      x(f = +1/(2 Delta),
                               -1/(2 Delta))
.        ........        ..................
N-3      z(t = N-3)      x(f = -3/(N Delta))
N-2      z(t = N-2)      x(f = -2/(N Delta))
N-1      z(t = N-1)      x(f = -1/(N Delta))
@end example

@noindent
When @math{N} is even the location @math{N/2} contains the most positive
and negative frequencies (@math{+1/(2 \Delta)}, @math{-1/(2 \Delta)})
which are equivalent.  If @math{N} is odd then general structure of the
table above still applies, but @math{N/2} does not appear.


@node Radix-2 FFT routines for complex data
@section Radix-2 FFT routines for complex data
@cindex FFT of complex data, radix-2 algorithm
@cindex Radix-2 FFT, complex data

The radix-2 algorithms described in this section are simple and compact,
although not necessarily the most efficient.  They use the Cooley-Tukey
algorithm to compute in-place complex FFTs for lengths which are a power
of 2---no additional storage is required.  The corresponding
self-sorting mixed-radix routines offer better performance at the
expense of requiring additional working space.

All the functions described in this section are declared in the header file @file{gsl_fft_complex.h}.

@deftypefun int gsl_fft_complex_radix2_forward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

@deftypefunx int gsl_fft_complex_radix2_transform (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, gsl_fft_direction @var{sign})

@deftypefunx int gsl_fft_complex_radix2_backward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

@deftypefunx int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

These functions compute forward, backward and inverse FFTs of length
@var{n} with stride @var{stride}, on the packed complex array @var{data}
using an in-place radix-2 decimation-in-time algorithm.  The length of
the transform @var{n} is restricted to powers of two.  For the
@code{transform} version of the function the @var{sign} argument can be
either @code{forward} (@math{-1}) or @code{backward} (@math{+1}).

The functions return a value of @code{GSL_SUCCESS} if no errors were
detected, or @code{GSL_EDOM} if the length of the data @var{n} is not a
power of two.
@end deftypefun

@deftypefun int gsl_fft_complex_radix2_dif_forward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

@deftypefunx int gsl_fft_complex_radix2_dif_transform (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, gsl_fft_direction @var{sign})

@deftypefunx int gsl_fft_complex_radix2_dif_backward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

@deftypefunx int gsl_fft_complex_radix2_dif_inverse (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n})

These are decimation-in-frequency versions of the radix-2 FFT functions.

@end deftypefun


@comment @node Example of using radix-2 FFT routines for complex data
@comment @subsection Example of using radix-2 FFT routines for complex data

Here is an example program which computes the FFT of a short pulse in a
sample of length 128.  To make the resulting fourier transform real the
pulse is defined for equal positive and negative times (@math{-10}
@dots{} @math{10}), where the negative times wrap around the end of the
array.

@example
@verbatiminclude examples/fft.c
@end example

@noindent
Note that we have assumed that the program is using the default error
handler (which calls @code{abort} for any errors).  If you are not using
a safe error handler you would need to check the return status of
@code{gsl_fft_complex_radix2_forward}.

The transformed data is rescaled by @math{1/\sqrt N} so that it fits on
the same plot as the input.  Only the real part is shown, by the choice
of the input data the imaginary part is zero.  Allowing for the
wrap-around of negative times at @math{t=128}, and working in units of
@math{k/N}, the DFT approximates the continuum fourier transform, giving
a modulated sine function.
@iftex
@tex
\beforedisplay
$$
\int_{-a}^{+a} e^{-2 \pi i k x} dx = {\sin(2\pi k a) \over\pi k}
$$
\afterdisplay
@end tex

@sp 1
@center @image{fft-complex-radix2-t,2.8in} 
@center @image{fft-complex-radix2-f,2.8in}
@quotation
A pulse and its discrete fourier transform, output from
the example program.
@end quotation
@end iftex

@node Mixed-radix FFT routines for complex data
@section Mixed-radix FFT routines for complex data
@cindex FFT of complex data, mixed-radix algorithm
@cindex Mixed-radix FFT, complex data

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 Paul Swarztrauber's Fortran @sc{fftpack} library.
The theory is explained in the review article @cite{Self-sorting
Mixed-radix FFTs} by Clive Temperton.  The routines here use the same
indexing scheme and basic algorithms as @sc{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 @cite{GSL FFT Algorithms} included in the GSL distribution
if you encounter this problem).

The mixed-radix initialization function @code{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
@cite{GSL FFT Algorithms} for details on adding support for other
factors.

@comment First, the space for the trigonometric lookup tables and scratch area is
@comment allocated by a call to one of the @code{alloc} functions.  We
@comment call the combination of factorization, scratch space and trigonometric
@comment lookup arrays a @dfn{wavetable}.  It contains the sine and cosine
@comment waveforms for the all the frequencies that will be used in the FFT.

@comment The wavetable is initialized by a call to the corresponding @code{init}
@comment function.  It factorizes the data length, using the implemented
@comment subtransforms as preferred factors wherever possible.  The trigonometric
@comment lookup table for the chosen factorization is also computed.

@comment An FFT is computed by a call to one of the @code{forward},
@comment @code{backward} or @code{inverse} functions, with the data, length and
@comment wavetable as arguments.

All the functions described in this section are declared in the header
file @file{gsl_fft_complex.h}.

@deftypefun {gsl_fft_complex_wavetable *} gsl_fft_complex_wavetable_alloc (size_t @var{n})
This function prepares a trigonometric lookup table for a complex FFT of
length @var{n}. The function returns a pointer to the newly allocated
@code{gsl_fft_complex_wavetable} if no errors were detected, and a null
pointer in the case of error.  The length @var{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 @code{sin} and @code{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.
@end deftypefun

@deftypefun void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * @var{wavetable})
This function frees the memory associated with the wavetable
@var{wavetable}.  The wavetable can be freed if no further FFTs of the
same length will be needed.
@end deftypefun

@noindent
These functions operate on a @code{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
@file{gsl_fft_complex.h}.

@deftp {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:

@table @code
@item size_t n
This is the number of complex data points

@item size_t nf
This is the number of factors that the length @code{n} was decomposed into.

@item size_t factor[64]
This is the array of factors.  Only the first @code{nf} elements are
used. 

@comment (FIXME: This is a fixed length array and therefore probably in
@comment violation of the GNU Coding Standards).

@item gsl_complex * trig
This is a pointer to a preallocated trigonometric lookup table of
@code{n} complex elements.

@item gsl_complex * twiddle[64]
This is an array of pointers into @code{trig}, giving the twiddle
factors for each pass.
@end table
@end deftp

@noindent
The mixed radix algorithms require additional working space to hold
the intermediate steps of the transform.

@deftypefun {gsl_fft_complex_workspace *} gsl_fft_complex_workspace_alloc (size_t @var{n})
This function allocates a workspace for a complex transform of length
@var{n}.
@end deftypefun

@deftypefun void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * @var{workspace})
This function frees the memory associated with the workspace
@var{workspace}. The workspace can be freed if no further FFTs of the
same length will be needed.
@end deftypefun

@comment @deftp {Data Type} gsl_fft_complex_workspace
@comment This is a structure that holds the workspace for the mixed radix fft
@comment algorithm.  It has the following components:
@comment
@comment @table @code
@comment @item gsl_complex * scratch
@comment This is a pointer to a workspace of @code{n} complex elements,
@comment capable of holding intermediate copies of the original data set.
@comment @end table
@comment @end deftp

@noindent
The following functions compute the transform,

@deftypefun int gsl_fft_complex_forward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, const gsl_fft_complex_wavetable * @var{wavetable}, gsl_fft_complex_workspace * @var{work})
@deftypefunx int gsl_fft_complex_transform (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, const gsl_fft_complex_wavetable * @var{wavetable}, gsl_fft_complex_workspace * @var{work}, gsl_fft_direction @var{sign})
@deftypefunx int gsl_fft_complex_backward (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, const gsl_fft_complex_wavetable * @var{wavetable}, gsl_fft_complex_workspace * @var{work})
@deftypefunx int gsl_fft_complex_inverse (gsl_complex_packed_array @var{data}, size_t @var{stride}, size_t @var{n}, const gsl_fft_complex_wavetable * @var{wavetable}, gsl_fft_complex_workspace * @var{work})

These functions compute forward, backward and inverse FFTs of length
@var{n} with stride @var{stride}, on the packed complex array
@var{data}, using a mixed radix decimation-in-frequency algorithm.
There is no restriction on the length @var{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 @var{wavetable} containing the
trigonometric lookup tables and a workspace @var{work}.  For the
@code{transform} version of the function the @var{sign} argument can be
either @code{forward} (@math{-1}) or @code{backward} (@math{+1}).

The functions return a value of @code{0} if no errors were detected. The
following @code{gsl_errno} conditions are defined for these functions:

@table @code
@item GSL_EDOM
The length of the data @var{n} is not a positive integer (i.e. @var{n}
is zero).

@item GSL_EINVAL
The length of the data @var{n} and the length used to compute the given
@var{wavetable} do not match.
@end table
@end deftypefun

@comment @node Example of using mixed-radix FFT routines for complex data
@comment @subsection Example of using mixed-radix FFT routines for complex data

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.

@example
@verbatiminclude examples/fftmr.c
@end example

@noindent
Note that we have assumed that the program is using the default
@code{gsl} error handler (which calls @code{abort} for any errors).  If
you are not using a safe error handler you would need to check the
return status of all the @code{gsl} routines.

@node Overview of real data FFTs
@section Overview of real data FFTs
@cindex FFT of real data
The functions for real data are similar to those for complex data.
However, there is an important difference between forward and inverse
transforms.  The fourier transform of a real sequence is not real.  It is
a complex sequence with a special symmetry:
@tex
\beforedisplay
$$
z_k = z_{N-k}^*
$$
\afterdisplay
@end tex
@ifinfo

@example
z_k = z_@{N-k@}^*
@end example

@end ifinfo
@noindent
A sequence with this symmetry is called @dfn{conjugate-complex} or
@dfn{half-complex}.  This different structure requires different
storage layouts for the forward transform (from real to half-complex)
and inverse transform (from half-complex back to real).  As a
consequence the routines are divided into two sets: functions in
@code{gsl_fft_real} which operate on real sequences and functions in
@code{gsl_fft_halfcomplex} which operate on half-complex sequences.

Functions in @code{gsl_fft_real} compute the frequency coefficients of a
real sequence.  The half-complex coefficients @math{c} of a real sequence
@math{x} are given by fourier analysis,
@tex
\beforedisplay
$$
c_k = \sum_{j=0}^{N-1} x_j \exp(-2 \pi i j k /N)
$$
\afterdisplay
@end tex
@ifinfo

@example
c_k = \sum_@{j=0@}^@{N-1@} x_j \exp(-2 \pi i j k /N)
@end example

@end ifinfo
@noindent
Functions in @code{gsl_fft_halfcomplex} compute inverse or backwards
transforms.  They reconstruct real sequences by fourier synthesis from
their half-complex frequency coefficients, @math{c},
@tex
\beforedisplay
$$
x_j = {1 \over N} \sum_{k=0}^{N-1} c_k \exp(2 \pi i j k /N)
$$
\afterdisplay
@end tex
@ifinfo

@example
x_j = @{1 \over N@} \sum_@{k=0@}^@{N-1@} c_k \exp(2 \pi i j k /N)
@end example

@end ifinfo
@noindent
The symmetry of the half-complex sequence implies that only half of the
complex numbers in the output need to be stored.  The remaining half can
be reconstructed using the half-complex symmetry condition. This works
for all lengths, even and odd---when the length is even the middle value
where @math{k=N/2} is also real.  Thus only @var{N} real numbers are
required to store the half-complex sequence, and the transform of a real
sequence can be stored in the same size array as the original data.

The precise storage arrangements depend on the algorithm, and are
different for radix-2 and mixed-radix routines.  The radix-2 function
operates in-place, which constrains the locations where each element can
be stored.  The restriction forces real and imaginary parts to be stored
far apart.  The mixed-radix algorithm does not have this restriction, and
it stores the real and imaginary parts of a given term in neighboring
locations (which is desirable for better locality of memory accesses).

@node Radix-2 FFT routines for real data
@section Radix-2 FFT routines for real data
@cindex FFT of real data, radix-2 algorithm
@cindex Radix-2 FFT 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
@file{gsl_fft_real.h} 

@deftypefun int gsl_fft_real_radix2_transform (double @var{data}[], size_t @var{stride}, size_t @var{n})

This function computes an in-place radix-2 FFT of length @var{n} and
stride @var{stride} on the real array @var{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 @math{k < N/2} the
real part of the @math{k}-th term is stored in location @math{k}, and
the corresponding imaginary part is stored in location @math{N-k}.  Terms
with @math{k > N/2} can be reconstructed using the symmetry 
@c{$z_k = z^*_{N-k}$}
@math{z_k = z^*_@{N-k@}}. 
The terms for @math{k=0} and @math{k=N/2} are both purely
real, and count as a special case.  Their real parts are stored in
locations @math{0} and @math{N/2} respectively, while their imaginary
parts which are zero are not stored.

The following table shows the correspondence between the output
@var{data} and the equivalent results obtained by considering the input
data as a complex sequence with zero imaginary part,

@example
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]
@end example
@noindent
Note that the output data can be converted into the full complex
sequence using the function @code{gsl_fft_halfcomplex_unpack}
described in the next section.

@end deftypefun
The radix-2 FFT functions for halfcomplex data are declared in the
header file @file{gsl_fft_halfcomplex.h}.

@deftypefun int gsl_fft_halfcomplex_radix2_inverse (double @var{data}[], size_t @var{stride}, size_t @var{n})
@deftypefunx int gsl_fft_halfcomplex_radix2_backward (double @var{data}[], size_t @var{stride}, size_t @var{n})

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

@end deftypefun



@node Mixed-radix FFT routines for real data
@section Mixed-radix FFT routines for real data
@cindex FFT of real data, mixed-radix algorithm
@cindex Mixed-radix FFT, real data

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 @sc{fftpack} library
by Paul Swarztrauber.  The theory behind the algorithm is explained in
the article @cite{Fast Mixed-Radix Real Fourier Transforms} by Clive
Temperton.  The routines here use the same indexing scheme and basic
algorithms as @sc{fftpack}.

The functions use the @sc{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 @math{n/2} is not stored
either, since the symmetry 
@c{$z_k = z_{N-k}^*$}
@math{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, @math{n=5}.  The two columns
give the correspondence between the 5 values in the half-complex
sequence returned by @code{gsl_fft_real_transform}, @var{halfcomplex}[] and the
values @var{complex}[] that would be returned if the same real input
sequence were passed to @code{gsl_fft_complex_backward} as a complex
sequence (with imaginary parts set to @code{0}),

@example
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]
@end example

@noindent
The upper elements of the @var{complex} array, @code{complex[3]} and
@code{complex[4]} are filled in using the symmetry condition.  The
imaginary part of the zero-frequency term @code{complex[0].imag} is
known to be zero by the symmetry.

The next table shows the output for an even-length sequence, @math{n=6}
In the even case there are two values which are purely real,

@example
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] 
@end example

@noindent
The upper elements of the @var{complex} array, @code{complex[4]} and
@code{complex[5]} are filled in using the symmetry condition.  Both
@code{complex[0].imag} and @code{complex[3].imag} are known to be zero.

All these functions are declared in the header files
@file{gsl_fft_real.h} and @file{gsl_fft_halfcomplex.h}.

@deftypefun {gsl_fft_real_wavetable *} gsl_fft_real_wavetable_alloc (size_t @var{n})
@deftypefunx {gsl_fft_halfcomplex_wavetable *} gsl_fft_halfcomplex_wavetable_alloc (size_t @var{n})
These functions prepare trigonometric lookup tables for an FFT of size
@math{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 @var{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 @code{sin} and @code{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.
@end deftypefun

@deftypefun void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * @var{wavetable})
@deftypefunx void gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * @var{wavetable})
These functions free the memory associated with the wavetable
@var{wavetable}. The wavetable can be freed if no further FFTs of the
same length will be needed.
@end deftypefun

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

@deftypefun {gsl_fft_real_workspace *} gsl_fft_real_workspace_alloc (size_t @var{n})
This function allocates a workspace for a real transform of length
@var{n}.  The same workspace can be used for both forward real and inverse
halfcomplex transforms.
@end deftypefun

@deftypefun void gsl_fft_real_workspace_free (gsl_fft_real_workspace * @var{workspace})
This function frees the memory associated with the workspace
@var{workspace}. The workspace can be freed if no further FFTs of the
same length will be needed.
@end deftypefun

@noindent
The following functions compute the transforms of real and half-complex
data,

@deftypefun int gsl_fft_real_transform (double @var{data}[], size_t @var{stride}, size_t @var{n}, const gsl_fft_real_wavetable * @var{wavetable}, gsl_fft_real_workspace * @var{work})
@deftypefunx int gsl_fft_halfcomplex_transform (double @var{data}[], size_t @var{stride}, size_t @var{n}, const gsl_fft_halfcomplex_wavetable * @var{wavetable}, gsl_fft_real_workspace * @var{work})
These functions compute the FFT of @var{data}, a real or half-complex
array of length @var{n}, using a mixed radix decimation-in-frequency
algorithm.  For @code{gsl_fft_real_transform} @var{data} is an array of
time-ordered real data.  For @code{gsl_fft_halfcomplex_transform}
@var{data} contains fourier coefficients in the half-complex ordering
described above.  There is no restriction on the length @var{n}.
Efficient modules are provided for subtransforms of length 2, 3, 4 and
5.  Any remaining factors are computed with a slow, @math{O(n^2)},
general-n module.  The caller must supply a @var{wavetable} containing
trigonometric lookup tables and a workspace @var{work}. 
@end deftypefun

@deftypefun int gsl_fft_real_unpack (const double @var{real_coefficient}[], gsl_complex_packed_array @var{complex_coefficient}, size_t @var{stride}, size_t @var{n})

This function converts a single real array, @var{real_coefficient} into
an equivalent complex array, @var{complex_coefficient}, (with imaginary
part set to zero), suitable for @code{gsl_fft_complex} routines.  The
algorithm for the conversion is simply,

@example
for (i = 0; i < n; i++)
  @{
    complex_coefficient[i].real 
      = real_coefficient[i];
    complex_coefficient[i].imag 
      = 0.0;
  @}
@end example
@end deftypefun

@deftypefun int gsl_fft_halfcomplex_unpack (const double @var{halfcomplex_coefficient}[], gsl_complex_packed_array @var{complex_coefficient}, size_t @var{stride}, size_t @var{n})

This function converts @var{halfcomplex_coefficient}, an array of
half-complex coefficients as returned by @code{gsl_fft_real_transform}, into an
ordinary complex array, @var{complex_coefficient}.  It fills in the
complex array using the symmetry 
@c{$z_k = z_{N-k}^*$}
@math{z_k = z_@{N-k@}^*}
to reconstruct the redundant elements.  The algorithm for the conversion
is,

@example  
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];
    double hc_imag 
      = halfcomplex_coefficient[2 * i];
    complex_coefficient[i].real = hc_real;
    complex_coefficient[i].imag = hc_imag;
    complex_coefficient[n - i].real = hc_real;
    complex_coefficient[n - i].imag = -hc_imag;
  @}

if (i == n - i)
  @{
    complex_coefficient[i].real 
      = halfcomplex_coefficient[n - 1];
    complex_coefficient[i].imag 
      = 0.0;
  @}
@end example
@end deftypefun

@comment @node Example of using mixed-radix FFT routines for real data
@comment @subsection Example of using mixed-radix FFT routines for real data

Here is an example program using @code{gsl_fft_real_transform} and
@code{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
@code{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.

@example
@verbatiminclude examples/fftreal.c
@end example

@iftex
@sp 1
@center @image{fft-real-mixedradix,3.4in}

@center Low-pass filtered version of a real pulse, 
@center output from the example program.
@end iftex

@node FFT References and Further Reading
@section References and Further Reading

A good starting point for learning more about the FFT is the review
article @cite{Fast Fourier Transforms: A Tutorial Review and A State of
the Art} by Duhamel and Vetterli,

@itemize @asis
@item
P. Duhamel and M. Vetterli.
Fast fourier transforms: A tutorial review and a state of the art.
@cite{Signal Processing}, 19:259--299, 1990.
@end itemize

@noindent
To find out about the algorithms used in the GSL routines you may want
to consult the document @cite{GSL FFT Algorithms} (it is included
in GSL, as @file{doc/fftalgorithms.tex}).  This has general information
on FFTs and explicit derivations of the implementation for each
routine.  There are also references to the relevant literature.  For
convenience some of the more important references are reproduced below.

@noindent
There are several introductory books on the FFT with example programs,
such as @cite{The Fast Fourier Transform} by Brigham and @cite{DFT/FFT
and Convolution Algorithms} by Burrus and Parks,

@itemize @asis 
@item
E. Oran Brigham.
@cite{The Fast Fourier Transform}.
Prentice Hall, 1974.

@item
C. S. Burrus and T. W. Parks.
@cite{DFT/FFT and Convolution Algorithms}.
Wiley, 1984.
@end itemize

@noindent
Both these introductory books cover the radix-2 FFT in some detail.
The mixed-radix algorithm at the heart of the @sc{fftpack} routines is
reviewed in Clive Temperton's paper,

@itemize @asis 
@item
Clive Temperton.
Self-sorting mixed-radix fast fourier transforms.
@cite{Journal of Computational Physics}, 52(1):1--23, 1983.
@end itemize

@noindent
The derivation of FFTs for real-valued data is explained in the
following two articles,

@itemize @asis
@item
Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney
Burrus.
Real-valued fast fourier transform algorithms.
@cite{IEEE Transactions on Acoustics, Speech, and Signal Processing},
ASSP-35(6):849--863, 1987.

@item
Clive Temperton.
Fast mixed-radix real fourier transforms.
@cite{Journal of Computational Physics}, 52:340--350, 1983.
@end itemize

@noindent
In 1979 the IEEE published a compendium of carefully-reviewed Fortran
FFT programs in @cite{Programs for Digital Signal Processing}.  It is a
useful reference for implementations of many different FFT
algorithms,

@itemize @asis 
@item
Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal
Processing Committee, editors.
@cite{Programs for Digital Signal Processing}.
IEEE Press, 1979.
@end itemize
@comment  @noindent
@comment  There is also an annotated bibliography of papers on the FFT and related
@comment  topics by Burrus,

@comment  @itemize @asis 
@comment  @item C. S. Burrus.  Notes on the FFT.
@comment  @end itemize
@comment  @noindent
@comment The notes are available from @url{http://www-dsp.rice.edu/res/fft/fftnote.asc}.

@noindent
For large-scale FFT work we recommend the use of the dedicated FFTW library
by Frigo and Johnson.  The FFTW library is self-optimizing---it
automatically tunes itself for each hardware platform in order to
achieve maximum performance.  It is available under the GNU GPL.

@itemize @asis
@item
FFTW Website, @uref{http://www.fftw.org/}
@end itemize

@noindent
The source code for @sc{fftpack} is available from Netlib,

@itemize @asis
@item
FFTPACK, @uref{http://www.netlib.org/fftpack/}
@end itemize