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
10#endif
11
12#ifndef BALL_DATATYPE_REGULARDATA1D_H
14#endif
15
16#include <cmath>
17#include <complex>
18#include <fftw3.h>
19
22
23namespace 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
51
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_;
251 double minPhys_;
252 double maxPhys_;
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());
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
#define BALL_CREATE(name)
Definition: create.h:62
BALL_EXPORT std::ostream & operator<<(std::ostream &os, const Exception::GeneralException &e)
Definition: constants.h:13
TFFT1D< BALL_FFTW_DEFAULT_TRAITS > FFT1D
Definition: FFT1D.h:267
std::complex< BALL_COMPLEX_PRECISION > Complex
Definition: complex.h:21
BALL_INDEX_TYPE Index
BALL_EXTERN_VARIABLE const double PI
PI.
Definition: constants.h:35
BALL_EXTERN_VARIABLE const double h
Definition: constants.h:102
T max(const T &a, const T &b)
Definition: MATHS/common.h:75
double rint(double x)
round to integral value in floating-point format
Definition: MATHS/common.h:327
T min(const T &a, const T &b)
Definition: MATHS/common.h:102
long floor(const T &t)
Definition: MATHS/common.h:284
BALL_INLINE size_type size() const
BALL_INLINE void setDimension(const CoordinateType &dimension)
BALL_INLINE void setOrigin(const CoordinateType &origin)
const ValueType & operator[](const IndexType &index) const
bool translate(double trans_origin)
Definition: FFT1D.h:330
Size length_
Definition: FFT1D.h:244
double getFourierSpaceMax() const
Definition: FFT1D.h:402
double getGridCoordinates(Position position) const
Definition: FFT1D.h:421
TFFT1D(const TFFT1D &data)
Copy constructor.
Size numFourierToPhys_
Definition: FFT1D.h:247
void setNumberOfFFTTransforms(Size num)
Definition: FFT1D.h:220
ComplexTraits::FftwPlan planForward_
Definition: FFT1D.h:256
std::complex< typename ComplexTraits::ComplexPrecision > Complex
Definition: FFT1D.h:42
virtual void clear()
void setData(double pos, Complex val)
Definition: FFT1D.h:516
ComplexTraits::FftwPlan planBackward_
Definition: FFT1D.h:257
const Complex & operator[](const Position &pos) const
Definition: FFT1D.h:215
bool inFourierSpace_
Definition: FFT1D.h:245
Complex & operator[](const Position &pos)
Definition: FFT1D.h:206
bool isInFourierSpace() const
Definition: FFT1D.h:603
double origin_
Definition: FFT1D.h:248
Complex getInterpolatedValue(const double pos) const
Definition: FFT1D.h:483
virtual ~TFFT1D()
Destructor.
Complex phase(const double pos) const
Definition: FFT1D.h:592
double getPhysStepWidth() const
Definition: FFT1D.h:372
Size getNumberOfInverseTransforms() const
Definition: FFT1D.h:414
double stepPhys_
Definition: FFT1D.h:249
TRegularData1D< std::complex< typename ComplexTraits::ComplexPrecision > > ComplexVector
Definition: FFT1D.h:43
Complex getData(const double pos) const
Definition: FFT1D.h:460
bool setPhysStepWidth(double new_width)
Definition: FFT1D.h:351
Complex * dataAdress_
Definition: FFT1D.h:260
double minPhys_
Definition: FFT1D.h:251
void setNumberOfiFFTTransforms(Size num)
Definition: FFT1D.h:225
double maxFourier_
Definition: FFT1D.h:254
TFFT1D()
Default constructor.
Definition: FFT1D.h:285
bool operator==(const TFFT1D &fft1d) const
Definition: FFT1D.h:296
double stepFourier_
Definition: FFT1D.h:250
Size getMaxIndex() const
Definition: FFT1D.h:408
double getFourierSpaceMin() const
Definition: FFT1D.h:396
double minFourier_
Definition: FFT1D.h:253
Complex & operator[](const double pos)
Definition: FFT1D.h:538
bool planCalculated_
Definition: FFT1D.h:261
double getPhysSpaceMin() const
Definition: FFT1D.h:384
double getFourierStepWidth() const
Definition: FFT1D.h:378
double getPhysSpaceMax() const
Definition: FFT1D.h:390
double maxPhys_
Definition: FFT1D.h:252
virtual void destroy()
const TFFT1D & operator=(const TFFT1D &fft1d)
Assignment operator.
TFFT1D(Size ldn, double stepPhys=1., double origin=0., bool inFourierSpace=false)
Size numPhysToFourier_
Definition: FFT1D.h:246