BALL  1.4.1
FFT1D.h
Go to the documentation of this file.
00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 
00005 #ifndef BALL_MATHS_TFFT1D_H
00006 #define BALL_MATHS_TFFT1D_H
00007 
00008 #ifndef BALL_COMMON_EXCEPTION_H
00009 # include <BALL/COMMON/exception.h>
00010 #endif
00011 
00012 #ifndef BALL_DATATYPE_REGULARDATA1D_H
00013 # include <BALL/DATATYPE/regularData1D.h>
00014 #endif
00015 
00016 #include <math.h>
00017 #include <complex>
00018 #include <fftw3.h>
00019 
00020 #include <BALL/MATHS/fftwCommon.h>
00021 
00022 namespace BALL
00023 {
00026 
00027 
00035   template <typename ComplexTraits>
00036   class TFFT1D 
00037     : public TRegularData1D<std::complex<typename ComplexTraits::ComplexPrecision> >
00038   {
00039     public:
00040 
00041     typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
00042     typedef TRegularData1D<std::complex<typename ComplexTraits::ComplexPrecision> > ComplexVector;
00043 
00044     BALL_CREATE(TFFT1D)
00045 
00046     
00049 
00051     TFFT1D();
00052 
00054     TFFT1D(const TFFT1D &data);
00055 
00065      // AR: ldn is not any longer the binary logarithm but the absolute number of grid points
00066     TFFT1D(Size ldn, double stepPhys = 1., double origin = 0., bool inFourierSpace = false);
00067     
00069     virtual ~TFFT1D();
00070     
00072 
00076 
00078     const TFFT1D& operator = (const TFFT1D& fft1d);
00079     
00082     virtual void clear();
00083     
00086     virtual void destroy();
00087 
00089 
00093 
00096     bool operator == (const TFFT1D& fft1d) const;
00098     
00099     // @name Accessors
00100     
00103     void doFFT();
00104 
00107     void doiFFT();
00108 
00113     bool translate(double trans_origin);
00114 
00120     bool setPhysStepWidth(double new_width);
00121 
00124     double getPhysStepWidth() const;
00125 
00128     double getFourierStepWidth() const;
00129 
00132     double getPhysSpaceMin() const;
00133 
00136     double getPhysSpaceMax() const;
00137 
00140     double getFourierSpaceMin() const;
00141 
00144     double getFourierSpaceMax() const;
00145       
00151     Size getMaxIndex() const;
00152 
00156     Size getNumberOfInverseTransforms() const;
00157   
00160     double getGridCoordinates(Position position) const;
00161 
00168     Complex getData(const double pos) const;
00169 
00177     Complex getInterpolatedValue(const double pos) const;
00178 
00185     void setData(double pos, Complex val);
00186 
00192     Complex& operator [] (const double pos);
00193 
00199     const Complex& operator [] (const double pos) const;
00200     
00205     Complex& operator[](const Position& pos)
00206     {
00207       return TRegularData1D<Complex>::operator[](pos);
00208     }
00209       
00214     const Complex& operator[](const Position& pos) const
00215     {
00216       return TRegularData1D<Complex>::operator[](pos);
00217     } 
00218     
00219     void setNumberOfFFTTransforms(Size num)
00220     {
00221       numPhysToFourier_ = num;
00222     }
00223       
00224     void setNumberOfiFFTTransforms(Size num)
00225     {
00226       numFourierToPhys_ = num;
00227     }
00228     
00232     bool isInFourierSpace() const;
00233     
00239     Complex phase(const double pos) const;
00240 
00241     protected:
00242 
00243     Size length_;
00244     bool inFourierSpace_;
00245     Size numPhysToFourier_;
00246     Size numFourierToPhys_;
00247     double origin_;
00248     double stepPhys_;
00249     double stepFourier_;
00250     double minPhys_;
00251     double maxPhys_;
00252     double minFourier_;
00253     double maxFourier_;
00254 
00255     typename ComplexTraits::FftwPlan planForward_;
00256     typename ComplexTraits::FftwPlan planBackward_;
00257     
00258     // AR: to control plan calculation with new fftw3
00259     Complex *dataAdress_;
00260     bool planCalculated_;
00261   };
00263   
00266   typedef TFFT1D<BALL_FFTW_DEFAULT_TRAITS> FFT1D;
00267   
00268   // AR:
00271 //  const TRegularData1D<Complex>& operator << (TRegularData1D<Complex>& to, const TFFT1D& from)
00272   //; ?????
00273   
00278 //  const RegularData1D& operator << (RegularData1D& to, const TFFT1D& from)
00279 //; ???????
00280 
00281 
00282 
00283   template <typename ComplexTraits>
00284   TFFT1D<ComplexTraits>::TFFT1D()
00285     : TRegularData1D<Complex>(0, 0, 1.),  // AR: old case: This is necessary because FFTW_COMPLEX has no default constructor
00286       length_(0),
00287       inFourierSpace_(false),
00288       dataAdress_(0),
00289       planCalculated_(false)
00290   {
00291   }
00292     
00293 
00294   template <typename ComplexTraits>
00295   bool TFFT1D<ComplexTraits>::operator == (const TFFT1D<ComplexTraits>& fft1d) const
00296   {
00297     if (ComplexVector::size() == fft1d.size() &&
00298         origin_ == fft1d.origin_ &&
00299         stepPhys_ == fft1d.stepPhys_ &&
00300         stepFourier_ == fft1d.stepFourier_ &&
00301         minPhys_ == fft1d.minPhys_ &&
00302         maxPhys_ == fft1d.maxPhys_ &&
00303         minFourier_ == fft1d.minFourier_ &&
00304         maxFourier_ == fft1d.maxFourier_ &&
00305         inFourierSpace_ == fft1d.inFourierSpace_ &&
00306         numPhysToFourier_ == fft1d.numPhysToFourier_ &&
00307         numFourierToPhys_ == fft1d.numFourierToPhys_ &&
00308         planCalculated_ == fft1d.planCalculated_)
00309     {
00310       double min  = inFourierSpace_ ?  minFourier_ :  minPhys_;
00311       double max  = inFourierSpace_ ?  maxFourier_ :  maxPhys_;
00312       double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
00313         
00314       for (double pos=min; pos<=max; pos+=step)
00315       {
00316         if (getData(pos) != fft1d.getData(pos))
00317         {
00318           return false;
00319         }
00320       }
00321       
00322       return true;
00323     }
00324   
00325     return false;
00326   } 
00327 
00328   template <typename ComplexTraits>
00329   bool TFFT1D<ComplexTraits>::translate(double trans_origin)
00330   {
00331     Position internalOrigin = (int) rint(trans_origin/stepPhys_);
00332     if (internalOrigin <= length_)
00333     {
00334       origin_ = trans_origin;
00335 
00336       minPhys_ = ((-1.)*origin_);
00337       maxPhys_ = (((length_-1)*stepPhys_)-origin_);
00338       minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
00339       maxFourier_ = ((length_/2.)*stepFourier_);
00340       
00341       return true;
00342     }
00343     else
00344     {
00345       return false;
00346     }
00347   }
00348 
00349   template <typename ComplexTraits>
00350   bool TFFT1D<ComplexTraits>::setPhysStepWidth(double new_width)
00351   {
00352     if (new_width < 0)
00353     {
00354       return false;
00355     }
00356     else
00357     {
00358       stepPhys_ = new_width;
00359       stepFourier_ = 2.*M_PI/(stepPhys_*length_);
00360 
00361       minPhys_ = ((-1.)*origin_);
00362       maxPhys_ = (((length_-1)*stepPhys_)-origin_);
00363       minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
00364       maxFourier_ = ((length_/2.)*stepFourier_);
00365 
00366       return true;
00367     }
00368   }
00369 
00370   template <typename ComplexTraits> 
00371   double TFFT1D<ComplexTraits>::getPhysStepWidth() const
00372   {
00373     return stepPhys_;
00374   }
00375 
00376   template <typename ComplexTraits>
00377   double TFFT1D<ComplexTraits>::getFourierStepWidth() const
00378   {
00379     return stepFourier_;
00380   }
00381 
00382   template <typename ComplexTraits>
00383   double TFFT1D<ComplexTraits>::getPhysSpaceMin() const
00384   {
00385     return minPhys_;
00386   }
00387 
00388   template <typename ComplexTraits>
00389   double TFFT1D<ComplexTraits>::getPhysSpaceMax() const
00390   {
00391     return maxPhys_;
00392   }
00393 
00394   template <typename ComplexTraits>
00395   double TFFT1D<ComplexTraits>::getFourierSpaceMin() const
00396   {
00397     return minFourier_;
00398   }
00399 
00400   template <typename ComplexTraits>
00401   double TFFT1D<ComplexTraits>::getFourierSpaceMax() const
00402   {
00403     return maxFourier_;
00404   }
00405 
00406   template <typename ComplexTraits> 
00407   Size TFFT1D<ComplexTraits>::getMaxIndex() const
00408   {
00409     return (length_ - 1);
00410   }
00411 
00412   template <typename ComplexTraits> 
00413   Size TFFT1D<ComplexTraits>::getNumberOfInverseTransforms() const
00414   {
00415     return numFourierToPhys_;
00416   }
00417   
00418   // AR: new 
00419   template <typename ComplexTraits>
00420   double TFFT1D<ComplexTraits>::getGridCoordinates(Position position) const
00421   {
00422     if (!inFourierSpace_)
00423     {
00424       if (position >= ComplexVector::size())
00425       {
00426         throw Exception::OutOfGrid(__FILE__, __LINE__);
00427       }
00428     
00429       double r;
00430       
00431       r = -origin_ + (float)position * stepPhys_;
00432 
00433       return r;
00434     }
00435     else
00436     {
00437       if (position >= ComplexVector::size())
00438       {
00439         throw Exception::OutOfGrid(__FILE__, __LINE__);
00440       }
00441     
00442       double r;
00443       Index x;
00444   
00445       x = position;
00446 
00447       if (x>=length_/2.)
00448       {
00449         x-=length_;
00450       }
00451       
00452       r = (float)x * stepFourier_;
00453 
00454       return r;
00455     }
00456   }
00457 
00458   template <typename ComplexTraits>
00459   typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::getData(const double pos) const
00460   {
00461     Complex result;
00462     double normalization=1.;
00463 
00464     if (!inFourierSpace_)
00465     {
00466       result = (*this)[pos];//Complex((*this)[pos].real(), (*this)[pos].imag());
00467       normalization=1./pow((double)length_,(int)numFourierToPhys_);
00468     }
00469     else
00470     {
00471       result = (*this)[pos]*phase(pos);//Complex((*this)[pos].real(),(*this)[pos].imag()) * phase(pos);
00472       //Log.error() << pos <<  " " << phase(pos).real() <<  std::endl;
00473       normalization=1./(sqrt(2.*M_PI))/pow((double)length_,(int)numFourierToPhys_);
00474     }
00475 
00476     result *= normalization;
00477     
00478     return result;
00479   }
00480 
00481   template <typename ComplexTraits>
00482   typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::getInterpolatedValue(const double pos) const
00483   {
00484     Complex result;
00485     
00486     double min  = inFourierSpace_ ? minFourier_  :  minPhys_;
00487     double max  = inFourierSpace_ ? maxFourier_  :  maxPhys_;
00488     double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
00489     
00490     if ((pos < min) || (pos > max))
00491     {
00492       throw Exception::OutOfGrid(__FILE__, __LINE__);
00493     }
00494 
00495     double h = pos - min;
00496     double mod = fmod(h,step);
00497 
00498     if (mod ==0) // we are on the grid
00499     {
00500       return getData(pos);
00501     }
00502 
00503     double before = floor(h/step)*step+ min;
00504     double after  = ceil(h/step)*step+ min;
00505       
00506     double t = (pos - before)/step;
00507 
00508     result  =  getData(before)*(typename ComplexTraits::ComplexPrecision)(1.-t);
00509     result +=  getData(after)*(typename ComplexTraits::ComplexPrecision)t; 
00510 
00511     return result;
00512   }
00513 
00514   template <typename ComplexTraits>
00515   void TFFT1D<ComplexTraits>::setData(double pos, Complex val)
00516   {
00517     Complex dummy;
00518     if (!inFourierSpace_)
00519     {
00520       dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)),
00521                        val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)));
00522   
00523       (*this)[pos]=dummy;
00524     }
00525     else
00526     {
00527       val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*M_PI)/stepPhys_))
00528                      *(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)length_,(int)numFourierToPhys_);
00529       
00530       dummy = val;
00531       
00532       (*this)[pos]=dummy;
00533     }
00534   }
00535 
00536   template <typename ComplexTraits>
00537   typename TFFT1D<ComplexTraits>::Complex& TFFT1D<ComplexTraits>::operator[](const double pos) 
00538   {
00539     Index internalPos;
00540 
00541     if (!inFourierSpace_)
00542     {
00543       internalPos = (Index) rint((pos+origin_)/stepPhys_);
00544     }
00545     else
00546     {
00547       internalPos =  (Index) rint(pos/stepFourier_);
00548 
00549       if (internalPos < 0)
00550       {
00551         internalPos+=length_;
00552       }
00553     }
00554 
00555     if ((internalPos < 0) || (internalPos>=(Index) length_))
00556     {
00557       throw Exception::OutOfGrid(__FILE__, __LINE__);
00558     }
00559     
00560     return operator [] ((Position)internalPos);
00561   }
00562 
00563   template <typename ComplexTraits>
00564   const typename TFFT1D<ComplexTraits>::Complex& TFFT1D<ComplexTraits>::operator[](const double pos) const
00565   {
00566     Index internalPos;
00567 
00568     if (!inFourierSpace_)
00569     {
00570       internalPos = (Index) rint((pos+origin_)/stepPhys_);
00571     }
00572     else
00573     {
00574       internalPos =  (Index) rint(pos/stepFourier_);
00575 
00576       if (internalPos < 0)
00577       {
00578         internalPos+=length_;
00579       }
00580     }
00581 
00582     if ((internalPos < 0) || (internalPos>=(Index) length_))
00583     {
00584       throw Exception::OutOfGrid(__FILE__, __LINE__);
00585     }
00586     
00587     return operator [] ((Position)internalPos);
00588   }
00589 
00590   template <typename ComplexTraits>
00591   typename TFFT1D<ComplexTraits>::Complex TFFT1D<ComplexTraits>::phase(const double pos) const
00592   {
00593     double phase = 2.*M_PI*(rint(pos/stepFourier_))
00594                          *(rint(origin_/stepPhys_))
00595                          /length_;
00596     Complex result = Complex(cos(phase), sin(phase));
00597             
00598     return result;
00599   }
00600 
00601   template <typename ComplexTraits>
00602   bool TFFT1D<ComplexTraits>::isInFourierSpace() const
00603   {
00604     return inFourierSpace_;
00605   }
00606 /*  
00607   const TRegularData1D<Complex >& operator << (TRegularData1D<Complex >& to, const TFFT1D<ComplexTraits>& from)
00608   {
00609     // first decide if the TFFT1D data is in Fourier space.
00610     if (!from.isInFourierSpace())
00611     {
00612       // create a new grid
00613       Size lengthX = from.getMaxIndex()+1;
00614       
00615       TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX), from.getPhysSpaceMin(), from.getPhysSpaceMax());
00616 
00617       // and fill it
00618       double normalization=1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
00619       
00620       for (Position i = 0; i < from.size(); i++)
00621       {
00622         newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization;
00623       }
00624 
00625       to = newGrid;
00626 
00627       return to;
00628     }
00629     else
00630     {
00631       // we are in Fourier space
00632       
00633       // create a new grid
00634       Size lengthX = from.getMaxIndex()+1;
00635       //float stepPhysX = from.getPhysStepWidthX();
00636       float stepFourierX = from.getFourierStepWidth();
00637 
00638       
00639       TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX),
00640                                       from.getFourierSpaceMin(), 
00641                                       from.getFourierSpaceMax());
00642 
00643       // and fill it
00644       // AR: old double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
00645       double normalization=1./sqrt(2.*M_PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
00646       
00647       
00648       Index x;
00649       float r;
00650   
00651       for (Position i = 0; i < from.size(); i++)
00652       {
00653         x = i;
00654 
00655         if (x>lengthX/2.)
00656         {
00657           x-=lengthX;
00658         }
00659 
00660         r = (float)x * stepFourierX;
00661 
00662         newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization*from.phase(r);
00663       }
00664 
00665       to = newGrid;
00666 
00667       return to;
00668     }
00669   }
00670   */
00671   
00672   template <typename ComplexTraits>
00673   const RegularData1D& operator << (RegularData1D& to, const TFFT1D<ComplexTraits>& from)
00674   {
00675     // first decide if the FFT1D data is in Fourier space.
00676     if (!from.isInFourierSpace())
00677     {
00678       // create a new grid
00679       Size lengthX = from.getMaxIndex()+1;
00680       
00681       RegularData1D newGrid(lengthX);
00682       newGrid.setOrigin(from.getPhysSpaceMin());
00683       newGrid.setDimension(from.getPhysSpaceMax()-from.getPhysSpaceMin());
00684 
00685       // and fill it
00686       double normalization = 1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
00687       
00688       for (Position i = 0; i < from.size(); i++)
00689       {
00690         newGrid[i] = from[i].real()*normalization;
00691       }
00692 
00693       to = newGrid;
00694 
00695       return to;
00696     }
00697     else
00698     {
00699       // we are in Fourier space
00700       
00701       // create a new grid
00702       Size lengthX = from.getMaxIndex()+1;
00703       //float stepPhysX = from.getPhysStepWidth();
00704       float stepFourierX = from.getFourierStepWidth();
00705 
00706       RegularData1D newGrid(lengthX);
00707       newGrid.setOrigin(from.getFourierSpaceMin());
00708       newGrid.setDimension(from.getFourierSpaceMax()-from.getFourierSpaceMin());
00709 
00710       // and fill it
00711       // AR: old version double normalization=1./(sqrt(2.*M_PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
00712       double normalization=1./sqrt(2.*M_PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
00713       
00714       Index x;
00715       signed int xp;
00716       float r;
00717   
00718       for (Position i = 0; i < from.size(); i++)
00719       {
00720         x =  i;
00721         
00722         xp = x;
00723 
00724         if (xp>=lengthX/2.)
00725         {
00726           xp-=(int)lengthX;
00727         }
00728 
00729         if (x>=lengthX/2.)
00730         {
00731           x-=(int)(lengthX/2.);
00732         }
00733         else
00734         {
00735           x+=(int)(lengthX/2.);
00736         }
00737 
00738 
00739         r = ((float)xp * stepFourierX);
00740 
00741         newGrid[i] = (from[i]*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
00742       }
00743 
00744       to = newGrid;
00745 
00746       return to;
00747     }
00748   }
00749 }
00750 #endif // BALL_MATHS_TFFT1D_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines