libmoldeo (Moldeo 1.0 Core)  1.0
libmoldeo is the group of objects and functions that executes the basic operations of Moldeo 1.0 Platform.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
moMathFFT.h
Go to the documentation of this file.
1 /*******************************************************************************
2 
3  moMathFFT.h
4 
5  ****************************************************************************
6  * *
7  * This source is free software; you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation; either version 2 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This code is distributed in the hope that it will be useful, but *
13  * WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
15  * General Public License for more details. *
16  * *
17  * A copy of the GNU General Public License is available on the World *
18  * Wide Web at <http://www.gnu.org/copyleft/gpl.html>. You can also *
19  * obtain it by writing to the Free Software Foundation, *
20  * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
21  * *
22  ****************************************************************************
23 
24  Copyright(C) 2006 Fabricio Costa
25 
26  Authors:
27  Fabricio Costa
28  Andrés Colubri
29 
30  Portions taken from
31  Fast Fourier transform C++ header class for the FFTW3 Library
32  Copyright (C) 2004 John C. Bowman, University of Alberta
33 
34 *******************************************************************************/
35 
36 /*
37 FFTW++ is a C++ header class for Version 3 of the highly optimized FFTW
38 (http://www.fftw.org) Fourier Transform library. It provides a simple
39 interface for 1d, 2d, and 3d complex-to-complex, real-to-complex, and
40 complex-to-real Fast Fourier Transforms that takes care of the technical
41 aspects of memory allocation, alignment, planning, and wisdom. Wrappers for
42 multiple 1d transforms are also provided. As with the FFTW3 library itself,
43 both in-place and out-of-place transforms of arbitrary size are supported.
44 
45 Convenient optional shift routines that place the Fourier origin in the logical
46 center of the domain are provided for 2d and 3d complex-to-real convolution
47 applications; see fftw++.h for details.
48 
49 FFTW++ can also exploit the high-performance Array class available at
50 http://www.math.ualberta.ca/~bowman/Array (version 1.39 or higher). The
51 arrays in that package do memory bounds checking in debugging mode, but can
52 be optimized by specifying the -DNDEBUG compilation option (1d arrays
53 optimize completely to pointer operations).
54 
55 Detailed documentation is provided before each class in the fftw++.h header
56 file. The included examples illustrate how easy it is to use FFTW in C++
57 with the FFTW++ header class. Use of the Array class is optional, but
58 encouraged. If for some reason the Array class is not used, memory should
59 be allocated with FFTWComplex (or FFTWdouble) to ensure that the data is
60 optimally aligned to sizeof(Complex), to enable the SIMD extensions.
61 The optional alignment check in fftw++.h can be disabled with the
62 -DNO_CHECK_ALIGN compiler option.
63 */
64 
65 #ifndef __MO_MATH_FFT_H__
66 #define __MO_MATH_FFT_H__
67 
68 #define __FFTWPP_H_VERSION__ 1.02
69 
70 #include <fstream>
71 #include <iostream>
72 #include <fftw3.h>
73 
74 #ifndef __Complex_h__
75 #include <complex>
76 typedef std::complex<double> Complex;
77 #endif
78 
79 using std::ifstream;
80 using std::ofstream;
81 using std::cerr;
82 using std::endl;
83 
84 #ifdef __Array_h__
85 using Array::array1;
86 using Array::array2;
87 using Array::array3;
88 
89 static array1<Complex> NULL1;
90 static array2<Complex> NULL2;
91 static array3<Complex> NULL3;
92 #endif
93 
94 inline Complex *FFTWComplex(size_t size)
95 {
96  static const size_t offset = sizeof(size_t)/sizeof(Complex)+
97  (sizeof(size_t) % sizeof(Complex) > 0);
98  void *alloc=fftw_malloc((size+offset)*sizeof(Complex));
99  if(size && !alloc) cerr << endl << "Memory limits exceeded" << endl;
100  *(size_t *) alloc=size;
101  Complex*p=(Complex *)alloc+offset;
102  for(size_t i=0; i < size; i++) new(p+i) Complex;
103  return p;
104 }
105 
106 inline double *FFTWdouble(size_t size)
107 {
108  static const size_t offset = sizeof(size_t)/sizeof(Complex)+
109  (sizeof(size_t) % sizeof(Complex) > 0);
110  void *alloc=fftw_malloc(size*sizeof(double)+offset*sizeof(Complex));
111  if(size && !alloc) cerr << endl << "Memory limits exceeded" << endl;
112  *(size_t *) alloc=size;
113  double*p=(double*)((Complex *)alloc+offset);
114  for(size_t i=0; i < size; i++) new(p+i) double;
115  return p;
116 }
117 
118 template<class T>
119 inline void FFTWdelete(T *p)
120 {
121  static const size_t offset = sizeof(size_t)/sizeof(Complex)+
122  (sizeof(size_t) % sizeof(Complex) > 0);
123  void *alloc=(Complex *)p-offset;
124  size_t size=*(size_t *) alloc;
125  for(size_t i=size-1; i > 0; i--) ((T *) p)[i].~T();
126  ((T *) p)[0].~T();
127  fftw_free(alloc);
128 }
129 
130 inline void fftw_export_wisdom(void (*emitter)(char c, ofstream& s),
131  ofstream& s)
132 {
133  fftw_export_wisdom((void (*) (char, void *)) emitter,(void *) &s);
134 }
135 
136 inline int fftw_import_wisdom(int (*g)(ifstream& s), ifstream &s)
137 {
138  return fftw_import_wisdom((int (*) (void *)) g,(void *) &s);
139 }
140 
141 inline void PutWisdom(char c, ofstream& s) {s.put(c);}
142 inline int GetWisdom(ifstream& s) {return s.get();}
143 
144 // Base clase for fft routines
145 //
146 class fftw {
147 protected:
148  unsigned int size;
149  int sign;
150  double norm;
151  bool shift;
152 
153  bool inplace;
154  fftw_plan plan;
155 
156  static unsigned int effort;
157  static bool Wise;
158  static const char *WisdomName;
159  static ifstream ifWisdom;
160  static ofstream ofWisdom;
161 
162  unsigned int Dist(unsigned int n, unsigned int stride, unsigned int dist) {
163  return dist ? dist : ((stride == 1) ? n : 1);
164  }
165 
166  unsigned int realsize(unsigned int n, Complex *in, Complex *out) {
167  return n/2+(!out || in == out);
168  }
169 
170  unsigned int realsize(unsigned int n, Complex *in, double *out) {
171  return realsize(n,in,(Complex *) out);
172  }
173 
174  // Shift the Fourier origin to (nx/2,0)
175  void Shift(Complex *data, unsigned int nx, unsigned int ny) {
176  const unsigned int nyp=ny/2+1;
177  Complex *pstop=data+nx*nyp;
178  int pinc=2*nyp;
179  for(Complex *p=data+nyp; p < pstop; p += pinc) {
180  //#pragma ivdep
181  for(unsigned int j=0; j < nyp; j++) p[j]=-p[j];
182  }
183  }
184 
185  // Shift the Fourier origin to (nx/2,ny/2,0)
186  void Shift(Complex *data, unsigned int nx, unsigned int ny,
187  unsigned int nz) {
188  const unsigned int nzp=nz/2+1;
189  const unsigned int nyzp=ny*nzp;
190  const unsigned int pinc=2*nzp;
191  Complex *p,*pstop;
192  p=pstop=data;
193  for(unsigned i=0; i < nx; i++) {
194  if(i % 2) p -= nzp;
195  else p += nzp;
196  pstop += nyzp;
197  for(; p < pstop; p += pinc) {
198  //#pragma ivdep
199  for(unsigned int k=0; k < nzp; k++) p[k]=-p[k];
200  }
201  }
202  }
203 
204 public:
205  fftw(unsigned int size, int sign, unsigned int n=0) :
206  size(size), sign(sign), norm(1.0/(n ? n : size)), shift(false) {}
207 
208  virtual ~fftw() {}
209 
210  virtual fftw_plan Plan(Complex *in, Complex *out)=0;
211 
212  inline void CheckAlign(Complex *p, const char *s) {
213  if((size_t) p % sizeof(Complex) == 0) return;
214  cerr << "WARNING: " << s << " array is not " << sizeof(Complex)
215  << "-byte aligned: address " << p << endl;
216  }
217 
218  void Setup(Complex *in, Complex *out=NULL) {
219  if(!Wise) LoadWisdom();
220  bool alloc=!in;
221  if(alloc) in=FFTWComplex(size);
222 #ifndef NO_CHECK_ALIGN
223  CheckAlign(in,"constructor input");
224  if(out) CheckAlign(out,"constructor output");
225  else out=in;
226 #else
227  if(!out) out=in;
228 #endif
229  inplace=(out==in);
230  plan=Plan(in,out);
231  if(!plan) {
232  cerr << "Unable to construct FFTW plan" << endl;
233  exit(1);
234  }
235 
236  if(alloc) FFTWdelete(in);
237  SaveWisdom();
238  }
239 
240  void Setup(Complex *in, double *out) {Setup(in,(Complex *) out);}
241  void Setup(double *in, Complex *out) {Setup((Complex *) in,out);}
242 
243  void LoadWisdom() {
244  ifWisdom.open(WisdomName);
245  fftw_import_wisdom(GetWisdom,ifWisdom);
246  ifWisdom.close();
247  Wise=true;
248  }
249 
250  void SaveWisdom() {
251  ofWisdom.open(WisdomName);
252  fftw_export_wisdom(PutWisdom,ofWisdom);
253  ofWisdom.close();
254  }
255 
256  virtual void Execute(Complex *in, Complex *out) {
257  fftw_execute_dft(plan,(fftw_complex *) in,(fftw_complex *) out);
258  }
259 
260  void Setout(Complex *in, Complex *&out) {
261 #ifndef NO_CHECK_ALIGN
262  CheckAlign(in,"input");
263  if(out) CheckAlign(out,"output");
264  else out=in;
265 #else
266  if(!out) out=in;
267 #endif
268  if(inplace ^ (out == in)) {
269  cerr << "ERROR: fft constructor and call must be either both in place or out of place" << endl;
270  exit(1);
271  }
272  }
273 
274  void fft(Complex *in, Complex *out=NULL) {
275  Setout(in,out);
276  Execute(in,out);
277  }
278 
279  void fft(double *in, Complex *out) {
280  fft((Complex *) in,out);
281  }
282 
283  void fft(Complex *in, double *out) {
284  fft(in,(Complex *) out);
285  }
286 
287  void fft0(Complex *in, Complex *out=NULL) {
288  Setout(in,out);
289  shift=true;
290  Execute(in,out);
291  shift=false;
292  }
293 
294  void fft0(double *in, Complex *out) {
295  fft0((Complex *) in,out);
296  }
297 
298  void fft0(Complex *in, double *out) {
299  fft0(in,(Complex *) out);
300  }
301 
302  void Normalize(Complex *out) {
303  for(unsigned int i=0; i < size; i++) out[i] *= norm;
304  }
305 
306  virtual void fftNormalized(Complex *in, Complex *out=NULL) {
307  Setout(in,out);
308  Execute(in,out);
309  Normalize(out);
310  }
311 
312  void fftNormalized(Complex *in, double *out) {
313  fftNormalized(in,(Complex *) out);
314  }
315 
316  void fftNormalized(double *in, Complex *out) {
317  fftNormalized((Complex *) in,out);
318  }
319 
320  void fft0Normalized(Complex *in, Complex *out=NULL) {
321  Setout(in,out);
322  shift=true;
323  Execute(in,out);
324  shift=false;
325  Normalize(out);
326  }
327 
328  void fft0Normalized(Complex *in, double *out) {
329  fft0Normalized(in,(Complex *) out);
330  }
331 
332  void fft0Normalized(double *in, Complex *out) {
333  fft0Normalized((Complex *) in,out);
334  }
335 
336  void fftNormalized(Complex *in, Complex *out,
337  unsigned int nx, unsigned int m,
338  unsigned int stride, unsigned int dist) {
339  if(stride == 1 && dist == nx) fftw::fftNormalized(in,out);
340  else if(stride == nx && dist == 1) fftw::fftNormalized(in,out);
341  else {
342  Setout(in,out);
343  Execute(in,out);
344  for(unsigned int k=0; k < m; k++) {
345  for(unsigned int j=0; j < nx; j++) {
346  out[j*stride+k*dist] *= norm;
347  }
348  }
349  }
350  }
351 
352 };
353 
354 // Compute the complex Fourier transform of n complex values.
355 // Before calling fft(), the arrays in and out (which may coincide) must be
356 // allocated as Complex[n].
357 //
358 // Out-of-place usage:
359 //
360 // fft1d Forward(n,-1,in,out);
361 // Forward.fft(in,out);
362 //
363 // fft1d Backward(n,1,out,in);
364 // Backward.fft(out,in);
365 //
366 // fft1d Backward(n,1,out,in);
367 // Backward.fftNormalized(out,in); // True inverse of Forward.fft(in,out);
368 //
369 // In-place usage:
370 //
371 // fft1d Forward(n,-1);
372 // Forward.fft(in);
373 //
374 // fft1d Backward(n,1);
375 // Backward.fft(in);
376 //
377 class fft1d : public fftw {
378  unsigned int nx;
379 public:
380  fft1d(unsigned int nx, int sign, Complex *in=NULL, Complex *out=NULL)
381  : fftw(nx,sign), nx(nx) {Setup(in,out);}
382 
383 #ifdef __Array_h__
384  fft1d(int sign, const array1<Complex>& in, const array1<Complex>& out=NULL1)
385  : fftw(in.Nx(),sign), nx(in.Nx()) {Setup(in,out);}
386 #endif
387 
388  fftw_plan Plan(Complex *in, Complex *out) {
389  return fftw_plan_dft_1d(nx,(fftw_complex *) in,(fftw_complex *) out,
390  sign,effort);
391  }
392 };
393 
394 // Compute the complex Fourier transform of m complex vectors, each of
395 // length n.
396 // Before calling fft(), the arrays in and out (which may coincide) must be
397 // allocated as Complex[m*n].
398 //
399 // Out-of-place usage:
400 //
401 // mfft1d Forward(n,-1,m,stride,dist,in,out);
402 // Forward.fft(in,out);
403 //
404 // In-place usage:
405 //
406 // mfft1d Forward(n,-1,m,stride,dist);
407 // Forward.fft(in);
408 //
409 // Notes:
410 // stride is the spacing between the elements of each Complex vector;
411 // dist is the spacing between the first elements of the vectors;
412 //
413 //
414 class mfft1d : public fftw {
415  unsigned int nx;
416  unsigned int m;
417  unsigned int stride;
418  unsigned int dist;
419 public:
420  mfft1d(unsigned int nx, int sign, unsigned int m=1, unsigned int stride=1,
421  unsigned int dist=0, Complex *in=NULL, Complex *out=NULL)
422  : fftw((nx-1)*stride+(m-1)*Dist(nx,stride,dist)+1,sign,nx),
423  nx(nx), m(m), stride(stride), dist(Dist(nx,stride,dist))
424  {Setup(in,out);}
425 
426  fftw_plan Plan(Complex *in, Complex *out) {
427  int n[1]={nx};
428  return fftw_plan_many_dft(1,n,m,
429  (fftw_complex *) in,NULL,stride,dist,
430  (fftw_complex *) out,NULL,stride,dist,
431  sign,effort);
432  }
433 
434  void fftNormalized(Complex *in, Complex *out=NULL) {
435  fftw::fftNormalized(in,out,nx,m,stride,dist);
436  }
437 };
438 
439 // Compute the complex Fourier transform of n real values, using phase sign -1.
440 // Before calling fft(), the array in must be allocated as double[n] and
441 // the array out must be allocated as Complex[n/2+1]. The arrays in and out
442 // may coincide, in which case they must both be allocated as Complex[n/2+1].
443 //
444 // Out-of-place usage:
445 //
446 // rcfft1d Forward(n,in,out);
447 // Forward.fft(in,out);
448 //
449 // In-place usage:
450 //
451 // rcfft1d Forward(n);
452 // Forward.fft(out);
453 //
454 // Notes:
455 // in contains the n real values stored as a Complex array;
456 // out contains the first n/2+1 Complex Fourier values.
457 //
458 class rcfft1d : public fftw {
459  unsigned int nx;
460 public:
461  rcfft1d(unsigned int nx, Complex *out=NULL)
462  : fftw(nx/2+1,-1,nx), nx(nx) {Setup(out);}
463 
464  rcfft1d(unsigned int nx, double *in, Complex *out=NULL)
465  : fftw(nx/2+1,-1,nx), nx(nx) {Setup(in,out);}
466 
467 #ifdef __Array_h__
468  rcfft1d(const array1<Complex>& in)
469  : fftw(in.Size(),-1,2*(in.Nx()-1)), nx(2*(in.Nx()-1)) {Setup(in);}
470 
471  rcfft1d(const array1<double>& in, const array1<Complex>& out)
472  : fftw(out.Size(),-1,in.Size()), nx(in.Nx()) {Setup(in,out);}
473 #endif
474 
475  fftw_plan Plan(Complex *in, Complex *out) {
476  return fftw_plan_dft_r2c_1d(nx,(double *) in,(fftw_complex *) out, effort);
477  }
478 
479  void Execute(Complex *in, Complex *out) {
480  fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
481  }
482 };
483 
484 // Compute the real inverse Fourier transform of the n/2+1 Complex values
485 // corresponding to the non-negative part of the frequency spectrum, using
486 // phase sign +1.
487 // Before calling fft(), the array in must be allocated as Complex[n/2+1]
488 // and the array out must be allocated as double[n]. The arrays in and out
489 // may coincide, in which case they must both be allocated as Complex[n/2+1].
490 //
491 // Out-of-place usage:
492 //
493 // crfft1d Backward(n,in,out);
494 // Backward.fft(in,out);
495 //
496 // In-place usage:
497 //
498 // crfft1d Backward(n);
499 // Backward.fft(in);
500 //
501 // Notes:
502 // in contains the first n/2+1 Complex Fourier values.
503 // out contains the n real values stored as a Complex array;
504 //
505 class crfft1d : public fftw {
506  unsigned int nx;
507 public:
508  crfft1d(unsigned int nx, Complex *in=NULL)
509  : fftw(nx/2+1,1,nx), nx(nx) {Setup(in);}
510 
511  crfft1d(unsigned int nx, Complex *in, double *out)
512  : fftw(realsize(nx,in,out),1,nx), nx(nx) {Setup(in,out);}
513 
514 #ifdef __Array_h__
515  crfft1d(const array1<Complex>& in)
516  : fftw(in.Size(),1,2*(in.Nx()-1)), nx(2*(in.Nx()-1)) {Setup(in);}
517 
518  crfft1d(const array1<Complex>& in, const array1<double>& out)
519  : fftw(out.Size()/2,1,out.Size()), nx(out.Nx()) {Setup(in,out);}
520 #endif
521 
522  fftw_plan Plan(Complex *in, Complex *out) {
523  return fftw_plan_dft_c2r_1d(nx,(fftw_complex *) in,(double *) out,effort);
524  }
525 
526  void Execute(Complex *in, Complex *out) {
527  fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
528  }
529 };
530 
531 // Compute the real Fourier transform of m real vectors, each of length n,
532 // using phase sign -1. Before calling fft(), the array in must be
533 // allocated as double[m*n] and the array out must be allocated as
534 // Complex[m*(n/2+1)]. The arrays in and out may coincide, in which case
535 // they must both be allocated as Complex[m*(n/2+1)].
536 //
537 // Out-of-place usage:
538 //
539 // mrcfft1d Forward(n,m,stride,dist,in,out);
540 // Forward.fft(in,out);
541 //
542 // In-place usage:
543 //
544 // mrcfft1d Forward(n,m,stride,dist);
545 // Forward.fft(out);
546 //
547 // Notes:
548 // stride is the spacing between the elements of each Complex vector;
549 // dist is the spacing between the first elements of the vectors;
550 // in contains the n real values stored as a Complex array;
551 // out contains the first n/2+1 Complex Fourier values.
552 //
553 class mrcfft1d : public fftw {
554  unsigned int nx;
555  unsigned int m;
556  unsigned int stride;
557  unsigned int dist;
558 public:
559  mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
560  unsigned int dist=0, Complex *out=NULL)
561  : fftw(nx/2*stride+(m-1)*Dist(nx,stride,dist)+1,-1,nx), nx(nx), m(m),
562  stride(stride), dist(Dist(nx,stride,dist)) {Setup(out);}
563 
564  mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
565  unsigned int dist=0, double *in=NULL, Complex *out=NULL)
566  : fftw(nx/2*stride+(m-1)*Dist(nx,stride,dist)+1,-1,nx), nx(nx), m(m),
567  stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
568 
569  fftw_plan Plan(Complex *in, Complex *out) {
570  const int n[1]={nx};
571  return fftw_plan_many_dft_r2c(1,n,m,
572  (double *) in,NULL,stride,2*dist,
573  (fftw_complex *) out,NULL,stride,dist,
574  effort);
575  }
576 
577  void Execute(Complex *in, Complex *out) {
578  fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
579  }
580 
581  void fftNormalized(Complex *in, Complex *out=NULL) {
582  fftw::fftNormalized(in,out,nx/2+1,m,stride,dist);
583  }
584 };
585 
586 // Compute the real inverse Fourier transform of m complex vectors, each of
587 // length n/2+1, corresponding to the non-negative parts of the frequency
588 // spectra, using phase sign +1. Before calling fft(), the array in must be
589 // allocated as Complex[m*(n/2+1)] and the array out must be allocated as
590 // double[m*n]. The arrays in and out may coincide, in which case they
591 // must both be allocated as Complex[m*(n/2+1)].
592 //
593 // Out-of-place usage:
594 //
595 // mcrfft1d Backward(n,m,stride,dist,in,out);
596 // Backward.fft(in,out);
597 //
598 // In-place usage:
599 //
600 // mcrfft1d Backward(n,m,stride,dist);
601 // Backward.fft(out);
602 //
603 // Notes:
604 // stride is the spacing between the elements of each Complex vector;
605 // dist is the spacing between the first elements of the vectors;
606 // in contains the first n/2+1 Complex Fourier values.
607 // out contains the n real values stored as a Complex array;
608 //
609 class mcrfft1d : public fftw {
610  unsigned int nx;
611  unsigned int m;
612  unsigned int stride;
613  unsigned int dist;
614 public:
615  mcrfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
616  unsigned int dist=0, Complex *in=NULL)
617  : fftw(nx/2*stride+(m-1)*Dist(nx,stride,dist)+1,1,nx),
618  nx(nx), m(m), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in);}
619 
620  mcrfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1,
621  unsigned int dist=0, Complex *in=NULL, double *out=NULL)
622  : fftw((realsize(nx,in,out)-1)*stride+(m-1)*Dist(nx,stride,dist)+1,1,nx),
623  nx(nx), m(m), stride(stride), dist(Dist(nx,stride,dist)) {Setup(in,out);}
624 
625  fftw_plan Plan(Complex *in, Complex *out) {
626  const int n[1]={nx};
627  return fftw_plan_many_dft_c2r(1,n,m,
628  (fftw_complex *) in,NULL,stride,dist,
629  (double *) out,NULL,stride,2*dist,
630  effort);
631  }
632 
633  void Execute(Complex *in, Complex *out) {
634  fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
635  }
636 
637  void fftNormalized(Complex *in, Complex *out=NULL) {
638  fftw::fftNormalized(in,out,(nx/2+1),m,stride,dist);
639  }
640 };
641 
642 // Compute the complex two-dimensional Fourier transform of nx times ny
643 // complex values. Before calling fft(), the arrays in and out (which may
644 // coincide) must be allocated as Complex[nx*ny].
645 //
646 // Out-of-place usage:
647 //
648 // fft2d Forward(nx,ny,-1,in,out);
649 // Forward.fft(in,out);
650 //
651 // fft2d Backward(nx,ny,1,out,in);
652 // Backward.fft(out,in);
653 //
654 // fft2d Backward(nx,ny,1,out,in);
655 // Backward.fftNormalized(out,in); // True inverse of Forward.fft(in,out);
656 //
657 // In-place usage:
658 //
659 // fft2d Forward(nx,ny,-1);
660 // Forward.fft(in);
661 //
662 // fft2d Backward(nx,ny,1);
663 // Backward.fft(in);
664 //
665 // Note:
666 // in[ny*i+j] contains the ny Complex values for each i=0,...,nx-1.
667 //
668 class fft2d : public fftw {
669  unsigned int nx;
670  unsigned int ny;
671 public:
672  fft2d(unsigned int nx, unsigned int ny, int sign, Complex *in=NULL,
673  Complex *out=NULL)
674  : fftw(nx*ny,sign), nx(nx), ny(ny) {Setup(in,out);}
675 
676 #ifdef __Array_h__
677  fft2d(int sign, const array2<Complex>& in, const array2<Complex>& out=NULL2)
678  : fftw(in.Size(),sign), nx(in.Nx()), ny(in.Ny()) {Setup(in,out);}
679 #endif
680 
681  fftw_plan Plan(Complex *in, Complex *out) {
682  return fftw_plan_dft_2d(nx,ny,(fftw_complex *) in,(fftw_complex *) out,
683  sign,effort);
684  }
685 };
686 
687 // Compute the complex two-dimensional Fourier transform of nx times ny real
688 // values, using phase sign -1.
689 // Before calling fft(), the array in must be allocated as double[nx*ny] and
690 // the array out must be allocated as Complex[nx*(ny/2+1)]. The arrays in
691 // and out may coincide, in which case they must both be allocated as
692 // Complex[nx*(ny/2+1)].
693 //
694 // Out-of-place usage:
695 //
696 // rcfft2d Forward(nx,ny,in,out);
697 // Forward.fft(in,out); // Origin of Fourier domain at (0,0)
698 // Forward.fft0(in,out); // Origin of Fourier domain at (nx/2,0)
699 //
700 //
701 // In-place usage:
702 //
703 // rcfft2d Forward(nx,ny);
704 // Forward.fft(in); // Origin of Fourier domain at (0,0)
705 // Forward.fft0(in); // Origin of Fourier domain at (nx/2,0)
706 //
707 // Notes:
708 // in contains the nx*ny real values stored as a Complex array;
709 // out contains the upper-half portion (ky >= 0) of the Complex transform.
710 //
711 class rcfft2d : public fftw {
712  unsigned int nx;
713  unsigned int ny;
714 public:
715  rcfft2d(unsigned int nx, unsigned int ny, Complex *out=NULL) :
716  fftw(nx*(ny/2+1),-1,nx*ny), nx(nx), ny(ny) {Setup(out);}
717 
718  rcfft2d(unsigned int nx, unsigned int ny, double *in, Complex *out=NULL) :
719  fftw(nx*(ny/2+1),-1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
720 
721 #ifdef __Array_h__
722  rcfft2d(const array2<Complex>& in) :
723  fftw(in.Size(),-1,in.Nx()*2*(in.Ny()-1)), nx(in.Nx()), ny(2*(in.Ny()-1))
724  {Setup(in);}
725 
726  rcfft2d(const array2<double>& in, const array2<Complex>& out) :
727  fftw(out.Size(),-1,in.Size()), nx(in.Nx()), ny(in.Ny()) {Setup(in,out);}
728 #endif
729 
730  fftw_plan Plan(Complex *in, Complex *out) {
731  return fftw_plan_dft_r2c_2d(nx,ny,(double *) in,(fftw_complex *) out,
732  effort);
733  }
734 
735  void Execute(Complex *in, Complex *out) {
736  if(shift) Shift(in,nx,ny);
737  fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
738  }
739 };
740 
741 // Compute the real two-dimensional inverse Fourier transform of the
742 // nx*(ny/2+1) Complex values corresponding to the spectral values in the
743 // half-plane ky >= 0, using phase sign +1.
744 // Before calling fft(), the array in must be allocated as
745 // Complex[nx*(ny+1)/2] and the array out must be allocated as
746 // double[nx*ny]. The arrays in and out may coincide, in which case they
747 // must both be allocated as Complex[nx*(ny/2+1)].
748 //
749 // Out-of-place usage:
750 //
751 // crfft2d Backward(nx,ny,in,out);
752 // Backward.fft(in,out); // Origin of Fourier domain at (0,0)
753 // Backward.fft0(in,out); // Origin of Fourier domain at (nx/2,0)
754 //
755 //
756 // In-place usage:
757 //
758 // crfft2d Backward(nx,ny);
759 // Backward.fft(in); // Origin of Fourier domain at (0,0)
760 // Backward.fft0(in); // Origin of Fourier domain at (nx/2,0)
761 //
762 // Notes:
763 // in contains the upper-half portion (ky >= 0) of the Complex transform;
764 // out contains the nx*ny real values stored as a Complex array.
765 //
766 class crfft2d : public fftw {
767  unsigned int nx;
768  unsigned int ny;
769 public:
770  crfft2d(unsigned int nx, unsigned int ny, Complex *in=NULL)
771  : fftw(nx*(ny/2+1),1,nx*ny), nx(nx), ny(ny) {Setup(in);}
772 
773  crfft2d(unsigned int nx, unsigned int ny, Complex *in, double *out)
774  : fftw(nx*(realsize(ny,in,out)),1,nx*ny), nx(nx), ny(ny) {Setup(in,out);}
775 
776 #ifdef __Array_h__
777  crfft2d(const array2<Complex>& in) :
778  fftw(in.Size(),1,in.Nx()*2*(in.Ny()-1)), nx(in.Nx()), ny(2*(in.Ny()-1))
779  {Setup(in);}
780 
781  crfft2d(const array2<Complex>& in, const array2<double>& out) :
782  fftw(out.Size()/2,1,out.Size()), nx(out.Nx()), ny(out.Ny())
783  {Setup(in,out);}
784 #endif
785 
786  fftw_plan Plan(Complex *in, Complex *out) {
787  return fftw_plan_dft_c2r_2d(nx,ny,(fftw_complex *) in,(double *) out,
788  effort);
789  }
790 
791  void Execute(Complex *in, Complex *out) {
792  fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
793  if(shift) Shift(out,nx,ny);
794  }
795 };
796 
797 // Compute the complex three-dimensional Fourier transform of
798 // nx times ny times nz complex values. Before calling fft(), the arrays in
799 // and out (which may coincide) must be allocated as Complex[nx*ny*nz].
800 //
801 // Out-of-place usage:
802 //
803 // fft3d Forward(nx,ny,nz,-1,in,out);
804 // Forward.fft(in,out);
805 //
806 // fft3d Backward(nx,ny,nz,1,out,in);
807 // Backward.fft(out,in);
808 //
809 // fft3d Backward(nx,ny,nz,1,out,in);
810 // Backward.fftNormalized(out,in); // True inverse of Forward.fft(in,out);
811 //
812 // In-place usage:
813 //
814 // fft3d Forward(nx,ny,nz,-1);
815 // Forward.fft(in);
816 //
817 // fft3d Backward(nx,ny,nz,1);
818 // Backward.fft(in);
819 //
820 // Note:
821 // in[nz*(ny*i+j)+k] contains the (i,j,k)th Complex value,
822 // indexed by i=0,...,nx-1, j=0,...,ny-1, and k=0,...,nz-1.
823 //
824 class fft3d : public fftw {
825  unsigned int nx;
826  unsigned int ny;
827  unsigned int nz;
828 public:
829  fft3d(unsigned int nx, unsigned int ny, unsigned int nz,
830  int sign, Complex *in=NULL, Complex *out=NULL)
831  : fftw(nx*ny*nz,sign), nx(nx), ny(ny), nz(nz) {Setup(in,out);}
832 
833 #ifdef __Array_h__
834  fft3d(int sign, const array3<Complex>& in, const array3<Complex>& out=NULL3)
835  : fftw(in.Size(),sign), nx(in.Nx()), ny(in.Ny()), nz(in.Nz())
836  {Setup(in,out);}
837 #endif
838 
839  fftw_plan Plan(Complex *in, Complex *out) {
840  return fftw_plan_dft_3d(nx,ny,nz,(fftw_complex *) in,
841  (fftw_complex *) out, sign, effort);
842  }
843 };
844 
845 // Compute the complex two-dimensional Fourier transform of
846 // nx times ny times nz real values, using phase sign -1.
847 // Before calling fft(), the array in must be allocated as double[nx*ny*nz]
848 // and the array out must be allocated as Complex[nx*ny*(nz/2+1)]. The
849 // arrays in and out may coincide, in which case they must both be allocated as
850 // Complex[nx*ny*(nz/2+1)].
851 //
852 // Out-of-place usage:
853 //
854 // rcfft3d Forward(nx,ny,nz,in,out);
855 // Forward.fft(in,out); // Origin of Fourier domain at (0,0)
856 // Forward.fft0(in,out); // Origin of Fourier domain at (nx/2,ny/2,0)
857 //
858 //
859 // In-place usage:
860 //
861 // rcfft3d Forward(nx,ny,nz);
862 // Forward.fft(in); // Origin of Fourier domain at (0,0)
863 // Forward.fft0(in); // Origin of Fourier domain at (nx/2,ny/2,0)
864 //
865 // Notes:
866 // in contains the nx*ny*nz real values stored as a Complex array;
867 // out contains the upper-half portion (kz >= 0) of the Complex transform.
868 //
869 class rcfft3d : public fftw {
870  unsigned int nx;
871  unsigned int ny;
872  unsigned int nz;
873 public:
874  rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL)
875  : fftw(nx*ny*(nz/2+1),-1,nx*ny*nz), nx(nx), ny(ny), nz(nz) {Setup(out);}
876 
877  rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *in,
878  Complex *out=NULL) : fftw(nx*ny*(nz/2+1),-1,nx*ny*nz),
879  nx(nx), ny(ny), nz(nz) {Setup(in,out);}
880 
881 #ifdef __Array_h__
882  rcfft3d(const array3<Complex>& in) :
883  fftw(in.Size(),-1,in.Nx()*in.Ny()*2*(in.Nz()-1)),
884  nx(in.Nx()), ny(in.Ny()), nz(2*(in.Nz()-1)) {Setup(in);}
885 
886  rcfft3d(const array3<double>& in, const array3<Complex>& out) :
887  fftw(out.Size(),-1,in.Size()),
888  nx(in.Nx()), ny(in.Ny()), nz(in.Nz()) {Setup(in,out);}
889 #endif
890 
891  fftw_plan Plan(Complex *in, Complex *out) {
892  return fftw_plan_dft_r2c_3d(nx,ny,nz,(double *) in,(fftw_complex *) out,
893  effort);
894  }
895 
896  void Execute(Complex *in, Complex *out) {
897  if(shift) Shift(in,nx,ny,nz);
898  fftw_execute_dft_r2c(plan,(double *) in,(fftw_complex *) out);
899  }
900 };
901 
902 // Compute the real two-dimensional inverse Fourier transform of the
903 // nx*ny*(nz/2+1) Complex values corresponding to the spectral values in the
904 // half-plane kz >= 0, using phase sign +1.
905 // Before calling fft(), the array in must be allocated as
906 // Complex[nx*ny*(nz+1)/2] and the array out must be allocated as
907 // double[nx*ny*nz]. The arrays in and out may coincide, in which case they
908 // must both be allocated as Complex[nx*ny*(nz/2+1)].
909 //
910 // Out-of-place usage:
911 //
912 // crfft3d Backward(nx,ny,nz,in,out);
913 // Backward.fft(in,out); // Origin of Fourier domain at (0,0)
914 // Backward.fft0(in,out); // Origin of Fourier domain at (nx/2,ny/2,0)
915 //
916 // In-place usage:
917 //
918 // crfft3d Backward(nx,ny,nz);
919 // Backward.fft(in); // Origin of Fourier domain at (0,0)
920 // Backward.fft0(in); // Origin of Fourier domain at (nx/2,ny/2,0)
921 //
922 // Notes:
923 // in contains the upper-half portion (kz >= 0) of the Complex transform;
924 // out contains the nx*ny*nz real values stored as a Complex array.
925 //
926 class crfft3d : public fftw {
927  unsigned int nx;
928  unsigned int ny;
929  unsigned int nz;
930 public:
931  crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in=NULL)
932  : fftw(nx*ny*(nz/2+1),1,nx*ny*nz), nx(nx), ny(ny), nz(nz)
933  {Setup(in);}
934  crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in,
935  double *out) : fftw(nx*ny*(realsize(nz,in,out)),1,nx*ny*nz),
936  nx(nx), ny(ny), nz(nz) {Setup(in,out);}
937 #ifdef __Array_h__
938  crfft3d(const array3<Complex>& in) :
939  fftw(in.Size(),1,in.Nx()*in.Ny()*2*(in.Nz()-1)),
940  nx(in.Nx()), ny(in.Ny()), nz(2*(in.Nz()-1)) {Setup(in);}
941 
942  crfft3d(const array3<Complex>& in, const array3<double>& out) :
943  fftw(out.Size()/2,1,out.Size()),
944  nx(out.Nx()), ny(out.Ny()), nz(out.Nz()) {Setup(in,out);}
945 #endif
946 
947  fftw_plan Plan(Complex *in, Complex *out) {
948  return fftw_plan_dft_c2r_3d(nx,ny,nz,(fftw_complex *) in,(double *) out,
949  effort);
950  }
951 
952  void Execute(Complex *in, Complex *out) {
953  fftw_execute_dft_c2r(plan,(fftw_complex *) in,(double *) out);
954  if(shift) Shift(out,nx,ny,nz);
955  }
956 };
957 
958 #endif
959 
void Shift(Complex *data, unsigned int nx, unsigned int ny, unsigned int nz)
Definition: moMathFFT.h:186
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:475
void Setout(Complex *in, Complex *&out)
Definition: moMathFFT.h:260
int GetWisdom(ifstream &s)
Definition: moMathFFT.h:142
crfft1d(unsigned int nx, Complex *in, double *out)
Definition: moMathFFT.h:511
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:791
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:633
mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1, unsigned int dist=0, Complex *out=NULL)
Definition: moMathFFT.h:559
static bool Wise
Definition: moMathFFT.h:157
crfft2d(unsigned int nx, unsigned int ny, Complex *in=NULL)
Definition: moMathFFT.h:770
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:952
void fftNormalized(double *in, Complex *out)
Definition: moMathFFT.h:316
void fftNormalized(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:434
bool inplace
Definition: moMathFFT.h:153
fft3d(unsigned int nx, unsigned int ny, unsigned int nz, int sign, Complex *in=NULL, Complex *out=NULL)
Definition: moMathFFT.h:829
crfft1d(unsigned int nx, Complex *in=NULL)
Definition: moMathFFT.h:508
fft1d(unsigned int nx, int sign, Complex *in=NULL, Complex *out=NULL)
Definition: moMathFFT.h:380
unsigned int realsize(unsigned int n, Complex *in, double *out)
Definition: moMathFFT.h:170
void Setup(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:218
std::complex< double > Complex
Definition: moMathFFT.h:76
void fftNormalized(Complex *in, double *out)
Definition: moMathFFT.h:312
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:896
mcrfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1, unsigned int dist=0, Complex *in=NULL)
Definition: moMathFFT.h:615
double norm
Definition: moMathFFT.h:150
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:526
unsigned int size
Definition: moMathFFT.h:148
void fft(double *in, Complex *out)
Definition: moMathFFT.h:279
void SaveWisdom()
Definition: moMathFFT.h:250
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:786
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:569
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:730
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:388
crfft2d(unsigned int nx, unsigned int ny, Complex *in, double *out)
Definition: moMathFFT.h:773
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:891
unsigned int Dist(unsigned int n, unsigned int stride, unsigned int dist)
Definition: moMathFFT.h:162
virtual fftw_plan Plan(Complex *in, Complex *out)=0
static unsigned int effort
Definition: moMathFFT.h:156
void fft(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:274
static const char * WisdomName
Definition: moMathFFT.h:158
void fft0Normalized(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:320
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:577
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:681
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:735
rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *out=NULL)
Definition: moMathFFT.h:874
mrcfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1, unsigned int dist=0, double *in=NULL, Complex *out=NULL)
Definition: moMathFFT.h:564
static ifstream ifWisdom
Definition: moMathFFT.h:159
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:426
virtual void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:256
void fft0Normalized(double *in, Complex *out)
Definition: moMathFFT.h:332
void fftNormalized(Complex *in, Complex *out, unsigned int nx, unsigned int m, unsigned int stride, unsigned int dist)
Definition: moMathFFT.h:336
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:522
void Execute(Complex *in, Complex *out)
Definition: moMathFFT.h:479
void fftw_export_wisdom(void(*emitter)(char c, ofstream &s), ofstream &s)
Definition: moMathFFT.h:130
void PutWisdom(char c, ofstream &s)
Definition: moMathFFT.h:141
rcfft2d(unsigned int nx, unsigned int ny, double *in, Complex *out=NULL)
Definition: moMathFFT.h:718
void fft0(double *in, Complex *out)
Definition: moMathFFT.h:294
int sign
Definition: moMathFFT.h:149
crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in=NULL)
Definition: moMathFFT.h:931
unsigned int realsize(unsigned int n, Complex *in, Complex *out)
Definition: moMathFFT.h:166
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:839
virtual ~fftw()
Definition: moMathFFT.h:208
rcfft3d(unsigned int nx, unsigned int ny, unsigned int nz, double *in, Complex *out=NULL)
Definition: moMathFFT.h:877
bool shift
Definition: moMathFFT.h:151
int fftw_import_wisdom(int(*g)(ifstream &s), ifstream &s)
Definition: moMathFFT.h:136
void Setup(Complex *in, double *out)
Definition: moMathFFT.h:240
Complex * FFTWComplex(size_t size)
Definition: moMathFFT.h:94
void FFTWdelete(T *p)
Definition: moMathFFT.h:119
mfft1d(unsigned int nx, int sign, unsigned int m=1, unsigned int stride=1, unsigned int dist=0, Complex *in=NULL, Complex *out=NULL)
Definition: moMathFFT.h:420
rcfft1d(unsigned int nx, double *in, Complex *out=NULL)
Definition: moMathFFT.h:464
void Setup(double *in, Complex *out)
Definition: moMathFFT.h:241
void fft0(Complex *in, double *out)
Definition: moMathFFT.h:298
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:947
static ofstream ofWisdom
Definition: moMathFFT.h:160
rcfft1d(unsigned int nx, Complex *out=NULL)
Definition: moMathFFT.h:461
fftw_plan plan
Definition: moMathFFT.h:154
fftw(unsigned int size, int sign, unsigned int n=0)
Definition: moMathFFT.h:205
rcfft2d(unsigned int nx, unsigned int ny, Complex *out=NULL)
Definition: moMathFFT.h:715
void Normalize(Complex *out)
Definition: moMathFFT.h:302
void CheckAlign(Complex *p, const char *s)
Definition: moMathFFT.h:212
virtual void fftNormalized(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:306
void fftNormalized(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:581
void Shift(Complex *data, unsigned int nx, unsigned int ny)
Definition: moMathFFT.h:175
fft2d(unsigned int nx, unsigned int ny, int sign, Complex *in=NULL, Complex *out=NULL)
Definition: moMathFFT.h:672
crfft3d(unsigned int nx, unsigned int ny, unsigned int nz, Complex *in, double *out)
Definition: moMathFFT.h:934
double * FFTWdouble(size_t size)
Definition: moMathFFT.h:106
mcrfft1d(unsigned int nx, unsigned int m=1, unsigned int stride=1, unsigned int dist=0, Complex *in=NULL, double *out=NULL)
Definition: moMathFFT.h:620
void fft0(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:287
void fftNormalized(Complex *in, Complex *out=NULL)
Definition: moMathFFT.h:637
fftw_plan Plan(Complex *in, Complex *out)
Definition: moMathFFT.h:625
void fft0Normalized(Complex *in, double *out)
Definition: moMathFFT.h:328
void LoadWisdom()
Definition: moMathFFT.h:243
void fft(Complex *in, double *out)
Definition: moMathFFT.h:283