Loading...
Searching...
No Matches
FFT1D.h
Go to the documentation of this file.
// -*- Mode: C++; tab-width: 2; -*-
// vi: set ts=2:
//
#ifndef BALL_MATHS_TFFT1D_H
#define BALL_MATHS_TFFT1D_H
#ifndef BALL_COMMON_EXCEPTION_H
# include <BALL/COMMON/exception.h>
#endif
#ifndef BALL_DATATYPE_REGULARDATA1D_H
# include <BALL/DATATYPE/regularData1D.h>
#endif
#include <cmath>
#include <complex>
#include <fftw3.h>
#include <BALL/COMMON/constants.h>
#include <BALL/MATHS/fftwCommon.h>
namespace BALL
{
template <typename ComplexTraits>
class TFFT1D
: public TRegularData1D<std::complex<typename ComplexTraits::ComplexPrecision> >
{
public:
typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
typedef TRegularData1D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
BALL_CREATE(TFFT1D)
TFFT1D();
TFFT1D(const TFFT1D &data);
// AR: ldn is not any longer the binary logarithm but the absolute number of grid points
TFFT1D(Size ldn, double stepPhys = 1., double origin = 0., bool inFourierSpace = false);
virtual ~TFFT1D();
const TFFT1D& operator = (const TFFT1D& fft1d);
virtual void clear();
virtual void destroy();
bool operator == (const TFFT1D& fft1d) const;
// @name Accessors
void doFFT();
void doiFFT();
bool translate(double trans_origin);
bool setPhysStepWidth(double new_width);
double getPhysStepWidth() const;
double getFourierStepWidth() const;
double getPhysSpaceMin() const;
double getPhysSpaceMax() const;
double getFourierSpaceMin() const;
double getFourierSpaceMax() const;
Size getMaxIndex() const;
Size getNumberOfInverseTransforms() const;
double getGridCoordinates(Position position) const;
Complex getData(const double pos) const;
Complex getInterpolatedValue(const double pos) const;
void setData(double pos, Complex val);
Complex& operator [] (const double pos);
const Complex& operator [] (const double pos) const;
Complex& operator[](const Position& pos)
{
return TRegularData1D<Complex>::operator[](pos);
}
const Complex& operator[](const Position& pos) const
{
return TRegularData1D<Complex>::operator[](pos);
}
void setNumberOfFFTTransforms(Size num)
{
numPhysToFourier_ = num;
}
void setNumberOfiFFTTransforms(Size num)
{
numFourierToPhys_ = num;
}
bool isInFourierSpace() const;
Complex phase(const double pos) const;
protected:
Size length_;
bool inFourierSpace_;
Size numPhysToFourier_;
Size numFourierToPhys_;
double origin_;
double stepPhys_;
double stepFourier_;
double minPhys_;
double maxPhys_;
double minFourier_;
double maxFourier_;
typename ComplexTraits::FftwPlan planForward_;
typename ComplexTraits::FftwPlan planBackward_;
// AR: to control plan calculation with new fftw3
Complex *dataAdress_;
bool planCalculated_;
};
typedef TFFT1D<BALL_FFTW_DEFAULT_TRAITS> FFT1D;
// AR:
// const TRegularData1D<Complex>& operator << (TRegularData1D<Complex>& to, const TFFT1D& from)
//; ?????
// const RegularData1D& operator << (RegularData1D& to, const TFFT1D& from)
//; ???????
template <typename ComplexTraits>
TFFT1D<ComplexTraits>::TFFT1D()
: TRegularData1D<Complex>(0, 0, 1.), // AR: old case: This is necessary because FFTW_COMPLEX has no default constructor
length_(0),
inFourierSpace_(false),
dataAdress_(0),
planCalculated_(false)
{
}
template <typename ComplexTraits>
bool TFFT1D<ComplexTraits>::operator == (const TFFT1D<ComplexTraits>& fft1d) const
{
if (ComplexVector::size() == fft1d.size() &&
origin_ == fft1d.origin_ &&
stepPhys_ == fft1d.stepPhys_ &&
stepFourier_ == fft1d.stepFourier_ &&
minPhys_ == fft1d.minPhys_ &&
maxPhys_ == fft1d.maxPhys_ &&
minFourier_ == fft1d.minFourier_ &&
maxFourier_ == fft1d.maxFourier_ &&
inFourierSpace_ == fft1d.inFourierSpace_ &&
numPhysToFourier_ == fft1d.numPhysToFourier_ &&
numFourierToPhys_ == fft1d.numFourierToPhys_ &&
planCalculated_ == fft1d.planCalculated_)
{
double min = inFourierSpace_ ? minFourier_ : minPhys_;
double max = inFourierSpace_ ? maxFourier_ : maxPhys_;
double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
for (double pos=min; pos<=max; pos+=step)
{
if (getData(pos) != fft1d.getData(pos))
{
return false;
}
}
return true;
}
return false;
}
template <typename ComplexTraits>
bool TFFT1D<ComplexTraits>::translate(double trans_origin)
{
Position internalOrigin = (int) rint(trans_origin/stepPhys_);
if (internalOrigin <= length_)
{
origin_ = trans_origin;
minPhys_ = ((-1.)*origin_);
maxPhys_ = (((length_-1)*stepPhys_)-origin_);
minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
maxFourier_ = ((length_/2.)*stepFourier_);
return true;
}
else
{
return false;
}
}
template <typename ComplexTraits>
bool TFFT1D<ComplexTraits>::setPhysStepWidth(double new_width)
{
if (new_width < 0)
{
return false;
}
else
{
stepPhys_ = new_width;
stepFourier_ = 2.*Constants::PI/(stepPhys_*length_);
minPhys_ = ((-1.)*origin_);
maxPhys_ = (((length_-1)*stepPhys_)-origin_);
minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
maxFourier_ = ((length_/2.)*stepFourier_);
return true;
}
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getPhysStepWidth() const
{
return stepPhys_;
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getFourierStepWidth() const
{
return stepFourier_;
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getPhysSpaceMin() const
{
return minPhys_;
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getPhysSpaceMax() const
{
return maxPhys_;
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getFourierSpaceMin() const
{
return minFourier_;
}
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getFourierSpaceMax() const
{
return maxFourier_;
}
template <typename ComplexTraits>
Size TFFT1D<ComplexTraits>::getMaxIndex() const
{
return (length_ - 1);
}
template <typename ComplexTraits>
Size TFFT1D<ComplexTraits>::getNumberOfInverseTransforms() const
{
return numFourierToPhys_;
}
// AR: new
template <typename ComplexTraits>
double TFFT1D<ComplexTraits>::getGridCoordinates(Position position) const
{
if (!inFourierSpace_)
{
if (position >= ComplexVector::size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
double r;
r = -origin_ + (float)position * stepPhys_;
return r;
}
else
{
if (position >= ComplexVector::size())
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
double r;
Index x;
x = position;
if (x>=length_/2.)
{
x-=length_;
}
r = (float)x * stepFourier_;
return r;
}
}
template <typename ComplexTraits>
typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::getData(const double pos) const
{
Complex result;
double normalization=1.;
if (!inFourierSpace_)
{
result = (*this)[pos];//Complex((*this)[pos].real(), (*this)[pos].imag());
normalization=1./pow((double)length_,(int)numFourierToPhys_);
}
else
{
result = (*this)[pos]*phase(pos);//Complex((*this)[pos].real(),(*this)[pos].imag()) * phase(pos);
//Log.error() << pos << " " << phase(pos).real() << std::endl;
normalization=1./(sqrt(2.*Constants::PI))/pow((double)length_,(int)numFourierToPhys_);
}
result *= normalization;
return result;
}
template <typename ComplexTraits>
typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::getInterpolatedValue(const double pos) const
{
Complex result;
double min = inFourierSpace_ ? minFourier_ : minPhys_;
double max = inFourierSpace_ ? maxFourier_ : maxPhys_;
double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
if ((pos < min) || (pos > max))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
double h = pos - min;
double mod = fmod(h,step);
if (mod ==0) // we are on the grid
{
return getData(pos);
}
double before = floor(h/step)*step+ min;
double after = ceil(h/step)*step+ min;
double t = (pos - before)/step;
result = getData(before)*(typename ComplexTraits::ComplexPrecision)(1.-t);
result += getData(after)*(typename ComplexTraits::ComplexPrecision)t;
return result;
}
template <typename ComplexTraits>
void TFFT1D<ComplexTraits>::setData(double pos, Complex val)
{
Complex dummy;
if (!inFourierSpace_)
{
dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)),
val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)));
(*this)[pos]=dummy;
}
else
{
val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*Constants::PI)/stepPhys_))
*(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)length_,(int)numFourierToPhys_);
dummy = val;
(*this)[pos]=dummy;
}
}
template <typename ComplexTraits>
typename TFFT1D<ComplexTraits>::Complex& TFFT1D<ComplexTraits>::operator[](const double pos)
{
Index internalPos;
if (!inFourierSpace_)
{
internalPos = (Index) rint((pos+origin_)/stepPhys_);
}
else
{
internalPos = (Index) rint(pos/stepFourier_);
if (internalPos < 0)
{
internalPos+=length_;
}
}
if ((internalPos < 0) || (internalPos>=(Index) length_))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return operator [] ((Position)internalPos);
}
template <typename ComplexTraits>
const typename TFFT1D<ComplexTraits>::Complex& TFFT1D<ComplexTraits>::operator[](const double pos) const
{
Index internalPos;
if (!inFourierSpace_)
{
internalPos = (Index) rint((pos+origin_)/stepPhys_);
}
else
{
internalPos = (Index) rint(pos/stepFourier_);
if (internalPos < 0)
{
internalPos+=length_;
}
}
if ((internalPos < 0) || (internalPos>=(Index) length_))
{
throw Exception::OutOfGrid(__FILE__, __LINE__);
}
return operator [] ((Position)internalPos);
}
template <typename ComplexTraits>
typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::phase(const double pos) const
{
double phase = 2.*Constants::PI*(rint(pos/stepFourier_))
*(rint(origin_/stepPhys_))
/length_;
Complex result = Complex(cos(phase), sin(phase));
return result;
}
template <typename ComplexTraits>
bool TFFT1D<ComplexTraits>::isInFourierSpace() const
{
return inFourierSpace_;
}
/*
const TRegularData1D<Complex >& operator << (TRegularData1D<Complex >& to, const TFFT1D<ComplexTraits>& from)
{
// first decide if the TFFT1D data is in Fourier space.
if (!from.isInFourierSpace())
{
// create a new grid
Size lengthX = from.getMaxIndex()+1;
TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX), from.getPhysSpaceMin(), from.getPhysSpaceMax());
// and fill it
double normalization=1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
for (Position i = 0; i < from.size(); i++)
{
newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization;
}
to = newGrid;
return to;
}
else
{
// we are in Fourier space
// create a new grid
Size lengthX = from.getMaxIndex()+1;
//float stepPhysX = from.getPhysStepWidthX();
float stepFourierX = from.getFourierStepWidth();
TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX),
from.getFourierSpaceMin(),
from.getFourierSpaceMax());
// and fill it
// AR: old double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
double normalization=1./sqrt(2.*Constants::PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
Index x;
float r;
for (Position i = 0; i < from.size(); i++)
{
x = i;
if (x>lengthX/2.)
{
x-=lengthX;
}
r = (float)x * stepFourierX;
newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization*from.phase(r);
}
to = newGrid;
return to;
}
}
*/
template <typename ComplexTraits>
const RegularData1D& operator << (RegularData1D& to, const TFFT1D<ComplexTraits>& from)
{
// first decide if the FFT1D data is in Fourier space.
if (!from.isInFourierSpace())
{
// create a new grid
Size lengthX = from.getMaxIndex()+1;
RegularData1D newGrid(lengthX);
newGrid.setOrigin(from.getPhysSpaceMin());
newGrid.setDimension(from.getPhysSpaceMax()-from.getPhysSpaceMin());
// and fill it
double normalization = 1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
for (Position i = 0; i < from.size(); i++)
{
newGrid[i] = from[i].real()*normalization;
}
to = newGrid;
return to;
}
else
{
// we are in Fourier space
// create a new grid
Size lengthX = from.getMaxIndex()+1;
//float stepPhysX = from.getPhysStepWidth();
float stepFourierX = from.getFourierStepWidth();
RegularData1D newGrid(lengthX);
newGrid.setOrigin(from.getFourierSpaceMin());
newGrid.setDimension(from.getFourierSpaceMax()-from.getFourierSpaceMin());
// and fill it
// AR: old version double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
double normalization=1./sqrt(2.*Constants::PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
Index x;
signed int xp;
float r;
for (Position i = 0; i < from.size(); i++)
{
x = i;
xp = x;
if (xp>=lengthX/2.)
{
xp-=(int)lengthX;
}
if (x>=lengthX/2.)
{
x-=(int)(lengthX/2.);
}
else
{
x+=(int)(lengthX/2.);
}
r = ((float)xp * stepFourierX);
newGrid[i] = (from[i]*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
}
to = newGrid;
return to;
}
}
}
#endif // BALL_MATHS_TFFT1D_H
exception.h
BALL::Exception::OutOfGrid
Definition COMMON/exception.h:351
BALL::TFFT1D
Definition FFT1D.h:39
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::translate
bool translate(double trans_origin)
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::length_
Size length_
Definition FFT1D.h:244
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMax
double getFourierSpaceMax() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getGridCoordinates
double getGridCoordinates(Position position) const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::numFourierToPhys_
Size numFourierToPhys_
Definition FFT1D.h:247
BALL::TFFT1D::setNumberOfFFTTransforms
void setNumberOfFFTTransforms(Size num)
Definition FFT1D.h:220
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::planForward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planForward_
Definition FFT1D.h:256
BALL::TFFT1D::Complex
std::complex< typename ComplexTraits::ComplexPrecision > Complex
Definition FFT1D.h:42
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::clear
virtual void clear()
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::setData
void setData(double pos, Complex val)
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::planBackward_
BALL_FFTW_DEFAULT_TRAITS::FftwPlan planBackward_
Definition FFT1D.h:257
BALL::TFFT1D::operator[]
const Complex & operator[](const Position &pos) const
Definition FFT1D.h:215
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::inFourierSpace_
bool inFourierSpace_
Definition FFT1D.h:245
BALL::TFFT1D::isInFourierSpace
bool isInFourierSpace() const
Definition FFT1D.h:603
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::doFFT
void doFFT()
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::origin_
double origin_
Definition FFT1D.h:248
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getInterpolatedValue
Complex getInterpolatedValue(const double pos) const
BALL::TFFT1D::phase
Complex phase(const double pos) const
Definition FFT1D.h:592
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getPhysStepWidth
double getPhysStepWidth() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getNumberOfInverseTransforms
Size getNumberOfInverseTransforms() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::stepPhys_
double stepPhys_
Definition FFT1D.h:249
BALL::TFFT1D::ComplexVector
TRegularData1D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
Definition FFT1D.h:43
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getData
Complex getData(const double pos) const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::setPhysStepWidth
bool setPhysStepWidth(double new_width)
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::dataAdress_
Complex * dataAdress_
Definition FFT1D.h:260
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::minPhys_
double minPhys_
Definition FFT1D.h:251
BALL::TFFT1D::setNumberOfiFFTTransforms
void setNumberOfiFFTTransforms(Size num)
Definition FFT1D.h:225
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::maxFourier_
double maxFourier_
Definition FFT1D.h:254
BALL::TFFT1D::TFFT1D
TFFT1D()
Default constructor.
Definition FFT1D.h:285
BALL::TFFT1D::operator==
bool operator==(const TFFT1D &fft1d) const
Definition FFT1D.h:296
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::stepFourier_
double stepFourier_
Definition FFT1D.h:250
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getMaxIndex
Size getMaxIndex() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getFourierSpaceMin
double getFourierSpaceMin() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::minFourier_
double minFourier_
Definition FFT1D.h:253
BALL::TFFT1D::operator[]
Complex & operator[](const double pos)
Definition FFT1D.h:538
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::planCalculated_
bool planCalculated_
Definition FFT1D.h:261
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMin
double getPhysSpaceMin() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getFourierStepWidth
double getFourierStepWidth() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::getPhysSpaceMax
double getPhysSpaceMax() const
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::maxPhys_
double maxPhys_
Definition FFT1D.h:252
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::doiFFT
void doiFFT()
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::destroy
virtual void destroy()
BALL::TFFT1D< BALL_FFTW_DEFAULT_TRAITS >::numPhysToFourier_
Size numPhysToFourier_
Definition FFT1D.h:246
BALL::TRegularData1D< std::complex< typename ComplexTraits::ComplexPrecision > >::size
BALL_INLINE size_type size() const
Definition regularData1D.h:161
BALL::TRegularData1D::setDimension
BALL_INLINE void setDimension(const CoordinateType &dimension)
Definition regularData1D.h:290
BALL::TRegularData1D::setOrigin
BALL_INLINE void setOrigin(const CoordinateType &origin)
Definition regularData1D.h:276
BALL::TRegularData1D< std::complex< ComplexTraits::ComplexPrecision > >::TRegularData1D
TRegularData1D()
Definition regularData1D.h:361
BALL::TRegularData1D::operator[]
const ValueType & operator[](const IndexType &index) const
Definition regularData1D.h:181
constants.h
BALL_CREATE
#define BALL_CREATE(name)
Definition create.h:62
fftwCommon.h
BALL::operator<<
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
BALL::Constants::PI
BALL_EXTERN_VARIABLE const double PI
PI.
Definition constants.h:35
BALL
Definition constants.h:13
BALL::RegularData1D
TRegularData1D< float > RegularData1D
Definition regularData1D.h:358
BALL::Size
BALL_SIZE_TYPE Size
Definition COMMON/global.h:114
BALL::Position
BALL_SIZE_TYPE Position
Definition COMMON/global.h:140
BALL::FFT1D
TFFT1D< BALL_FFTW_DEFAULT_TRAITS > FFT1D
Definition FFT1D.h:267
BALL::Complex
std::complex< BALL_COMPLEX_PRECISION > Complex
Definition complex.h:21
BALL::Index
BALL_INDEX_TYPE Index
Definition COMMON/global.h:105
regularData1D.h