fft.texi

开放gsl矩阵运算

@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 by 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
gsl_complex_packed_array data[6];
@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
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.

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

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

@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.

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})

@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
#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;
  double data[2*128];

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

  REAL(data,0) = 1.0;

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

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

  gsl_fft_complex_radix2_forward (data, 1, 128);

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

  return 0;
@}
@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 @math{\sin} 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,3.4in} 
@center @image{fft-complex-radix2-f,3.4in}

@display
A pulse and its discrete fourier transform, output from
the example program.
@end display
@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 the Fortran @sc{fftpack} library by Paul
Swarztrauber.  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