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