This is Info file gsl-ref.info, produced by Makeinfo version 1.67 from
the input file gsl-ref.texi.
INFO-DIR-SECTION Scientific software
START-INFO-DIR-ENTRY
* gsl-ref: (gsl-ref). GNU Scientific Library - Reference
END-INFO-DIR-ENTRY
This file documents the *GNU Scientific Library*.
Copyright (C) 1996, 1997, 1998 The GSL Project.
Permission is granted to make and distribute verbatim copies of this
manual provided the copyright notice and this permission notice are
preserved on all copies.
Permission is granted to copy and distribute modified versions of
this manual under the conditions for verbatim copying, provided that
the entire resulting derived work is distributed under the terms of a
permission notice identical to this one.
Permission is granted to copy and distribute translations of this
manual into another language, under the above conditions for modified
versions, except that this permission notice may be stated in a
translation approved by the Foundation.
File: gsl-ref.info, Node: Example of using mixed-radix FFT routines for complex data, Prev: Mixed-radix FFT routines for complex data, Up: FFTs of complex data
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 (=2*3*3*5*7) using the mixed-radix algorithm.
#include <stdio.h>
#include <math.h>
#include <gsl_errno.h>
#include <gsl_fft_complex.h>
int main ()
{
int i;
const int n = 630 ;
complex data[n];
gsl_fft_complex_wavetable wavetable;
for (i = 0; i < n; i++)
{
data[i].real = 0.0 ;
data[i].imag = 0.0 ;
}
data[0].real = 1.0 ;
for (i = 1; i <= 10; i++)
{
data[i].real = data[n-i].real = 1.0 ;
}
for (i = 0; i < n; i++)
{
printf ("%d: %e %e\n", i, data[i].real, data[i].imag);
}
printf ("\n");
gsl_fft_complex_wavetable_alloc (n, &wavetable);
gsl_fft_complex_init (n, &wavetable);
for (i = 0; i < wavetable.nf; i++)
{
printf("# factor %d: %d\n", i, wavetable.factor[i]);
}
gsl_fft_complex_forward (data, n, &wavetable);
for (i = 0; i < n; i++)
{
printf ("%d: %e %e\n", i, data[i].real, data[i].imag);
}
gsl_fft_complex_wavetable_free (&wavetable);
}
Note that we have assumed that the program is using the default `gsl'
error handler (which calls `abort' for any errors). If you are not
using a safe error handler you would need to check the return status of
of all the `gsl' routines.
File: gsl-ref.info, Node: FFTs of real data, Next: Test Signals, Prev: FFTs of complex data, Up: FFTs
FFTs 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:
z_k = z_{N-k}^*
A sequence with this symmetry is called "conjugate-complex" or
"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
`gsl_fft_real' which operate on real sequences and functions in
`gsl_fft_halfcomplex' which operate on half-complex sequences.
Functions in `gsl_fft_real' compute the frequency coefficients of a
real sequence. The half-complex coefficients c of a real sequence x are
given by fourier analysis,
c_k = \sum_{j=0}^{N-1} x_k \exp(-2 \pi i j k /N)
Functions in `gsl_fft_halfcomplex' compute inverse or backwards
transforms. They reconstruct real sequences by fourier synthesis from
their half-complex frequency coefficients, c,
x_j = {1 \over N} \sum_{k=0}^{N-1} c_k \exp(2 \pi i j k /N)
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 k=N/2, is also real). Thus only 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.
* Menu:
* Radix-2 FFT routines for real data::
* Mixed-radix FFT routines for real data::
* Example of using mixed-radix FFT routines for real data::
File: gsl-ref.info, Node: Radix-2 FFT routines for real data, Next: Mixed-radix FFT routines for real data, Up: FFTs of real data
Radix-2 FFT routines 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 `gsl_fft_real.h'
- Function: int gsl_fft_real_radix2 (double data[], const unsigned int
n)
This function computes an in-place radix-2 FFT of DATA, a real
array of length N. The output is a half-complex sequence, which is
stored in-place. The arrangement of the half-complex terms uses the
following scheme: for k < N/2 the real part of the k-th term is
stored in location k, and the corresponding imaginary part is
stored in location N-k. Terms with k > N/2 can be reconstructed
using the symmetry z_k = z^*_{N-k}. The terms for k=0 and k=N/2
are both purely real, and count as a special case. Their real
parts are stored in locations 0 and N/2 respectively, while their
imaginary parts which are zero are not stored.
The following table shows the correspondence between the output
DATA and the equivalent results obtained by considering the input
data as a complex sequence with zero imaginary part,
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]
The radix-2 FFT functions for halfcomplex data are declared in the
header file `gsl_fft_halfcomplex.h'.
- Function: int gsl_fft_halfcomplex_radix2_inverse (double data[],
const unsigned int n)
- Function: int gsl_fft_halfcomplex_radix2_backwards (double data[],
const unsigned int n)
These functions compute the inverse or backwards in-place radix-2
FFT of the half-complex sequence DATA, a real of length N stored
according the output scheme used by `gsl_fft_real_radix2'. The
result is a real array stored in natural order.
File: gsl-ref.info, Node: Mixed-radix FFT routines for real data, Next: Example of using mixed-radix FFT routines for real data, Prev: Radix-2 FFT routines for real data, Up: FFTs of real data
Mixed-radix FFT routines for 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 FFTPACK library
by Paul Swarztrauber. The theory behind the algorithm is explained in
the article `Fast Mixed-Radix Real Fourier Transforms' by Clive
Temperton. The routines here use the same indexing scheme and basic
algorithms as FFTPACK.
The functions use the FFTPACK storage convention for half-complex
sequences. In this convention the half-complex transform of a real
sequence is stored in 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 n/2 is not stored either,
since the symmetry 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, n=5. The two columns give
the correspondence between the 5 values in the half-complex sequence
returned by `gsl_fft_real', HALFCOMPLEX[] and the values COMPLEX[] that
would be returned if the same real input sequence were passed to
`gsl_fft_complex_backward' as a complex sequence (with imaginary parts
set to `0'),
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]
The upper elements of the COMPLEX array, `complex[3]' and `complex[4]'
are filled in using the symmetry condition. The imaginary part of the
zero-frequency term `complex[0].imag' is known to be zero by the
symmetry.
The next table shows the output for an even-length sequence, n=5 In
the even case both the there are two values which are purely real,
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]
The upper elements of the COMPLEX array, `complex[4]' and `complex[5]'
are be filled in using the symmetry condition. Both `complex[0].imag'
and `complex[3].imag' are known to be zero.
All these functions are declared in the header files
`gsl_fft_real.h' and `gsl_fft_halfcomplex.h'.
- Function: int gsl_fft_real_wavetable_alloc (unsigned int N,
gsl_fft_real_wavetable * WAVETABLE);
- Function: int gsl_fft_halfcomplex_wavetable_alloc (unsigned int N,
gsl_fft_halfcomplex_wavetable * WAVETABLE);
These functions allocate space for a scratch area of size N real
elements and a trigonometric lookup table, of size n/2 complex
elements. The relevant components of WAVETABLE are modified to
point to the newly allocated memory.
These functions return a value of `0' if no errors were detected,
and `-1' in the case of error. The following `gsl_errno'
conditions are defined for these functions:
`GSL_EDOM'
The length of the data N is not a positive integer (i.e. N is
zero).
`GSL_ENOMEM'
The requested memory could not be allocated (`malloc' returned
a null pointer).
- Function: int gsl_fft_real_init (const unsigned int N,
gsl_fft_real_wavetable * WAVETABLE);
- Function: int gsl_fft_halfcomplex_init (const unsigned int N,
gsl_fft_halfcomplex_wavetable * WAVETABLE);
These functions initialize WAVETABLE. They first select a
factorization of the length N into the implemented subtransforms,
storing the details of the factorization in WAVETABLE.
Using this factorization they then prepare a trigonometric lookup
table in the memory previously allocated by
`gsl_fft_real_wavetable_alloc' or
`gsl_fft_halfcomplex_wavetable_alloc'. The wavetable is computed
using direct calls to `sin' and `cos', for accuracy. It could be
computed faster using recursion relations if necessary. 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 same wavetable can be used for both
forward and backward (or inverse) transforms of a given length.
The functions return a value of `0' if no errors were detected,
and `-1' in the case of error. The following `gsl_errno'
conditions are defined for these functions:
`GSL_EDOM'
The length of the data N is not a positive integer (i.e. N is
zero).
`GSL_EFACTOR'
The length N could not be factorized (this shouldn't happen).
`GSL_EFAILED'
A failure was detected in the wavetable generation. This
could be caused by an inconsistency in a user-supplied
WAVETABLE structure.
- Function: int gsl_fft_real_wavetable_free (gsl_fft_real_wavetable *
WAVETABLE);
- Function: int gsl_fft_halfcomplex_wavetable_free
(gsl_fft_halfcomplex_wavetable * WAVETABLE);
These functions free the blocks of memory previously allocated by
`gsl_fft_real_wavetable_alloc' or
`gsl_fft_halfcomplex_wavetable_alloc' and pointed to by the
components of WAVETABLE.
The wavetable should be freed if no further FFTs of the same
length will be needed.
- Function: int gsl_fft_real (double DATA[], const unsigned int N,
const gsl_fft_real_wavetable * WAVETABLE)
- Function: int gsl_fft_halfcomplex (double DATA[], const unsigned int
N, const gsl_fft_halfcomplex_wavetable * WAVETABLE)
These functions compute the FFT of DATA, a real or half-complex
array of length N, using a mixed radix decimation-in-frequency
algorithm. For `gsl_fft_real' DATA is an array of time-ordered
real data. For `gsl_fft_halfcomplex' DATA contains fourier
coefficients in the half-complex ordering described above.
The caller must supply a WAVETABLE containing the chosen
factorization, trigonometric lookup tables and scratch area. The
wavetable can be easily prepared using the functions
`gsl_fft_real_alloc' and `gsl_fft_real_init' or
`gsl_fft_halfcomplex_alloc' and `gsl_fft_halfcomplex_init'.
There is no restriction on the length N. Efficient modules are
provided for subtransforms of length 2, 3, 4 and 5. Any remaining
factors are computed with a slow, O(n^2), general-n module.
The functions return a value of `0' if no errors were detected,
and `-1' in the case of error. The following `gsl_errno'
conditions are defined for these functions:
`GSL_EDOM'
The length of the data N is not a positive integer (i.e. N is
zero).
`GSL_EINVAL'
The length of the data N and the length used to compute the
given WAVETABLE do not match.
- Function: int gsl_fft_real_unpack (const double real_coefficient[],
complex complex_coefficient[], const unsigned int n)
This function converts a single real array, REAL_COEFFICIENT into
an equivalent complex array, COMPLEX_COEFFICIENT, (with imaginary
part set to zero), suitable for `gsl_fft_complex' routines. The
This chapter describes functions for finding a root of an arbitrary
function which you provide. It discusses proper use and possible
pitfalls of each function and gives an overview of the algorithms
involved.
The *GNU Scientific Library* provides two types of root finding
functions: high-level functions which take few arguments and make most
decisions for you, and low-level functions which take many arguments
and can be precisely controlled. For most tasks, the simpler high-level
functions will be sufficient, but the greater degree of control
provided by the low-level functions is sometimes useful.
For examples of GSL root finding functions in action, look at the
source code for the testing program which verifies that the they work
correctly (it is run when you do a `make check'). It is in the file
`roots/test.c' and implements one way of properly handling errors.
The header file `gsl_roots.h' contains prototypes for GSL root
finding functions and some other related declarations.
*FIXME: add a note about performance.*
* Menu:
* Root Finding Overview::
* High Level vs. Low Level::
* High Level Root Finding Functions::
* Low Level Root Finding Functions::
* Root Finder Error Handling::
File: gsl-ref.info, Node: Root Finding Overview, Next: High Level vs. Low Level, Up: Root finding
Root Finding Overview
=====================
*FIXME: Insert discussion of numerical root finding here.*
GSL root finding functions operate on continous functions only.
While it is not absolutely required that f have a root within the
search region, numerical root finding functions should not be used
haphazardly to check for the *existence* of roots. There are better
ways to do this! Because it is so easy to create situations where
numerical root finders go awry, it is a bad idea to throw a root finder
at a function you do not know much about.
Note that GSL root finding functions can only search for one root at
a time. While their behavior in cases of multiple roots (several roots
in the search area, not one root of degree greater than 1) is
deterministic, it depends entirely on the algorithm and the function
under search and is rather complex; it is therefore difficult to
predict. *In most cases, no error will be reported if you try to find a
root in an area where there is more than one.*
File: gsl-ref.info, Node: High Level vs. Low Level, Next: High Level Root Finding Functions, Prev: Root Finding Overview, Up: Root finding
High Level vs. Low Level
========================
*FIXME: fill in this section.*
In general, you should use high level functions unless there is a
definite reason to use low level functions.
File: gsl-ref.info, Node: High Level Root Finding Functions, Next: Low Level Root Finding Functions, Prev: High Level vs. Low Level, Up: Root finding
High Level Root Finding Functions
=================================
GSL provides three high level root finding functions. They are
simple to use, robust, and adequate for most tasks, but certain
situations require more control than the high level root finding
functions provide.
* Menu:
* Root Finder Exit Values:: How to tell if an error occured and
how the root location is returned
* Providing the Function to Search:: How to provide a function for the
root finder to operate on
* Automatic Control Decisions:: How high level root finders set
parameters such as tolerance
* High Level Functions:: Names and arguments of the high
level root finding functions
File: gsl-ref.info, Node: Root Finder Exit Values, Next: Providing the Function to Search, Up: High Level Root Finding Functions
Root Finder Exit Values
-----------------------
Since the return value of GSL root finding functions is reserved for
the error status, you must provide storage for the location of the found
root.
- Function Argument: double * root
A pointer to a place for the GSL root finder to store the location
of the found root. This must be a valid pointer; GSL root finders
will not allocate any memory for you.
If a GSL root finder succeeds, it will return `0' and store the
location of the found root in `*root'.
If a GSL root finder fails, it will return `-1' and set `gsl_errno'
to a diagnostic value. *Note Root Finder Error Handling::, for a
discussion of possible error codes. Nothing useful will be stored in
`*root' if the function failed.
File: gsl-ref.info, Node: Providing the Function to Search, Next: Automatic Control Decisions, Prev: Root Finder Exit Values, Up: High Level Root Finding Functions
Providing the Function to Search
--------------------------------
You must provide a continous function of one variable for GSL root
finder(s) to operate on, and, sometimes, its first derivative.
Recall that when passing pointers to functions, you give the name of
the function you are passing. For example:
foo = i_take_a_function_pointer(my_function);
- Function Argument: double (* f)(double)
A pointer to the function whose root you are searching for. It is
called by the root finding function many times during its search.
It must be continous within the region of interest.
Here is an example function which you could pass to a GSL root
finder:
double
my_f (double x) {
return sin (2 * x) + 2 * cos (x);
}
- Function Argument: double (* df)(double)
A pointer to the first derivative of the function whose root you
are searching for.
If we were looking for a root of the function in the previous
File: gsl-ref.info, Node: Automatic Control Decisions, Next: High Level Functions, Prev: Providing the Function to Search, Up: High Level Root Finding Functions
Automatic Control Decisions
---------------------------
*FIXME: Coming soon...*
File: gsl-ref.info, Node: High Level Functions, Prev: Automatic Control Decisions, Up: High Level Root Finding Functions
(.txt)