fft.texi

开放gsl矩阵运算

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 these functions 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 an 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})
@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}.

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
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
@{
  int i;
  const int n = 630;
  double data[2*n];

  gsl_fft_complex_wavetable * wavetable;
  gsl_fft_complex_workspace * workspace;

  for (i = 0; i < n; i++)
    @{
      REAL(data,i) = 0.0;
      IMAG(data,i) = 0.0;
    @}

  data[0].real = 1.0;

  for (i = 1; i <= 10; i++)
    @{
      REAL(data,i) = REAL(data,n-i) = 1.0;
    @}

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

  wavetable = gsl_fft_complex_wavetable_alloc (n);
  workspace = gsl_fft_complex_workspace_alloc (n);

  for (i = 0; i < wavetable->nf; i++)
    @{
       printf("# factor %d: %d\n", i, 
              wavetable->factor[i]);
    @}

  gsl_fft_complex_forward (data, 1, n, 
                           wavetable, workspace);

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

  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);
  return 0;
@}
@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_k \exp(-2 \pi i j k /N)
$$
\afterdisplay
@end tex
@ifinfo
@example
c_k = \sum_@{j=0@}^@{N-1@} x_k \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 constrain 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.  This 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].real  =    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

@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