fft.texi

开放gsl矩阵运算

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=5}
In the even case both the 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 an 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 is 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
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>

int
main (void)
@{
  int i, n = 100;
  double data[n];

  gsl_fft_real_wavetable * real;
  gsl_fft_halfcomplex_wavetable * hc;
  gsl_fft_real_workspace * work;

  for (i = 0; i < n; i++)
    @{
      data[i] = 0.0;
    @}

  for (i = n / 3; i < 2 * n / 3; i++)
    @{
      data[i] = 1.0;
    @}

  for (i = 0; i < n; i++)
    @{
      printf ("%d: %e\n", i, data[i]);
    @}
  printf ("\n");

  work = gsl_fft_real_workspace_alloc (n);
  real = gsl_fft_real_wavetable_alloc (n);

  gsl_fft_real_transform (data, 1, n, 
                          real, work);

  gsl_fft_real_wavetable_free (real);

  for (i = 11; i < n; i++)
    @{
      data[i] = 0;
    @}

  hc = gsl_fft_halfcomplex_wavetable_alloc (n);

  gsl_fft_halfcomplex_inverse (data, 1, n, 
                               hc, work);
  gsl_fft_halfcomplex_wavetable_free (hc);

  for (i = 0; i < n; i++)
    @{
      printf ("%d: %e\n", i, data[i]);
    @}

  gsl_fft_real_workspace_free (work);
  return 0;
@}
@end example

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

@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

@noindent
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 latex 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 serious 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, @url{http://www.fftw.org/}
@end itemize