casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
FFTServer.h
Go to the documentation of this file.
00001 //# FFTServer.h: A class with methods for Fast Fourier Transforms
00002 //# Copyright (C) 1994,1995,1996,1997,1999,2003
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id: FFTServer.h 21051 2011-04-20 11:46:29Z gervandiepen $
00027 
00028 #ifndef SCIMATH_FFTSERVER_H
00029 #define SCIMATH_FFTSERVER_H
00030 
00031 //# Includes
00032 #include <casa/aips.h>
00033 #include <scimath/Mathematics/FFTW.h>
00034 #include <casa/Arrays/IPosition.h>
00035 #include <casa/Containers/Block.h>
00036 #include <vector>
00037 
00038 namespace casa { //# NAMESPACE CASA - BEGIN
00039 
00040 //# Forward Declarations
00041 template <class T> class Array;
00042 template <class S> class Matrix;
00043 // <summary>Lists the different types of FFT's that can be done</summary>
00044 // <synopsis>This enumerator is brought out as a separate class because g++
00045 // currently cannot handle enumerators in a templated class. When it can this
00046 // class will go away and this enumerator moved into the FFTServer
00047 // class</synopsis>
00048 class FFTEnums {
00049 public:
00050   enum TransformType {
00051     // Forward Complex to Complex transforms.
00052     COMPLEX,
00053     // Inverse Complex to Complex transforms.
00054     INVCOMPLEX,
00055     // Real to Complex or Complex to Real transforms.
00056     REALTOCOMPLEX,
00057     // Real to Complex or Complex to Real transforms.
00058     COMPLEXTOREAL,
00059     // Real to Real transforms with symmetric Arrays (not used)
00060     REALSYMMETRIC
00061   };
00062 };
00063 
00064 // <summary>A class with methods for Fast Fourier Transforms</summary>
00065 
00066 // <use visibility=export>
00067 
00068 // <reviewed reviewer="wbrouw" date="1997/10/29" tests="tFFTServer">
00069 // </reviewed>
00070 
00071 // <prerequisite> 
00072 // <li> Basic concepts of Fast Fourier Transforms.
00073 // <li> <linkto module=Arrays>The Arrays module</linkto>
00074 // </prerequisite>
00075 
00076 // <etymology> The FFTServer class, can do Fast Fourier Transforms of
00077 // any length and dimensionality.
00078 // </etymology>
00079 
00080 
00081 // <synopsis>
00082 
00083 // The FFTServer class provides methods for performing n-dimensional Fast
00084 // Fourier Transforms with real and complex Array's of arbitrary size and
00085 // dimensionality. It can do either real to complex, complex to real, or
00086 // complex to complex transforms with the "origin" of the transform either at
00087 // the centre of the Array or at the first element.
00088 
00089 // Because the output from a real to complex transform is Hermitian only half
00090 // of the complex result is returned. Similarly with a complex to real
00091 // transform only half of the complex plane is required, the other half is
00092 // implicitly assumed to be the complex conjugate of the supplied half-plane.
00093 // <note role=warning> The complex to real transform does not check that the
00094 // imaginary component of the values where u=0 are zero</note>
00095 
00096 // This class can be initialised with a shape that indicates the length of the
00097 // transforms that will be performed, and whether they are going to be
00098 // real<->complex transforms or complex<->complex ones. This initialisation
00099 // sets up a variety of internal buffers and computes factorizations and
00100 // twiddle factors used during the transform. The initialised transform shape
00101 // is always compared with the shape of the supplied arguments when a transform
00102 // is done and the FFTServer class will automatically resize itself if
00103 // necessary. So the default constructor is perfectly safe to use.
00104 
00105 // With any transform the output Array must either be the correct shape for the
00106 // desired output or zero length (ie not contain any elements). If it is zero
00107 // length then it will be resized to the correct shape. For a complex->complex
00108 // transform the output Array will be the same shape as the input Array. For a
00109 // real->complex transform the output Array will be the same size as the input
00110 // Array except in the first dimension which will have a length of (nx+2)/2. So
00111 // if nx=7 the output length will be 4 and if nx=8 the output length will be 5,
00112 // on the first axis. nx is the length of the first axis on the <em>real</em>
00113 // input Array and cx (which is used later) is the length of the first axis on
00114 // the <em>complex</em> input Array.
00115 
00116 // <strong>For complex to real transforms the output length on the first axis
00117 // is not uniquely defined by the shape of the complex input
00118 // Array</strong>. This class uses the following algorithm to work out the
00119 // length of the first axis on the output Array.
00120 // <ul> 
00121 // <li> If the size of the output Array is non-zero then its shape must match
00122 // the size of the input Array except for the first axis. The length of the
00123 // first axis must either be 2*cx-2 or 2*cx-1 and this determines the length of
00124 // the transform on the first axis.
00125 // <li> If the size of the output Array is zero then scan the imaginary
00126 // components of the values at the end of the first axis on the input Array (ie
00127 // at <src>[cx-1,....]</src> If any of these are non-zero the output Array
00128 // will have an odd length.
00129 // <li> Otherwise if all the imaginary components described above are zero then
00130 // look at the current size of the FFTServer object (either defined at
00131 // construction time or with the resize function). If it matches the size of
00132 // the input Array except for the first axis and if the length on this axis is
00133 // either 2*cx-2 or 2*cx-1 then use that to determine the size of the output
00134 // Array.
00135 // <li> Otherwise assume the output Array will an even length of 2*cx-2 on its
00136 // first axis.
00137 // </ul>
00138 
00139 // This class does transforms using the widely used FORTRAN fftpack
00140 // package or the highly optimized FFTW package (to be chosen at build time).
00141 // <br>
00142 // <em> P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
00143 // (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83. </em><br>
00144 // The fftpack package only does one dimensional transforms and this class
00145 // decomposes multi-dimensional transforms into a series of 1-dimensional ones.
00146 // <br>If at build time it is chosen to use FFTW in a multi-threaded way,
00147 // it will try to use as many cores as possible.
00148 
00149 // In this class a forward transform is defined as one that goes from the real
00150 // to the complex (or the time to frequency) domain. In a forward transform the
00151 // sign of the exponent is negative and no scaling is done on the output.  The
00152 // backward transform goes from the complex to the real (or the frequency to
00153 // the time) domain. The sign of the exponent is positive and the result is
00154 // always scaled by 1/N were N is the total number of elements in the Array.
00155 
00156 // The origin of the transform is defined as the point where setting only that
00157 // element to one, and then doing a forward transform results in an Array that
00158 // is all one. The <src>fft</src> member functions in this class all assume
00159 // that the origin of the Transform is at the centre of the Array ie. at
00160 // <src>[nx/2,ny/2,...]</src> were the indexing begins at zero. Because the
00161 // fftpack software assumes the origin of the transform is at the first element
00162 // ie.,<src>[0,0,...]</src> this class flips the data in the Array around to
00163 // compensate. For fftpack this flipping takes about one 20% of the total
00164 // transform time, while for FFTW it can easily exceed the transform time.
00165 // Flipping can be avoided by using the <src>fft0</src> member
00166 // functions which do not flip the data.
00167 
00168 // Some of the member functions in this class scramble the input Array,
00169 // possibly by flipping the quandrants of the data although this is not
00170 // guaranteed. Modification of the input Array can be avoided, at the expense
00171 // of copying the data to temporary storage, by either:
00172 // <ul> <li> Ensuring the input Array is a const Array.
00173 //      <li> Setting the constInput Flag to True.
00174 // </ul>
00175 // The latter option is provided to avoid users having to cast non-const
00176 // Arrays to const ones in order to prevent there input Array from being
00177 // scrambled.
00178 
00179 // <note role=warning> This class assumes that a Complex array is stored as
00180 // pairs of floating point numbers, with no intervening gaps, and with the real
00181 // component first ie., <src>[re0,im0,re1,im1, ...]</src>. This means that the
00182 // following type casts work,
00183 // <srcblock>
00184 // S * complexPtr;
00185 // T * realPtr = (T * ) complexPtr;
00186 // </srcblock>
00187 // and allow a Complex number to be accessed as a pair of real numbers. If this
00188 // assumption is bad then real Arrays will have to generated by copying the
00189 // complex ones. Ultimately this assumption about Complex<->Real Array
00190 // conversion should be put somewhere central like Array2Math.cc.
00191 // </note>
00192 // </synopsis>
00193 
00194 // <templating arg=T>
00195 // <li> The T argument must be of type Float or Double. These are the only 
00196 // possible instantiations of this class.
00197 // </templating>
00198 
00199 // <templating arg=S>
00200 // <li> The S argument must be of type Complex, if T is Float, or DComplex, if T is
00201 // Double. These are the only possible instantiations of this class.
00202 // </templating>
00203 //
00204 // <example>
00205 // Do a real to complex Transform of a 1-Dimensional Vector.  The following
00206 // example can trivially be extended to any number of dimensions.
00207 // <srcblock>
00208 // FFTServer<Float,Complex> server;
00209 // Vector<Float> input(32);
00210 // Vector<Complex> output(17);
00211 // input = 0.0f;
00212 // input(16) = 1.0f;
00213 // cout << "Input:" << input << endl;
00214 // server.fft(output, input);
00215 // cout << "Output:" << output << endl;
00216 // </srcblock>
00217 // </example>
00218 //
00219 // <thrown>
00220 // <li> AipsError: If the input and output Array have bad or incompatible
00221 // shapes. See the individual function descriptions for what Array shapes are
00222 // required.
00223 // </thrown>
00224 //
00225 // <todo asof="1997/10/22">
00226 //   <li> The time taken to flip the Array can be reduced, if all the Array
00227 //   dimensions are even, by pre-multiplying the every other element on the
00228 //   input Array by -1. Then no flipping needs to be done on the output Array.
00229 // </todo>
00230 
00231 template<class T, class S> class FFTServer
00232 {
00233 public:
00234 
00235   // The default constructor. The server will automatically resize to do
00236   // transforms of the appropriate length when necessary.
00237   FFTServer();
00238 
00239   // Initialise the server to do transforms on Arrays of the specified
00240   // shape. The server will, however, resize to do transforms of other lengths
00241   // if necessary. See the resize function for a description of the
00242   // TransformType enumerator.
00243   FFTServer(const IPosition & fftSize, 
00244             const FFTEnums::TransformType transformType 
00245             = FFTEnums::REALTOCOMPLEX);
00246   
00247   // copy constructor. The copied server is initialised to do transforms of the
00248   // same length as the other server. Uses copy (and not reference) semantics
00249   // so that changing the transform length of one server does not affect the
00250   // other server.
00251   FFTServer(const FFTServer<T,S> & other);
00252 
00253   // destructor
00254   ~FFTServer();
00255   
00256   // The assignment operator which does the same thing as the copy
00257   // constructor.
00258   FFTServer<T,S> & operator=(const FFTServer<T,S> & other);
00259 
00260   // Modify the FFTServer object to do transforms of the supplied shape. The
00261   // amount of internal storage, and the initialisation, depends on the type of
00262   // transform that will be done. The transform type is specified with the
00263   // TransformTypes enumerator. Currently there is no difference in
00264   // initialisation for the COMPLEXTOREAL and REALTOCOMPLEX transforms. The
00265   // shape argument is the shape of the real array (or complex one if complex
00266   // to complex transforms are being done). In general it is not necessary to
00267   // use this function as all the fft & fft0 functions will automatically
00268   // resize the server, if necessary, to match their input arguments.
00269   void resize(const IPosition & fftSize,
00270               const FFTEnums::TransformType transformType
00271               = FFTEnums::REALTOCOMPLEX);
00272 
00273   // Real to complex fft. The origin of the transform is in the centre of the
00274   // Array. Because of the Hermitian property the output Array only contains
00275   // half of the complex result. The output Array must either have no elements
00276   // or be a size that is appropriate to the input Array size,
00277   // ie. <src>shape = [(nx+2)/2, ny, nz,...]</src>.  Otherwise an AipsError is
00278   // thrown. See the synopsis for a description of the constInput flag.
00279   // <group>
00280   void fft(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
00281   void fft(Array<S> & cResult, const Array<T> & rData);
00282   // </group>
00283 
00284   // Complex to real fft. The origin of the transform is in the centre of the
00285   // Array. Because of the Hermitian property the input Array only contains
00286   // half of the complex values. The output Array must either have no elements,
00287   // or be a size that is appropriate to the input Array size ie.,<br>
00288   // <src>shape = [2*cx-2, cy, cz,...]</src> or <br>
00289   // <src>shape = [2*cx-1, cy, cz,...]</src>.  <br>
00290   // Otherwise an AipsError is thrown. See the description in the synopsis for
00291   // the algorithm used to choose between the two possible output shapes and a
00292   // description of the constInput Flag.  
00293   // <group>
00294   void fft(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
00295   void fft(Array<T> & rResult, const Array<S> & cData);
00296   // </group>
00297 
00298   // Complex to complex in-place fft. The origin of the transform is in the
00299   // centre of the Array. The direction of the transform is controlled by the
00300   // toFrequency variable. If True then a forward, or time to frequency,
00301   // transform is performed. If False a backward or frequency to time transform
00302   // is done. Scaling is always done on the backward transform.
00303   void fft(Array<S> & cValues, const Bool toFrequency=True);
00304 
00305   // Complex to complex fft. The origin of the transform is in the centre of
00306   // the Array. The direction of the transform is controlled by the toFrequency
00307   // variable. If True then a forward, or time to frequency, transform is
00308   // performed. If False a backward or frequency to time transform is
00309   // done. Scaling is always done on the backward transform. The output Array
00310   // must either either contain no elements or be the same as the input Array,
00311   // ie. <src>shape = [cx, cy, cz,...]</src>.  Otherwise an AipsError is
00312   // thrown.
00313   void fft(Array<S> & cResult, const Array<S> & cData,
00314            const Bool toFrequency=True);
00315 
00316   // The <src>fft0</src> functions are equivalent to the <src>fft</src>
00317   // functions described above except that the origin of the transform is the
00318   // first element of the Array, ie. [0,0,0...], rather than the centre
00319   // element, ie [nx/2, ny/2, nz/2, ...]. As the underlying functions
00320   // assume that the origin of the transform is the first element these
00321   // routines are in general faster than the equivalent ones with the origin
00322   // at the centre of the Array.
00323   // <group>
00324   void fft0(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
00325   void fft0(Array<S> & cResult, const Array<T> & rData);
00326   void fft0(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
00327   void fft0(Array<T> & rResult, const Array<S> & cData);
00328   void fft0(Array<S> & cValues, const Bool toFrequency=True);
00329   void fft0(Array<S> & cResult, const Array<S> & cData,
00330             const Bool toFrequency=True);
00331   //# void fft0(Array<T> & rValues, const Bool toFrequency=True);
00332 
00333   // </group>
00334   //# Flips the quadrants in a complex Array so that the point at
00335   //# cData.shape()/2 moves to the origin. This moves, for example, the point
00336   //# at [8,3] to the origin ([0,0]) in an array of shape [16,7]. Usually two
00337   //# flips will restore an Array to its original state.  But for Array's
00338   //# where one or more dimension is an odd length two flips do NOT restore
00339   //# the data to its original state.  So the when toZero=False this routine
00340   //# does an unflip operation (ie moves the data at [0,0] to the centre) and
00341   //# restores the data to its original state for odd length arrays.  When
00342   //# passed a Hermitian Array where half the complex plane is implicit (eg as
00343   //# produced by a real->complex Transform) it is not necessary to flip the
00344   //# first dimension of the Array. In this case the isHermitian flag should
00345   //# be set to True.  For complex<->complex transforms this should be False.
00346   // <group>
00347   void flip(Array<T> & rData, const Bool toZero, const Bool isHermitian);
00348   void flip(Array<S> & cData, const Bool toZero, const Bool isHermitian);
00349   // </group>
00350 
00351   // N-D in-place complex->complex FFT shift (FFT - phase-mult - inverse FFT)
00352   // If toFrequency is true, the first FFT will be from time to frequency. 
00353   // relshift is the freq shift normalised to the bandwidth.
00354   // Only transform over selected dimension. Iterate over the others. 
00355   void fftshift(Array<S> & cValues, const uInt& whichAxis, 
00356                 const Double& relshift, const Bool toFrequency=True);
00357 
00358   // N-D complex->complex FFT shift (FFT - phase-mult - inverse FFT)
00359   // with flagging.
00360   // If toFrequency is true, the first FFT will be from time to frequency. 
00361   // relshift is the freq shift normalised to the bandwidth.
00362   // Only transform over selected dimension. Iterate over the others. 
00363   void fftshift(Array<S> & outValues, Array<Bool> & outFlags,
00364                 const Array<S> & cValues, const Array<Bool>& inFlags,
00365                 const uInt& whichAxis, 
00366                 const Double& relshift, 
00367                 const Bool goodIsTrue=False,
00368                 const Bool toFrequency=True);
00369 
00370   // N-D real->real FFT shift (FFT to complex - phase-mult - inverse FFT)
00371   // with flagging.
00372   // relshift is the freq shift normalised to the bandwidth.
00373   // Only transform over selected dimension. Iterate over the others. 
00374   void fftshift(Array<T> & outValues, Array<Bool> & outFlags,
00375                 const Array<T> & rValues, const Array<Bool>& inFlags,
00376                 const uInt& whichAxis, 
00377                 const Double& relshift, 
00378                 const Bool goodIsTrue=False);
00379 
00380 private:
00381   //# finds the shape of the output array when doing complex->real transforms
00382   IPosition determineShape(const IPosition & rShape, const Array<S> & cData);
00383 
00384   //# Data members.
00385   // The size of the last FFT done by this object
00386   IPosition itsSize;
00387   // Whether the last FFT was complex<->complex or not
00388   FFTEnums::TransformType itsTransformType;
00389   //# twiddle factors and factorisations used by fftpack
00390   PtrBlock<Block<T> *> itsWork;
00391   // buffer for copying non-contigious arrays to contigious ones. This is done
00392   // so that the FFT's have a better chance of fitting into cache and hence
00393   // going faster. 
00394   // This buffer is also used as temporary storage when flipping the data.
00395   Block<S> itsBuffer;
00396   // FFTW specific members. Do not harm if FFTPack is used.
00397   FFTW           itsFFTW;
00398   std::vector<T> itsWorkIn;
00399   std::vector<S> itsWorkOut;
00400   std::vector<S> itsWorkC2C;
00401 };
00402 
00403 
00404 } //# NAMESPACE CASA - END
00405 
00406 //# Do NOT include the .tcc file here like done for other templated classes.
00407 //# The instantiations are done explicitly.
00408 //# In this way the HAVE_FFTW ifdef is only used in .cc files and does
00409 //# not appear in headers, so other packages using FFTServer do not need
00410 //# to (un)set HAVE_FFTW.
00411 
00412 #endif