A class with methods for Fast Fourier Transforms. More...

#include <FFTServer.h>

Public Member Functions
	FFTServer ()
	The default constructor.
	FFTServer (const IPosition &fftSize, const FFTEnums::TransformType transformType=FFTEnums::REALTOCOMPLEX)
	Initialise the server to do transforms on Arrays of the specified shape.
	FFTServer (const FFTServer< T, S > &other)
	copy constructor.
	~FFTServer ()
	destructor
FFTServer< T, S > &	operator= (const FFTServer< T, S > &other)
	The assignment operator which does the same thing as the copy constructor.
void	resize (const IPosition &fftSize, const FFTEnums::TransformType transformType=FFTEnums::REALTOCOMPLEX)
	Modify the FFTServer object to do transforms of the supplied shape.
void	fft (Array< S > &cResult, Array< T > &rData, const Bool constInput=False)
	Real to complex fft.
void	fft (Array< S > &cResult, const Array< T > &rData)
void	fft (Array< T > &rResult, Array< S > &cData, const Bool constInput=False)
	Complex to real fft.
void	fft (Array< T > &rResult, const Array< S > &cData)
void	fft (Array< S > &cValues, const Bool toFrequency=True)
	Complex to complex in-place fft.
void	fft (Array< S > &cResult, const Array< S > &cData, const Bool toFrequency=True)
	Complex to complex fft.
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 transform is the first element of the Array, ie.
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	flip (Array< T > &rData, const Bool toZero, const Bool isHermitian)

void	flip (Array< S > &cData, const Bool toZero, const Bool isHermitian)
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, the first FFT will be from time to frequency.
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< 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.
Private Member Functions
IPosition	determineShape (const IPosition &rShape, const Array< S > &cData)
Private Attributes
IPosition	itsSize
	The size of the last FFT done by this object.
FFTEnums::TransformType	itsTransformType
	Whether the last FFT was complex<->complex or not.
PtrBlock< Block< T > * >	itsWork
Block< S >	itsBuffer
	buffer for copying non-contigious arrays to contigious ones.
FFTW	itsFFTW
	FFTW specific members.
std::vector< T >	itsWorkIn
std::vector< S >	itsWorkOut
std::vector< S >	itsWorkC2C

Detailed Description

template<class T, class S>
class casa::FFTServer< T, S >

A class with methods for Fast Fourier Transforms.

Intended use:

Public interface

 <h3>Review Status</h3><dl><dt>Reviewed By:<dd>wbrouw<dt>Date Reviewed:<dd>1997/10/29<dt>Test programs:<dd>tFFTServer</dl>

Prerequisite

Basic concepts of Fast Fourier Transforms.
The Arrays module

Etymology

The FFTServer class, can do Fast Fourier Transforms of any length and dimensionality.

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.
Warning: The complex to real transform does not check that the imaginary component of the values where u=0 are zero

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 real input Array and cx (which is used later) is the length of the first axis on the complex input Array.

For complex to real transforms the output length on the first axis is not uniquely defined by the shape of the complex input Array. This class uses the following algorithm to work out the length of the first axis on the output Array.

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.
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 [cx-1,....] If any of these are non-zero the output Array will have an odd length.
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.
Otherwise assume the output Array will an even length of 2*cx-2 on its first axis.

This class does transforms using the widely used FORTRAN fftpack package or the highly optimized FFTW package (to be chosen at build time).
P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
The fftpack package only does one dimensional transforms and this class decomposes multi-dimensional transforms into a series of 1-dimensional ones.
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 fft member functions in this class all assume that the origin of the Transform is at the centre of the Array ie. at [nx/2,ny/2,...] were the indexing begins at zero. Because the fftpack software assumes the origin of the transform is at the first element ie.,[0,0,...] 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 fft0 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:

Ensuring the input Array is a const Array.
Setting the constInput Flag to True.

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.

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;, [re0,im0,re1,im1, ;;;]; This means that the following type casts work,

S * complexPtr; T * realPtr = (T * ) complexPtr;

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;

Template Type Argument Requirements (T)

The T argument must be of type Float or Double. These are the only possible instantiations of this class.

Template Type Argument Requirements (S)

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.

Example

Do a real to complex Transform of a 1-Dimensional Vector. The following example can trivially be extended to any number of dimensions.

    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;

Thrown Exceptions

AipsError: If the input and output Array have bad or incompatible shapes. See the individual function descriptions for what Array shapes are required.

To Do

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.

Definition at line 231 of file FFTServer.h.

Constructor & Destructor Documentation

template<class T, class S>

casa::FFTServer< T, S >::FFTServer ( )

The default constructor.

The server will automatically resize to do transforms of the appropriate length when necessary.

template<class T, class S>

casa::FFTServer< T, S >::FFTServer	(	const IPosition &	fftSize,
		const FFTEnums::TransformType	transformType = `FFTEnums::REALTOCOMPLEX`
	)

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.

