BALL  1.5.0
FFT1D.h
Go to the documentation of this file.
1 // -*- Mode: C++; tab-width: 2; -*-
2 // vi: set ts=2:
3 //
4 
5 #ifndef BALL_MATHS_TFFT1D_H
6 #define BALL_MATHS_TFFT1D_H
7 
8 #ifndef BALL_COMMON_EXCEPTION_H
9 # include <BALL/COMMON/exception.h>
10 #endif
11 
12 #ifndef BALL_DATATYPE_REGULARDATA1D_H
14 #endif
15 
16 #include <cmath>
17 #include <complex>
18 #include <fftw3.h>
19 
20 #include <BALL/COMMON/constants.h>
21 #include <BALL/MATHS/fftwCommon.h>
22 
23 namespace BALL
24 {
27 
28 
36  template <typename ComplexTraits>
37  class TFFT1D
38  : public TRegularData1D<std::complex<typename ComplexTraits::ComplexPrecision> >
39  {
40  public:
41 
42  typedef std::complex<typename ComplexTraits::ComplexPrecision> Complex;
44 
46 
47 
50 
52  TFFT1D();
53 
55  TFFT1D(const TFFT1D &data);
56 
66  // AR: ldn is not any longer the binary logarithm but the absolute number of grid points
67  TFFT1D(Size ldn, double stepPhys = 1., double origin = 0., bool inFourierSpace = false);
68 
70  virtual ~TFFT1D();
71 
73 
77 
79  const TFFT1D& operator = (const TFFT1D& fft1d);
80 
83  virtual void clear();
84 
87  virtual void destroy();
88 
90 
94 
97  bool operator == (const TFFT1D& fft1d) const;
99 
100  // @name Accessors
101 
104  void doFFT();
105 
108  void doiFFT();
109 
114  bool translate(double trans_origin);
115 
121  bool setPhysStepWidth(double new_width);
122 
125  double getPhysStepWidth() const;
126 
129  double getFourierStepWidth() const;
130 
133  double getPhysSpaceMin() const;
134 
137  double getPhysSpaceMax() const;
138 
141  double getFourierSpaceMin() const;
142 
145  double getFourierSpaceMax() const;
146 
152  Size getMaxIndex() const;
153 
158 
161  double getGridCoordinates(Position position) const;
162 
169  Complex getData(const double pos) const;
170 
178  Complex getInterpolatedValue(const double pos) const;
179 
186  void setData(double pos, Complex val);
187 
193  Complex& operator [] (const double pos);
194 
200  const Complex& operator [] (const double pos) const;
201 
207  {
209  }
210 
215  const Complex& operator[](const Position& pos) const
216  {
218  }
219 
221  {
222  numPhysToFourier_ = num;
223  }
224 
226  {
227  numFourierToPhys_ = num;
228  }
229 
233  bool isInFourierSpace() const;
234 
240  Complex phase(const double pos) const;
241 
242  protected:
243 
248  double origin_;
249  double stepPhys_;
250  double stepFourier_;
251  double minPhys_;
252  double maxPhys_;
253  double minFourier_;
254  double maxFourier_;
255 
256  typename ComplexTraits::FftwPlan planForward_;
257  typename ComplexTraits::FftwPlan planBackward_;
258 
259  // AR: to control plan calculation with new fftw3
262  };
264 
268 
269  // AR:
272 // const TRegularData1D<Complex>& operator << (TRegularData1D<Complex>& to, const TFFT1D& from)
273  //; ?????
274 
279 // const RegularData1D& operator << (RegularData1D& to, const TFFT1D& from)
280 //; ???????
281 
282 
283 
284  template <typename ComplexTraits>
286  : TRegularData1D<Complex>(0, 0, 1.), // AR: old case: This is necessary because FFTW_COMPLEX has no default constructor
287  length_(0),
288  inFourierSpace_(false),
289  dataAdress_(0),
290  planCalculated_(false)
291  {
292  }
293 
294 
295  template <typename ComplexTraits>
297  {
298  if (ComplexVector::size() == fft1d.size() &&
299  origin_ == fft1d.origin_ &&
300  stepPhys_ == fft1d.stepPhys_ &&
301  stepFourier_ == fft1d.stepFourier_ &&
302  minPhys_ == fft1d.minPhys_ &&
303  maxPhys_ == fft1d.maxPhys_ &&
304  minFourier_ == fft1d.minFourier_ &&
305  maxFourier_ == fft1d.maxFourier_ &&
306  inFourierSpace_ == fft1d.inFourierSpace_ &&
307  numPhysToFourier_ == fft1d.numPhysToFourier_ &&
308  numFourierToPhys_ == fft1d.numFourierToPhys_ &&
309  planCalculated_ == fft1d.planCalculated_)
310  {
311  double min = inFourierSpace_ ? minFourier_ : minPhys_;
312  double max = inFourierSpace_ ? maxFourier_ : maxPhys_;
313  double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
314 
315  for (double pos=min; pos<=max; pos+=step)
316  {
317  if (getData(pos) != fft1d.getData(pos))
318  {
319  return false;
320  }
321  }
322 
323  return true;
324  }
325 
326  return false;
327  }
328 
329  template <typename ComplexTraits>
330  bool TFFT1D<ComplexTraits>::translate(double trans_origin)
331  {
332  Position internalOrigin = (int) rint(trans_origin/stepPhys_);
333  if (internalOrigin <= length_)
334  {
335  origin_ = trans_origin;
336 
337  minPhys_ = ((-1.)*origin_);
338  maxPhys_ = (((length_-1)*stepPhys_)-origin_);
339  minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
340  maxFourier_ = ((length_/2.)*stepFourier_);
341 
342  return true;
343  }
344  else
345  {
346  return false;
347  }
348  }
349 
350  template <typename ComplexTraits>
352  {
353  if (new_width < 0)
354  {
355  return false;
356  }
357  else
358  {
359  stepPhys_ = new_width;
360  stepFourier_ = 2.*Constants::PI/(stepPhys_*length_);
361 
362  minPhys_ = ((-1.)*origin_);
363  maxPhys_ = (((length_-1)*stepPhys_)-origin_);
364  minFourier_ = ((-1.)*(length_/2.-1)*stepFourier_);
365  maxFourier_ = ((length_/2.)*stepFourier_);
366 
367  return true;
368  }
369  }
370 
371  template <typename ComplexTraits>
373  {
374  return stepPhys_;
375  }
376 
377  template <typename ComplexTraits>
379  {
380  return stepFourier_;
381  }
382 
383  template <typename ComplexTraits>
385  {
386  return minPhys_;
387  }
388 
389  template <typename ComplexTraits>
391  {
392  return maxPhys_;
393  }
394 
395  template <typename ComplexTraits>
397  {
398  return minFourier_;
399  }
400 
401  template <typename ComplexTraits>
403  {
404  return maxFourier_;
405  }
406 
407  template <typename ComplexTraits>
409  {
410  return (length_ - 1);
411  }
412 
413  template <typename ComplexTraits>
415  {
416  return numFourierToPhys_;
417  }
418 
419  // AR: new
420  template <typename ComplexTraits>
422  {
423  if (!inFourierSpace_)
424  {
425  if (position >= ComplexVector::size())
426  {
427  throw Exception::OutOfGrid(__FILE__, __LINE__);
428  }
429 
430  double r;
431 
432  r = -origin_ + (float)position * stepPhys_;
433 
434  return r;
435  }
436  else
437  {
438  if (position >= ComplexVector::size())
439  {
440  throw Exception::OutOfGrid(__FILE__, __LINE__);
441  }
442 
443  double r;
444  Index x;
445 
446  x = position;
447 
448  if (x>=length_/2.)
449  {
450  x-=length_;
451  }
452 
453  r = (float)x * stepFourier_;
454 
455  return r;
456  }
457  }
458 
459  template <typename ComplexTraits>
461  {
462  Complex result;
463  double normalization=1.;
464 
465  if (!inFourierSpace_)
466  {
467  result = (*this)[pos];//Complex((*this)[pos].real(), (*this)[pos].imag());
468  normalization=1./pow((double)length_,(int)numFourierToPhys_);
469  }
470  else
471  {
472  result = (*this)[pos]*phase(pos);//Complex((*this)[pos].real(),(*this)[pos].imag()) * phase(pos);
473  //Log.error() << pos << " " << phase(pos).real() << std::endl;
474  normalization=1./(sqrt(2.*Constants::PI))/pow((double)length_,(int)numFourierToPhys_);
475  }
476 
477  result *= normalization;
478 
479  return result;
480  }
481 
482  template <typename ComplexTraits>
484  {
485  Complex result;
486 
487  double min = inFourierSpace_ ? minFourier_ : minPhys_;
488  double max = inFourierSpace_ ? maxFourier_ : maxPhys_;
489  double step = inFourierSpace_ ? stepFourier_ : stepPhys_;
490 
491  if ((pos < min) || (pos > max))
492  {
493  throw Exception::OutOfGrid(__FILE__, __LINE__);
494  }
495 
496  double h = pos - min;
497  double mod = fmod(h,step);
498 
499  if (mod ==0) // we are on the grid
500  {
501  return getData(pos);
502  }
503 
504  double before = floor(h/step)*step+ min;
505  double after = ceil(h/step)*step+ min;
506 
507  double t = (pos - before)/step;
508 
509  result = getData(before)*(typename ComplexTraits::ComplexPrecision)(1.-t);
510  result += getData(after)*(typename ComplexTraits::ComplexPrecision)t;
511 
512  return result;
513  }
514 
515  template <typename ComplexTraits>
517  {
518  Complex dummy;
519  if (!inFourierSpace_)
520  {
521  dummy = Complex(val.real()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)),
522  val.imag()*((typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)(length_),(int)numFourierToPhys_)));
523 
524  (*this)[pos]=dummy;
525  }
526  else
527  {
528  val*=phase(pos)*(typename ComplexTraits::ComplexPrecision)((sqrt(2*Constants::PI)/stepPhys_))
529  *(typename ComplexTraits::ComplexPrecision)pow((typename ComplexTraits::ComplexPrecision)length_,(int)numFourierToPhys_);
530 
531  dummy = val;
532 
533  (*this)[pos]=dummy;
534  }
535  }
536 
537  template <typename ComplexTraits>
539  {
540  Index internalPos;
541 
542  if (!inFourierSpace_)
543  {
544  internalPos = (Index) rint((pos+origin_)/stepPhys_);
545  }
546  else
547  {
548  internalPos = (Index) rint(pos/stepFourier_);
549 
550  if (internalPos < 0)
551  {
552  internalPos+=length_;
553  }
554  }
555 
556  if ((internalPos < 0) || (internalPos>=(Index) length_))
557  {
558  throw Exception::OutOfGrid(__FILE__, __LINE__);
559  }
560 
561  return operator [] ((Position)internalPos);
562  }
563 
564  template <typename ComplexTraits>
566  {
567  Index internalPos;
568 
569  if (!inFourierSpace_)
570  {
571  internalPos = (Index) rint((pos+origin_)/stepPhys_);
572  }
573  else
574  {
575  internalPos = (Index) rint(pos/stepFourier_);
576 
577  if (internalPos < 0)
578  {
579  internalPos+=length_;
580  }
581  }
582 
583  if ((internalPos < 0) || (internalPos>=(Index) length_))
584  {
585  throw Exception::OutOfGrid(__FILE__, __LINE__);
586  }
587 
588  return operator [] ((Position)internalPos);
589  }
590 
591  template <typename ComplexTraits>
593  {
594  double phase = 2.*Constants::PI*(rint(pos/stepFourier_))
595  *(rint(origin_/stepPhys_))
596  /length_;
597  Complex result = Complex(cos(phase), sin(phase));
598 
599  return result;
600  }
601 
602  template <typename ComplexTraits>
604  {
605  return inFourierSpace_;
606  }
607 /*
608  const TRegularData1D<Complex >& operator << (TRegularData1D<Complex >& to, const TFFT1D<ComplexTraits>& from)
609  {
610  // first decide if the TFFT1D data is in Fourier space.
611  if (!from.isInFourierSpace())
612  {
613  // create a new grid
614  Size lengthX = from.getMaxIndex()+1;
615 
616  TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX), from.getPhysSpaceMin(), from.getPhysSpaceMax());
617 
618  // and fill it
619  double normalization=1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
620 
621  for (Position i = 0; i < from.size(); i++)
622  {
623  newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization;
624  }
625 
626  to = newGrid;
627 
628  return to;
629  }
630  else
631  {
632  // we are in Fourier space
633 
634  // create a new grid
635  Size lengthX = from.getMaxIndex()+1;
636  //float stepPhysX = from.getPhysStepWidthX();
637  float stepFourierX = from.getFourierStepWidth();
638 
639 
640  TRegularData1D<Complex > newGrid(TRegularData1D<Complex >::IndexType(lengthX),
641  from.getFourierSpaceMin(),
642  from.getFourierSpaceMax());
643 
644  // and fill it
645  // AR: old double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
646  double normalization=1./sqrt(2.*Constants::PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
647 
648 
649  Index x;
650  float r;
651 
652  for (Position i = 0; i < from.size(); i++)
653  {
654  x = i;
655 
656  if (x>lengthX/2.)
657  {
658  x-=lengthX;
659  }
660 
661  r = (float)x * stepFourierX;
662 
663  newGrid[i] = from[i]*(ComplexTraits::ComplexPrecision)normalization*from.phase(r);
664  }
665 
666  to = newGrid;
667 
668  return to;
669  }
670  }
671  */
672 
673  template <typename ComplexTraits>
675  {
676  // first decide if the FFT1D data is in Fourier space.
677  if (!from.isInFourierSpace())
678  {
679  // create a new grid
680  Size lengthX = from.getMaxIndex()+1;
681 
682  RegularData1D newGrid(lengthX);
683  newGrid.setOrigin(from.getPhysSpaceMin());
684  newGrid.setDimension(from.getPhysSpaceMax()-from.getPhysSpaceMin());
685 
686  // and fill it
687  double normalization = 1./(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
688 
689  for (Position i = 0; i < from.size(); i++)
690  {
691  newGrid[i] = from[i].real()*normalization;
692  }
693 
694  to = newGrid;
695 
696  return to;
697  }
698  else
699  {
700  // we are in Fourier space
701 
702  // create a new grid
703  Size lengthX = from.getMaxIndex()+1;
704  //float stepPhysX = from.getPhysStepWidth();
705  float stepFourierX = from.getFourierStepWidth();
706 
707  RegularData1D newGrid(lengthX);
708  newGrid.setOrigin(from.getFourierSpaceMin());
709  newGrid.setDimension(from.getFourierSpaceMax()-from.getFourierSpaceMin());
710 
711  // and fill it
712  // AR: old version double normalization=1./(sqrt(2.*Constants::PI))*(stepPhysX*stepPhysY*stepPhysZ)/(pow((float)(lengthX*lengthY*lengthZ),from.getNumberOfInverseTransforms()));
713  double normalization=1./sqrt(2.*Constants::PI)/(pow((float)(lengthX),from.getNumberOfInverseTransforms()));
714 
715  Index x;
716  signed int xp;
717  float r;
718 
719  for (Position i = 0; i < from.size(); i++)
720  {
721  x = i;
722 
723  xp = x;
724 
725  if (xp>=lengthX/2.)
726  {
727  xp-=(int)lengthX;
728  }
729 
730  if (x>=lengthX/2.)
731  {
732  x-=(int)(lengthX/2.);
733  }
734  else
735  {
736  x+=(int)(lengthX/2.);
737  }
738 
739 
740  r = ((float)xp * stepFourierX);
741 
742  newGrid[i] = (from[i]*(typename ComplexTraits::ComplexPrecision)normalization*from.phase(r)).real();
743  }
744 
745  to = newGrid;
746 
747  return to;
748  }
749  }
750 }
751 #endif // BALL_MATHS_TFFT1D_H
BALL::TFFT1D::getFourierSpaceMin
double getFourierSpaceMin() const
Definition: FFT1D.h:396
BALL::TFFT1D::clear
virtual void clear()
BALL::TFFT1D::stepPhys_
double stepPhys_
Definition: FFT1D.h:249
BALL::TFFT1D::getPhysStepWidth
double getPhysStepWidth() const
Definition: FFT1D.h:372
BALL_INDEX_TYPE
BALL::TFFT1D::getGridCoordinates
double getGridCoordinates(Position position) const
Definition: FFT1D.h:421
BALL::TFFT1D::getInterpolatedValue
Complex getInterpolatedValue(const double pos) const
Definition: FFT1D.h:483
BALL::TFFT1D::stepFourier_
double stepFourier_
Definition: FFT1D.h:250
BALL::TFFT1D::setNumberOfiFFTTransforms
void setNumberOfiFFTTransforms(Size num)
Definition: FFT1D.h:225
BALL::Constants::PI
const BALL_EXTERN_VARIABLE double PI
PI.
Definition: constants.h:35
BALL::TFFT1D::setPhysStepWidth
bool setPhysStepWidth(double new_width)
Definition: FFT1D.h:351
BALL::TFFT1D::operator[]
Complex & operator[](const Position &pos)
Definition: FFT1D.h:206
BALL::TFFT1D::getMaxIndex
Size getMaxIndex() const
Definition: FFT1D.h:408
BALL::TFFT1D::doiFFT
void doiFFT()
BALL::TFFT1D::TFFT1D
TFFT1D()
Default constructor.
Definition: FFT1D.h:285
BALL::TFFT1D::setNumberOfFFTTransforms
void setNumberOfFFTTransforms(Size num)
Definition: FFT1D.h:220
BALL::Exception::OutOfGrid
Definition: COMMON/exception.h:349
BALL::TFFT1D::maxFourier_
double maxFourier_
Definition: FFT1D.h:254
BALL::TFFT1D::planBackward_
ComplexTraits::FftwPlan planBackward_
Definition: FFT1D.h:257
BALL::TFFT1D::getPhysSpaceMin
double getPhysSpaceMin() const
Definition: FFT1D.h:384
BALL::TFFT1D::origin_
double origin_
Definition: FFT1D.h:248
BALL::Maths::rint
double rint(double x)
round to integral value in floating-point format
Definition: MATHS/common.h:327
BALL::TFFT1D::numPhysToFourier_
Size numPhysToFourier_
Definition: FFT1D.h:246
BALL::TFFT1D::phase
Complex phase(const double pos) const
Definition: FFT1D.h:592
BALL::TRegularData1D< std::complex< ComplexTraits::ComplexPrecision > >::size
BALL_INLINE size_type size() const
Definition: regularData1D.h:161
BALL_SIZE_TYPE
BALL::Constants::h
const BALL_EXTERN_VARIABLE double h
Definition: constants.h:102
BALL::TRegularData1D
Definition: regularData1D.h:41
BALL::TFFT1D::numFourierToPhys_
Size numFourierToPhys_
Definition: FFT1D.h:247
BALL::TFFT1D::minFourier_
double minFourier_
Definition: FFT1D.h:253
BALL
Definition: constants.h:12
BALL::Maths::floor
long floor(const T &t)
Definition: MATHS/common.h:284
BALL::TFFT1D::dataAdress_
Complex * dataAdress_
Definition: FFT1D.h:260
fftwCommon.h
exception.h
BALL::Index
BALL_INDEX_TYPE Index
Definition: COMMON/global.h:105
BALL::TFFT1D::translate
bool translate(double trans_origin)
Definition: FFT1D.h:330
BALL::TFFT1D::inFourierSpace_
bool inFourierSpace_
Definition: FFT1D.h:245
BALL::FFT1D
TFFT1D< BALL_FFTW_DEFAULT_TRAITS > FFT1D
Definition: FFT1D.h:267
BALL::TFFT1D::ComplexVector
TRegularData1D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
Definition: FFT1D.h:43
BALL::TFFT1D::Complex
std::complex< typename ComplexTraits::ComplexPrecision > Complex
Definition: FFT1D.h:42
BALL::TRegularData1D::operator[]
const ValueType & operator[](const IndexType &index) const
Definition: regularData1D.h:181
BALL::TFFT1D::planCalculated_
bool planCalculated_
Definition: FFT1D.h:261
BALL::TFFT1D::isInFourierSpace
bool isInFourierSpace() const
Definition: FFT1D.h:603
BALL::TFFT1D::minPhys_
double minPhys_
Definition: FFT1D.h:251
constants.h
BALL::TFFT1D
Definition: FFT1D.h:37
BALL::TFFT1D::operator[]
const Complex & operator[](const Position &pos) const
Definition: FFT1D.h:215
BALL::Maths::min
T min(const T &a, const T &b)
Definition: MATHS/common.h:102
BALL::TFFT1D::length_
Size length_
Definition: FFT1D.h:244
BALL::TFFT1D::destroy
virtual void destroy()
BALL::TFFT1D::planForward_
ComplexTraits::FftwPlan planForward_
Definition: FFT1D.h:256
BALL::TFFT1D::getFourierSpaceMax
double getFourierSpaceMax() const
Definition: FFT1D.h:402
BALL::TFFT1D::getFourierStepWidth
double getFourierStepWidth() const
Definition: FFT1D.h:378
BALL::TFFT1D::getData
Complex getData(const double pos) const
Definition: FFT1D.h:460
BALL::TFFT1D::~TFFT1D
virtual ~TFFT1D()
Destructor.
BALL::operator<<
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
BALL::Complex
std::complex< BALL_COMPLEX_PRECISION > Complex
Definition: complex.h:21
BALL::TFFT1D::operator[]
Complex & operator[](const double pos)
Definition: FFT1D.h:538
BALL::TFFT1D::operator=
const TFFT1D & operator=(const TFFT1D &fft1d)
Assignment operator.
BALL::TFFT1D::maxPhys_
double maxPhys_
Definition: FFT1D.h:252
BALL::Maths::max
T max(const T &a, const T &b)
Definition: MATHS/common.h:75
BALL_CREATE
#define BALL_CREATE(name)
Definition: create.h:62
BALL::TFFT1D::getNumberOfInverseTransforms
Size getNumberOfInverseTransforms() const
Definition: FFT1D.h:414
BALL::TFFT1D::doFFT
void doFFT()
BALL::TFFT1D::operator==
bool operator==(const TFFT1D &fft1d) const
Definition: FFT1D.h:296
BALL::TFFT1D::getPhysSpaceMax
double getPhysSpaceMax() const
Definition: FFT1D.h:390
BALL::TFFT1D::setData
void setData(double pos, Complex val)
Definition: FFT1D.h:516
BALL::TRegularData1D::setDimension
BALL_INLINE void setDimension(const CoordinateType &dimension)
Definition: regularData1D.h:290
float
regularData1D.h
BALL::TRegularData1D::setOrigin
BALL_INLINE void setOrigin(const CoordinateType &origin)
Definition: regularData1D.h:276