|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT
|
|
Packit |
67cb25 |
single: Fast Fourier Transforms, see FFT
|
|
Packit |
67cb25 |
single: Fourier Transforms, see FFT
|
|
Packit |
67cb25 |
single: Discrete Fourier Transforms, see FFT
|
|
Packit |
67cb25 |
single: DFTs, see FFT
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
******************************
|
|
Packit |
67cb25 |
Fast Fourier Transforms (FFTs)
|
|
Packit |
67cb25 |
******************************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: include.rst
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This chapter describes functions for performing Fast Fourier Transforms
|
|
Packit |
67cb25 |
(FFTs). The library includes radix-2 routines (for lengths which are a
|
|
Packit |
67cb25 |
power of two) and mixed-radix routines (which work for any length). For
|
|
Packit |
67cb25 |
efficiency there are separate versions of the routines for real data and
|
|
Packit |
67cb25 |
for complex data. The mixed-radix routines are a reimplementation of the
|
|
Packit |
67cb25 |
|fftpack| library of Paul Swarztrauber. Fortran code for |fftpack| is
|
|
Packit |
67cb25 |
available on Netlib (|fftpack| also includes some routines for sine and
|
|
Packit |
67cb25 |
cosine transforms but these are currently not available in GSL). For
|
|
Packit |
67cb25 |
details and derivations of the underlying algorithms consult the
|
|
Packit |
67cb25 |
document "GSL FFT Algorithms" (see :ref:`References and Further Reading <fft-references>`)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: FFT mathematical definition
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Mathematical Definitions
|
|
Packit |
67cb25 |
========================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Fast Fourier Transforms are efficient algorithms for
|
|
Packit |
67cb25 |
calculating the discrete Fourier transform (DFT),
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_j = \sum_{k=0}^{n-1} z_k \exp(-2 \pi i j k / n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The DFT usually arises as an approximation to the continuous Fourier
|
|
Packit |
67cb25 |
transform when functions are sampled at discrete intervals in space or
|
|
Packit |
67cb25 |
time. The naive evaluation of the discrete Fourier transform is a
|
|
Packit |
67cb25 |
matrix-vector multiplication :math:`W\vec{z}`.
|
|
Packit |
67cb25 |
A general matrix-vector multiplication takes
|
|
Packit |
67cb25 |
:math:`O(n^2)` operations for :math:`n` data-points. Fast Fourier
|
|
Packit |
67cb25 |
transform algorithms use a divide-and-conquer strategy to factorize the
|
|
Packit |
67cb25 |
matrix :math:`W` into smaller sub-matrices, corresponding to the integer
|
|
Packit |
67cb25 |
factors of the length :math:`n`. If :math:`n` can be factorized into a
|
|
Packit |
67cb25 |
product of integers :math:`f_1 f_2 \ldots f_m`
|
|
Packit |
67cb25 |
then the DFT can be computed in :math:`O(n \sum f_i)`
|
|
Packit |
67cb25 |
operations. For a radix-2 FFT this gives an operation count of
|
|
Packit |
67cb25 |
:math:`O(n \log_2 n)`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
All the FFT functions offer three types of transform: forwards, inverse
|
|
Packit |
67cb25 |
and backwards, based on the same mathematical definitions. The
|
|
Packit |
67cb25 |
definition of the *forward Fourier transform*,
|
|
Packit |
67cb25 |
:math:`x = \hbox{FFT}(z)`, is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_j = \sum_{k=0}^{n-1} z_k \exp(-2 \pi i j k / n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
and the definition of the *inverse Fourier transform*,
|
|
Packit |
67cb25 |
:math:`x = \hbox{IFFT}(z)`, is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: z_j = {1 \over n} \sum_{k=0}^{n-1} x_k \exp(2 \pi i j k / n).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The factor of :math:`1/n` makes this a true inverse. For example, a call
|
|
Packit |
67cb25 |
to :func:`gsl_fft_complex_forward` followed by a call to
|
|
Packit |
67cb25 |
:func:`gsl_fft_complex_inverse` should return the original data (within
|
|
Packit |
67cb25 |
numerical errors).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In general there are two possible choices for the sign of the
|
|
Packit |
67cb25 |
exponential in the transform/ inverse-transform pair. GSL follows the
|
|
Packit |
67cb25 |
same convention as |fftpack|, using a negative exponential for the forward
|
|
Packit |
67cb25 |
transform. The advantage of this convention is that the inverse
|
|
Packit |
67cb25 |
transform recreates the original function with simple Fourier
|
|
Packit |
67cb25 |
synthesis. Numerical Recipes uses the opposite convention, a positive
|
|
Packit |
67cb25 |
exponential in the forward transform.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The *backwards FFT* is simply our terminology for an unscaled
|
|
Packit |
67cb25 |
version of the inverse FFT,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: z^{backwards}_j = \sum_{k=0}^{n-1} x_k \exp(2 \pi i j k / n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
When the overall scale of the result is unimportant it is often
|
|
Packit |
67cb25 |
convenient to use the backwards FFT instead of the inverse to save
|
|
Packit |
67cb25 |
unnecessary divisions.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT, complex data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Overview of complex data FFTs
|
|
Packit |
67cb25 |
=============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The inputs and outputs for the complex FFT routines are *packed arrays*
|
|
Packit |
67cb25 |
of floating point numbers. In a packed array the real and
|
|
Packit |
67cb25 |
imaginary parts of each complex number are placed in alternate
|
|
Packit |
67cb25 |
neighboring elements. For example, the following definition of a packed
|
|
Packit |
67cb25 |
array of length 6::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double x[3*2];
|
|
Packit |
67cb25 |
gsl_complex_packed_array data = x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
can be used to hold an array of three complex numbers, :code:`z[3]`, in
|
|
Packit |
67cb25 |
the following way::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
data[0] = Re(z[0])
|
|
Packit |
67cb25 |
data[1] = Im(z[0])
|
|
Packit |
67cb25 |
data[2] = Re(z[1])
|
|
Packit |
67cb25 |
data[3] = Im(z[1])
|
|
Packit |
67cb25 |
data[4] = Re(z[2])
|
|
Packit |
67cb25 |
data[5] = Im(z[2])
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The array indices for the data have the same ordering as those
|
|
Packit |
67cb25 |
in the definition of the DFT---i.e. there are no index transformations
|
|
Packit |
67cb25 |
or permutations of the data.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A *stride* parameter allows the user to perform transforms on the
|
|
Packit |
67cb25 |
elements :code:`z[stride*i]` instead of :code:`z[i]`. A stride greater
|
|
Packit |
67cb25 |
than 1 can be used to take an in-place FFT of the column of a matrix. A
|
|
Packit |
67cb25 |
stride of 1 accesses the array without any additional spacing between
|
|
Packit |
67cb25 |
elements.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
To perform an FFT on a vector argument, such as :code:`gsl_vector_complex * v`,
|
|
Packit |
67cb25 |
use the following definitions (or their equivalents) when calling
|
|
Packit |
67cb25 |
the functions described in this chapter::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_complex_packed_array data = v->data;
|
|
Packit |
67cb25 |
size_t stride = v->stride;
|
|
Packit |
67cb25 |
size_t n = v->size;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For physical applications it is important to remember that the index
|
|
Packit |
67cb25 |
appearing in the DFT does not correspond directly to a physical
|
|
Packit |
67cb25 |
frequency. If the time-step of the DFT is :math:`\Delta` then the
|
|
Packit |
67cb25 |
frequency-domain includes both positive and negative frequencies,
|
|
Packit |
67cb25 |
ranging from :math:`-1/(2\Delta)` through 0 to :math:`+1/(2\Delta)`. The
|
|
Packit |
67cb25 |
positive frequencies are stored from the beginning of the array up to
|
|
Packit |
67cb25 |
the middle, and the negative frequencies are stored backwards from the
|
|
Packit |
67cb25 |
end of the array.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is a table which shows the layout of the array :data:`data`, and the
|
|
Packit |
67cb25 |
correspondence between the time-domain data :math:`z`, and the
|
|
Packit |
67cb25 |
frequency-domain data :math:`x`::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
index z x = FFT(z)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
0 z(t = 0) x(f = 0)
|
|
Packit |
67cb25 |
1 z(t = 1) x(f = 1/(n Delta))
|
|
Packit |
67cb25 |
2 z(t = 2) x(f = 2/(n Delta))
|
|
Packit |
67cb25 |
. ........ ..................
|
|
Packit |
67cb25 |
n/2 z(t = n/2) x(f = +1/(2 Delta),
|
|
Packit |
67cb25 |
-1/(2 Delta))
|
|
Packit |
67cb25 |
. ........ ..................
|
|
Packit |
67cb25 |
n-3 z(t = n-3) x(f = -3/(n Delta))
|
|
Packit |
67cb25 |
n-2 z(t = n-2) x(f = -2/(n Delta))
|
|
Packit |
67cb25 |
n-1 z(t = n-1) x(f = -1/(n Delta))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
When :math:`n` is even the location :math:`n/2` contains the most positive
|
|
Packit |
67cb25 |
and negative frequencies (:math:`+1/(2 \Delta)`, :math:`-1/(2 \Delta)`)
|
|
Packit |
67cb25 |
which are equivalent. If :math:`n` is odd then general structure of the
|
|
Packit |
67cb25 |
table above still applies, but :math:`n/2` does not appear.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT of complex data, radix-2 algorithm
|
|
Packit |
67cb25 |
single: Radix-2 FFT, complex data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Radix-2 FFT routines for complex data
|
|
Packit |
67cb25 |
=====================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The radix-2 algorithms described in this section are simple and compact,
|
|
Packit |
67cb25 |
although not necessarily the most efficient. They use the Cooley-Tukey
|
|
Packit |
67cb25 |
algorithm to compute in-place complex FFTs for lengths which are a power
|
|
Packit |
67cb25 |
of 2---no additional storage is required. The corresponding
|
|
Packit |
67cb25 |
self-sorting mixed-radix routines offer better performance at the
|
|
Packit |
67cb25 |
expense of requiring additional working space.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
All the functions described in this section are declared in the header file :file:`gsl_fft_complex.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_transform (gsl_complex_packed_array data, size_t stride, size_t n, gsl_fft_direction sign)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_backward (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute forward, backward and inverse FFTs of length
|
|
Packit |
67cb25 |
:data:`n` with stride :data:`stride`, on the packed complex array :data:`data`
|
|
Packit |
67cb25 |
using an in-place radix-2 decimation-in-time algorithm. The length of
|
|
Packit |
67cb25 |
the transform :data:`n` is restricted to powers of two. For the
|
|
Packit |
67cb25 |
:code:`transform` version of the function the :data:`sign` argument can be
|
|
Packit |
67cb25 |
either :code:`forward` (:math:`-1`) or :code:`backward` (:math:`+1`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions return a value of :macro:`GSL_SUCCESS` if no errors were
|
|
Packit |
67cb25 |
detected, or :macro:`GSL_EDOM` if the length of the data :data:`n` is not a
|
|
Packit |
67cb25 |
power of two.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_complex_radix2_dif_forward (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_dif_transform (gsl_complex_packed_array data, size_t stride, size_t n, gsl_fft_direction sign)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_dif_backward (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
int gsl_fft_complex_radix2_dif_inverse (gsl_complex_packed_array data, size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These are decimation-in-frequency versions of the radix-2 FFT functions.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is an example program which computes the FFT of a short pulse in a
|
|
Packit |
67cb25 |
sample of length 128. To make the resulting Fourier transform real the
|
|
Packit |
67cb25 |
pulse is defined for equal positive and negative times (:math:`-10 \dots 10`),
|
|
Packit |
67cb25 |
where the negative times wrap around the end of the array.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/fft.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that we have assumed that the program is using the default error
|
|
Packit |
67cb25 |
handler (which calls :func:`abort` for any errors). If you are not using
|
|
Packit |
67cb25 |
a safe error handler you would need to check the return status of
|
|
Packit |
67cb25 |
:func:`gsl_fft_complex_radix2_forward`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The transformed data is rescaled by :math:`1/\sqrt n` so that it fits on
|
|
Packit |
67cb25 |
the same plot as the input. Only the real part is shown, by the choice
|
|
Packit |
67cb25 |
of the input data the imaginary part is zero. Allowing for the
|
|
Packit |
67cb25 |
wrap-around of negative times at :math:`t=128`, and working in units of
|
|
Packit |
67cb25 |
:math:`k/n`, the DFT approximates the continuum Fourier transform, giving
|
|
Packit |
67cb25 |
a modulated sine function.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: \int_{-a}^{+a} e^{-2 \pi i k x} dx = {\sin(2\pi k a) \over\pi k}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The output of the example program is plotted in :numref:`fig_fft-complex-radix2`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _fig_fft-complex-radix2:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. figure:: /images/fft-complex-radix2.png
|
|
Packit |
67cb25 |
:scale: 60%
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A pulse and its discrete Fourier transform, output from
|
|
Packit |
67cb25 |
the example program.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT of complex data, mixed-radix algorithm
|
|
Packit |
67cb25 |
single: Mixed-radix FFT, complex data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Mixed-radix FFT routines for complex data
|
|
Packit |
67cb25 |
=========================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This section describes mixed-radix FFT algorithms for complex data. The
|
|
Packit |
67cb25 |
mixed-radix functions work for FFTs of any length. They are a
|
|
Packit |
67cb25 |
reimplementation of Paul Swarztrauber's Fortran |fftpack| library.
|
|
Packit |
67cb25 |
The theory is explained in the review article "Self-sorting
|
|
Packit |
67cb25 |
Mixed-radix FFTs" by Clive Temperton. The routines here use the same
|
|
Packit |
67cb25 |
indexing scheme and basic algorithms as |fftpack|.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The mixed-radix algorithm is based on sub-transform modules---highly
|
|
Packit |
67cb25 |
optimized small length FFTs which are combined to create larger FFTs.
|
|
Packit |
67cb25 |
There are efficient modules for factors of 2, 3, 4, 5, 6 and 7. The
|
|
Packit |
67cb25 |
modules for the composite factors of 4 and 6 are faster than combining
|
|
Packit |
67cb25 |
the modules for :math:`2*2` and :math:`2*3`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For factors which are not implemented as modules there is a fall-back to
|
|
Packit |
67cb25 |
a general length-:math:`n` module which uses Singleton's method for
|
|
Packit |
67cb25 |
efficiently computing a DFT. This module is :math:`O(n^2)`, and slower
|
|
Packit |
67cb25 |
than a dedicated module would be but works for any length :math:`n`. Of
|
|
Packit |
67cb25 |
course, lengths which use the general length-:math:`n` module will still
|
|
Packit |
67cb25 |
be factorized as much as possible. For example, a length of 143 will be
|
|
Packit |
67cb25 |
factorized into :math:`11*13`. Large prime factors are the worst case
|
|
Packit |
67cb25 |
scenario, e.g. as found in :math:`n=2*3*99991`, and should be avoided
|
|
Packit |
67cb25 |
because their :math:`O(n^2)` scaling will dominate the run-time (consult
|
|
Packit |
67cb25 |
the document "GSL FFT Algorithms" included in the GSL distribution
|
|
Packit |
67cb25 |
if you encounter this problem).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The mixed-radix initialization function :func:`gsl_fft_complex_wavetable_alloc`
|
|
Packit |
67cb25 |
returns the list of factors chosen by the library for a given length
|
|
Packit |
67cb25 |
:math:`n`. It can be used to check how well the length has been
|
|
Packit |
67cb25 |
factorized, and estimate the run-time. To a first approximation the
|
|
Packit |
67cb25 |
run-time scales as :math:`n \sum f_i`, where the :math:`f_i` are the
|
|
Packit |
67cb25 |
factors of :math:`n`. For programs under user control you may wish to
|
|
Packit |
67cb25 |
issue a warning that the transform will be slow when the length is
|
|
Packit |
67cb25 |
poorly factorized. If you frequently encounter data lengths which
|
|
Packit |
67cb25 |
cannot be factorized using the existing small-prime modules consult
|
|
Packit |
67cb25 |
"GSL FFT Algorithms" for details on adding support for other
|
|
Packit |
67cb25 |
factors.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. First, the space for the trigonometric lookup tables and scratch area is
|
|
Packit |
67cb25 |
.. allocated by a call to one of the :code:`alloc` functions. We
|
|
Packit |
67cb25 |
.. call the combination of factorization, scratch space and trigonometric
|
|
Packit |
67cb25 |
.. lookup arrays a *wavetable*. It contains the sine and cosine
|
|
Packit |
67cb25 |
.. waveforms for the all the frequencies that will be used in the FFT.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. The wavetable is initialized by a call to the corresponding :code:`init`
|
|
Packit |
67cb25 |
.. function. It factorizes the data length, using the implemented
|
|
Packit |
67cb25 |
.. subtransforms as preferred factors wherever possible. The trigonometric
|
|
Packit |
67cb25 |
.. lookup table for the chosen factorization is also computed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. An FFT is computed by a call to one of the :code:`forward`,
|
|
Packit |
67cb25 |
.. :code:`backward` or :code:`inverse` functions, with the data, length and
|
|
Packit |
67cb25 |
.. wavetable as arguments.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
All the functions described in this section are declared in the header
|
|
Packit |
67cb25 |
file :file:`gsl_fft_complex.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_fft_complex_wavetable * gsl_fft_complex_wavetable_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function prepares a trigonometric lookup table for a complex FFT of
|
|
Packit |
67cb25 |
length :data:`n`. The function returns a pointer to the newly allocated
|
|
Packit |
67cb25 |
:type:`gsl_fft_complex_wavetable` if no errors were detected, and a null
|
|
Packit |
67cb25 |
pointer in the case of error. The length :data:`n` is factorized into a
|
|
Packit |
67cb25 |
product of subtransforms, and the factors and their trigonometric
|
|
Packit |
67cb25 |
coefficients are stored in the wavetable. The trigonometric coefficients
|
|
Packit |
67cb25 |
are computed using direct calls to :code:`sin` and :code:`cos`, for
|
|
Packit |
67cb25 |
accuracy. Recursion relations could be used to compute the lookup table
|
|
Packit |
67cb25 |
faster, but if an application performs many FFTs of the same length then
|
|
Packit |
67cb25 |
this computation is a one-off overhead which does not affect the final
|
|
Packit |
67cb25 |
throughput.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The wavetable structure can be used repeatedly for any transform of the
|
|
Packit |
67cb25 |
same length. The table is not modified by calls to any of the other FFT
|
|
Packit |
67cb25 |
functions. The same wavetable can be used for both forward and backward
|
|
Packit |
67cb25 |
(or inverse) transforms of a given length.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * wavetable)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the wavetable
|
|
Packit |
67cb25 |
:data:`wavetable`. The wavetable can be freed if no further FFTs of the
|
|
Packit |
67cb25 |
same length will be needed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions operate on a :type:`gsl_fft_complex_wavetable` structure
|
|
Packit |
67cb25 |
which contains internal parameters for the FFT. It is not necessary to
|
|
Packit |
67cb25 |
set any of the components directly but it can sometimes be useful to
|
|
Packit |
67cb25 |
examine them. For example, the chosen factorization of the FFT length
|
|
Packit |
67cb25 |
is given and can be used to provide an estimate of the run-time or
|
|
Packit |
67cb25 |
numerical error. The wavetable structure is declared in the header file
|
|
Packit |
67cb25 |
:file:`gsl_fft_complex.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_fft_complex_wavetable
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a structure that holds the factorization and trigonometric
|
|
Packit |
67cb25 |
lookup tables for the mixed radix fft algorithm. It has the following
|
|
Packit |
67cb25 |
components:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
================================= ==============================================================================================
|
|
Packit |
67cb25 |
:code:`size_t n` This is the number of complex data points
|
|
Packit |
67cb25 |
:code:`size_t nf` This is the number of factors that the length :code:`n` was decomposed into.
|
|
Packit |
67cb25 |
:code:`size_t factor[64]` This is the array of factors. Only the first :code:`nf` elements are used.
|
|
Packit |
67cb25 |
:code:`gsl_complex * trig` This is a pointer to a preallocated trigonometric lookup table of :code:`n` complex elements.
|
|
Packit |
67cb25 |
:code:`gsl_complex * twiddle[64]` This is an array of pointers into :code:`trig`, giving the twiddle factors for each pass.
|
|
Packit |
67cb25 |
================================= ==============================================================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. (FIXME: factor[64] is a fixed length array and therefore probably in
|
|
Packit |
67cb25 |
.. violation of the GNU Coding Standards).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_fft_complex_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The mixed radix algorithms require additional working space to hold
|
|
Packit |
67cb25 |
the intermediate steps of the transform.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_fft_complex_workspace * gsl_fft_complex_workspace_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for a complex transform of length
|
|
Packit |
67cb25 |
:data:`n`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * workspace)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace
|
|
Packit |
67cb25 |
:data:`workspace`. The workspace can be freed if no further FFTs of the
|
|
Packit |
67cb25 |
same length will be needed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions compute the transform,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_complex_forward (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
|
|
Packit |
67cb25 |
int gsl_fft_complex_transform (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work, gsl_fft_direction sign)
|
|
Packit |
67cb25 |
int gsl_fft_complex_backward (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
|
|
Packit |
67cb25 |
int gsl_fft_complex_inverse (gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute forward, backward and inverse FFTs of length
|
|
Packit |
67cb25 |
:data:`n` with stride :data:`stride`, on the packed complex array
|
|
Packit |
67cb25 |
:data:`data`, using a mixed radix decimation-in-frequency algorithm.
|
|
Packit |
67cb25 |
There is no restriction on the length :data:`n`. Efficient modules are
|
|
Packit |
67cb25 |
provided for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining
|
|
Packit |
67cb25 |
factors are computed with a slow, :math:`O(n^2)`, general-:math:`n`
|
|
Packit |
67cb25 |
module. The caller must supply a :data:`wavetable` containing the
|
|
Packit |
67cb25 |
trigonometric lookup tables and a workspace :data:`work`. For the
|
|
Packit |
67cb25 |
:code:`transform` version of the function the :data:`sign` argument can be
|
|
Packit |
67cb25 |
either :code:`forward` (:math:`-1`) or :code:`backward` (:math:`+1`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions return a value of :code:`0` if no errors were detected. The
|
|
Packit |
67cb25 |
following :data:`gsl_errno` conditions are defined for these functions:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
=================================== =========================================================================================================
|
|
Packit |
67cb25 |
:macro:`GSL_EDOM` The length of the data :data:`n` is not a positive integer (i.e. :data:`n` is zero).
|
|
Packit |
67cb25 |
:macro:`GSL_EINVAL` The length of the data :data:`n` and the length used to compute the given :data:`wavetable` do not match.
|
|
Packit |
67cb25 |
=================================== =========================================================================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is an example program which computes the FFT of a short pulse in a
|
|
Packit |
67cb25 |
sample of length 630 (:math:`=2*3*3*5*7`) using the mixed-radix
|
|
Packit |
67cb25 |
algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/fftmr.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that we have assumed that the program is using the default
|
|
Packit |
67cb25 |
:code:`gsl` error handler (which calls :func:`abort` for any errors). If
|
|
Packit |
67cb25 |
you are not using a safe error handler you would need to check the
|
|
Packit |
67cb25 |
return status of all the :code:`gsl` routines.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: FFT of real data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Overview of real data FFTs
|
|
Packit |
67cb25 |
==========================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions for real data are similar to those for complex data.
|
|
Packit |
67cb25 |
However, there is an important difference between forward and inverse
|
|
Packit |
67cb25 |
transforms. The Fourier transform of a real sequence is not real. It is
|
|
Packit |
67cb25 |
a complex sequence with a special symmetry:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: z_k = z_{n-k}^*
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A sequence with this symmetry is called *conjugate-complex* or
|
|
Packit |
67cb25 |
*half-complex*. This different structure requires different
|
|
Packit |
67cb25 |
storage layouts for the forward transform (from real to half-complex)
|
|
Packit |
67cb25 |
and inverse transform (from half-complex back to real). As a
|
|
Packit |
67cb25 |
consequence the routines are divided into two sets: functions in
|
|
Packit |
67cb25 |
:code:`gsl_fft_real` which operate on real sequences and functions in
|
|
Packit |
67cb25 |
:code:`gsl_fft_halfcomplex` which operate on half-complex sequences.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Functions in :code:`gsl_fft_real` compute the frequency coefficients of a
|
|
Packit |
67cb25 |
real sequence. The half-complex coefficients :math:`c` of a real sequence
|
|
Packit |
67cb25 |
:math:`x` are given by Fourier analysis,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: c_k = \sum_{j=0}^{n-1} x_j \exp(-2 \pi i j k /n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Functions in :code:`gsl_fft_halfcomplex` compute inverse or backwards
|
|
Packit |
67cb25 |
transforms. They reconstruct real sequences by Fourier synthesis from
|
|
Packit |
67cb25 |
their half-complex frequency coefficients, :math:`c`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_j = {1 \over n} \sum_{k=0}^{n-1} c_k \exp(2 \pi i j k /n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The symmetry of the half-complex sequence implies that only half of the
|
|
Packit |
67cb25 |
complex numbers in the output need to be stored. The remaining half can
|
|
Packit |
67cb25 |
be reconstructed using the half-complex symmetry condition. This works
|
|
Packit |
67cb25 |
for all lengths, even and odd---when the length is even the middle value
|
|
Packit |
67cb25 |
where :math:`k=n/2` is also real. Thus only :data:`n` real numbers are
|
|
Packit |
67cb25 |
required to store the half-complex sequence, and the transform of a real
|
|
Packit |
67cb25 |
sequence can be stored in the same size array as the original data.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The precise storage arrangements depend on the algorithm, and are
|
|
Packit |
67cb25 |
different for radix-2 and mixed-radix routines. The radix-2 function
|
|
Packit |
67cb25 |
operates in-place, which constrains the locations where each element can
|
|
Packit |
67cb25 |
be stored. The restriction forces real and imaginary parts to be stored
|
|
Packit |
67cb25 |
far apart. The mixed-radix algorithm does not have this restriction, and
|
|
Packit |
67cb25 |
it stores the real and imaginary parts of a given term in neighboring
|
|
Packit |
67cb25 |
locations (which is desirable for better locality of memory accesses).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT of real data, radix-2 algorithm
|
|
Packit |
67cb25 |
single: Radix-2 FFT for real data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Radix-2 FFT routines for real data
|
|
Packit |
67cb25 |
==================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This section describes radix-2 FFT algorithms for real data. They use
|
|
Packit |
67cb25 |
the Cooley-Tukey algorithm to compute in-place FFTs for lengths which
|
|
Packit |
67cb25 |
are a power of 2.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The radix-2 FFT functions for real data are declared in the header files
|
|
Packit |
67cb25 |
:file:`gsl_fft_real.h`
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_real_radix2_transform (double data[], size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function computes an in-place radix-2 FFT of length :data:`n` and
|
|
Packit |
67cb25 |
stride :data:`stride` on the real array :data:`data`. The output is a
|
|
Packit |
67cb25 |
half-complex sequence, which is stored in-place. The arrangement of the
|
|
Packit |
67cb25 |
half-complex terms uses the following scheme: for :math:`k < n/2` the
|
|
Packit |
67cb25 |
real part of the :math:`k`-th term is stored in location :math:`k`, and
|
|
Packit |
67cb25 |
the corresponding imaginary part is stored in location :math:`n-k`. Terms
|
|
Packit |
67cb25 |
with :math:`k > n/2` can be reconstructed using the symmetry
|
|
Packit |
67cb25 |
:math:`z_k = z^*_{n-k}`.
|
|
Packit |
67cb25 |
The terms for :math:`k=0` and :math:`k=n/2` are both purely
|
|
Packit |
67cb25 |
real, and count as a special case. Their real parts are stored in
|
|
Packit |
67cb25 |
locations :math:`0` and :math:`n/2` respectively, while their imaginary
|
|
Packit |
67cb25 |
parts which are zero are not stored.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following table shows the correspondence between the output
|
|
Packit |
67cb25 |
:data:`data` and the equivalent results obtained by considering the input
|
|
Packit |
67cb25 |
data as a complex sequence with zero imaginary part (assuming :data:`stride` = 1})::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
complex[0].real = data[0]
|
|
Packit |
67cb25 |
complex[0].imag = 0
|
|
Packit |
67cb25 |
complex[1].real = data[1]
|
|
Packit |
67cb25 |
complex[1].imag = data[n-1]
|
|
Packit |
67cb25 |
............... ................
|
|
Packit |
67cb25 |
complex[k].real = data[k]
|
|
Packit |
67cb25 |
complex[k].imag = data[n-k]
|
|
Packit |
67cb25 |
............... ................
|
|
Packit |
67cb25 |
complex[n/2].real = data[n/2]
|
|
Packit |
67cb25 |
complex[n/2].imag = 0
|
|
Packit |
67cb25 |
............... ................
|
|
Packit |
67cb25 |
complex[k'].real = data[k] k' = n - k
|
|
Packit |
67cb25 |
complex[k'].imag = -data[n-k]
|
|
Packit |
67cb25 |
............... ................
|
|
Packit |
67cb25 |
complex[n-1].real = data[1]
|
|
Packit |
67cb25 |
complex[n-1].imag = -data[n-1]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that the output data can be converted into the full complex
|
|
Packit |
67cb25 |
sequence using the function :func:`gsl_fft_halfcomplex_radix2_unpack`
|
|
Packit |
67cb25 |
described below.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The radix-2 FFT functions for halfcomplex data are declared in the
|
|
Packit |
67cb25 |
header file :file:`gsl_fft_halfcomplex.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_halfcomplex_radix2_inverse (double data[], size_t stride, size_t n)
|
|
Packit |
67cb25 |
int gsl_fft_halfcomplex_radix2_backward (double data[], size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the inverse or backwards in-place radix-2 FFT of
|
|
Packit |
67cb25 |
length :data:`n` and stride :data:`stride` on the half-complex sequence
|
|
Packit |
67cb25 |
:data:`data` stored according the output scheme used by
|
|
Packit |
67cb25 |
:func:`gsl_fft_real_radix2`. The result is a real array stored in natural
|
|
Packit |
67cb25 |
order.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_halfcomplex_radix2_unpack (const double halfcomplex_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function converts :data:`halfcomplex_coefficient`, an array of
|
|
Packit |
67cb25 |
half-complex coefficients as returned by :func:`gsl_fft_real_radix2_transform`, into an ordinary complex array, :data:`complex_coefficient`. It fills in the
|
|
Packit |
67cb25 |
complex array using the symmetry :math:`z_k = z_{n-k}^*`
|
|
Packit |
67cb25 |
to reconstruct the redundant elements. The algorithm for the conversion
|
|
Packit |
67cb25 |
is::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
complex_coefficient[0].real = halfcomplex_coefficient[0];
|
|
Packit |
67cb25 |
complex_coefficient[0].imag = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < n - i; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double hc_real = halfcomplex_coefficient[i*stride];
|
|
Packit |
67cb25 |
double hc_imag = halfcomplex_coefficient[(n-i)*stride];
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].real = hc_real;
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].imag = hc_imag;
|
|
Packit |
67cb25 |
complex_coefficient[(n - i)*stride].real = hc_real;
|
|
Packit |
67cb25 |
complex_coefficient[(n - i)*stride].imag = -hc_imag;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == n - i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].real = halfcomplex_coefficient[(n - 1)*stride];
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].imag = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: FFT of real data, mixed-radix algorithm
|
|
Packit |
67cb25 |
single: Mixed-radix FFT, real data
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Mixed-radix FFT routines for real data
|
|
Packit |
67cb25 |
======================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This section describes mixed-radix FFT algorithms for real data. The
|
|
Packit |
67cb25 |
mixed-radix functions work for FFTs of any length. They are a
|
|
Packit |
67cb25 |
reimplementation of the real-FFT routines in the Fortran |fftpack| library
|
|
Packit |
67cb25 |
by Paul Swarztrauber. The theory behind the algorithm is explained in
|
|
Packit |
67cb25 |
the article "Fast Mixed-Radix Real Fourier Transforms" by Clive
|
|
Packit |
67cb25 |
Temperton. The routines here use the same indexing scheme and basic
|
|
Packit |
67cb25 |
algorithms as |fftpack|.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions use the |fftpack| storage convention for half-complex
|
|
Packit |
67cb25 |
sequences. In this convention the half-complex transform of a real
|
|
Packit |
67cb25 |
sequence is stored with frequencies in increasing order, starting at
|
|
Packit |
67cb25 |
zero, with the real and imaginary parts of each frequency in neighboring
|
|
Packit |
67cb25 |
locations. When a value is known to be real the imaginary part is not
|
|
Packit |
67cb25 |
stored. The imaginary part of the zero-frequency component is never
|
|
Packit |
67cb25 |
stored. It is known to be zero (since the zero frequency component is
|
|
Packit |
67cb25 |
simply the sum of the input data (all real)). For a sequence of even
|
|
Packit |
67cb25 |
length the imaginary part of the frequency :math:`n/2` is not stored
|
|
Packit |
67cb25 |
either, since the symmetry :math:`z_k = z_{n-k}^*`
|
|
Packit |
67cb25 |
implies that this is purely real too.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The storage scheme is best shown by some examples. The table below
|
|
Packit |
67cb25 |
shows the output for an odd-length sequence, :math:`n=5`. The two columns
|
|
Packit |
67cb25 |
give the correspondence between the 5 values in the half-complex
|
|
Packit |
67cb25 |
sequence returned by :func:`gsl_fft_real_transform`, :data:`halfcomplex[]` and the
|
|
Packit |
67cb25 |
values :data:`complex[]` that would be returned if the same real input
|
|
Packit |
67cb25 |
sequence were passed to :func:`gsl_fft_complex_backward` as a complex
|
|
Packit |
67cb25 |
sequence (with imaginary parts set to :code:`0`)::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
complex[0].real = halfcomplex[0]
|
|
Packit |
67cb25 |
complex[0].imag = 0
|
|
Packit |
67cb25 |
complex[1].real = halfcomplex[1]
|
|
Packit |
67cb25 |
complex[1].imag = halfcomplex[2]
|
|
Packit |
67cb25 |
complex[2].real = halfcomplex[3]
|
|
Packit |
67cb25 |
complex[2].imag = halfcomplex[4]
|
|
Packit |
67cb25 |
complex[3].real = halfcomplex[3]
|
|
Packit |
67cb25 |
complex[3].imag = -halfcomplex[4]
|
|
Packit |
67cb25 |
complex[4].real = halfcomplex[1]
|
|
Packit |
67cb25 |
complex[4].imag = -halfcomplex[2]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The upper elements of the :data:`complex` array, :code:`complex[3]` and
|
|
Packit |
67cb25 |
:code:`complex[4]` are filled in using the symmetry condition. The
|
|
Packit |
67cb25 |
imaginary part of the zero-frequency term :code:`complex[0].imag` is
|
|
Packit |
67cb25 |
known to be zero by the symmetry.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The next table shows the output for an even-length sequence, :math:`n=6`.
|
|
Packit |
67cb25 |
In the even case there are two values which are purely real::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
complex[0].real = halfcomplex[0]
|
|
Packit |
67cb25 |
complex[0].imag = 0
|
|
Packit |
67cb25 |
complex[1].real = halfcomplex[1]
|
|
Packit |
67cb25 |
complex[1].imag = halfcomplex[2]
|
|
Packit |
67cb25 |
complex[2].real = halfcomplex[3]
|
|
Packit |
67cb25 |
complex[2].imag = halfcomplex[4]
|
|
Packit |
67cb25 |
complex[3].real = halfcomplex[5]
|
|
Packit |
67cb25 |
complex[3].imag = 0
|
|
Packit |
67cb25 |
complex[4].real = halfcomplex[3]
|
|
Packit |
67cb25 |
complex[4].imag = -halfcomplex[4]
|
|
Packit |
67cb25 |
complex[5].real = halfcomplex[1]
|
|
Packit |
67cb25 |
complex[5].imag = -halfcomplex[2]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The upper elements of the :data:`complex` array, :code:`complex[4]` and
|
|
Packit |
67cb25 |
:code:`complex[5]` are filled in using the symmetry condition. Both
|
|
Packit |
67cb25 |
:code:`complex[0].imag` and :code:`complex[3].imag` are known to be zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
All these functions are declared in the header files
|
|
Packit |
67cb25 |
:file:`gsl_fft_real.h` and :file:`gsl_fft_halfcomplex.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_fft_real_wavetable
|
|
Packit |
67cb25 |
gsl_fft_halfcomplex_wavetable
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These data structures contain lookup tables for an FFT of a fixed size.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_fft_real_wavetable * gsl_fft_real_wavetable_alloc (size_t n)
|
|
Packit |
67cb25 |
gsl_fft_halfcomplex_wavetable * gsl_fft_halfcomplex_wavetable_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions prepare trigonometric lookup tables for an FFT of size
|
|
Packit |
67cb25 |
:math:`n` real elements. The functions return a pointer to the newly
|
|
Packit |
67cb25 |
allocated struct if no errors were detected, and a null pointer in the
|
|
Packit |
67cb25 |
case of error. The length :data:`n` is factorized into a product of
|
|
Packit |
67cb25 |
subtransforms, and the factors and their trigonometric coefficients are
|
|
Packit |
67cb25 |
stored in the wavetable. The trigonometric coefficients are computed
|
|
Packit |
67cb25 |
using direct calls to :code:`sin` and :code:`cos`, for accuracy.
|
|
Packit |
67cb25 |
Recursion relations could be used to compute the lookup table faster,
|
|
Packit |
67cb25 |
but if an application performs many FFTs of the same length then
|
|
Packit |
67cb25 |
computing the wavetable is a one-off overhead which does not affect the
|
|
Packit |
67cb25 |
final throughput.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The wavetable structure can be used repeatedly for any transform of the
|
|
Packit |
67cb25 |
same length. The table is not modified by calls to any of the other FFT
|
|
Packit |
67cb25 |
functions. The appropriate type of wavetable must be used for forward
|
|
Packit |
67cb25 |
real or inverse half-complex transforms.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * wavetable)
|
|
Packit |
67cb25 |
void gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * wavetable)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions free the memory associated with the wavetable
|
|
Packit |
67cb25 |
:data:`wavetable`. The wavetable can be freed if no further FFTs of the
|
|
Packit |
67cb25 |
same length will be needed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The mixed radix algorithms require additional working space to hold
|
|
Packit |
67cb25 |
the intermediate steps of the transform,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_fft_real_workspace
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This workspace contains parameters needed to compute a real FFT.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_fft_real_workspace * gsl_fft_real_workspace_alloc (size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function allocates a workspace for a real transform of length
|
|
Packit |
67cb25 |
:data:`n`. The same workspace can be used for both forward real and inverse
|
|
Packit |
67cb25 |
halfcomplex transforms.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_fft_real_workspace_free (gsl_fft_real_workspace * workspace)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees the memory associated with the workspace
|
|
Packit |
67cb25 |
:data:`workspace`. The workspace can be freed if no further FFTs of the
|
|
Packit |
67cb25 |
same length will be needed.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions compute the transforms of real and half-complex
|
|
Packit |
67cb25 |
data,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. 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)
|
|
Packit |
67cb25 |
int gsl_fft_halfcomplex_transform (double data[], size_t stride, size_t n, const gsl_fft_halfcomplex_wavetable * wavetable, gsl_fft_real_workspace * work)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions compute the FFT of :data:`data`, a real or half-complex
|
|
Packit |
67cb25 |
array of length :data:`n`, using a mixed radix decimation-in-frequency
|
|
Packit |
67cb25 |
algorithm. For :func:`gsl_fft_real_transform` :data:`data` is an array of
|
|
Packit |
67cb25 |
time-ordered real data. For :func:`gsl_fft_halfcomplex_transform`
|
|
Packit |
67cb25 |
:data:`data` contains Fourier coefficients in the half-complex ordering
|
|
Packit |
67cb25 |
described above. There is no restriction on the length :data:`n`.
|
|
Packit |
67cb25 |
Efficient modules are provided for subtransforms of length 2, 3, 4 and
|
|
Packit |
67cb25 |
5. Any remaining factors are computed with a slow, :math:`O(n^2)`,
|
|
Packit |
67cb25 |
general-n module. The caller must supply a :data:`wavetable` containing
|
|
Packit |
67cb25 |
trigonometric lookup tables and a workspace :data:`work`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_real_unpack (const double real_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function converts a single real array, :data:`real_coefficient` into
|
|
Packit |
67cb25 |
an equivalent complex array, :data:`complex_coefficient`, (with imaginary
|
|
Packit |
67cb25 |
part set to zero), suitable for :code:`gsl_fft_complex` routines. The
|
|
Packit |
67cb25 |
algorithm for the conversion is simply::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i < n; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].real = real_coefficient[i*stride];
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].imag = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[], gsl_complex_packed_array complex_coefficient, size_t stride, size_t n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function converts :data:`halfcomplex_coefficient`, an array of
|
|
Packit |
67cb25 |
half-complex coefficients as returned by :func:`gsl_fft_real_transform`, into an
|
|
Packit |
67cb25 |
ordinary complex array, :data:`complex_coefficient`. It fills in the
|
|
Packit |
67cb25 |
complex array using the symmetry :math:`z_k = z_{n-k}^*`
|
|
Packit |
67cb25 |
to reconstruct the redundant elements. The algorithm for the conversion
|
|
Packit |
67cb25 |
is::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
complex_coefficient[0].real = halfcomplex_coefficient[0];
|
|
Packit |
67cb25 |
complex_coefficient[0].imag = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 1; i < n - i; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double hc_real = halfcomplex_coefficient[(2 * i - 1)*stride];
|
|
Packit |
67cb25 |
double hc_imag = halfcomplex_coefficient[(2 * i)*stride];
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].real = hc_real;
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].imag = hc_imag;
|
|
Packit |
67cb25 |
complex_coefficient[(n - i)*stride].real = hc_real;
|
|
Packit |
67cb25 |
complex_coefficient[(n - i)*stride].imag = -hc_imag;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (i == n - i)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].real = halfcomplex_coefficient[(n - 1)*stride];
|
|
Packit |
67cb25 |
complex_coefficient[i*stride].imag = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is an example program using :func:`gsl_fft_real_transform` and
|
|
Packit |
67cb25 |
:func:`gsl_fft_halfcomplex_inverse`. It generates a real signal in the
|
|
Packit |
67cb25 |
shape of a square pulse. The pulse is Fourier transformed to frequency
|
|
Packit |
67cb25 |
space, and all but the lowest ten frequency components are removed from
|
|
Packit |
67cb25 |
the array of Fourier coefficients returned by
|
|
Packit |
67cb25 |
:func:`gsl_fft_real_transform`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The remaining Fourier coefficients are transformed back to the
|
|
Packit |
67cb25 |
time-domain, to give a filtered version of the square pulse. Since
|
|
Packit |
67cb25 |
Fourier coefficients are stored using the half-complex symmetry both
|
|
Packit |
67cb25 |
positive and negative frequencies are removed and the final filtered
|
|
Packit |
67cb25 |
signal is also real.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/fftreal.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The program output is shown in :numref:`fig_fft-real-mixedradix`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _fig_fft-real-mixedradix:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. figure:: /images/fft-real-mixedradix.png
|
|
Packit |
67cb25 |
:scale: 100%
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Low-pass filtered version of a real pulse, output from the example program.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. _fft-references:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A good starting point for learning more about the FFT is the following review
|
|
Packit |
67cb25 |
article,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. Duhamel and M. Vetterli.
|
|
Packit |
67cb25 |
Fast Fourier transforms: A tutorial review and a state of the art.
|
|
Packit |
67cb25 |
Signal Processing, 19:259--299, 1990.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
To find out about the algorithms used in the GSL routines you may want
|
|
Packit |
67cb25 |
to consult the document "GSL FFT Algorithms" (it is included
|
|
Packit |
67cb25 |
in GSL, as :file:`doc/fftalgorithms.tex`). This has general information
|
|
Packit |
67cb25 |
on FFTs and explicit derivations of the implementation for each
|
|
Packit |
67cb25 |
routine. There are also references to the relevant literature. For
|
|
Packit |
67cb25 |
convenience some of the more important references are reproduced below.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
There are several introductory books on the FFT with example programs,
|
|
Packit |
67cb25 |
such as "The Fast Fourier Transform" by Brigham and "DFT/FFT
|
|
Packit |
67cb25 |
and Convolution Algorithms" by Burrus and Parks,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* E. Oran Brigham. "The Fast Fourier Transform". Prentice Hall, 1974.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* C. S. Burrus and T. W. Parks. "DFT/FFT and Convolution Algorithms",
|
|
Packit |
67cb25 |
Wiley, 1984.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Both these introductory books cover the radix-2 FFT in some detail.
|
|
Packit |
67cb25 |
The mixed-radix algorithm at the heart of the |fftpack| routines is
|
|
Packit |
67cb25 |
reviewed in Clive Temperton's paper,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Clive Temperton, Self-sorting mixed-radix fast Fourier transforms,
|
|
Packit |
67cb25 |
Journal of Computational Physics, 52(1):1--23, 1983.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The derivation of FFTs for real-valued data is explained in the
|
|
Packit |
67cb25 |
following two articles,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney
|
|
Packit |
67cb25 |
Burrus. Real-valued fast Fourier transform algorithms.
|
|
Packit |
67cb25 |
"IEEE Transactions on Acoustics, Speech, and Signal Processing",
|
|
Packit |
67cb25 |
ASSP-35(6):849--863, 1987.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Clive Temperton. Fast mixed-radix real Fourier transforms.
|
|
Packit |
67cb25 |
"Journal of Computational Physics", 52:340--350, 1983.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In 1979 the IEEE published a compendium of carefully-reviewed Fortran
|
|
Packit |
67cb25 |
FFT programs in "Programs for Digital Signal Processing". It is a
|
|
Packit |
67cb25 |
useful reference for implementations of many different FFT
|
|
Packit |
67cb25 |
algorithms,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal
|
|
Packit |
67cb25 |
Processing Committee, editors.
|
|
Packit |
67cb25 |
Programs for Digital Signal Processing. IEEE Press, 1979.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For large-scale FFT work we recommend the use of the dedicated FFTW library
|
|
Packit |
67cb25 |
by Frigo and Johnson. The FFTW library is self-optimizing---it
|
|
Packit |
67cb25 |
automatically tunes itself for each hardware platform in order to
|
|
Packit |
67cb25 |
achieve maximum performance. It is available under the GNU GPL.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* FFTW Website, http://www.fftw.org/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The source code for |fftpack| is available from http://www.netlib.org/fftpack/
|