template<class T, class S>

casa::FFTServer< T, S >::FFTServer ( const FFTServer< T, S > & other )

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.

template<class T, class S>

casa::FFTServer< T, S >::~FFTServer ( )

destructor

Member Function Documentation

template<class T, class S>

IPosition casa::FFTServer< T, S >::determineShape	(	const IPosition &	rShape,
		const Array< S > &	cData
	)		`[private]`

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< S > &	cResult,
		Array< T > &	rData,
		const Bool	constInput = `False`
	)

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. shape = [(nx+2)/2, ny, nz,...]. Otherwise an AipsError is thrown. See the synopsis for a description of the constInput flag.

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< S > &	cResult,
		const Array< T > &	rData
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< T > &	rResult,
		Array< S > &	cData,
		const Bool	constInput = `False`
	)

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.,
shape = [2*cx-2, cy, cz,...] or
shape = [2*cx-1, cy, cz,...].
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.

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< T > &	rResult,
		const Array< S > &	cData
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< S > &	cValues,
		const Bool	toFrequency = `True`
	)

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.

template<class T, class S>

void casa::FFTServer< T, S >::fft	(	Array< S > &	cResult,
		const Array< S > &	cData,
		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. shape = [cx, cy, cz,...]. Otherwise an AipsError is thrown.

template<class T, class S>

void casa::FFTServer< T, S >::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 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.

template<class T, class S>

void casa::FFTServer< T, S >::fft0	(	Array< S > &	cResult,
		const Array< T > &	rData
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft0	(	Array< T > &	rResult,
		Array< S > &	cData,
		const Bool	constInput = `False`
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft0	(	Array< T > &	rResult,
		const Array< S > &	cData
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft0	(	Array< S > &	cValues,
		const Bool	toFrequency = `True`
	)

template<class T, class S>

void casa::FFTServer< T, S >::fft0	(	Array< S > &	cResult,
		const Array< S > &	cData,
		const Bool	toFrequency = `True`
	)

template<class T, class S>

void casa::FFTServer< T, S >::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, 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.

template<class T, class S>

void casa::FFTServer< T, S >::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.

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.

template<class T, class S>

void casa::FFTServer< T, S >::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.

relshift is the freq shift normalised to the bandwidth. Only transform over selected dimension. Iterate over the others.

template<class T, class S>

void casa::FFTServer< T, S >::flip	(	Array< T > &	rData,
		const Bool	toZero,
		const Bool	isHermitian
	)

template<class T, class S>

void casa::FFTServer< T, S >::flip	(	Array< S > &	cData,
		const Bool	toZero,
		const Bool	isHermitian
	)

template<class T, class S>

FFTServer<T,S>& casa::FFTServer< T, S >::operator= ( const FFTServer< T, S > & other )

The assignment operator which does the same thing as the copy constructor.

template<class T, class S>

void casa::FFTServer< T, S >::resize	(	const IPosition &	fftSize,
		const FFTEnums::TransformType	transformType = `FFTEnums::REALTOCOMPLEX`
	)

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.

Member Data Documentation

template<class T, class S>

Block<S> casa::FFTServer< T, S >::itsBuffer [private]

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.

Definition at line 395 of file FFTServer.h.

template<class T, class S>

FFTW casa::FFTServer< T, S >::itsFFTW [private]

FFTW specific members.

Do not harm if FFTPack is used.

Definition at line 397 of file FFTServer.h.

template<class T, class S>

IPosition casa::FFTServer< T, S >::itsSize [private]

The size of the last FFT done by this object.

Definition at line 386 of file FFTServer.h.

template<class T, class S>

FFTEnums::TransformType casa::FFTServer< T, S >::itsTransformType [private]

Whether the last FFT was complex<->complex or not.

Definition at line 388 of file FFTServer.h.

template<class T, class S>

PtrBlock<Block<T> *> casa::FFTServer< T, S >::itsWork [private]

Definition at line 390 of file FFTServer.h.

template<class T, class S>

std::vector<S> casa::FFTServer< T, S >::itsWorkC2C [private]

Definition at line 400 of file FFTServer.h.

template<class T, class S>

std::vector<T> casa::FFTServer< T, S >::itsWorkIn [private]

Definition at line 398 of file FFTServer.h.

template<class T, class S>

std::vector<S> casa::FFTServer< T, S >::itsWorkOut [private]

Definition at line 399 of file FFTServer.h.

The documentation for this class was generated from the following file:

casacore/scimath/Mathematics/FFTServer.h

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

template<class T, class S> class casa::FFTServer< T, S >

Intended use:

Prerequisite

Etymology

Synopsis

Template Type Argument Requirements (T)

Template Type Argument Requirements (S)

Example

Thrown Exceptions

To Do

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<class T, class S>
class casa::FFTServer< T, S >