This is gsl-ref.info, produced by makeinfo version 4.0 from 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, 1999, 2000, 2001 The GSL Team. Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.1 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled "GNU Free Documentation License".  File: gsl-ref.info, Node: Complex Hermitian Matrices, Next: Sorting Eigenvalues and Eigenvectors, Prev: Real Symmetric Matrices, Up: Eigensystems Complex Hermitian Matrices ========================== - Function: gsl_eigen_herm_workspace * gsl_eigen_herm_alloc (const size_t N) This function allocates a workspace for computing eigenvalues of N-by-N complex hermitian matrices. The size of the workspace is O(3n). - Function: void gsl_eigen_herm_free (gsl_eigen_herm_workspace * W) This function frees the memory associated with the workspace W. - Function: int gsl_eigen_herm (gsl_matrix_complex * A, gsl_vector * EVAL, gsl_eigen_herm_workspace * W) This function computes the eigenvalues of the complex hermitian matrix A. Additional workspace of the appropriate size must be provided in W. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The imaginary parts of the diagonal are assumed to be zero and are not referenced. The eigenvalues are stored in the vector EVAL and are unordered. - Function: gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t N) This function allocates a workspace for computing eigenvalues and eigenvectors of N-by-N complex hermitian matrices. The size of the workspace is O(5n). - Function: void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * W) This function frees the memory associated with the workspace W. - Function: int gsl_eigen_hermv (gsl_matrix_complex * A, gsl_vector * EVAL, gsl_matrix_complex * EVEC, gsl_eigen_hermv_workspace * W) This function computes the eigenvalues and eigenvectors of the complex hermitian matrix A. Additional workspace of the appropriate size must be provided in W. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The imaginary parts of the diagonal are assumed to be zero and are not referenced. The eigenvalues are stored in the vector EVAL and are unordered. The corresponding complex eigenvectors are stored in the columns of the matrix EVEC. For example, the eigenvector in the first column corresponds to the first eigenvalue. The eigenvectors are guaranteed to be mutually orthogonal and normalised to unit magnitude.  File: gsl-ref.info, Node: Sorting Eigenvalues and Eigenvectors, Next: Eigenvalue and Eigenvector Examples, Prev: Complex Hermitian Matrices, Up: Eigensystems Sorting Eigenvalues and Eigenvectors ==================================== - Function: int gsl_eigen_symmv_sort (gsl_vector * EVAL, gsl_matrix * EVEC, gsl_eigen_sort_t SORT_TYPE) This function simultaneously sorts the eigenvalues stored in the vector EVAL and the corresponding real eigenvectors stored in the columns of the matrix EVEC into ascending or descending order according to the value of the parameter SORT_TYPE, `GSL_EIGEN_SORT_VAL_ASC' ascending order in numerical value `GSL_EIGEN_SORT_VAL_DESC' descending order in numerical value `GSL_EIGEN_SORT_ABS_ASC' ascending order in magnitude `GSL_EIGEN_SORT_ABS_DESC' descending order in magnitude - Function: int gsl_eigen_hermv_sort (gsl_vector * EVAL, gsl_matrix_complex * EVEC, gsl_eigen_sort_t SORT_TYPE) This function simultaneously sorts the eigenvalues stored in the vector EVAL and the corresponding complex eigenvectors stored in the columns of the matrix EVEC into ascending or descending order according to the value of the parameter SORT_TYPE as shown above.  File: gsl-ref.info, Node: Eigenvalue and Eigenvector Examples, Next: Eigenvalue and Eigenvector References, Prev: Sorting Eigenvalues and Eigenvectors, Up: Eigensystems Examples ======== The following program computes the eigenvalues and eigenvectors of the 4-th order Hilbert matrix, H(i,j) = 1/(i + j + 1). #include #include #include int main (void) { double data[] = { 1.0 , 1/2.0, 1/3.0, 1/4.0, 1/2.0, 1/3.0, 1/4.0, 1/5.0, 1/3.0, 1/4.0, 1/5.0, 1/6.0, 1/4.0, 1/5.0, 1/6.0, 1/7.0 }; gsl_matrix_view m = gsl_matrix_view_array(data, 4, 4); gsl_vector *eval = gsl_vector_alloc (4); gsl_matrix *evec = gsl_matrix_alloc (4, 4); gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (4); gsl_eigen_symmv (&m.matrix, eval, evec, w); gsl_eigen_symmv_free(w); gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_ASC); { int i; for (i = 0; i < 4; i++) { double eval_i = gsl_vector_get(eval, i); gsl_vector_view evec_i = gsl_matrix_column(evec, i); printf("eigenvalue = %g\n", eval_i); printf("eigenvector = \n"); gsl_vector_fprintf(stdout, &evec_i.vector, "%g"); } } return 0; } Here is the beginning of the output from the program, $ ./a.out eigenvalue = 9.67023e-05 eigenvector = -0.0291933 0.328712 -0.791411 0.514553 ... This can be compared with the corresponding output from GNU OCTAVE, octave> [v,d] = eig(hilb(4)); octave> diag(d) ans = 9.6702e-05 6.7383e-03 1.6914e-01 1.5002e+00 octave> v v = 0.029193 0.179186 -0.582076 0.792608 -0.328712 -0.741918 0.370502 0.451923 0.791411 0.100228 0.509579 0.322416 -0.514553 0.638283 0.514048 0.252161 Note that the eigenvectors can differ by a change of sign, since the sign of an eigenvector is arbitrary.  File: gsl-ref.info, Node: Eigenvalue and Eigenvector References, Prev: Eigenvalue and Eigenvector Examples, Up: Eigensystems References and Further Reading ============================== Further information on the algorithms described in this section can be found in the following book, G. H. Golub, C. F. Van Loan, `Matrix Computations' (3rd Ed, 1996), Johns Hopkins University Press, ISBN 0-8018-5414-8. The LAPACK library is described in, `LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM, ISBN 0-89871-447-8. The LAPACK source code can be found at the website above along with an online copy of the users guide.  File: gsl-ref.info, Node: Fast Fourier Transforms, Next: Numerical Integration, Prev: Eigensystems, Up: Top Fast Fourier Transforms (FFTs) ****************************** 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 FFTPACK library by Paul Swarztrauber. Fortran code for FFTPACK is available on Netlib (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 `GSL FFT Algorithms' (*note 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::  File: gsl-ref.info, Node: Mathematical Definitions, Next: Overview of complex data FFTs, Up: Fast Fourier Transforms Mathematical Definitions ======================== Fast Fourier Transforms are efficient algorithms for calculating the discrete fourier transform (DFT), x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N) 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 W\vec{z}. A general matrix-vector multiplication takes O(N^2) operations for N data-points. Fast fourier transform algorithms use a divide-and-conquer strategy to factorize the matrix W into smaller sub-matrices, corresponding to the integer factors of the length N. If N can be factorized into a product of integers f_1 f_2 ... f_n then the DFT can be computed in O(N \sum f_i) operations. For a radix-2 FFT this gives an operation count of 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 "forward fourier transform", x = FFT(z), is, x_j = \sum_{k=0}^{N-1} z_k \exp(-2\pi i j k / N) and the definition of the "inverse fourier transform", x = IFFT(z), is, z_j = {1 \over N} \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N). The factor of 1/N makes this a true inverse. For example, a call to `gsl_fft_complex_forward' followed by a call to `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 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 "backwards FFT" is simply our terminology for an unscaled version of the inverse FFT, z^{backwards}_j = \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N). 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.  File: gsl-ref.info, Node: Overview of complex data FFTs, Next: Radix-2 FFT routines for complex data, Prev: Mathematical Definitions, Up: Fast Fourier Transforms Overview of complex data FFTs ============================= The inputs and outputs for the complex FFT routines are "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, gsl_complex_packed_array data[6]; can be used to hold an array of three complex numbers, `z[3]', in the following way, 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]) A "stride" parameter allows the user to perform transforms on the elements `z[stride*i]' instead of `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 \Delta then the frequency-domain includes both positive and negative frequencies, ranging from -1/(2\Delta) through 0 to +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 DATA, and the correspondence between the time-domain data z, and the frequency-domain data x. 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)) When N is even the location N/2 contains the most positive and negative frequencies +1/(2 \Delta), -1/(2 \Delta)) which are equivalent. If N is odd then general structure of the table above still applies, but N/2 does not appear.  File: gsl-ref.info, Node: Radix-2 FFT routines for complex data, Next: Mixed-radix FFT routines for complex data, Prev: Overview of complex data FFTs, Up: Fast Fourier Transforms Radix-2 FFT routines for 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 `gsl_fft_complex.h'. - Function: int gsl_fft_complex_radix2_forward (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_transform (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_backward (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) These functions compute forward, backward and inverse FFTs of length N with stride STRIDE, on the packed complex array DATA using an in-place radix-2 decimation-in-time algorithm. The length of the transform N is restricted to powers of two. The functions return a value of `GSL_SUCCESS' if no errors were detected, or `GSL_EDOM' if the length of the data N is not a power of two. - Function: int gsl_fft_complex_radix2_dif_forward (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_dif_transform (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_dif_backward (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) - Function: int gsl_fft_complex_radix2_dif_inverse (gsl_complex_packed_array DATA[], size_t STRIDE, size_t N) These are decimation-in-frequency versions of the radix-2 FFT functions. 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 (-10 ... 10), where the negative times wrap around the end of the array. #include #include #include #include #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; } Note that we have assumed that the program is using the default 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 `gsl_fft_complex_radix2_forward'. The transformed data is rescaled by 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 t=128, and working in units of k/N, the DFT approximates the continuum fourier transform, giving a modulated \sin function.  File: gsl-ref.info, Node: Mixed-radix FFT routines for complex data, Next: Overview of real data FFTs, Prev: Radix-2 FFT routines for complex data, Up: Fast Fourier Transforms Mixed-radix FFT routines for 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 FFTPACK library by Paul Swarztrauber. The theory is explained in the review article `Self-sorting Mixed-radix FFTs' by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as 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 2*2 and 2*3. For factors which are not implemented as modules there is a fall-back to a general length-n module which uses Singleton's method for efficiently computing a DFT. This module is O(n^2), and slower than a dedicated module would be but works for any length n. Of course, lengths which use the general length-n module will still be factorized as much as possible. For example, a length of 143 will be factorized into 11*13. Large prime factors are the worst case scenario, e.g. as found in n=2*3*99991, and should be avoided because their O(n^2) scaling will dominate the run-time (consult the document `GSL FFT Algorithms' included in the GSL distribution if you encounter this problem). The mixed-radix initialization function `gsl_fft_complex_wavetable_alloc' returns the list of factors chosen by the library for a given length 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 N \sum f_i, where the f_i are the factors of N. For programs under user control you may wish to 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 `GSL FFT Algorithms' for details on adding support for other factors. All these functions are declared in the header file `gsl_fft_complex.h'. - Function: gsl_fft_complex_wavetable * gsl_fft_complex_wavetable_alloc (size_t N) This function prepares a trigonometric lookup table for a complex FFT of length N. The function returns a pointer to the newly allocated `gsl_fft_complex_wavetable' if no errors were detected, and a null pointer in the case of error. The length 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 `sin' and `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. - Function: void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * WAVETABLE) This function frees the memory associated with the wavetable WAVETABLE. The wavetable can be freed if no further FFTs of the same length will be needed. These functions operate on a `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 `gsl_fft_complex.h'. - 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: `size_t n' This is the number of complex data points `size_t nf' This is the number of factors that the length `n' was decomposed into. `size_t factor[64]' This is the array of factors. Only the first `nf' elements are used. `gsl_complex * trig' This is a pointer to a preallocated trigonometric lookup table of `n' complex elements. `gsl_complex * twiddle[64]' This is an array of pointers into `trig', giving the twiddle factors for each pass. The mixed radix algorithms require an additional working space to hold the intermediate steps of the transform. - Function: gsl_fft_complex_workspace * gsl_fft_complex_workspace_alloc (size_t N) This function allocates a workspace for a complex transform of length N. - Function: void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * WORKSPACE) This function frees the memory associated with the workspace WORKSPACE. The workspace can be freed if no further FFTs of the same length will be needed. The following functions compute the transform, - 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) - Function: 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) - Function: 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) - Function: 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) These functions compute forward, backward and inverse FFTs of length N with stride STRIDE, on the packed complex array DATA, using a mixed radix decimation-in-frequency algorithm. There is no restriction on the length N. Efficient modules are provided for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining factors are computed with a slow, O(n^2), general-n module. The caller must supply a WAVETABLE containing the trigonometric lookup tables and a workspace WORK. The functions return a value of `0' if no errors were detected. 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. 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 #include #include #include #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; } 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 all the `gsl' routines.  File: gsl-ref.info, Node: Overview of real data FFTs, Next: Radix-2 FFT routines for real data, Prev: Mixed-radix FFT routines for complex data, Up: Fast Fourier Transforms Overview of real data FFTs ========================== 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.  File: gsl-ref.info, Node: Radix-2 FFT routines for real data, Next: Mixed-radix FFT routines for real data, Prev: Overview of real data FFTs, Up: Fast Fourier Transforms 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_transform (double DATA[], size_t STRIDE, size_t N) This function computes an in-place radix-2 FFT of length N and stride STRIDE on the real array 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 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[], size_t STRIDE, size_t N) - Function: int gsl_fft_halfcomplex_radix2_backward (double DATA[], size_t STRIDE, size_t N) These functions compute the inverse or backwards in-place radix-2 FFT of length N and stride STRIDE on the half-complex sequence DATA 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: FFT References and Further Reading, Prev: Radix-2 FFT routines for real data, Up: Fast Fourier Transforms 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 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_transform', 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 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: gsl_fft_real_wavetable * gsl_fft_real_wavetable_alloc (size_t N) - Function: gsl_fft_halfcomplex_wavetable * gsl_fft_halfcomplex_wavetable_alloc (size_t N) These functions prepare trigonometric lookup tables for an FFT of size 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 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 `sin' and `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. - Function: void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * WAVETABLE) - Function: void gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * WAVETABLE) These functions free the memory associated with the wavetable WAVETABLE. The wavetable can be freed if no further FFTs of the same length will be needed. The mixed radix algorithms require an additional working space to hold the intermediate steps of the transform, - Function: gsl_fft_real_workspace * gsl_fft_real_workspace_alloc (size_t N) This function allocates a workspace for a real transform of length N. The same workspace is used for both forward real and inverse halfcomplex transforms. - Function: void gsl_fft_real_workspace_free (gsl_fft_real_workspace * WORKSPACE) This function frees the memory associated with the workspace WORKSPACE. The workspace can be freed if no further FFTs of the same length will be needed. The following functions compute the transforms of real and half-complex data, - 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) - Function: int gsl_fft_halfcomplex_transform (double DATA[], size_t STRIDE, size_t N, const gsl_fft_halfcomplex_wavetable * WAVETABLE, gsl_fft_real_workspace * WORK) 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_transform' DATA is an array of time-ordered real data. For `gsl_fft_halfcomplex_transform' DATA contains fourier coefficients in the half-complex ordering described above. 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 caller must supply a WAVETABLE containing trigonometric lookup tables and a workspace WORK. - Function: int gsl_fft_real_unpack (const double REAL_COEFFICIENT[], gsl_complex_packed_array COMPLEX_COEFFICIENT[], size_t STRIDE, size_t 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 algorithm for the conversion is simply, for (i = 0; i < n; i++) { complex_coefficient[i].real = real_coefficient[i]; complex_coefficient[i].imag = 0.0; } - Function: int gsl_fft_halfcomplex_unpack (const double HALFCOMPLEX_COEFFICIENT[], gsl_complex_packed_array COMPLEX_COEFFICIENT[], size_t STRIDE, size_t N) This function converts HALFCOMPLEX_COEFFICIENT, an array of half-complex coefficients as returned by `gsl_fft_real_transform', into an ordinary complex array, COMPLEX_COEFFICIENT. It fills in the complex array using the symmetry z_k = z_{N-k}^* to reconstruct the redundant elements. The algorithm for the conversion is, 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; } Here is an example program using `gsl_fft_real_transform' and `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 `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. #include #include #include #include #include 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; }  File: gsl-ref.info, Node: FFT References and Further Reading, Prev: Mixed-radix FFT routines for real data, Up: Fast Fourier Transforms References and Further Reading ============================== A good starting point for learning more about the FFT is the review article `Fast Fourier Transforms: A Tutorial Review and A State of the Art' by Duhamel and Vetterli, P. Duhamel and M. Vetterli. Fast fourier transforms: A tutorial review and a state of the art. `Signal Processing', 19:259-299, 1990. To find out about the algorithms used in the GSL routines you may want to consult the latex document `GSL FFT Algorithms' (it is included in GSL, as `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. There are several introductory books on the FFT with example programs, such as `The Fast Fourier Transform' by Brigham and `DFT/FFT and Convolution Algorithms' by Burrus and Parks, E. Oran Brigham. `The Fast Fourier Transform'. Prentice Hall, 1974. C. S. Burrus and T. W. Parks. `DFT/FFT and Convolution Algorithms'. Wiley, 1984. Both these introductory books cover the radix-2 FFT in some detail. The mixed-radix algorithm at the heart of the FFTPACK routines is reviewed in Clive Temperton's paper, Clive Temperton. Self-sorting mixed-radix fast fourier transforms. `Journal of Computational Physics', 52(1):1-23, 1983. The derivation of FFTs for real-valued data is explained in the following two articles, Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney Burrus. Real-valued fast fourier transform algorithms. `IEEE Transactions on Acoustics, Speech, and Signal Processing', ASSP-35(6):849-863, 1987. Clive Temperton. Fast mixed-radix real fourier transforms. `Journal of Computational Physics', 52:340-350, 1983. In 1979 the IEEE published a compendium of carefully-reviewed Fortran FFT programs in `Programs for Digital Signal Processing'. It is a useful reference for implementations of many different FFT algorithms, Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal Processing Committee, editors. `Programs for Digital Signal Processing'. IEEE Press, 1979. 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. FFTW Website,  File: gsl-ref.info, Node: Numerical Integration, Next: Random Number Generation, Prev: Fast Fourier Transforms, Up: Top Numerical Integration ********************* This chapter describes routines for performing numerical integration (quadrature) of a function in one dimension. There are routines for adaptive and non-adaptive integration of general functions, with specialised routines for specific cases. These include integration over infinite and semi-infinite ranges, singular integrals, including logarithmic singularities, computation of Cauchy principal values and oscillatory integrals. The library reimplements the algorithms used in QUADPACK, a numerical integration package written by Piessens, Doncker-Kapenga, Uberhuber and Kahaner. Fortran code for QUADPACK is available on Netlib. The functions described in this chapter are declared in the header file `gsl_integration.h'. * Menu: * Numerical Integration Introduction:: * QNG non-adaptive Gauss-Kronrod integration:: * QAG adaptive integration:: * QAGS adaptive integration with singularities:: * QAGP adaptive integration with known singular points:: * QAGI adaptive integration on infinite intervals:: * QAWC adaptive integration for Cauchy principal values:: * QAWS adaptive integration for singular functions:: * QAWO adaptive integration for oscillatory functions:: * QAWF adaptive integration for Fourier integrals:: * Numerical integration error codes:: * Numerical integration examples:: * Numerical integration References and Further Reading::