Loading...
Searching...
No Matches
FFTServer.h
Go to the documentation of this file.
//# FFTServer.h: A class with methods for Fast Fourier Transforms
//# Copyright (C) 1994,1995,1996,1997,1999,2003
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//# Internet email: casa-feedback@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
#ifndef SCIMATH_FFTSERVER_H
#define SCIMATH_FFTSERVER_H
//# Includes
#include <casacore/casa/aips.h>
#include <casacore/scimath/Mathematics/FFTW.h>
#include <casacore/casa/Arrays/IPosition.h>
#include <casacore/casa/Arrays/ArrayFwd.h>
#include <casacore/casa/Containers/Block.h>
#include <vector>
namespace casacore { //# NAMESPACE CASACORE - BEGIN
// <summary>Lists the different types of FFT's that can be done</summary>
// <synopsis>This enumerator is brought out as a separate class because g++
// currently cannot handle enumerators in a templated class. When it can this
// class will go away and this enumerator moved into the FFTServer
// class</synopsis>
class FFTEnums {
public:
enum TransformType {
// Forward Complex to Complex transforms.
COMPLEX,
// Inverse Complex to Complex transforms.
INVCOMPLEX,
// Real to Complex or Complex to Real transforms.
REALTOCOMPLEX,
// Real to Complex or Complex to Real transforms.
COMPLEXTOREAL,
// Real to Real transforms with symmetric Arrays (not used)
REALSYMMETRIC
};
};
// <summary>A class with methods for Fast Fourier Transforms</summary>
// <use visibility=export>
// <reviewed reviewer="wbrouw" date="1997/10/29" tests="tFFTServer">
// </reviewed>
// <prerequisite>
// <li> Basic concepts of Fast Fourier Transforms.
// <li> <linkto module=Arrays>The Arrays module</linkto>
// </prerequisite>
// <etymology> The FFTServer class, can do Fast Fourier Transforms of
// any length and dimensionality.
// </etymology>
// <synopsis>
// The FFTServer class provides methods for performing n-dimensional Fast
// Fourier Transforms with real and complex Array's of arbitrary size and
// dimensionality. It can do either real to complex, complex to real, or
// complex to complex transforms with the "origin" of the transform either at
// the centre of the Array or at the first element.
// Because the output from a real to complex transform is Hermitian only half
// of the complex result is returned. Similarly with a complex to real
// transform only half of the complex plane is required, the other half is
// implicitly assumed to be the complex conjugate of the supplied half-plane.
// <note role=warning> The complex to real transform does not check that the
// imaginary component of the values where u=0 are zero</note>
// This class can be initialised with a shape that indicates the length of the
// transforms that will be performed, and whether they are going to be
// real<->complex transforms or complex<->complex ones. This initialisation
// sets up a variety of internal buffers and computes factorizations and
// twiddle factors used during the transform. The initialised transform shape
// is always compared with the shape of the supplied arguments when a transform
// is done and the FFTServer class will automatically resize itself if
// necessary. So the default constructor is perfectly safe to use.
// With any transform the output Array must either be the correct shape for the
// desired output or zero length (ie not contain any elements). If it is zero
// length then it will be resized to the correct shape. For a complex->complex
// transform the output Array will be the same shape as the input Array. For a
// real->complex transform the output Array will be the same size as the input
// Array except in the first dimension which will have a length of (nx+2)/2. So
// if nx=7 the output length will be 4 and if nx=8 the output length will be 5,
// on the first axis. nx is the length of the first axis on the <em>real</em>
// input Array and cx (which is used later) is the length of the first axis on
// the <em>complex</em> input Array.
// <strong>For complex to real transforms the output length on the first axis
// is not uniquely defined by the shape of the complex input
// Array</strong>. This class uses the following algorithm to work out the
// length of the first axis on the output Array.
// <ul>
// <li> If the size of the output Array is non-zero then its shape must match
// the size of the input Array except for the first axis. The length of the
// first axis must either be 2*cx-2 or 2*cx-1 and this determines the length of
// the transform on the first axis.
// <li> If the size of the output Array is zero then scan the imaginary
// components of the values at the end of the first axis on the input Array (ie
// at <src>[cx-1,....]</src> If any of these are non-zero the output Array
// will have an odd length.
// <li> Otherwise if all the imaginary components described above are zero then
// look at the current size of the FFTServer object (either defined at
// construction time or with the resize function). If it matches the size of
// the input Array except for the first axis and if the length on this axis is
// either 2*cx-2 or 2*cx-1 then use that to determine the size of the output
// Array.
// <li> Otherwise assume the output Array will an even length of 2*cx-2 on its
// first axis.
// </ul>
// This class does transforms using
// the highly optimized FFTW package.
// <br>
// <em> P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
// (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83. </em><br>
// <br>If at build time it is chosen to use FFTW in a multi-threaded way,
// it will try to use as many cores as possible.
// In this class a forward transform is defined as one that goes from the real
// to the complex (or the time to frequency) domain. In a forward transform the
// sign of the exponent is negative and no scaling is done on the output. The
// backward transform goes from the complex to the real (or the frequency to
// the time) domain. The sign of the exponent is positive and the result is
// always scaled by 1/N were N is the total number of elements in the Array.
// The origin of the transform is defined as the point where setting only that
// element to one, and then doing a forward transform results in an Array that
// is all one. The <src>fft</src> member functions in this class all assume
// that the origin of the Transform is at the centre of the Array ie. at
// <src>[nx/2,ny/2,...]</src> were the indexing begins at zero. Because the
// fftpack software assumes the origin of the transform is at the first element
// ie.,<src>[0,0,...]</src> this class flips the data in the Array around to
// compensate. For fftpack this flipping takes about one 20% of the total
// transform time, while for FFTW it can easily exceed the transform time.
// Flipping can be avoided by using the <src>fft0</src> member
// functions which do not flip the data.
// Some of the member functions in this class scramble the input Array,
// possibly by flipping the quandrants of the data although this is not
// guaranteed. Modification of the input Array can be avoided, at the expense
// of copying the data to temporary storage, by either:
// <ul> <li> Ensuring the input Array is a const Array.
// <li> Setting the constInput Flag to True.
// </ul>
// The latter option is provided to avoid users having to cast non-const
// Arrays to const ones in order to prevent there input Array from being
// scrambled.
// <note role=warning> This class assumes that a Complex array is stored as
// pairs of floating point numbers, with no intervening gaps, and with the real
// component first ie., <src>[re0,im0,re1,im1, ...]</src>. This means that the
// following type casts work,
// <srcblock>
// S * complexPtr;
// T * realPtr = (T * ) complexPtr;
// </srcblock>
// and allow a Complex number to be accessed as a pair of real numbers. If this
// assumption is bad then real Arrays will have to generated by copying the
// complex ones. Ultimately this assumption about Complex<->Real Array
// conversion should be put somewhere central like Array2Math.cc.
// </note>
// </synopsis>
// <templating arg=T>
// <li> The T argument must be of type Float or Double. These are the only
// possible instantiations of this class.
// </templating>
// <templating arg=S>
// <li> The S argument must be of type Complex, if T is Float, or DComplex, if T is
// Double. These are the only possible instantiations of this class.
// </templating>
//
// <example>
// Do a real to complex Transform of a 1-Dimensional Vector. The following
// example can trivially be extended to any number of dimensions.
// <srcblock>
// FFTServer<Float,Complex> server;
// Vector<Float> input(32);
// Vector<Complex> output(17);
// input = 0.0f;
// input(16) = 1.0f;
// cout << "Input:" << input << endl;
// server.fft(output, input);
// cout << "Output:" << output << endl;
// </srcblock>
// </example>
//
// <thrown>
// <li> AipsError: If the input and output Array have bad or incompatible
// shapes. See the individual function descriptions for what Array shapes are
// required.
// </thrown>
//
// <todo asof="1997/10/22">
// <li> The time taken to flip the Array can be reduced, if all the Array
// dimensions are even, by pre-multiplying the every other element on the
// input Array by -1. Then no flipping needs to be done on the output Array.
// </todo>
template<class T, class S> class FFTServer
{
public:
// The default constructor. The server will automatically resize to do
// transforms of the appropriate length when necessary.
FFTServer();
// Initialise the server to do transforms on Arrays of the specified
// shape. The server will, however, resize to do transforms of other lengths
// if necessary. See the resize function for a description of the
// TransformType enumerator.
FFTServer(const IPosition & fftSize,
const FFTEnums::TransformType transformType
= FFTEnums::REALTOCOMPLEX);
// copy constructor. The copied server is initialised to do transforms of the
// same length as the other server. Uses copy (and not reference) semantics
// so that changing the transform length of one server does not affect the
// other server.
FFTServer(const FFTServer<T,S> & other);
// destructor
~FFTServer();
// The assignment operator which does the same thing as the copy
// constructor.
FFTServer<T,S> & operator=(const FFTServer<T,S> & other);
// Modify the FFTServer object to do transforms of the supplied shape. The
// amount of internal storage, and the initialisation, depends on the type of
// transform that will be done. The transform type is specified with the
// TransformTypes enumerator. Currently there is no difference in
// initialisation for the COMPLEXTOREAL and REALTOCOMPLEX transforms. The
// shape argument is the shape of the real array (or complex one if complex
// to complex transforms are being done). In general it is not necessary to
// use this function as all the fft & fft0 functions will automatically
// resize the server, if necessary, to match their input arguments.
void resize(const IPosition & fftSize,
const FFTEnums::TransformType transformType
= FFTEnums::REALTOCOMPLEX);
// Real to complex fft. The origin of the transform is in the centre of the
// Array. Because of the Hermitian property the output Array only contains
// half of the complex result. The output Array must either have no elements
// or be a size that is appropriate to the input Array size,
// ie. <src>shape = [(nx+2)/2, ny, nz,...]</src>. Otherwise an AipsError is
// thrown. See the synopsis for a description of the constInput flag.
// <group>
void fft(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
void fft(Array<S> & cResult, const Array<T> & rData);
// </group>
// Complex to real fft. The origin of the transform is in the centre of the
// Array. Because of the Hermitian property the input Array only contains
// half of the complex values. The output Array must either have no elements,
// or be a size that is appropriate to the input Array size ie.,<br>
// <src>shape = [2*cx-2, cy, cz,...]</src> or <br>
// <src>shape = [2*cx-1, cy, cz,...]</src>. <br>
// Otherwise an AipsError is thrown. See the description in the synopsis for
// the algorithm used to choose between the two possible output shapes and a
// description of the constInput Flag.
// <group>
void fft(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
void fft(Array<T> & rResult, const Array<S> & cData);
// </group>
// Complex to complex in-place fft. The origin of the transform is in the
// centre of the Array. The direction of the transform is controlled by the
// toFrequency variable. If True then a forward, or time to frequency,
// transform is performed. If False a backward or frequency to time transform
// is done. Scaling is always done on the backward transform.
void fft(Array<S> & cValues, const Bool toFrequency=True);
// Complex to complex fft. The origin of the transform is in the centre of
// the Array. The direction of the transform is controlled by the toFrequency
// variable. If True then a forward, or time to frequency, transform is
// performed. If False a backward or frequency to time transform is
// done. Scaling is always done on the backward transform. The output Array
// must either either contain no elements or be the same as the input Array,
// ie. <src>shape = [cx, cy, cz,...]</src>. Otherwise an AipsError is
// thrown.
void fft(Array<S> & cResult, const Array<S> & cData,
const Bool toFrequency=True);
// The <src>fft0</src> functions are equivalent to the <src>fft</src>
// functions described above except that the origin of the transform is the
// first element of the Array, ie. [0,0,0...], rather than the centre
// element, ie [nx/2, ny/2, nz/2, ...]. As the underlying functions
// assume that the origin of the transform is the first element these
// routines are in general faster than the equivalent ones with the origin
// at the centre of the Array.
// <group>
void fft0(Array<S> & cResult, Array<T> & rData, const Bool constInput=False);
void fft0(Array<S> & cResult, const Array<T> & rData);
void fft0(Array<T> & rResult, Array<S> & cData, const Bool constInput=False);
void fft0(Array<T> & rResult, const Array<S> & cData);
void fft0(Array<S> & cValues, const Bool toFrequency=True);
void fft0(Array<S> & cResult, const Array<S> & cData,
const Bool toFrequency=True);
//# void fft0(Array<T> & rValues, const Bool toFrequency=True);
// </group>
//# Flips the quadrants in a complex Array so that the point at
//# cData.shape()/2 moves to the origin. This moves, for example, the point
//# at [8,3] to the origin ([0,0]) in an array of shape [16,7]. Usually two
//# flips will restore an Array to its original state. But for Array's
//# where one or more dimension is an odd length two flips do NOT restore
//# the data to its original state. So the when toZero=False this routine
//# does an unflip operation (ie moves the data at [0,0] to the centre) and
//# restores the data to its original state for odd length arrays. When
//# passed a Hermitian Array where half the complex plane is implicit (eg as
//# produced by a real->complex Transform) it is not necessary to flip the
//# first dimension of the Array. In this case the isHermitian flag should
//# be set to True. For complex<->complex transforms this should be False.
// <group>
void flip(Array<T> & rData, const Bool toZero, const Bool isHermitian);
void flip(Array<S> & cData, const Bool toZero, const Bool isHermitian);
// </group>
// N-D in-place complex->complex FFT shift (FFT - phase-mult - inverse FFT)
// If toFrequency is true, the first FFT will be from time to frequency.
// relshift is the freq shift normalised to the bandwidth.
// Only transform over selected dimension. Iterate over the others.
void fftshift(Array<S> & cValues, const uInt& whichAxis,
const Double& relshift, const Bool toFrequency=True);
// N-D complex->complex FFT shift (FFT - phase-mult - inverse FFT)
// with flagging.
// If toFrequency is true, the first FFT will be from time to frequency.
// relshift is the freq shift normalised to the bandwidth.
// Only transform over selected dimension. Iterate over the others.
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 real->real FFT shift (FFT to complex - phase-mult - inverse FFT)
// with flagging.
// relshift is the freq shift normalised to the bandwidth.
// Only transform over selected dimension. Iterate over the others.
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);
private:
//# finds the shape of the output array when doing complex->real transforms
IPosition determineShape(const IPosition & rShape, const Array<S> & cData);
//# Data members.
// The size of the last FFT done by this object
IPosition itsSize;
// Whether the last FFT was complex<->complex or not
FFTEnums::TransformType itsTransformType;
// buffer for copying non-contigious arrays to contigious ones. This is done
// so that the FFT's have a better chance of fitting into cache and hence
// going faster.
// This buffer is also used as temporary storage when flipping the data.
Block<S> itsBuffer;
// FFTW specific members.
FFTW itsFFTW;
std::vector<T> itsWorkIn;
std::vector<S> itsWorkOut;
std::vector<S> itsWorkC2C;
};
} //# NAMESPACE CASACORE - END
//# Do NOT include the .tcc file here like done for other templated classes.
//# The instantiations are done explicitly.
//# In this way the HAVE_FFTW ifdef is only used in .cc files and does
//# not appear in headers, so other packages using FFTServer do not need
//# to (un)set HAVE_FFTW.
#endif
casacore::Array
Definition Array.h:155
casacore::Block
simple 1-D array
Definition Block.h:198
casacore::FFTEnums
Definition FFTServer.h:45
casacore::FFTEnums::TransformType
TransformType
Definition FFTServer.h:47
casacore::FFTEnums::REALTOCOMPLEX
@ REALTOCOMPLEX
Real to Complex or Complex to Real transforms.
Definition FFTServer.h:53
casacore::FFTEnums::REALSYMMETRIC
@ REALSYMMETRIC
Real to Real transforms with symmetric Arrays (not used).
Definition FFTServer.h:57
casacore::FFTEnums::COMPLEXTOREAL
@ COMPLEXTOREAL
Real to Complex or Complex to Real transforms.
Definition FFTServer.h:55
casacore::FFTEnums::COMPLEX
@ COMPLEX
Forward Complex to Complex transforms.
Definition FFTServer.h:49
casacore::FFTEnums::INVCOMPLEX
@ INVCOMPLEX
Inverse Complex to Complex transforms.
Definition FFTServer.h:51
casacore::FFTServer::flip
void flip(Array< S > &cData, const Bool toZero, const Bool isHermitian)
casacore::FFTServer::fft
void fft(Array< S > &cValues, const Bool toFrequency=True)
Complex to complex in-place fft.
casacore::FFTServer::itsWorkOut
std::vector< S > itsWorkOut
Definition FFTServer.h:391
casacore::FFTServer::determineShape
IPosition determineShape(const IPosition &rShape, const Array< S > &cData)
casacore::FFTServer::fftshift
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.
casacore::FFTServer::itsBuffer
Block< S > itsBuffer
buffer for copying non-contigious arrays to contigious ones.
Definition FFTServer.h:387
casacore::FFTServer::fft
void fft(Array< T > &rResult, const Array< S > &cData)
casacore::FFTServer::FFTServer
FFTServer(const IPosition &fftSize, const FFTEnums::TransformType transformType=FFTEnums::REALTOCOMPLEX)
Initialise the server to do transforms on Arrays of the specified shape.
casacore::FFTServer::fft
void fft(Array< T > &rResult, Array< S > &cData, const Bool constInput=False)
Complex to real fft.
casacore::FFTServer::itsWorkIn
std::vector< T > itsWorkIn
Definition FFTServer.h:390
casacore::FFTServer::itsSize
IPosition itsSize
The size of the last FFT done by this object.
Definition FFTServer.h:380
casacore::FFTServer::fft0
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 ...
casacore::FFTServer::itsTransformType
FFTEnums::TransformType itsTransformType
Whether the last FFT was complex<->complex or not.
Definition FFTServer.h:382
casacore::FFTServer::fft
void fft(Array< S > &cResult, const Array< T > &rData)
casacore::FFTServer::resize
void resize(const IPosition &fftSize, const FFTEnums::TransformType transformType=FFTEnums::REALTOCOMPLEX)
Modify the FFTServer object to do transforms of the supplied shape.
casacore::FFTServer::fft0
void fft0(Array< S > &cResult, const Array< S > &cData, const Bool toFrequency=True)
casacore::FFTServer::operator=
FFTServer< T, S > & operator=(const FFTServer< T, S > &other)
The assignment operator which does the same thing as the copy constructor.
casacore::FFTServer::flip
void flip(Array< T > &rData, const Bool toZero, const Bool isHermitian)
casacore::FFTServer::fft
void fft(Array< S > &cResult, const Array< S > &cData, const Bool toFrequency=True)
Complex to complex fft.
casacore::FFTServer::fft0
void fft0(Array< T > &rResult, const Array< S > &cData)
casacore::FFTServer::fft0
void fft0(Array< T > &rResult, Array< S > &cData, const Bool constInput=False)
casacore::FFTServer::fftshift
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.
casacore::FFTServer::fftshift
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,...
casacore::FFTServer::FFTServer
FFTServer(const FFTServer< T, S > &other)
copy constructor.
casacore::FFTServer::itsWorkC2C
std::vector< S > itsWorkC2C
Definition FFTServer.h:392
casacore::FFTServer::FFTServer
FFTServer()
The default constructor.
casacore::FFTServer::fft0
void fft0(Array< S > &cResult, const Array< T > &rData)
casacore::FFTServer::fft0
void fft0(Array< S > &cValues, const Bool toFrequency=True)
casacore::FFTServer::fft
void fft(Array< S > &cResult, Array< T > &rData, const Bool constInput=False)
Real to complex fft.
casacore::FFTServer::~FFTServer
~FFTServer()
destructor
casacore::FFTServer::itsFFTW
FFTW itsFFTW
FFTW specific members.
Definition FFTServer.h:389
casacore::FFTW
Definition FFTW.h:55
casacore::IPosition
Definition IPosition.h:117
casacore
this file contains all the compiler specific defines
Definition mainpage.dox:28
casacore::False
const Bool False
Definition aipstype.h:42
casacore::uInt
unsigned int uInt
Definition aipstype.h:49
casacore::Bool
bool Bool
Define the standard types used by Casacore.
Definition aipstype.h:40
casacore::True
const Bool True
Definition aipstype.h:41
casacore::Double
double Double
Definition aipstype.h:53