ls_fft.h

00001 /*
00002  *  This file is part of libfftpack.
00003  *
00004  *  libfftpack is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  libfftpack is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with libfftpack; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  */
00018 
00019 /*
00020  *  libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
00021  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00022  *  (DLR).
00023  */
00024 
00025 /*! \file ls_fft.h
00026  *  Interface for the LevelS FFT package.
00027  *
00028  *  Copyright (C) 2004 Max-Planck-Society
00029  *  \author Martin Reinecke
00030  */
00031 
00032 #ifndef PLANCK_LS_FFT_H
00033 #define PLANCK_LS_FFT_H
00034 
00035 #ifdef __cplusplus
00036 extern "C" {
00037 #endif
00038 
00039 /*!\defgroup fftgroup FFT interface
00040 This package is intended to calculate one-dimensional real or complex FFTs
00041 with high accuracy and good efficiency even for lengths containing large
00042 prime factors.
00043 The code is written in C, but a Fortran wrapper exists as well.
00044 
00045 Before any FFT is executed, a plan must be generated for it. Plan creation
00046 is designed to be fast, so that there is no significant overhead if the
00047 plan is only used once or a few times.
00048 
00049 The main component of the code is a C port of Swarztrauber's FFTPACK
00050 (http://www.netlib.org/fftpack/), which was originally done by Pekka Janhunen
00051 and reformatted by Joerg Arndt.
00052 
00053 I added a few digits to the floating-point constants to achieve higher
00054 precision, split the complex transforms into separate routines for forward
00055 and backward transforms (to increase performance a bit), and replaced
00056 the iterative twiddling factor calculation in radfg() and radbg() by an exact
00057 calculation of the factors.
00058 
00059 Since FFTPACK becomes quite slow for FFT lengths with large prime factors
00060 (in the worst case of prime lengths it reaches O(n*n) complexity), I
00061 implemented Bluestein's algorithm, which computes a FFT of length n by
00062 several FFTs of length n2>=2*n and a convolution. Since n2 can be chosen
00063 to be highly composite, this algorithm is more efficient if n has large
00064 prime factors. The longer FFTs themselves are then computed using the FFTPACK
00065 routines.
00066 Bluestein's algorithm was implemented according to the description at
00067 <a href="http://en.wikipedia.org/wiki/Bluestein%27s_FFT_algorithm">Wikipedia</a>.
00068 
00069 \b Thread-safety:
00070 All routines can be called concurrently; all information needed by ls_fft
00071 is stored in the plan variable. However, using the same plan variable on
00072 multiple threads simultaneously is not supported and will lead to data
00073 corruption.
00074 */
00075 /*! \{ */
00076 
00077 typedef struct
00078   {
00079   double *work;
00080   int length;
00081   int bluestein;
00082   } complex_plan_i;
00083 
00084 /*! The opaque handle type for complex-FFT plans. */
00085 typedef complex_plan_i * complex_plan;
00086 
00087 /*! Returns a plan for a complex FFT with \a length elements. */
00088 complex_plan make_complex_plan (int length);
00089 /*! Destroys a plan for a complex FFT. */
00090 void kill_complex_plan (complex_plan plan);
00091 /*! Computes a complex forward FFT on \a data, using \a plan.
00092     \a Data has the form r0, i0, r1, i1, ..., r[length-1], i[length-1]. */
00093 void complex_plan_forward (complex_plan plan, double *data);
00094 /*! Computes a complex backward FFT on \a data, using \a plan.
00095     \a Data has the form r0, i0, r1, i1, ..., r[length-1], i[length-1]. */
00096 void complex_plan_backward (complex_plan plan, double *data);
00097 
00098 typedef struct
00099   {
00100   double *work;
00101   int length;
00102   int bluestein;
00103   } real_plan_i;
00104 
00105 /*! The opaque handle type for real-FFT plans. */
00106 typedef real_plan_i * real_plan;
00107 
00108 /*! Returns a plan for a real FFT with \a length elements. */
00109 real_plan make_real_plan (int length);
00110 /*! Destroys a plan for a real FFT. */
00111 void kill_real_plan (real_plan plan);
00112 /*! Computes a real forward FFT on \a data, using \a plan
00113     and assuming the FFTPACK storage scheme:
00114     - on entry, \a data has the form r0, r1, ..., r[length-1];
00115     - on exit, it has the form r0, r1, i1, r2, i2, ...
00116       (a total of \a length values). */
00117 void real_plan_forward_fftpack (real_plan plan, double *data);
00118 /*! Computes a real forward FFT on \a data, using \a plan
00119     and assuming the FFTPACK storage scheme:
00120     - on entry, \a data has the form r0, r1, i1, r2, i2, ...
00121     (a total of \a length values);
00122     - on exit, it has the form r0, r1, ..., r[length-1]. */
00123 void real_plan_backward_fftpack (real_plan plan, double *data);
00124 /*! Computes a real forward FFT on \a data, using \a plan
00125     and assuming the FFTW halfcomplex storage scheme:
00126     - on entry, \a data has the form r0, r1, ..., r[length-1];
00127     - on exit, it has the form r0, r1, r2, ..., i2, i1. */
00128 void real_plan_forward_fftw (real_plan plan, double *data);
00129 /*! Computes a real backward FFT on \a data, using \a plan
00130     and assuming the FFTW halfcomplex storage scheme:
00131     - on entry, \a data has the form r0, r1, r2, ..., i2, i1.
00132     - on exit, it has the form r0, r1, ..., r[length-1]. */
00133 void real_plan_backward_fftw (real_plan plan, double *data);
00134 /*! Computes a real forward FFT on \a data, using \a plan
00135     and assuming a full-complex storage scheme:
00136     - on entry, \a data has the form r0, [ignored], r1, [ignored], ...,
00137       r[length-1], [ignored];
00138     - on exit, it has the form r0, i0, r1, i1, ..., r[length-1], i[length-1].
00139     */
00140 void real_plan_forward_c (real_plan plan, double *data);
00141 /*! Computes a real backward FFT on \a data, using \a plan
00142     and assuming a full-complex storage scheme:
00143     - on entry, \a data has the form r0, i0, r1, i1, ...,
00144       r[length-1], i[length-1];
00145     - on exit, it has the form r0, 0, r1, 0, ..., r[length-1], 0. */
00146 void real_plan_backward_c (real_plan plan, double *data);
00147 
00148 /*! \} */
00149 
00150 #ifdef __cplusplus
00151 }
00152 #endif
00153 
00154 #endif