BALL
1.4.1
|
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