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/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.*M_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.*M_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*M_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.*M_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.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
      double normalization=1./sqrt(2.*M_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.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
      double normalization=1./sqrt(2.*M_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