casacore
FFTServer.h
Go to the documentation of this file.
1 //# FFTServer.h: A class with methods for Fast Fourier Transforms
2 //# Copyright (C) 1994,1995,1996,1997,1999,2003
3 //# Associated Universities, Inc. Washington DC, USA.
4 //#
5 //# This library is free software; you can redistribute it and/or modify it
6 //# under the terms of the GNU Library General Public License as published by
7 //# the Free Software Foundation; either version 2 of the License, or (at your
8 //# option) any later version.
9 //#
10 //# This library is distributed in the hope that it will be useful, but WITHOUT
11 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 //# License for more details.
14 //#
15 //# You should have received a copy of the GNU Library General Public License
16 //# along with this library; if not, write to the Free Software Foundation,
17 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 //#
19 //# Correspondence concerning AIPS++ should be addressed as follows:
20 //# Internet email: aips2-request@nrao.edu.
21 //# Postal address: AIPS++ Project Office
22 //# National Radio Astronomy Observatory
23 //# 520 Edgemont Road
24 //# Charlottesville, VA 22903-2475 USA
25 //#
26 //# $Id$
27 
28 #ifndef SCIMATH_FFTSERVER_H
29 #define SCIMATH_FFTSERVER_H
30 
31 //# Includes
32 #include <casacore/casa/aips.h>
33 #include <casacore/scimath/Mathematics/FFTW.h>
34 #include <casacore/casa/Arrays/IPosition.h>
35 #include <casacore/casa/Containers/Block.h>
36 #include <vector>
37 
38 namespace casacore { //# NAMESPACE CASACORE - BEGIN
39 
40 //# Forward Declarations
41 template <class T> class Array;
42 template <class S> class Matrix;
43 // <summary>Lists the different types of FFT's that can be done</summary>
44 // <synopsis>This enumerator is brought out as a separate class because g++
45 // currently cannot handle enumerators in a templated class. When it can this
46 // class will go away and this enumerator moved into the FFTServer
47 // class</synopsis>
48 class FFTEnums {
49 public:
50  enum TransformType {
51  // Forward Complex to Complex transforms.
52  COMPLEX,
53  // Inverse Complex to Complex transforms.
54  INVCOMPLEX,
55  // Real to Complex or Complex to Real transforms.
57  // Real to Complex or Complex to Real transforms.
59  // Real to Real transforms with symmetric Arrays (not used)
61  };
62 };
63 
64 // <summary>A class with methods for Fast Fourier Transforms</summary>
65 
66 // <use visibility=export>
67 
68 // <reviewed reviewer="wbrouw" date="1997/10/29" tests="tFFTServer">
69 // </reviewed>
70 
71 // <prerequisite>
72 // <li> Basic concepts of Fast Fourier Transforms.
73 // <li> <linkto module=Arrays>The Arrays module</linkto>
74 // </prerequisite>
75 
76 // <etymology> The FFTServer class, can do Fast Fourier Transforms of
77 // any length and dimensionality.
78 // </etymology>
79 
80 
81 // <synopsis>
82 
83 // The FFTServer class provides methods for performing n-dimensional Fast
84 // Fourier Transforms with real and complex Array's of arbitrary size and
85 // dimensionality. It can do either real to complex, complex to real, or
86 // complex to complex transforms with the "origin" of the transform either at
87 // the centre of the Array or at the first element.
88 
89 // Because the output from a real to complex transform is Hermitian only half
90 // of the complex result is returned. Similarly with a complex to real
91 // transform only half of the complex plane is required, the other half is
92 // implicitly assumed to be the complex conjugate of the supplied half-plane.
93 // <note role=warning> The complex to real transform does not check that the
94 // imaginary component of the values where u=0 are zero</note>
95 
96 // This class can be initialised with a shape that indicates the length of the
97 // transforms that will be performed, and whether they are going to be
98 // real<->complex transforms or complex<->complex ones. This initialisation
99 // sets up a variety of internal buffers and computes factorizations and
100 // twiddle factors used during the transform. The initialised transform shape
101 // is always compared with the shape of the supplied arguments when a transform
102 // is done and the FFTServer class will automatically resize itself if
103 // necessary. So the default constructor is perfectly safe to use.
104 
105 // With any transform the output Array must either be the correct shape for the
106 // desired output or zero length (ie not contain any elements). If it is zero
107 // length then it will be resized to the correct shape. For a complex->complex
108 // transform the output Array will be the same shape as the input Array. For a
109 // real->complex transform the output Array will be the same size as the input
110 // Array except in the first dimension which will have a length of (nx+2)/2. So
111 // if nx=7 the output length will be 4 and if nx=8 the output length will be 5,
112 // on the first axis. nx is the length of the first axis on the <em>real</em>
113 // input Array and cx (which is used later) is the length of the first axis on
114 // the <em>complex</em> input Array.
115 
116 // <strong>For complex to real transforms the output length on the first axis
117 // is not uniquely defined by the shape of the complex input
118 // Array</strong>. This class uses the following algorithm to work out the
119 // length of the first axis on the output Array.
120 // <ul>
121 // <li> If the size of the output Array is non-zero then its shape must match
122 // the size of the input Array except for the first axis. The length of the
123 // first axis must either be 2*cx-2 or 2*cx-1 and this determines the length of
124 // the transform on the first axis.
125 // <li> If the size of the output Array is zero then scan the imaginary
126 // components of the values at the end of the first axis on the input Array (ie
127 // at <src>[cx-1,....]</src> If any of these are non-zero the output Array
128 // will have an odd length.
129 // <li> Otherwise if all the imaginary components described above are zero then
130 // look at the current size of the FFTServer object (either defined at
131 // construction time or with the resize function). If it matches the size of
132 // the input Array except for the first axis and if the length on this axis is
133 // either 2*cx-2 or 2*cx-1 then use that to determine the size of the output
134 // Array.
135 // <li> Otherwise assume the output Array will an even length of 2*cx-2 on its
136 // first axis.
137 // </ul>
138 
139 // This class does transforms using the widely used FORTRAN fftpack
140 // package or the highly optimized FFTW package (to be chosen at build time).
141 // <br>
142 // <em> P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
143 // (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83. </em><br>
144 // The fftpack package only does one dimensional transforms and this class
145 // decomposes multi-dimensional transforms into a series of 1-dimensional ones.
146 // <br>If at build time it is chosen to use FFTW in a multi-threaded way,
147 // it will try to use as many cores as possible.
148 
149 // In this class a forward transform is defined as one that goes from the real
150 // to the complex (or the time to frequency) domain. In a forward transform the
151 // sign of the exponent is negative and no scaling is done on the output. The
152 // backward transform goes from the complex to the real (or the frequency to
153 // the time) domain. The sign of the exponent is positive and the result is
154 // always scaled by 1/N were N is the total number of elements in the Array.
155 
156 // The origin of the transform is defined as the point where setting only that
157 // element to one, and then doing a forward transform results in an Array that
158 // is all one. The <src>fft</src> member functions in this class all assume
159 // that the origin of the Transform is at the centre of the Array ie. at
160 // <src>[nx/2,ny/2,...]</src> were the indexing begins at zero. Because the
161 // fftpack software assumes the origin of the transform is at the first element
162 // ie.,<src>[0,0,...]</src> this class flips the data in the Array around to
163 // compensate. For fftpack this flipping takes about one 20% of the total
164 // transform time, while for FFTW it can easily exceed the transform time.
165 // Flipping can be avoided by using the <src>fft0</src> member
166 // functions which do not flip the data.
167 
168 // Some of the member functions in this class scramble the input Array,
169 // possibly by flipping the quandrants of the data although this is not
170 // guaranteed. Modification of the input Array can be avoided, at the expense
171 // of copying the data to temporary storage, by either:
172 // <ul> <li> Ensuring the input Array is a const Array.
173 // <li> Setting the constInput Flag to True.
174 // </ul>
175 // The latter option is provided to avoid users having to cast non-const
176 // Arrays to const ones in order to prevent there input Array from being
177 // scrambled.
178 
179 // <note role=warning> This class assumes that a Complex array is stored as
180 // pairs of floating point numbers, with no intervening gaps, and with the real
181 // component first ie., <src>[re0,im0,re1,im1, ...]</src>. This means that the
182 // following type casts work,
183 // <srcblock>
184 // S * complexPtr;
185 // T * realPtr = (T * ) complexPtr;
186 // </srcblock>
187 // and allow a Complex number to be accessed as a pair of real numbers. If this
188 // assumption is bad then real Arrays will have to generated by copying the
189 // complex ones. Ultimately this assumption about Complex<->Real Array
190 // conversion should be put somewhere central like Array2Math.cc.
191 // </note>
192 // </synopsis>
193 
194 // <templating arg=T>
195 // <li> The T argument must be of type Float or Double. These are the only
196 // possible instantiations of this class.
197 // </templating>
198 
199 // <templating arg=S>
200 // <li> The S argument must be of type Complex, if T is Float, or DComplex, if T is
201 // Double. These are the only possible instantiations of this class.
202 // </templating>
203 //
204 // <example>
205 // Do a real to complex Transform of a 1-Dimensional Vector. The following
206 // example can trivially be extended to any number of dimensions.
207 // <srcblock>
208 // FFTServer<Float,Complex> server;
209 // Vector<Float> input(32);
210 // Vector<Complex> output(17);
211 // input = 0.0f;
212 // input(16) = 1.0f;
213 // cout << "Input:" << input << endl;
214 // server.fft(output, input);
215 // cout << "Output:" << output << endl;
216 // </srcblock>
217 // </example>
218 //
219 // <thrown>
220 // <li> AipsError: If the input and output Array have bad or incompatible
221 // shapes. See the individual function descriptions for what Array shapes are
222 // required.
223 // </thrown>
224 //
225 // <todo asof="1997/10/22">
226 // <li> The time taken to flip the Array can be reduced, if all the Array
227 // dimensions are even, by pre-multiplying the every other element on the
228 // input Array by -1. Then no flipping needs to be done on the output Array.
229 // </todo>
230 
231 template<class T, class S> class FFTServer
232 {
233 public:
234 
235  // The default constructor. The server will automatically resize to do
236  // transforms of the appropriate length when necessary.
237  FFTServer();
238 
239  // Initialise the server to do transforms on Arrays of the specified
240  // shape. The server will, however, resize to do transforms of other lengths
241  // if necessary. See the resize function for a description of the
242  // TransformType enumerator.
243  FFTServer(const IPosition & fftSize,
244  const FFTEnums::TransformType transformType
246 
247  // copy constructor. The copied server is initialised to do transforms of the
248  // same length as the other server. Uses copy (and not reference) semantics
249  // so that changing the transform length of one server does not affect the
250  // other server.
251  FFTServer(const FFTServer<T,S> & other);
252 
253  // destructor
254  ~FFTServer();
255 
256  // The assignment operator which does the same thing as the copy
257  // constructor.
258  FFTServer<T,S> & operator=(const FFTServer<T,S> & other);
259 
260  // Modify the FFTServer object to do transforms of the supplied shape. The
261  // amount of internal storage, and the initialisation, depends on the type of
262  // transform that will be done. The transform type is specified with the
263  // TransformTypes enumerator. Currently there is no difference in
264  // initialisation for the COMPLEXTOREAL and REALTOCOMPLEX transforms. The
265  // shape argument is the shape of the real array (or complex one if complex
266  // to complex transforms are being done). In general it is not necessary to
267  // use this function as all the fft & fft0 functions will automatically
268  // resize the server, if necessary, to match their input arguments.
269  void resize(const IPosition & fftSize,
270  const FFTEnums::TransformType transformType
272 
273  // Real to complex fft. The origin of the transform is in the centre of the
274  // Array. Because of the Hermitian property the output Array only contains
275  // half of the complex result. The output Array must either have no elements
276  // or be a size that is appropriate to the input Array size,
277  // ie. <src>shape = [(nx+2)/2, ny, nz,...]</src>. Otherwise an AipsError is
278  // thrown. See the synopsis for a description of the constInput flag.
279  // <group>
280  void fft(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
281  void fft(Array<S> & cResult, const Array<T> & rData);
282  // </group>
283 
284  // Complex to real fft. The origin of the transform is in the centre of the
285  // Array. Because of the Hermitian property the input Array only contains
286  // half of the complex values. The output Array must either have no elements,
287  // or be a size that is appropriate to the input Array size ie.,<br>
288  // <src>shape = [2*cx-2, cy, cz,...]</src> or <br>
289  // <src>shape = [2*cx-1, cy, cz,...]</src>. <br>
290  // Otherwise an AipsError is thrown. See the description in the synopsis for
291  // the algorithm used to choose between the two possible output shapes and a
292  // description of the constInput Flag.
293  // <group>
294  void fft(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
295  void fft(Array<T> & rResult, const Array<S> & cData);
296  // </group>
297 
298  // Complex to complex in-place fft. The origin of the transform is in the
299  // centre of the Array. The direction of the transform is controlled by the
300  // toFrequency variable. If True then a forward, or time to frequency,
301  // transform is performed. If False a backward or frequency to time transform
302  // is done. Scaling is always done on the backward transform.
303  void fft(Array<S> & cValues, const Bool toFrequency=True);
304 
305  // Complex to complex fft. The origin of the transform is in the centre of
306  // the Array. The direction of the transform is controlled by the toFrequency
307  // variable. If True then a forward, or time to frequency, transform is
308  // performed. If False a backward or frequency to time transform is
309  // done. Scaling is always done on the backward transform. The output Array
310  // must either either contain no elements or be the same as the input Array,
311  // ie. <src>shape = [cx, cy, cz,...]</src>. Otherwise an AipsError is
312  // thrown.
313  void fft(Array<S> & cResult, const Array<S> & cData,
314  const Bool toFrequency=True);
315 
316  // The <src>fft0</src> functions are equivalent to the <src>fft</src>
317  // functions described above except that the origin of the transform is the
318  // first element of the Array, ie. [0,0,0...], rather than the centre
319  // element, ie [nx/2, ny/2, nz/2, ...]. As the underlying functions
320  // assume that the origin of the transform is the first element these
321  // routines are in general faster than the equivalent ones with the origin
322  // at the centre of the Array.
323  // <group>
324  void fft0(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
325  void fft0(Array<S> & cResult, const Array<T> & rData);
326  void fft0(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
327  void fft0(Array<T> & rResult, const Array<S> & cData);
328  void fft0(Array<S> & cValues, const Bool toFrequency=True);
329  void fft0(Array<S> & cResult, const Array<S> & cData,
330  const Bool toFrequency=True);
331  //# void fft0(Array<T> & rValues, const Bool toFrequency=True);
332 
333  // </group>
334  //# Flips the quadrants in a complex Array so that the point at
335  //# cData.shape()/2 moves to the origin. This moves, for example, the point
336  //# at [8,3] to the origin ([0,0]) in an array of shape [16,7]. Usually two
337  //# flips will restore an Array to its original state. But for Array's
338  //# where one or more dimension is an odd length two flips do NOT restore
339  //# the data to its original state. So the when toZero=False this routine
340  //# does an unflip operation (ie moves the data at [0,0] to the centre) and
341  //# restores the data to its original state for odd length arrays. When
342  //# passed a Hermitian Array where half the complex plane is implicit (eg as
343  //# produced by a real->complex Transform) it is not necessary to flip the
344  //# first dimension of the Array. In this case the isHermitian flag should
345  //# be set to True. For complex<->complex transforms this should be False.
346  // <group>
347  void flip(Array<T> & rData, const Bool toZero, const Bool isHermitian);
348  void flip(Array<S> & cData, const Bool toZero, const Bool isHermitian);
349  // </group>
350 
351  // N-D in-place complex->complex FFT shift (FFT - phase-mult - inverse FFT)
352  // If toFrequency is true, the first FFT will be from time to frequency.
353  // relshift is the freq shift normalised to the bandwidth.
354  // Only transform over selected dimension. Iterate over the others.
355  void fftshift(Array<S> & cValues, const uInt& whichAxis,
356  const Double& relshift, const Bool toFrequency=True);
357 
358  // N-D complex->complex FFT shift (FFT - phase-mult - inverse FFT)
359  // with flagging.
360  // If toFrequency is true, the first FFT will be from time to frequency.
361  // relshift is the freq shift normalised to the bandwidth.
362  // Only transform over selected dimension. Iterate over the others.
363  void fftshift(Array<S> & outValues, Array<Bool> & outFlags,
364  const Array<S> & cValues, const Array<Bool>& inFlags,
365  const uInt& whichAxis,
366  const Double& relshift,
367  const Bool goodIsTrue=False,
368  const Bool toFrequency=True);
369 
370  // N-D real->real FFT shift (FFT to complex - phase-mult - inverse FFT)
371  // with flagging.
372  // relshift is the freq shift normalised to the bandwidth.
373  // Only transform over selected dimension. Iterate over the others.
374  void fftshift(Array<T> & outValues, Array<Bool> & outFlags,
375  const Array<T> & rValues, const Array<Bool>& inFlags,
376  const uInt& whichAxis,
377  const Double& relshift,
378  const Bool goodIsTrue=False);
379 
380 private:
381  //# finds the shape of the output array when doing complex->real transforms
382  IPosition determineShape(const IPosition & rShape, const Array<S> & cData);
383 
384  //# Data members.
385  // The size of the last FFT done by this object
387  // Whether the last FFT was complex<->complex or not
389  //# twiddle factors and factorisations used by fftpack
391  // buffer for copying non-contigious arrays to contigious ones. This is done
392  // so that the FFT's have a better chance of fitting into cache and hence
393  // going faster.
394  // This buffer is also used as temporary storage when flipping the data.
396  // FFTW specific members. Do not harm if FFTPack is used.
398  std::vector<T> itsWorkIn;
399  std::vector<S> itsWorkOut;
400  std::vector<S> itsWorkC2C;
401 };
402 
403 
404 } //# NAMESPACE CASACORE - END
405 
406 //# Do NOT include the .tcc file here like done for other templated classes.
407 //# The instantiations are done explicitly.
408 //# In this way the HAVE_FFTW ifdef is only used in .cc files and does
409 //# not appear in headers, so other packages using FFTServer do not need
410 //# to (un)set HAVE_FFTW.
411 
412 #endif
A Vector of integers, for indexing into Array<T> objects.
Definition: IPosition.h:119
Inverse Complex to Complex transforms.
Definition: FFTServer.h:55
FFTW itsFFTW
FFTW specific members.
Definition: FFTServer.h:397
PtrHolder< T > & operator=(const PtrHolder< T > &other)
IPosition itsSize
The size of the last FFT done by this object.
Definition: FFTServer.h:386
Real to Real transforms with symmetric Arrays (not used)
Definition: FFTServer.h:61
C++ interface to the FFTWw library.
Definition: FFTW.h:57
double Double
Definition: aipstype.h:55
std::vector< S > itsWorkOut
Definition: FFTServer.h:399
A class with methods for Fast Fourier Transforms.
Definition: FFTServer.h:231
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
Real to Complex or Complex to Real transforms.
Definition: FFTServer.h:57
const Bool False
Definition: aipstype.h:44
A drop-in replacement for Block<T*>.
Definition: Block.h:814
template <class T, class U> class vector;
Definition: Array.h:166
std::vector< T > itsWorkIn
Definition: FFTServer.h:398
Real to Complex or Complex to Real transforms.
Definition: FFTServer.h:59
PtrBlock< Block< T > * > itsWork
Definition: FFTServer.h:390
FFTEnums::TransformType itsTransformType
Whether the last FFT was complex<->complex or not.
Definition: FFTServer.h:388
std::vector< S > itsWorkC2C
Definition: FFTServer.h:400
Block< S > itsBuffer
buffer for copying non-contigious arrays to contigious ones.
Definition: FFTServer.h:395
Forward Complex to Complex transforms.
Definition: FFTServer.h:53
const Bool True
Definition: aipstype.h:43
this file contains all the compiler specific defines
Definition: mainpage.dox:28
unsigned int uInt
Definition: aipstype.h:51