SLIP  1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FFT.hpp
Go to the documentation of this file.
1 /*
2  * Copyright(c):
3  * Signal Image and Communications (SIC) Department
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * - University of Poitiers, France http://www.univ-poitiers.fr
6  * - XLIM Institute UMR CNRS 7252 http://www.xlim.fr/
7  *
8  * and
9  *
10  * D2 Fluid, Thermic and Combustion
11  * - University of Poitiers, France http://www.univ-poitiers.fr
12  * - PPRIME Institute - UPR CNRS 3346 http://www.pprime.fr
13  * - ISAE-ENSMA http://www.ensma.fr
14  *
15  * Contributor(s):
16  * The SLIP team,
17  * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
18  * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>,
19  * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>,
20  * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>,
21  * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>,
22  * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
23  *
24  * Description:
25  * The Simple Library of Image Processing (SLIP) is a new image processing
26  * library. It is written in the C++ language following as much as possible
27  * the ISO/ANSI C++ standard. It is consequently compatible with any system
28  * satisfying the ANSI C++ complience. It works on different Unix , Linux ,
29  * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that
30  * was created by the Signal, Image and Communications (SIC) departement of
31  * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic
32  * and Combustion departement of the P', UPR 3346 CNRS Institute of the
33  * University of Poitiers.
34  *
35  * The SLIP Library source code has been registered to the APP (French Agency
36  * for the Protection of Programs) by the University of Poitiers and CNRS,
37  * under registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.
38 
39  * http://www.sic.sp2mi.univ-poitiers.fr/slip/
40  *
41  * This software is governed by the CeCILL-C license under French law and
42  * abiding by the rules of distribution of free software. You can use,
43  * modify and/ or redistribute the software under the terms of the CeCILL-C
44  * license as circulated by CEA, CNRS and INRIA at the following URL
45  * http://www.cecill.info.
46  * As a counterpart to the access to the source code and rights to copy,
47  * modify and redistribute granted by the license, users are provided only
48  * with a limited warranty and the software's author, the holder of the
49  * economic rights, and the successive licensors have only limited
50  * liability.
51  *
52  * In this respect, the user's attention is drawn to the risks associated
53  * with loading, using, modifying and/or developing or reproducing the
54  * software by the user in light of its specific status of free software,
55  * that may mean that it is complicated to manipulate, and that also
56  * therefore means that it is reserved for developers and experienced
57  * professionals having in-depth computer knowledge. Users are therefore
58  * encouraged to load and test the software's suitability as regards their
59  * requirements in conditions enabling the security of their systems and/or
60  * data to be ensured and, more generally, to use and operate it in the
61  * same conditions as regards security.
62  *
63  * The fact that you are presently reading this means that you have had
64  * knowledge of the CeCILL-C license and that you accept its terms.
65  */
66 
67 
68 
76 #ifndef FFT_HPP
77 #define FFT_HPP
78 
79 #include <iostream>
80 #include <complex>
81 #include <cmath>
82 #include <iterator>
83 #include <algorithm>
84 #include "Array.hpp"
85 #include "macros.hpp"
86 #include "Matrix.hpp"
87 #include "Vector.hpp"
88 #include "statistics.hpp"
89 #include "FFTPACK.hpp"
90 #include <vector>
91 
92 #ifdef HAVE_FFTW
93 #include <fftw3.h>
94 #endif
95 
96 
97 
98 
99 
100 namespace slip
101 {
102 
104 
109 
116  unsigned int inline bitReverse(unsigned int x,
117  const int log2n)
118  {
119  int n = 0;
120  for(int i = 0; i < log2n; i++)
121  {
122  n <<= 1;
123  n |= (x & 1);
124  x >>= 1;
125  }
126  return n;
127  }
128 
135  unsigned int inline log2(unsigned int x)
136  {
137  assert (x!= 0);
138  unsigned int n = 0;
139  while(x > 1)
140  {
141  x = x >> 1;
142  n++;
143  }
144  return n;
145  }
146 
164  template<class ForwardIterator>
165  inline
166  void fftshift(ForwardIterator first,
167  ForwardIterator last)
168  {
169  ForwardIterator middle = first + ((slip::cardinal<std::size_t>(first,last) + 1) / 2);
170  std::rotate(first,middle,last);
171  }
172 
189  template<class ForwardIterator2D>
190  inline
191  void fftshift2d(ForwardIterator2D upper_left,
192  ForwardIterator2D bottom_right)
193  {
194  //recuperer la taille du container
195  typename ForwardIterator2D::difference_type size2d = bottom_right - upper_left;
196 
197  std::size_t nb_lig = size2d[0];
198  std::size_t nb_col = size2d[1];
199 
200 
201  //centered
202  for(size_t i = 0; i < nb_lig; ++i)
203  {
204  slip::fftshift(upper_left.row_begin(i),upper_left.row_end(i));
205  }
206  for(size_t j = 0; j < nb_col; ++j)
207  {
208  slip::fftshift(upper_left.col_begin(j),upper_left.col_end(j));
209  }
210  }
211 
228  template<class ForwardIterator3D>
229  inline
230  void fftshift3d(ForwardIterator3D front_upper_left,
231  ForwardIterator3D back_bottom_right)
232  {
233  //recuperer la taille du container
234  typename ForwardIterator3D::difference_type size3d = back_bottom_right - front_upper_left;
235 
236  std::size_t nb_sli = size3d[0];
237  std::size_t nb_lig = size3d[1];
238  std::size_t nb_col = size3d[2];
239 
240  for(size_t k = 0; k < nb_sli; ++k)
241  for(size_t i = 0; i < nb_lig; ++i)
242  {
243  slip::fftshift(front_upper_left.row_begin(k,i),front_upper_left.row_end(k,i));
244  }
245  for(size_t k = 0; k < nb_sli; ++k)
246  for(size_t j = 0; j < nb_col; ++j)
247  {
248  slip::fftshift(front_upper_left.col_begin(k,j),front_upper_left.col_end(k,j));
249  }
250  for(size_t i = 0; i < nb_lig; ++i)
251  for(size_t j = 0; j < nb_col; ++j)
252  {
253  slip::fftshift(front_upper_left.slice_begin(i,j),front_upper_left.slice_end(i,j));
254  }
255  }
256 
261 
267 
290  template<class InputIter_T,class Iter_T>
291  inline
292  void radix2_fft(InputIter_T a, Iter_T b, const int log2n)
293  {
294  typedef typename std::iterator_traits<Iter_T>::value_type complex;
295  typedef typename complex::value_type real;
296 
297  // typedef typename std::complex<typename std::iterator_traits<Iter_T>::value_type> complex;
298  real pi = slip::constants<real>::pi();
299  const complex J(0,1);
300  unsigned int n = 1 << log2n;
301  for(unsigned int i = 0; i < n; ++i)
302  {
303  b[bitReverse(i,log2n)] = a[i];
304  }
305  for(int s = 1; s <= log2n; ++s)
306  {
307  unsigned int m = 1 << s;
308  unsigned int m2 = m >> 1;
309  complex w(1,0);
310  complex wm = exp(-J*(pi/m2));
311  for(unsigned int j = 0; j < m2; ++j)
312  {
313  for(unsigned int k = j ; k < n; k += m)
314  {
315  complex t = w * b[k + m2];
316  complex u = b[k];
317  b[k] = u + t;
318  b[k + m2] = u - t;
319  }
320  w *= wm;
321  }
322  }
323 
324  }
325 
349  template<class InputIter_T,class Iter_T>
350  inline
351  void radix2_ifft(InputIter_T a, Iter_T b, const int log2n)
352  {
353  typedef typename std::iterator_traits<InputIter_T>::value_type complex;
354  typedef typename complex::value_type real;
355  real pi = slip::constants<real>::pi();
356  const complex J(0,1);
357  unsigned int n = 1 << log2n;
358  real coef = 1.0 / (real)n;
359  for(unsigned int i = 0; i < n; ++i)
360  {
361  b[bitReverse(i,log2n)] = a[i];
362  }
363  for(int s = 1; s <= log2n; ++s)
364  {
365  unsigned int m = 1 << s;
366  unsigned int m2 = m >> 1;
367  complex w(1,0);
368  complex wm = exp(J*(pi/m2));
369  for(unsigned int j = 0; j < m2; ++j)
370  {
371  for(unsigned int k = j ; k < n; k += m)
372  {
373  complex t = w * b[k + m2];
374  complex u = b[k];
375  b[k] = (u + t);
376  b[k + m2] = (u - t);
377  }
378  w *= wm;
379  }
380  }
381  for(unsigned int i = 0; i < n; ++i)
382  {
383  b[i] = coef * b[i];
384  }
385 
386  }
387 
388 
405  template<class InputIter_T, class Iter_T>
406  inline
407  void real_radix2_fft(InputIter_T a, Iter_T b, const int log2n)
408  {
409  slip::radix2_fft(a,b,log2n);
410  }
411 
412 
448  template<typename InputIterator1d,typename OuputIterator1d>
449  inline
450  void radix2_fft1d2d(InputIterator1d begin1,OuputIterator1d begin2,std::size_t nb_lig, std::size_t nb_col)
451  {
452  std::size_t n1 = slip::log2(nb_lig);
453  std::size_t n2 = slip::log2(nb_col);
454  typename std::iterator_traits<InputIterator1d>::difference_type step_horiz1 = 1;
455  typename std::iterator_traits<OuputIterator1d>::difference_type step_horiz2 = 1;
456  typename std::iterator_traits<OuputIterator1d>::difference_type step_vert2 = nb_col;
457 
458  //fft2d
459  for(size_t i = 0; i < nb_lig; ++i)
460  {
461  typename slip::stride_iterator<typename std::iterator_traits<InputIterator1d>::pointer> it_row1(&begin1[i*nb_col],step_horiz1);
462  typename slip::stride_iterator<typename std::iterator_traits<OuputIterator1d>::pointer> it_row2(&begin2[i*nb_col],step_horiz2);
463  slip::radix2_fft(it_row1,it_row2,n2);
464  }
465  slip::Vector<std::complex<double> > Vtmp(nb_lig);
466  for(size_t j = 0; j < nb_col; ++j)
467  {
469  slip::radix2_fft(it_col2,Vtmp.begin(),n1);
470  std::copy(Vtmp.begin(),Vtmp.end(),it_col2);
471  }
472  }
473 
505  template<typename InputIterator1d,typename OuputIterator1d>
506  inline
507  void radix2_ifft1d2d(InputIterator1d begin1,OuputIterator1d begin2,std::size_t nb_lig, std::size_t nb_col)
508  {
509  std::size_t n1 = slip::log2(nb_lig);
510  std::size_t n2 = slip::log2(nb_col);
511  typename std::iterator_traits<InputIterator1d>::difference_type step_horiz1 = 1;
512  typename std::iterator_traits<OuputIterator1d>::difference_type step_horiz2 = 1;
513  typename std::iterator_traits<OuputIterator1d>::difference_type step_vert2 = nb_col;
514  //ifft2d
515  for(size_t i = 0; i < nb_lig; ++i)
516  {
517  typename slip::stride_iterator<typename std::iterator_traits<InputIterator1d>::pointer> it_row1(&begin1[i*nb_col],step_horiz1);
518  typename slip::stride_iterator<typename std::iterator_traits<OuputIterator1d>::pointer> it_row2(&begin2[i*nb_col],step_horiz2);
519  slip::radix2_ifft(it_row1,it_row2,n2);
520  }
521  slip::Vector<std::complex<double> > Vtmp(nb_lig);
522  for(size_t j = 0; j < nb_col; ++j)
523  {
525  slip::radix2_ifft(it_col2,Vtmp.begin(),n1);
526  std::copy(Vtmp.begin(),Vtmp.end(),it_col2);
527  }
528  }
529 
556  template<typename InputBidirectionalIterator2d,
557  typename OutputBidirectionalIterator2d>
558  inline
559  void radix2_fft2d(InputBidirectionalIterator2d in_upper_left,
560  InputBidirectionalIterator2d in_bottom_right,
561  OutputBidirectionalIterator2d out_upper_left)
562  {
563  //recuperer la taille du container
564  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
565 
566  std::size_t nb_lig = size2d[0];
567  std::size_t nb_col = size2d[1];
568 
569  // std::size_t n1 = std::size_t(std::log(nb_lig)/std::log(2) + 0.5);
570  // std::size_t n2 = std::size_t(std::log(nb_col)/std::log(2) + 0.5);
571  std::size_t n1 = slip::log2(nb_lig);
572  std::size_t n2 = slip::log2(nb_col);
573 
574  //fft2d
575  for(size_t i = 0; i < nb_lig; ++i)
576  {
577  slip::radix2_fft(in_upper_left.row_begin(i),out_upper_left.row_begin(i),n2);
578  }
579 
580  slip::Vector<std::complex<double> > Vtmp(nb_lig);
581  for(size_t j = 0; j < nb_col; ++j)
582  {
583  slip::radix2_fft(out_upper_left.col_begin(j),Vtmp.begin(),n1);
584  std::copy(Vtmp.begin(),Vtmp.end(),out_upper_left.col_begin(j));
585  }
586 
587 
588  //centered fft2d
589  // slip::fftshift2d(out_upper_left,out_upper_left+size2d);
590  }
591 
617  template<typename InputBidirectionalIterator2d,
618  typename OutputBidirectionalIterator2d>
619  inline
620  void radix2_ifft2d(InputBidirectionalIterator2d in_upper_left,
621  InputBidirectionalIterator2d in_bottom_right,
622  OutputBidirectionalIterator2d out_upper_left)
623  {
624  //recuperer la taille du container
625  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
626 
627  std::size_t nb_lig = size2d[0];
628  std::size_t nb_col = size2d[1];
629 
630 
631  // std::size_t n1 = size_t(std::log(nb_lig)/std::log(2) + 0.5);
632  // std::size_t n2 = size_t(std::log(nb_col)/std::log(2) + 0.5);
633  std::size_t n1 = slip::log2(nb_lig);
634  std::size_t n2 = slip::log2(nb_col);
635  //centered fft2d
636  // slip::fftshift2d(in_upper_left,in_bottom_right);
637 
638  //fft2d
639  for(size_t i = 0; i < nb_lig; ++i)
640  {
641  slip::radix2_ifft(in_upper_left.row_begin(i),out_upper_left.row_begin(i),n2);
642  }
643 
644  slip::Vector<std::complex<double> > Vtmp(nb_lig);
645  for(size_t j = 0; j < nb_col; ++j)
646  {
647  slip::radix2_ifft(out_upper_left.col_begin(j),Vtmp.begin(),n1);
648  std::copy(Vtmp.begin(),Vtmp.end(),out_upper_left.col_begin(j));
649  }
650  }
651 
684  template<typename Matrix1, typename Matrix2>
685  inline
686  void radix2_fft2d(Matrix1 &datain, Matrix2 &dataout)
687  {
688  assert(datain.cols() == dataout.cols());
689  assert(datain.rows() == dataout.rows());
690  slip::radix2_fft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
691  }
692 
723  template<typename Matrix1, typename Matrix2>
724  inline
725  void radix2_ifft2d(Matrix1 &datain, Matrix2 &dataout)
726  {
727  assert(datain.cols() == dataout.cols());
728  assert(datain.rows() == dataout.rows());
729  slip::radix2_ifft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
730  }
731 
734 
740 
741 
772  template<class InputIter,class OutputIter>
773  inline
774  void real_split_radix_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
775  {
776  typedef typename std::iterator_traits<InputIter>::value_type real;
777  typedef typename std::iterator_traits<OutputIter>::value_type complex;
778  const complex J(0,1);
779  size_t length = in_end - in_begin;
780  assert(length > 1);
781  //data copy
782  slip::Array<real> in(2*length,0.0);
784  std::copy(in_begin,in_end,itst_in_begin);
785  //FFT initialization
786  slip::Array<real> r(4*length+15,0.0);
787  slip::Array<size_t> ifac(4*length+15,(size_t)0);
788  size_t double_length = 2*length;
789  slip::rffti(&double_length,r.begin(),ifac.begin());
790 
791  //FFT
792  slip::rfftf(&double_length,in.begin(),r.begin(),ifac.begin());
793 
794  //Output copy in complex type
795  out_begin[0] = in[0];
796  for(size_t i = 1; i < length; i++)
797  {
798  out_begin[i] = in[2*i-1] + J * in[2*i];
799  }
800  if(in[0] != in[2*length - 1])
801  out_begin[0]+= J * in[2*length - 1];
802  }
803 
826  template<class InputIter,class OutputIter>
827  inline
828  void complex_split_radix_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
829  {
830  typedef typename std::iterator_traits<InputIter>::value_type complex;
831  typedef typename complex::value_type real;
832  const complex J(0,1);
833  size_t length = in_end - in_begin;
834  assert(length > 1);
835 
836  //data copy
837  slip::Array<real> in(2*length,0.0);
838  typename slip::stride_iterator<typename slip::Array<real>::iterator> itst_in_oddbegin(in.begin(),2);
839  typename slip::stride_iterator<typename slip::Array<real>::iterator> itst_in_evenbegin(in.begin()+1,2);
840  std::transform(in_begin,in_end,itst_in_oddbegin,slip::un_real<complex,real>());
841  std::transform(in_begin,in_end,itst_in_evenbegin,slip::un_imag<complex,real>());
842 
843  //FFT initialization
844  slip::Array<real> r(8*length+15,0.0);
845  slip::Array<size_t> ifac(8*length+15,(size_t)0);
846  slip::cffti(&length,r.begin(),ifac.begin());
847 
848  //FFT
849  slip::cfftf(&length,in.begin(),r.begin(),ifac.begin());
850 
851  //Output copy
852  for(size_t i = 0; i < length; i++)
853  {
854  out_begin[i] = in[2*i] + J * in[2*i+1];
855  }
856  }
857 
878  template<class InputIter,class OutputIter>
879  inline
880  void split_radix_ifft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
881  {
882  typedef typename std::iterator_traits<InputIter>::value_type complex;
883  typedef typename complex::value_type real;
884  const complex J(0,1);
885  size_t length = in_end - in_begin;
886  assert(length > 1);
887 
888  // data copy
889  slip::Array<real> in(2*length,0.0);
890  typename slip::stride_iterator<typename slip::Array<real>::iterator> itst_in_oddbegin(in.begin(),2);
891  typename slip::stride_iterator<typename slip::Array<real>::iterator> itst_in_evenbegin(in.begin()+1,2);
892  std::transform(in_begin,in_end,itst_in_oddbegin,slip::un_real<complex,real>());
893  std::transform(in_begin,in_end,itst_in_evenbegin,slip::un_imag<complex,real>());
894 
895  // FFT initialization
896  slip::Array<real> r(5*length,0.0);
897  slip::Array<size_t> ifac(5*length,(size_t)0);
898  slip::cffti(&length,r.begin(),ifac.begin());
899 
900  // FFT
901  slip::cfftb(&length,in.begin(),r.begin(),ifac.begin());
902 
903  // output data copy
904  for(size_t i = 0; i < length; i++)
905  {
906  out_begin[i] = (in[2*i] + J * in[2*i+1]) / (real)length;
907  }
908  }
909 
910 
911 
941  template<typename InputBidirectionalIterator2d,
942  typename OutputBidirectionalIterator2d>
943  inline
944  void real_split_radix_fft2d(InputBidirectionalIterator2d in_upper_left,
945  InputBidirectionalIterator2d in_bottom_right,
946  OutputBidirectionalIterator2d out_upper_left)
947  {
948  //recuperer la taille du container
949  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
950 
951  std::size_t nb_lig = size2d[0];
952  std::size_t nb_col = size2d[1];
953 
954  // fft2d
955  for(size_t i = 0; i < nb_lig; ++i)
956  {
957  slip::real_split_radix_fft(in_upper_left.row_begin(i),in_upper_left.row_end(i),out_upper_left.row_begin(i));
958  }
959 
960  slip::Vector<std::complex<double> > Vtmp(nb_lig);
961  for(size_t j = 0; j < nb_col; ++j)
962  {
963  slip::complex_split_radix_fft(out_upper_left.col_begin(j),out_upper_left.col_end(j),Vtmp.begin());
964  std::copy(Vtmp.begin(),Vtmp.end(),out_upper_left.col_begin(j));
965  }
966  }
967 
981  template<typename InputBidirectionalIterator2d,
982  typename OutputBidirectionalIterator2d>
983  inline
984  void complex_split_radix_fft2d(InputBidirectionalIterator2d in_upper_left,
985  InputBidirectionalIterator2d in_bottom_right,
986  OutputBidirectionalIterator2d out_upper_left)
987  {
988  //recuperer la taille du container
989  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
990 
991  std::size_t nb_lig = size2d[0];
992  std::size_t nb_col = size2d[1];
993 
994  // fft2d
995  for(size_t i = 0; i < nb_lig; ++i)
996  {
997  slip::complex_split_radix_fft(in_upper_left.row_begin(i),in_upper_left.row_end(i),out_upper_left.row_begin(i));
998  }
999 
1000  slip::Vector<std::complex<double> > Vtmp(nb_lig);
1001  for(size_t j = 0; j < nb_col; ++j)
1002  {
1003  slip::complex_split_radix_fft(out_upper_left.col_begin(j),out_upper_left.col_end(j),Vtmp.begin());
1004  std::copy(Vtmp.begin(),Vtmp.end(),out_upper_left.col_begin(j));
1005  }
1006  }
1007 
1036  template<typename InputBidirectionalIterator2d,
1037  typename OutputBidirectionalIterator2d>
1038  inline
1039  void split_radix_ifft2d(InputBidirectionalIterator2d in_upper_left,
1040  InputBidirectionalIterator2d in_bottom_right,
1041  OutputBidirectionalIterator2d out_upper_left)
1042  {
1043  //recuperer la taille du container
1044  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
1045 
1046  std::size_t nb_lig = size2d[0];
1047  std::size_t nb_col = size2d[1];
1048 
1049 
1050  //centered fft2d
1051  // slip::fftshift2d(in_upper_left,in_bottom_right);
1052 
1053  //fft2d
1054  for(size_t i = 0; i < nb_lig; ++i)
1055  {
1056  slip::split_radix_ifft(in_upper_left.row_begin(i),in_upper_left.row_end(i),out_upper_left.row_begin(i));
1057  }
1058 
1059  slip::Vector<std::complex<double> > Vtmp(nb_lig);
1060  for(size_t j = 0; j < nb_col; ++j)
1061  {
1062  slip::split_radix_ifft(out_upper_left.col_begin(j),out_upper_left.col_end(j),Vtmp.begin());
1063  std::copy(Vtmp.begin(),Vtmp.end(),out_upper_left.col_begin(j));
1064  }
1065  }
1066 
1099  template<typename InputMatrix, typename OutputMatrix>
1100  inline void real_split_radix_fft2d(InputMatrix &datain, OutputMatrix &dataout)
1101  {
1102  assert(datain.cols() == dataout.cols());
1103  assert(datain.rows() == dataout.rows());
1104  slip::real_split_radix_fft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
1105  }
1106 
1136  template<typename InputMatrix, typename OutputMatrix>
1137  inline
1138  void split_radix_ifft2d(InputMatrix &datain, OutputMatrix &dataout)
1139  {
1140  assert(datain.cols() == dataout.cols());
1141  assert(datain.rows() == dataout.rows());
1142  slip::split_radix_ifft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
1143  }
1144 
1177  template<typename InputBidirectionalIterator3d,
1178  typename OutputBidirectionalIterator3d>
1179  inline
1180  void real_split_radix_fft3d(InputBidirectionalIterator3d in_front_upper_left,
1181  InputBidirectionalIterator3d in_back_bottom_right,
1182  OutputBidirectionalIterator3d out_front_upper_left)
1183  {
1184  //recuperer la taille du container
1185  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
1186  typedef typename std::iterator_traits<InputBidirectionalIterator3d>::value_type real;
1187  typedef typename std::iterator_traits<OutputBidirectionalIterator3d>::value_type complex;
1188 
1189  std::size_t nb_sli = size3d[0];
1190  std::size_t nb_lig = size3d[1];
1191  std::size_t nb_col = size3d[2];
1192  //fft3d
1193  for(size_t k = 0; k < nb_sli; ++k)
1194  for(size_t i = 0; i < nb_lig; ++i)
1195  {
1196  slip::real_split_radix_fft(in_front_upper_left.row_begin(k,i),in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
1197  }
1198 
1199  slip::Vector<std::complex<double> > Vtmp(nb_lig);
1200  for(size_t k = 0; k < nb_sli; ++k)
1201  for(size_t j = 0; j < nb_col; ++j)
1202  {
1203  slip::complex_split_radix_fft(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
1204  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
1205  }
1206 
1207  Vtmp.resize(nb_sli);
1208  for(size_t i = 0; i < nb_lig; ++i)
1209  for(size_t j = 0; j < nb_col; ++j)
1210  {
1211  slip::complex_split_radix_fft(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
1212  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
1213  }
1214  }
1215 
1248  template<typename InputBidirectionalIterator3d,
1249  typename OutputBidirectionalIterator3d>
1250  inline
1251  void split_radix_ifft3d(InputBidirectionalIterator3d in_front_upper_left,
1252  InputBidirectionalIterator3d in_back_bottom_right,
1253  OutputBidirectionalIterator3d out_front_upper_left)
1254  {
1255  //recuperer la taille du container
1256  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
1257 
1258  std::size_t nb_sli = size3d[0];
1259  std::size_t nb_lig = size3d[1];
1260  std::size_t nb_col = size3d[2];
1261 
1262  //ifft3d
1263  for(size_t k = 0; k < nb_sli; ++k)
1264  for(size_t i = 0; i < nb_lig; ++i)
1265  {
1266  slip::split_radix_ifft(in_front_upper_left.row_begin(k,i),in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
1267  }
1268 
1269  slip::Vector<std::complex<double> > Vtmp(nb_lig);
1270  for(size_t k = 0; k < nb_sli; ++k)
1271  for(size_t j = 0; j < nb_col; ++j)
1272  {
1273  slip::split_radix_ifft(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
1274  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
1275  }
1276 
1277  Vtmp.resize(nb_sli);
1278  for(size_t i = 0; i < nb_lig; ++i)
1279  for(size_t j = 0; j < nb_col; ++j)
1280  {
1281  slip::split_radix_ifft(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
1282  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
1283  }
1284  }
1285 
1286 
1300 template<typename InputBidirectionalIterator3d,
1301 typename OutputBidirectionalIterator3d>
1302 inline
1303 void complex_split_radix_fft3d(InputBidirectionalIterator3d in_front_upper_left,
1304  InputBidirectionalIterator3d in_back_bottom_right,
1305  OutputBidirectionalIterator3d out_front_upper_left)
1306 {
1307  //recuperer la taille du container
1308  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
1309  typedef typename std::iterator_traits<InputBidirectionalIterator3d>::value_type complex;
1310  typedef typename complex::value_type real;
1311 
1312 
1313  std::size_t nb_sli = size3d[0];
1314  std::size_t nb_lig = size3d[1];
1315  std::size_t nb_col = size3d[2];
1316  //fft3d
1317  for(size_t k = 0; k < nb_sli; ++k)
1318  for(size_t i = 0; i < nb_lig; ++i)
1319  {
1320  slip::complex_split_radix_fft(in_front_upper_left.row_begin(k,i),
1321  in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
1322  }
1323 
1324  slip::Vector<std::complex<double> > Vtmp(nb_lig);
1325  for(size_t k = 0; k < nb_sli; ++k)
1326  for(size_t j = 0; j < nb_col; ++j)
1327  {
1328  slip::complex_split_radix_fft(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
1329  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
1330  }
1331 
1332  Vtmp.resize(nb_sli);
1333  for(size_t i = 0; i < nb_lig; ++i)
1334  for(size_t j = 0; j < nb_col; ++j)
1335  {
1336  slip::complex_split_radix_fft(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
1337  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
1338  }
1339 }
1340 
1341 
1356 template<typename InputMatrix3d,
1357 typename OutputMatrix3d>
1358 inline
1359 void complex_split_radix_fft3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
1360 {
1361  assert(datain.cols() == dataout.cols());
1362  assert(datain.rows() == dataout.rows());
1363  assert(datain.slices() == dataout.slices());
1364  slip::complex_split_radix_fft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
1365 }
1366 
1401  template<typename InputMatrix3d,
1402  typename OutputMatrix3d>
1403  inline
1404  void real_split_radix_fft3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
1405  {
1406  assert(datain.cols() == dataout.cols());
1407  assert(datain.rows() == dataout.rows());
1408  assert(datain.slices() == dataout.slices());
1409  slip::real_split_radix_fft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
1410  }
1411 
1412 
1413 
1414 
1415 
1416 
1452  template<typename InputMatrix3d,
1453  typename OutputMatrix3d>
1454  inline
1455  void split_radix_ifft3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
1456  {
1457  assert(datain.cols() == dataout.cols());
1458  assert(datain.rows() == dataout.rows());
1459  assert(datain.slices() == dataout.slices());
1460  slip::split_radix_ifft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
1461  }
1462 
1463 
1464  //Discrete Cosinus Transform
1465 
1490  template<class InputIter,class OutputIter>
1491  inline
1492  void split_radix_idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
1493  {
1494  typedef typename std::iterator_traits<InputIter>::value_type real;
1495  size_t length = in_end - in_begin;
1496  assert(length > 1);
1497 
1498  // data copy
1499  std::copy(in_begin,in_end,out_begin);
1500 
1501  //DCT initialization
1502  slip::Array<real> r(8*length+15,0.0);
1503  slip::Array<size_t> ifac(8*length+15,(size_t)0);
1504  slip::cosqi(&length,r.begin(),ifac.begin());
1505 
1506  //DCT
1507  slip::cosqf(&length,out_begin,r.begin(),ifac.begin());
1508  std::transform(out_begin, out_begin+length, out_begin,std::bind2nd(std::divides<real>(),std::sqrt(2*length)));
1509  }
1510 
1534  template<class InputIter,class OutputIter>
1535  inline
1536  void split_radix_dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
1537  {
1538  typedef typename std::iterator_traits<InputIter>::value_type real;
1539  size_t length = in_end - in_begin;
1540  assert(length > 1);
1541 
1542  // data copy
1543  std::copy(in_begin,in_end,out_begin);
1544 
1545  // DCT initialization
1546  slip::Array<real> r(5*length,0.0);
1547  slip::Array<size_t> ifac(5*length,(size_t)0);
1548  slip::cosqi(&length,r.begin(),ifac.begin());
1549 
1550  // DCT
1551  slip::cosqb(&length,out_begin,r.begin(),ifac.begin());
1552  std::transform(out_begin, out_begin+length, out_begin,std::bind2nd(std::divides<real>(),2*std::sqrt(2*length)));
1553  }
1554 
1557 #ifdef HAVE_FFTW
1558 
1560 
1565 
1580  template<class ForwardIterator1, class ForwardIterator2>
1581  inline
1582  void fftduplicate(ForwardIterator1 begin1, ForwardIterator2 begin2,
1583  ForwardIterator2 end2)
1584  {
1585  typedef typename std::iterator_traits<ForwardIterator1>::value_type complex;
1586  complex J(0,1);
1587  begin1++;
1588  begin2++;
1589  end2--;
1590  while (begin2 < end2)
1591  {
1592  *end2 = std::real(*begin1) - J*std::imag(*begin1);
1593  begin1++;
1594  begin2++;
1595  end2--;
1596  }
1597  }
1598 
1611  template<class ForwardIterator2D>
1612  inline
1613  void fftduplicate2d(ForwardIterator2D upper_left,
1614  ForwardIterator2D bottom_right)
1615  {
1616  typename ForwardIterator2D::difference_type size2d = bottom_right - upper_left;
1617  std::size_t nb_lig = size2d[0];
1618  slip::fftduplicate(upper_left.row_begin(0),upper_left.row_begin(0),upper_left.row_end(0));
1619  for(size_t i = 1; i < nb_lig; ++i)
1620  {
1621  slip::fftduplicate(upper_left.row_begin(i),upper_left.row_begin(nb_lig - i),upper_left.row_end(nb_lig - i));
1622  }
1623  }
1624 
1639  template<class InputIter,class OutputIter>
1640  inline
1641  void real_fftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
1642  {
1643  // typedef typename std::iterator_traits<InputIter>::value_type real;
1644  typedef typename std::iterator_traits<OutputIter>::value_type complex;
1645  typedef typename complex::value_type real;
1646  complex J(0,1);
1647  size_t length = in_end - in_begin;
1648  size_t new_length = length/2 + 1;
1649  assert(length > 1);
1650  slip::Array<double> in(length);
1651  slip::Array<std::complex<double> > out(new_length);
1652  fftw_plan plan;
1653 
1654  std::copy(in_begin,in_end,in.begin());
1655  //With the FFTW_ESTIMATE flag, the input/output arrays are not overwritten during planning. !!! Default flag
1656  //gave some unexplained segmentation faults with the odd dimensions greater than 40 !!!
1657  plan = fftw_plan_dft_r2c_1d(int(length),(double*)&(*in.begin()),(fftw_complex*)&(*out.begin()),FFTW_ESTIMATE);
1658  fftw_execute(plan);
1659 
1660  slip::Array<std::complex<double> >::iterator it = out.begin();
1661  OutputIter itout = out_begin;
1662  while(it < out.end())
1663  {
1664  *itout = (real)std::real(*it) + J * (real)std::imag(*it);
1665  itout++;
1666  it++;
1667  }
1668 
1669  slip::fftduplicate(out_begin,out_begin,out_begin+length);
1670 
1671  fftw_destroy_plan(plan);
1672  }
1673 
1688  template<class InputIter,class OutputIter>
1689  inline
1690  void complex_fftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
1691  {
1692  typedef typename std::iterator_traits<OutputIter>::value_type complex;
1693  typedef typename complex::value_type real;
1694  complex J(0,1);
1695  // typedef typename complex::value_type real;
1696  size_t length = in_end - in_begin;
1697  assert(length > 1);
1698 
1699  slip::Array<std::complex<double> > in(length);
1700  slip::Array<std::complex<double> > out(length);
1701  fftw_plan plan;
1702 
1703  std::copy(in_begin,in_end,in.begin());
1704  plan = fftw_plan_dft_1d(int(length),(fftw_complex*)&(*in.begin()),(fftw_complex*)&(*out.begin()),-1,FFTW_ESTIMATE);
1705  fftw_execute(plan);
1706  slip::Array<std::complex<double> >::iterator it = out.begin();
1707  OutputIter itout = out_begin;
1708  while(it < out.end())
1709  {
1710  *itout = ((real)std::real(*it) + J * (real)std::imag(*it));
1711  itout++;
1712  it++;
1713  }
1714 
1715  fftw_destroy_plan(plan);
1716  }
1717 
1731  template<class InputIter,class OutputIter>
1732  inline
1733  void ifftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
1734  {
1735  typedef typename std::iterator_traits<OutputIter>::value_type complex;
1736  typedef typename complex::value_type real;
1737  complex J(0,1);
1738  size_t length = in_end - in_begin;
1739  assert(length > 1);
1740 
1741  slip::Array<std::complex<double> > in(length);
1742  slip::Array<std::complex<double> > out(length);
1743  fftw_plan plan;
1744 
1745  std::copy(in_begin,in_end,in.begin());
1746  plan = fftw_plan_dft_1d(int(length),(fftw_complex*)&(*in.begin()),(fftw_complex*)&(*out.begin()),1,FFTW_ESTIMATE);
1747  fftw_execute(plan);
1748  slip::Array<std::complex<double> >::iterator it = out.begin();
1749  OutputIter itout = out_begin;
1750  while(it < out.end())
1751  {
1752  *itout = ((real)std::real(*it) + J * (real)std::imag(*it)) / (real)length;
1753  itout++;
1754  it++;
1755  }
1756 
1757  fftw_destroy_plan(plan);
1758  }
1759 
1775  template<typename InputBidirectionalIterator2d,
1776  typename OutputBidirectionalIterator2d>
1777  inline
1778  void real_fftw2d(InputBidirectionalIterator2d in_upper_left,
1779  InputBidirectionalIterator2d in_bottom_right,
1780  OutputBidirectionalIterator2d out_upper_left)
1781  {
1782  typedef typename std::iterator_traits<OutputBidirectionalIterator2d>::value_type complex;
1783  typedef typename complex::value_type real;
1784  complex J(0,1);
1785 
1786  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
1787  std::size_t nb_lig = size2d[0];
1788  std::size_t nb_col = size2d[1];
1789  int new_nb_col = (nb_col/2)+1;
1790  assert((nb_lig > 1)&&(nb_col>1));
1791 
1792  slip::Array2d<double> in(nb_lig,nb_col);
1793  slip::Array2d<std::complex<double> > out(nb_lig,new_nb_col);
1794  fftw_plan plan;
1795 
1796  std::copy(in_upper_left,in_bottom_right,in.upper_left());
1797  //With the FFTW_ESTIMATE flag, the input/output arrays are not overwritten during planning. !!! Default flag
1798  //gave some unexplained segmentation faults with the odd dimensions greater than 40 !!!
1799  plan = fftw_plan_dft_r2c_2d(int(nb_lig),int(nb_col),(double*)&(*in.upper_left()),(fftw_complex*)&(*out.upper_left()),FFTW_ESTIMATE);
1800  fftw_execute(plan);
1801 
1802  for(int i=0; i<int(nb_lig);i++)
1803  for(int j=0; j<new_nb_col;j++)
1804  {
1805  slip::DPoint2d<int> dp(i,j);
1806  *(out_upper_left + dp) = (real)std::real(out[i][j]) + J * (real)std::imag(out[i][j]);
1807  }
1808  slip::fftduplicate2d(out_upper_left, out_upper_left+size2d);
1809 
1810  fftw_destroy_plan(plan);
1811  }
1812 
1828  template<typename InputBidirectionalIterator2d,
1829  typename OutputBidirectionalIterator2d>
1830  inline
1831  void complex_fftw2d(InputBidirectionalIterator2d in_upper_left,
1832  InputBidirectionalIterator2d in_bottom_right,
1833  OutputBidirectionalIterator2d out_upper_left)
1834  {
1835  typedef typename std::iterator_traits<OutputBidirectionalIterator2d>::value_type complex;
1836  typedef typename complex::value_type real;
1837  complex J(0,1);
1838 
1839  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
1840  std::size_t nb_lig = size2d[0];
1841  std::size_t nb_col = size2d[1];
1842  assert((nb_lig > 1)&&(nb_col>1));
1843 
1844  slip::Array2d<std::complex<double> > in(nb_lig,nb_col);
1845  slip::Array2d<std::complex<double> > out(nb_lig,nb_col);
1846  fftw_plan plan;
1847 
1848  std::copy(in_upper_left,in_bottom_right,in.upper_left());
1849  plan = fftw_plan_dft_2d(int(nb_lig),int(nb_col),(fftw_complex*)&(*in.upper_left()),(fftw_complex*)&(*out.upper_left()),-1,FFTW_ESTIMATE);
1850  fftw_execute(plan);
1851  slip::Array2d<std::complex<double> >::iterator it = out.begin();
1852  OutputBidirectionalIterator2d itout = out_upper_left;
1853  while(it < out.end())
1854  {
1855  *itout = ((real)std::real(*it) + J * (real)std::imag(*it));
1856  itout++;
1857  it++;
1858  }
1859 
1860  fftw_destroy_plan(plan);
1861  }
1862 
1878  template<typename InputBidirectionalIterator2d,
1879  typename OutputBidirectionalIterator2d>
1880  inline
1881  void ifftw2d(InputBidirectionalIterator2d in_upper_left,
1882  InputBidirectionalIterator2d in_bottom_right,
1883  OutputBidirectionalIterator2d out_upper_left)
1884  {
1885  typedef typename std::iterator_traits<OutputBidirectionalIterator2d>::value_type complex;
1886  typedef typename complex::value_type real;
1887  complex J(0,1);
1888 
1889  typename InputBidirectionalIterator2d::difference_type size2d = in_bottom_right - in_upper_left;
1890  std::size_t nb_lig = size2d[0];
1891  std::size_t nb_col = size2d[1];
1892  int length(nb_lig*nb_col);
1893  assert((nb_lig > 1)&&(nb_col>1));
1894 
1895  slip::Array2d<std::complex<double> > in(nb_lig,nb_col);
1896  slip::Array2d<std::complex<double> > out(nb_lig,nb_col);
1897  fftw_plan plan;
1898 
1899  std::copy(in_upper_left,in_bottom_right,in.upper_left());
1900  plan = fftw_plan_dft_2d(int(nb_lig),int(nb_col),(fftw_complex*)&(*in.upper_left()),(fftw_complex*)&(*out.upper_left()),1,FFTW_ESTIMATE);
1901  fftw_execute(plan);
1902 
1903  slip::Array2d<std::complex<double> >::iterator it = out.begin();
1904  OutputBidirectionalIterator2d itout = out_upper_left;
1905  while(it < out.end())
1906  {
1907  *itout = ((real)std::real(*it) + J * (real)std::imag(*it)) / (real)length;
1908  itout++;
1909  it++;
1910  }
1911 
1912  fftw_destroy_plan(plan);
1913  }
1914 
1932  template<typename InputMatrix, typename OutputMatrix>
1933  inline void real_fftw2d(InputMatrix &datain, OutputMatrix &dataout)
1934  {
1935  assert(datain.cols() == dataout.cols());
1936  assert(datain.rows() == dataout.rows());
1937  slip::real_fftw2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
1938  }
1939 
1955  template<typename InputMatrix, typename OutputMatrix>
1956  inline
1957  void ifftw2d(InputMatrix &datain, OutputMatrix &dataout)
1958  {
1959  assert(datain.cols() == dataout.cols());
1960  assert(datain.rows() == dataout.rows());
1961  slip::ifftw2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
1962  }
1963 
1979  template<typename InputBidirectionalIterator3d,
1980  typename OutputBidirectionalIterator3d>
1981  inline
1982  void real_fftw3d(InputBidirectionalIterator3d in_front_upper_left,
1983  InputBidirectionalIterator3d in_back_bottom_right,
1984  OutputBidirectionalIterator3d out_front_upper_left)
1985  {
1986  //recuperer la taille du container
1987  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
1988  typedef typename std::iterator_traits<InputBidirectionalIterator3d>::value_type real;
1989  typedef typename std::iterator_traits<OutputBidirectionalIterator3d>::value_type complex;
1990 
1991  std::size_t nb_sli = size3d[0];
1992  std::size_t nb_lig = size3d[1];
1993  std::size_t nb_col = size3d[2];
1994  assert((nb_sli > 1)&&(nb_lig > 1)&&(nb_col>1));
1995  //fft3d
1996  for(size_t k = 0; k < nb_sli; ++k)
1997  for(size_t i = 0; i < nb_lig; ++i)
1998  {
1999  slip::real_fftw(in_front_upper_left.row_begin(k,i),in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
2000  }
2001 
2002  slip::Vector<std::complex<double> > Vtmp(nb_lig);
2003  for(size_t k = 0; k < nb_sli; ++k)
2004  for(size_t j = 0; j < nb_col; ++j)
2005  {
2006  slip::complex_fftw(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
2007  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
2008  }
2009 
2010  Vtmp.resize(nb_sli);
2011  for(size_t i = 0; i < nb_lig; ++i)
2012  for(size_t j = 0; j < nb_col; ++j)
2013  {
2014  slip::complex_fftw(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
2015  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
2016  }
2017  }
2018 
2034  template<typename InputBidirectionalIterator3d,
2035  typename OutputBidirectionalIterator3d>
2036  inline
2037  void complex_fftw3d(InputBidirectionalIterator3d in_front_upper_left,
2038  InputBidirectionalIterator3d in_back_bottom_right,
2039  OutputBidirectionalIterator3d out_front_upper_left)
2040  {
2041  //recuperer la taille du container
2042  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
2043  typedef typename std::iterator_traits<InputBidirectionalIterator3d>::value_type complex;
2044  typedef typename complex::value_type real;
2045 
2046  std::size_t nb_sli = size3d[0];
2047  std::size_t nb_lig = size3d[1];
2048  std::size_t nb_col = size3d[2];
2049  assert((nb_sli > 1)&&(nb_lig > 1)&&(nb_col>1));
2050 
2051  //fft3d
2052  for(size_t k = 0; k < nb_sli; ++k)
2053  for(size_t i = 0; i < nb_lig; ++i)
2054  {
2055  slip::complex_fftw(in_front_upper_left.row_begin(k,i),in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
2056  }
2057 
2058  slip::Vector<std::complex<double> > Vtmp(nb_lig);
2059  for(size_t k = 0; k < nb_sli; ++k)
2060  for(size_t j = 0; j < nb_col; ++j)
2061  {
2062  slip::complex_fftw(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
2063  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
2064  }
2065 
2066  Vtmp.resize(nb_sli);
2067  for(size_t i = 0; i < nb_lig; ++i)
2068  for(size_t j = 0; j < nb_col; ++j)
2069  {
2070  slip::complex_fftw(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
2071  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
2072  }
2073  }
2110 template<typename InputBidirectionalIterator3d,
2111 typename OutputBidirectionalIterator3d>
2112 inline
2113 void complex_fft3d(InputBidirectionalIterator3d in_front_upper_left,
2114  InputBidirectionalIterator3d in_back_bottom_right,
2115  OutputBidirectionalIterator3d out_front_upper_left)
2116 {
2117 #ifdef HAVE_FFTW
2118  slip::complex_fftw3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2119 #else
2120  slip::complex_split_radix_fft3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2121 #endif
2122 }
2123 
2124 
2162 template<typename Volume1, typename Volume2>
2163 inline
2164 void complex_fft3d(Volume1 &datain, Volume2 &dataout)
2165 {
2166  assert(datain.cols() == dataout.cols());
2167  assert(datain.rows() == dataout.rows());
2168  assert(datain.slices() == dataout.slices());
2169 #ifdef HAVE_FFTW
2170  slip::complex_fftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2171 #else
2172  slip::complex_split_radix_fft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2173 #endif
2174 }
2175 
2191  template<typename InputBidirectionalIterator3d,
2192  typename OutputBidirectionalIterator3d>
2193  inline
2194  void ifftw3d(InputBidirectionalIterator3d in_front_upper_left,
2195  InputBidirectionalIterator3d in_back_bottom_right,
2196  OutputBidirectionalIterator3d out_front_upper_left)
2197  {
2198  //recuperer la taille du container
2199  typename InputBidirectionalIterator3d::difference_type size3d = in_back_bottom_right - in_front_upper_left;
2200 
2201  std::size_t nb_sli = size3d[0];
2202  std::size_t nb_lig = size3d[1];
2203  std::size_t nb_col = size3d[2];
2204 
2205  //ifft3d
2206  for(size_t k = 0; k < nb_sli; ++k)
2207  for(size_t i = 0; i < nb_lig; ++i)
2208  {
2209  slip::ifftw(in_front_upper_left.row_begin(k,i),in_front_upper_left.row_end(k,i),out_front_upper_left.row_begin(k,i));
2210  }
2211 
2212  slip::Vector<std::complex<double> > Vtmp(nb_lig);
2213  for(size_t k = 0; k < nb_sli; ++k)
2214  for(size_t j = 0; j < nb_col; ++j)
2215  {
2216  slip::ifftw(out_front_upper_left.col_begin(k,j),out_front_upper_left.col_end(k,j),Vtmp.begin());
2217  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.col_begin(k,j));
2218  }
2219 
2220  Vtmp.resize(nb_sli);
2221  for(size_t i = 0; i < nb_lig; ++i)
2222  for(size_t j = 0; j < nb_col; ++j)
2223  {
2224  slip::ifftw(out_front_upper_left.slice_begin(i,j),out_front_upper_left.slice_end(i,j),Vtmp.begin());
2225  std::copy(Vtmp.begin(),Vtmp.end(),out_front_upper_left.slice_begin(i,j));
2226  }
2227  }
2228 
2245  template<typename InputMatrix3d,
2246  typename OutputMatrix3d>
2247  inline
2248  void real_fftw3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
2249  {
2250  assert(datain.cols() == dataout.cols());
2251  assert(datain.rows() == dataout.rows());
2252  assert(datain.slices() == dataout.slices());
2253  slip::real_fftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2254  }
2255 
2272  template<typename InputMatrix3d,
2273  typename OutputMatrix3d>
2274  inline
2275  void complex_fftw3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
2276  {
2277  assert(datain.cols() == dataout.cols());
2278  assert(datain.rows() == dataout.rows());
2279  assert(datain.slices() == dataout.slices());
2280  slip::complex_fftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2281  }
2282 
2301  template<typename InputMatrix3d,
2302  typename OutputMatrix3d>
2303  inline
2304  void ifftw3d(InputMatrix3d & datain,OutputMatrix3d & dataout)
2305  {
2306  assert(datain.cols() == dataout.cols());
2307  assert(datain.rows() == dataout.rows());
2308  assert(datain.slices() == dataout.slices());
2309  slip::ifftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2310  }
2311 
2325  template<class InputIter,class OutputIter>
2326  inline
2327  void fftw_dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2328  {
2329  size_t length = in_end - in_begin;
2330  assert(length > 1);
2331 
2332  double *in = (double*) fftw_malloc(sizeof(double) * int(length));
2333  double *out = (double*) fftw_malloc(sizeof(double) * int(length));
2334  fftw_plan plan;
2335  fftw_r2r_kind kind = FFTW_REDFT10;
2336 
2337  std::copy(in_begin,in_end,in);
2338  plan = fftw_plan_r2r_1d(int(length),in,out,kind,FFTW_ESTIMATE);
2339  fftw_execute(plan);
2340  std::transform(out,out+length,out_begin,std::bind2nd(std::divides<double>(),std::sqrt(2*length)));
2341 
2342  fftw_destroy_plan(plan);
2343  fftw_free(in);
2344  fftw_free(out);
2345  }
2346 
2360  template<class InputIter,class OutputIter>
2361  inline
2362  void fftw_idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2363  {
2364  size_t length = in_end - in_begin;
2365  assert(length > 1);
2366 
2367  double *in = (double*) fftw_malloc(sizeof(double) * int(length));
2368  double *out = (double*) fftw_malloc(sizeof(double) * int(length));
2369  fftw_plan plan;
2370  fftw_r2r_kind kind = FFTW_REDFT01;
2371 
2372  std::copy(in_begin,in_end,in);
2373  plan = fftw_plan_r2r_1d(int(length),in,out,kind,FFTW_ESTIMATE);
2374  fftw_execute(plan);
2375  std::transform(out,out+length,out_begin,std::bind2nd(std::divides<double>(),std::sqrt(2*length)));
2376 
2377  fftw_destroy_plan(plan);
2378  fftw_free(in);
2379  fftw_free(out);
2380  }
2381 
2384 #endif //HAVE_FFTW
2385 
2387 
2392 
2428  template<class InputIter, class OutputIter>
2429  inline
2430  void real_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2431  {
2432 #ifdef HAVE_FFTW
2433  slip::real_fftw(in_begin,in_end,out_begin);
2434 #else
2435  slip::real_split_radix_fft(in_begin,in_end,out_begin);
2436 #endif
2437  }
2438 
2439 
2468  template<class InputIter, class OutputIter>
2469  inline
2470  void complex_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2471  {
2472 #ifdef HAVE_FFTW
2473  slip::complex_fftw(in_begin,in_end,out_begin);
2474 #else
2475  slip::complex_split_radix_fft(in_begin,in_end,out_begin);
2476 #endif
2477 }
2478 
2505  template<class InputIter, class OutputIter>
2506  inline
2507  void ifft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2508  {
2509 #ifdef HAVE_FFTW
2510  slip::ifftw(in_begin,in_end,out_begin);
2511 #else
2512  slip::split_radix_ifft(in_begin,in_end,out_begin);
2513 #endif
2514  }
2515 
2550  template<typename InputBidirectionalIterator2d,
2551  typename OutputBidirectionalIterator2d>
2552  inline
2553  void real_fft2d(InputBidirectionalIterator2d in_upper_left,
2554  InputBidirectionalIterator2d in_bottom_right,
2555  OutputBidirectionalIterator2d out_upper_left)
2556  {
2557 #ifdef HAVE_FFTW
2558  slip::real_fftw2d(in_upper_left,in_bottom_right,out_upper_left);
2559 #else
2560  slip::real_split_radix_fft2d(in_upper_left,in_bottom_right,out_upper_left);
2561 #endif
2562  }
2563 
2582  template<typename InputBidirectionalIterator2d,
2583  typename OutputBidirectionalIterator2d>
2584  inline
2585  void complex_fft2d(InputBidirectionalIterator2d in_upper_left,
2586  InputBidirectionalIterator2d in_bottom_right,
2587  OutputBidirectionalIterator2d out_upper_left)
2588  {
2589 #ifdef HAVE_FFTW
2590  slip::complex_fftw2d(in_upper_left,in_bottom_right,out_upper_left);
2591 #else
2592  slip::complex_split_radix_fft2d(in_upper_left,in_bottom_right,out_upper_left);
2593 #endif
2594  }
2595 
2626  template<typename InputBidirectionalIterator2d,
2627  typename OutputBidirectionalIterator2d>
2628  inline
2629  void ifft2d(InputBidirectionalIterator2d in_upper_left,
2630  InputBidirectionalIterator2d in_bottom_right,
2631  OutputBidirectionalIterator2d out_upper_left)
2632  {
2633 #ifdef HAVE_FFTW
2634  slip::ifftw2d(in_upper_left,in_bottom_right,out_upper_left);
2635 #else
2636  slip::split_radix_ifft2d(in_upper_left,in_bottom_right,out_upper_left);
2637 #endif
2638  }
2639 
2674  template<typename Matrix1, typename Matrix2>
2675  inline
2676  void real_fft2d(Matrix1 &datain, Matrix2 &dataout)
2677  {
2678  assert(datain.cols() == dataout.cols());
2679  assert(datain.rows() == dataout.rows());
2680 #ifdef HAVE_FFTW
2681  slip::real_fftw2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2682 #else
2683  slip::real_split_radix_fft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2684 #endif
2685  }
2686 
2705  template<typename Matrix1, typename Matrix2>
2706  inline
2707  void complex_fft2d(Matrix1 &datain, Matrix2 &dataout)
2708  {
2709  assert(datain.cols() == dataout.cols());
2710  assert(datain.rows() == dataout.rows());
2711 #ifdef HAVE_FFTW
2712  slip::complex_fftw2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2713 #else
2714  slip::complex_split_radix_fft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2715 #endif
2716  }
2717 
2751  template<typename Matrix1, typename Matrix2>
2752  inline
2753  void ifft2d(Matrix1 &datain, Matrix2 &dataout)
2754  {
2755  assert(datain.cols() == dataout.cols());
2756  assert(datain.rows() == dataout.rows());
2757 #ifdef HAVE_FFTW
2758  slip::ifftw2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2759 #else
2760  slip::split_radix_ifft2d(datain.upper_left(),datain.bottom_right(),dataout.upper_left());
2761 #endif
2762  }
2763 
2801  template<typename InputBidirectionalIterator3d,
2802  typename OutputBidirectionalIterator3d>
2803  inline
2804  void real_fft3d(InputBidirectionalIterator3d in_front_upper_left,
2805  InputBidirectionalIterator3d in_back_bottom_right,
2806  OutputBidirectionalIterator3d out_front_upper_left)
2807  {
2808 #ifdef HAVE_FFTW
2809  slip::real_fftw3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2810 #else
2811  slip::real_split_radix_fft3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2812 #endif
2813  }
2814 
2849  template<typename InputBidirectionalIterator3d,
2850  typename OutputBidirectionalIterator3d>
2851  inline
2852  void ifft3d(InputBidirectionalIterator3d in_front_upper_left,
2853  InputBidirectionalIterator3d in_back_bottom_right,
2854  OutputBidirectionalIterator3d out_front_upper_left)
2855  {
2856 #ifdef HAVE_FFTW
2857  slip::ifftw3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2858 #else
2859  slip::split_radix_ifft3d(in_front_upper_left,in_back_bottom_right,out_front_upper_left);
2860 #endif
2861  }
2862 
2901  template<typename Volume1, typename Volume2>
2902  inline
2903  void real_fft3d(Volume1 &datain, Volume2 &dataout)
2904  {
2905  assert(datain.cols() == dataout.cols());
2906  assert(datain.rows() == dataout.rows());
2907  assert(datain.slices() == dataout.slices());
2908 #ifdef HAVE_FFTW
2909  slip::real_fftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2910 #else
2911  slip::real_split_radix_fft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2912 #endif
2913  }
2914 
2953  template<typename Volume1, typename Volume2>
2954  inline
2955  void ifft3d(Volume1 &datain, Volume2 &dataout)
2956  {
2957  assert(datain.cols() == dataout.cols());
2958  assert(datain.rows() == dataout.rows());
2959  assert(datain.slices() == dataout.slices());
2960 #ifdef HAVE_FFTW
2961  slip::ifftw3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2962 #else
2963  slip::split_radix_ifft3d(datain.front_upper_left(),datain.back_bottom_right(),dataout.front_upper_left());
2964 #endif
2965  }
2966 
2991  template<class InputIter,class OutputIter>
2992  inline
2993  void dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
2994  {
2995 #ifdef HAVE_FFTW
2996  slip::fftw_dct(in_begin,in_end,out_begin);
2997 #else
2998  slip::split_radix_dct(in_begin,in_end,out_begin);
2999 #endif
3000  }
3001 
3027  template<class InputIter,class OutputIter>
3028  inline
3029  void idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
3030  {
3031 #ifdef HAVE_FFTW
3032  slip::fftw_idct(in_begin,in_end,out_begin);
3033 #else
3034  slip::split_radix_idct(in_begin,in_end,out_begin);
3035 #endif
3036  }
3037 
3038 
3043 
3056 template<class ForwardIterator4D>
3057 inline
3058 void fftshift4d(ForwardIterator4D first_front_upper_left,
3059  ForwardIterator4D last_back_bottom_right)
3060 {
3061  typename ForwardIterator4D::difference_type size4d = last_back_bottom_right - first_front_upper_left;
3062 
3063  std::size_t nb_sla = size4d[0];
3064  std::size_t nb_sli = size4d[1];
3065  std::size_t nb_lig = size4d[2];
3066  std::size_t nb_col = size4d[3];
3067 
3068  for(size_t t = 0; t < nb_sla; ++t)
3069  for(size_t k = 0; k < nb_sli; ++k)
3070  for(size_t i = 0; i < nb_lig; ++i)
3071  {
3072  slip::fftshift(first_front_upper_left.row_begin(t,k,i),first_front_upper_left.row_end(t,k,i));
3073  }
3074 
3075  for(size_t t = 0; t < nb_sla; ++t)
3076  for(size_t k = 0; k < nb_sli; ++k)
3077  for(size_t j = 0; j < nb_col; ++j)
3078  {
3079  slip::fftshift(first_front_upper_left.col_begin(t,k,j),first_front_upper_left.col_end(t,k,j));
3080  }
3081  for(size_t t = 0; t < nb_sla; ++t)
3082  for(size_t i = 0; i < nb_lig; ++i)
3083  for(size_t j = 0; j < nb_col; ++j)
3084  {
3085  slip::fftshift(first_front_upper_left.slice_begin(t,i,j),first_front_upper_left.slice_end(t,i,j));
3086  }
3087  for(size_t k = 0; k < nb_sli; ++k)
3088  for(size_t i = 0; i < nb_lig; ++i)
3089  for(size_t j = 0; j < nb_col; ++j)
3090  {
3091  slip::fftshift(first_front_upper_left.slab_begin(k,i,j),first_front_upper_left.slab_end(k,i,j));
3092  }
3093 }
3094 
3101 
3102 
3134 template<typename InputBidirectionalIterator4d,
3135 typename OutputBidirectionalIterator4d>
3136 inline
3137 void real_split_radix_fft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3138  InputBidirectionalIterator4d in_last_back_bottom_right,
3139  OutputBidirectionalIterator4d out_first_front_upper_left)
3140 {
3141  //recuperer la taille du container
3142  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3143  typedef typename std::iterator_traits<InputBidirectionalIterator4d>::value_type real;
3144  typedef typename std::iterator_traits<OutputBidirectionalIterator4d>::value_type complex;
3145 
3146  std::size_t nb_sla = size4d[0];
3147  std::size_t nb_sli = size4d[1];
3148  std::size_t nb_lig = size4d[2];
3149  std::size_t nb_col = size4d[3];
3150  //fft4d
3151  for(size_t t = 0; t < nb_sla; ++t)
3152  for(size_t k = 0; k < nb_sli; ++k)
3153  for(size_t i = 0; i < nb_lig; ++i)
3154  {
3155  slip::real_split_radix_fft(in_first_front_upper_left.row_begin(t,k,i),
3156  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3157  }
3158 
3159  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3160  for(size_t t = 0; t < nb_sla; ++t)
3161  for(size_t k = 0; k < nb_sli; ++k)
3162  for(size_t j = 0; j < nb_col; ++j)
3163  {
3164  slip::complex_split_radix_fft(out_first_front_upper_left.col_begin(t,k,j),
3165  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3166  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3167  }
3168 
3169  Vtmp.resize(nb_sli);
3170  for(size_t t = 0; t < nb_sla; ++t)
3171  for(size_t i = 0; i < nb_lig; ++i)
3172  for(size_t j = 0; j < nb_col; ++j)
3173  {
3174  slip::complex_split_radix_fft(out_first_front_upper_left.slice_begin(t,i,j),
3175  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3176  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3177  }
3178 
3179  Vtmp.resize(nb_sla);
3180  for(size_t k = 0; k < nb_sli; ++k)
3181  for(size_t i = 0; i < nb_lig; ++i)
3182  for(size_t j = 0; j < nb_col; ++j)
3183  {
3184  slip::complex_split_radix_fft(out_first_front_upper_left.slab_begin(k,i,j),
3185  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3186  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3187  }
3188 }
3189 
3203 template<typename InputBidirectionalIterator4d,
3204 typename OutputBidirectionalIterator4d>
3205 inline
3206 void complex_split_radix_fft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3207  InputBidirectionalIterator4d in_last_back_bottom_right,
3208  OutputBidirectionalIterator4d out_first_front_upper_left)
3209 {
3210  //recuperer la taille du container
3211  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3212  typedef typename std::iterator_traits<InputBidirectionalIterator4d>::value_type complex;
3213 
3214  std::size_t nb_sla = size4d[0];
3215  std::size_t nb_sli = size4d[1];
3216  std::size_t nb_lig = size4d[2];
3217  std::size_t nb_col = size4d[3];
3218  //fft4d
3219  for(size_t t = 0; t < nb_sla; ++t)
3220  for(size_t k = 0; k < nb_sli; ++k)
3221  for(size_t i = 0; i < nb_lig; ++i)
3222  {
3223  slip::complex_split_radix_fft(in_first_front_upper_left.row_begin(t,k,i),
3224  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3225  }
3226 
3227  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3228  for(size_t t = 0; t < nb_sla; ++t)
3229  for(size_t k = 0; k < nb_sli; ++k)
3230  for(size_t j = 0; j < nb_col; ++j)
3231  {
3232  slip::complex_split_radix_fft(out_first_front_upper_left.col_begin(t,k,j),
3233  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3234  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3235  }
3236 
3237  Vtmp.resize(nb_sli);
3238  for(size_t t = 0; t < nb_sla; ++t)
3239  for(size_t i = 0; i < nb_lig; ++i)
3240  for(size_t j = 0; j < nb_col; ++j)
3241  {
3242  slip::complex_split_radix_fft(out_first_front_upper_left.slice_begin(t,i,j),
3243  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3244  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3245  }
3246 
3247  Vtmp.resize(nb_sla);
3248  for(size_t k = 0; k < nb_sli; ++k)
3249  for(size_t i = 0; i < nb_lig; ++i)
3250  for(size_t j = 0; j < nb_col; ++j)
3251  {
3252  slip::complex_split_radix_fft(out_first_front_upper_left.slab_begin(k,i,j),
3253  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3254  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3255  }
3256 }
3257 
3290 template<typename InputBidirectionalIterator4d,
3291 typename OutputBidirectionalIterator4d>
3292 inline
3293 void split_radix_ifft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3294  InputBidirectionalIterator4d in_last_back_bottom_right,
3295  OutputBidirectionalIterator4d out_first_front_upper_left)
3296 {
3297  //recuperer la taille du container
3298  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3299 
3300  std::size_t nb_sla = size4d[0];
3301  std::size_t nb_sli = size4d[1];
3302  std::size_t nb_lig = size4d[2];
3303  std::size_t nb_col = size4d[3];
3304 
3305  //ifft4d
3306  for(size_t t = 0; t < nb_sla; ++t)
3307  for(size_t k = 0; k < nb_sli; ++k)
3308  for(size_t i = 0; i < nb_lig; ++i)
3309  {
3310  slip::split_radix_ifft(in_first_front_upper_left.row_begin(t,k,i),
3311  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3312  }
3313 
3314  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3315  for(size_t t = 0; t < nb_sla; ++t)
3316  for(size_t k = 0; k < nb_sli; ++k)
3317  for(size_t j = 0; j < nb_col; ++j)
3318  {
3319  slip::split_radix_ifft(out_first_front_upper_left.col_begin(t,k,j),
3320  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3321  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3322  }
3323 
3324  Vtmp.resize(nb_sli);
3325  for(size_t t = 0; t < nb_sla; ++t)
3326  for(size_t i = 0; i < nb_lig; ++i)
3327  for(size_t j = 0; j < nb_col; ++j)
3328  {
3329  slip::split_radix_ifft(out_first_front_upper_left.slice_begin(t,i,j),
3330  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3331  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3332  }
3333 
3334  Vtmp.resize(nb_sla);
3335  for(size_t k = 0; k < nb_sli; ++k)
3336  for(size_t i = 0; i < nb_lig; ++i)
3337  for(size_t j = 0; j < nb_col; ++j)
3338  {
3339  slip::split_radix_ifft(out_first_front_upper_left.slab_begin(k,i,j),
3340  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3341  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3342  }
3343 }
3344 
3345 
3346 
3380 template<typename InputMatrix4d,
3381 typename OutputMatrix4d>
3382 inline
3383 void real_split_radix_fft4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3384 {
3385  assert(datain.cols() == dataout.cols());
3386  assert(datain.rows() == dataout.rows());
3387  assert(datain.slices() == dataout.slices());
3388  real_split_radix_fft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),
3389  dataout.first_front_upper_left());
3390 }
3391 
3406 template<typename InputMatrix4d,
3407 typename OutputMatrix4d>
3408 inline
3409 void complex_split_radix_fft4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3410 {
3411  assert(datain.cols() == dataout.cols());
3412  assert(datain.rows() == dataout.rows());
3413  assert(datain.slices() == dataout.slices());
3414  complex_split_radix_fft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),
3415  dataout.first_front_upper_left());
3416 }
3417 
3453 template<typename InputMatrix4d,
3454 typename OutputMatrix4d>
3455 inline
3456 void split_radix_ifft4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3457 {
3458  assert(datain.cols() == dataout.cols());
3459  assert(datain.rows() == dataout.rows());
3460  assert(datain.slices() == dataout.slices());
3461  split_radix_ifft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3462 }
3463 
3466 #ifdef HAVE_FFTW
3467 
3472 
3488 template<typename InputBidirectionalIterator4d,
3489 typename OutputBidirectionalIterator4d>
3490 inline
3491 void real_fftw4d(InputBidirectionalIterator4d in_first_front_upper_left,
3492  InputBidirectionalIterator4d in_last_back_bottom_right,
3493  OutputBidirectionalIterator4d out_first_front_upper_left)
3494 {
3495  //recuperer la taille du container
3496  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3497  typedef typename std::iterator_traits<InputBidirectionalIterator4d>::value_type real;
3498  typedef typename std::iterator_traits<OutputBidirectionalIterator4d>::value_type complex;
3499 
3500  std::size_t nb_sla = size4d[0];
3501  std::size_t nb_sli = size4d[1];
3502  std::size_t nb_lig = size4d[2];
3503  std::size_t nb_col = size4d[3];
3504  assert((nb_sla > 1)&&(nb_sli > 1)&&(nb_lig > 1)&&(nb_col>1));
3505  //fft4d
3506  for(size_t t = 0; t < nb_sla; ++t)
3507  for(size_t k = 0; k < nb_sli; ++k)
3508  for(size_t i = 0; i < nb_lig; ++i)
3509  {
3510  slip::real_fftw(in_first_front_upper_left.row_begin(t,k,i),
3511  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3512  }
3513 
3514  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3515  for(size_t t = 0; t < nb_sla; ++t)
3516  for(size_t k = 0; k < nb_sli; ++k)
3517  for(size_t j = 0; j < nb_col; ++j)
3518  {
3519  slip::complex_fftw(out_first_front_upper_left.col_begin(t,k,j),
3520  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3521  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3522  }
3523 
3524  Vtmp.resize(nb_sli);
3525  for(size_t t = 0; t < nb_sla; ++t)
3526  for(size_t i = 0; i < nb_lig; ++i)
3527  for(size_t j = 0; j < nb_col; ++j)
3528  {
3529  slip::complex_fftw(out_first_front_upper_left.slice_begin(t,i,j),
3530  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3531  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3532  }
3533 
3534  Vtmp.resize(nb_sla);
3535  for(size_t k = 0; k < nb_sli; ++k)
3536  for(size_t i = 0; i < nb_lig; ++i)
3537  for(size_t j = 0; j < nb_col; ++j)
3538  {
3539  slip::complex_fftw(out_first_front_upper_left.slab_begin(k,i,j),
3540  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3541  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3542  }
3543 }
3544 
3560 template<typename InputBidirectionalIterator4d,
3561 typename OutputBidirectionalIterator4d>
3562 inline
3563 void complex_fftw4d(InputBidirectionalIterator4d in_first_front_upper_left,
3564  InputBidirectionalIterator4d in_last_back_bottom_right,
3565  OutputBidirectionalIterator4d out_first_front_upper_left)
3566 {
3567  //recuperer la taille du container
3568  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3569  typedef typename std::iterator_traits<InputBidirectionalIterator4d>::value_type complex;
3570  typedef typename complex::value_type real;
3571 
3572  std::size_t nb_sla = size4d[0];
3573  std::size_t nb_sli = size4d[1];
3574  std::size_t nb_lig = size4d[2];
3575  std::size_t nb_col = size4d[3];
3576  assert((nb_sla > 1)&&(nb_sli > 1)&&(nb_lig > 1)&&(nb_col>1));
3577 
3578  //fft4d
3579  for(size_t t = 0; t < nb_sla; ++t)
3580  for(size_t k = 0; k < nb_sli; ++k)
3581  for(size_t i = 0; i < nb_lig; ++i)
3582  {
3583  slip::complex_fftw(in_first_front_upper_left.row_begin(t,k,i),
3584  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3585  }
3586 
3587  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3588  for(size_t t = 0; t < nb_sla; ++t)
3589  for(size_t k = 0; k < nb_sli; ++k)
3590  for(size_t j = 0; j < nb_col; ++j)
3591  {
3592  slip::complex_fftw(out_first_front_upper_left.col_begin(t,k,j),
3593  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3594  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3595  }
3596 
3597  Vtmp.resize(nb_sli);
3598  for(size_t t = 0; t < nb_sla; ++t)
3599  for(size_t i = 0; i < nb_lig; ++i)
3600  for(size_t j = 0; j < nb_col; ++j)
3601  {
3602  slip::complex_fftw(out_first_front_upper_left.slice_begin(t,i,j),
3603  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3604  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3605  }
3606 
3607  Vtmp.resize(nb_sla);
3608  for(size_t k = 0; k < nb_sli; ++k)
3609  for(size_t i = 0; i < nb_lig; ++i)
3610  for(size_t j = 0; j < nb_col; ++j)
3611  {
3612  slip::complex_fftw(out_first_front_upper_left.slab_begin(k,i,j),
3613  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3614  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3615  }
3616 }
3617 
3633 template<typename InputBidirectionalIterator4d,
3634 typename OutputBidirectionalIterator4d>
3635 inline
3636 void ifftw4d(InputBidirectionalIterator4d in_first_front_upper_left,
3637  InputBidirectionalIterator4d in_last_back_bottom_right,
3638  OutputBidirectionalIterator4d out_first_front_upper_left)
3639 {
3640  //recuperer la taille du container
3641  typename InputBidirectionalIterator4d::difference_type size4d = in_last_back_bottom_right - in_first_front_upper_left;
3642 
3643  std::size_t nb_sla = size4d[0];
3644  std::size_t nb_sli = size4d[1];
3645  std::size_t nb_lig = size4d[2];
3646  std::size_t nb_col = size4d[3];
3647 
3648  //ifft4d
3649  for(size_t t = 0; t < nb_sla; ++t)
3650  for(size_t k = 0; k < nb_sli; ++k)
3651  for(size_t i = 0; i < nb_lig; ++i)
3652  {
3653  slip::ifftw(in_first_front_upper_left.row_begin(t,k,i),
3654  in_first_front_upper_left.row_end(t,k,i),out_first_front_upper_left.row_begin(t,k,i));
3655  }
3656 
3657  slip::Vector<std::complex<double> > Vtmp(nb_lig);
3658  for(size_t t = 0; t < nb_sla; ++t)
3659  for(size_t k = 0; k < nb_sli; ++k)
3660  for(size_t j = 0; j < nb_col; ++j)
3661  {
3662  slip::ifftw(out_first_front_upper_left.col_begin(t,k,j),
3663  out_first_front_upper_left.col_end(t,k,j),Vtmp.begin());
3664  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.col_begin(t,k,j));
3665  }
3666 
3667  Vtmp.resize(nb_sli);
3668  for(size_t t = 0; t < nb_sla; ++t)
3669  for(size_t i = 0; i < nb_lig; ++i)
3670  for(size_t j = 0; j < nb_col; ++j)
3671  {
3672  slip::ifftw(out_first_front_upper_left.slice_begin(t,i,j),
3673  out_first_front_upper_left.slice_end(t,i,j),Vtmp.begin());
3674  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slice_begin(t,i,j));
3675  }
3676 
3677  Vtmp.resize(nb_sla);
3678  for(size_t k = 0; k < nb_sli; ++k)
3679  for(size_t i = 0; i < nb_lig; ++i)
3680  for(size_t j = 0; j < nb_col; ++j)
3681  {
3682  slip::ifftw(out_first_front_upper_left.slab_begin(k,i,j),
3683  out_first_front_upper_left.slab_end(k,i,j),Vtmp.begin());
3684  std::copy(Vtmp.begin(),Vtmp.end(),out_first_front_upper_left.slab_begin(k,i,j));
3685  }
3686 }
3687 
3703 template<typename InputMatrix4d,
3704 typename OutputMatrix4d>
3705 inline
3706 void real_fftw4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3707 {
3708  assert(datain.cols() == dataout.cols());
3709  assert(datain.rows() == dataout.rows());
3710  assert(datain.slices() == dataout.slices());
3711  assert(datain.slabs() == dataout.slabs());
3712  real_fftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3713 }
3714 
3731 template<typename InputMatrix4d,
3732 typename OutputMatrix4d>
3733 inline
3734 void complex_fftw4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3735 {
3736  assert(datain.cols() == dataout.cols());
3737  assert(datain.rows() == dataout.rows());
3738  assert(datain.slices() == dataout.slices());
3739  assert(datain.slabs() == dataout.slabs());
3740  complex_fftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3741 }
3742 
3762 template<typename InputMatrix4d,
3763 typename OutputMatrix4d>
3764 inline
3765 void ifftw4d(const InputMatrix4d & datain,OutputMatrix4d & dataout)
3766 {
3767  assert(datain.cols() == dataout.cols());
3768  assert(datain.rows() == dataout.rows());
3769  assert(datain.slices() == dataout.slices());
3770  assert(datain.slabs() == dataout.slabs());
3771  ifftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3772 }
3773 
3774 
3777 #endif //HAVE_FFTW
3778 
3783 
3820 template<typename InputBidirectionalIterator4d,
3821 typename OutputBidirectionalIterator4d>
3822 inline
3823 void real_fft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3824  InputBidirectionalIterator4d in_last_back_bottom_right,
3825  OutputBidirectionalIterator4d out_first_front_upper_left)
3826 {
3827 #ifdef HAVE_FFTW
3828  real_fftw4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3829 #else
3830  real_split_radix_fft4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3831 #endif
3832 }
3833 
3870 template<typename InputBidirectionalIterator4d,
3871 typename OutputBidirectionalIterator4d>
3872 inline
3873 void complex_fft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3874  InputBidirectionalIterator4d in_last_back_bottom_right,
3875  OutputBidirectionalIterator4d out_first_front_upper_left)
3876 {
3877 #ifdef HAVE_FFTW
3878  complex_fftw4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3879 #else
3880  complex_split_radix_fft4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3881 #endif
3882 }
3883 
3918 template<typename InputBidirectionalIterator4d,
3919 typename OutputBidirectionalIterator4d>
3920 inline
3921 void ifft4d(InputBidirectionalIterator4d in_first_front_upper_left,
3922  InputBidirectionalIterator4d in_last_back_bottom_right,
3923  OutputBidirectionalIterator4d out_first_front_upper_left)
3924 {
3925 #ifdef HAVE_FFTW
3926  ifftw4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3927 #else
3928  split_radix_ifft4d(in_first_front_upper_left,in_last_back_bottom_right,out_first_front_upper_left);
3929 #endif
3930 }
3931 
3932 
3970 template<typename HyperVolume1, typename HyperVolume2>
3971 inline
3972 void real_fft4d(const HyperVolume1 &datain, HyperVolume2 &dataout)
3973 {
3974  assert(datain.cols() == dataout.cols());
3975  assert(datain.rows() == dataout.rows());
3976  assert(datain.slices() == dataout.slices());
3977 #ifdef HAVE_FFTW
3978  real_fftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3979 #else
3980  real_split_radix_fft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
3981 #endif
3982 }
3983 
4021 template<typename HyperVolume1, typename HyperVolume2>
4022 inline
4023 void complex_fft4d(const HyperVolume1 &datain, HyperVolume2 &dataout)
4024 {
4025  assert(datain.cols() == dataout.cols());
4026  assert(datain.rows() == dataout.rows());
4027  assert(datain.slices() == dataout.slices());
4028 #ifdef HAVE_FFTW
4029  complex_fftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
4030 #else
4031  complex_split_radix_fft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
4032 #endif
4033 }
4034 
4073 template<typename HyperVolume1, typename HyperVolume2>
4074 inline
4075 void ifft4d(const HyperVolume1 &datain, HyperVolume2 &dataout)
4076 {
4077  assert(datain.cols() == dataout.cols());
4078  assert(datain.rows() == dataout.rows());
4079  assert(datain.slices() == dataout.slices());
4080 #ifdef HAVE_FFTW
4081  ifftw4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
4082 #else
4083  split_radix_ifft4d(datain.first_front_upper_left(),datain.last_back_bottom_right(),dataout.first_front_upper_left());
4084 #endif
4085 }
4086 
4087 
4088 
4089 
4090 
4091 
4094 }//slip::
4095 
4096 
4097 
4098 
4099 
4100 #endif //FFT_HPP
4101 
void complex_split_radix_fft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the complex split-radix fft4d of a container with 4d iterators.
Definition: FFT.hpp:3206
void split_radix_ifft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the split-radix ifft2d of a container with 2d iterators.
Definition: FFT.hpp:1039
void ifftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the fftw backward of a container (from fftw3.h).
Definition: FFT.hpp:1733
iterator2d upper_left()
Returns a read/write iterator2d that points to the first element of the Array2d. It points to the upp...
Definition: Array2d.hpp:2744
void complex_fft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the complex fft4d of a container with 4d iterators.
Definition: FFT.hpp:3873
void fftshift(ForwardIterator first, ForwardIterator last)
Performs a shift of a container, for use with the fft and ifft functions, in order to move the freque...
Definition: FFT.hpp:166
void fftshift2d(ForwardIterator2D upper_left, ForwardIterator2D bottom_right)
Performs a shift of a 2d container, for use with the fft2d and ifft2d functions, in order to move the...
Definition: FFT.hpp:191
void complex_fftw2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the complex fftw2d of a container with 2d iterators.
Definition: FFT.hpp:1831
iterator end()
Returns a read/write iterator that points one past the last element in the Vector. Iteration is done in ordinary element order.
Definition: Vector.hpp:1481
void real_split_radix_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the real split radix fft of a container (from FFTPACK.hpp).
Definition: FFT.hpp:774
iterator begin()
Returns a read/write iterator that points to the first element in the Array. Iteration is done in ord...
Definition: Array.hpp:1077
void radix2_fft(InputIter_T a, Iter_T b, const int log2n)
Computes the Cooley-Tukey (radix2) fft of a container.
Definition: FFT.hpp:292
void fftshift4d(ForwardIterator4D first_front_upper_left, ForwardIterator4D last_back_bottom_right)
Performs a shift of a 4d container, for use with the fft4d and ifft4d functions, in order to move the...
Definition: FFT.hpp:3058
void ifftw2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the ifftw2d of a container with 2d iterators.
Definition: FFT.hpp:1881
void real_fftw4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the real fftw4d of a container with 4d iterators.
Definition: FFT.hpp:3491
Real functor. Return the real part of x.
Definition: macros.hpp:98
void complex_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the complex fft2d of a container with 2d iterators.
Definition: FFT.hpp:2585
void ifft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the ifft4d of a container with 4d iterators.
Definition: FFT.hpp:3921
HyperVolume< T > exp(const HyperVolume< T > &M)
void real_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the real fft2d of a container with 2d iterators.
Definition: FFT.hpp:2553
void radix2_fft1d2d(InputIterator1d begin1, OuputIterator1d begin2, std::size_t nb_lig, std::size_t nb_col)
Computes the radix 2 fft2d of a container with 1d iterator (compatible with simple pointers)...
Definition: FFT.hpp:450
void ifftw4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the ifftw4d of a container with 4d iterators.
Definition: FFT.hpp:3636
void real(const Matrix1 &C, Matrix2 &R)
Extract the real Matrix of a complex Matrix. .
void complex_split_radix_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the complex split-radix fft2d of a container with 2d iterators.
Definition: FFT.hpp:984
static _Tp pi()
Constant .
Definition: macros.hpp:402
void fftduplicate(ForwardIterator1 begin1, ForwardIterator2 begin2, ForwardIterator2 end2)
Calculates the conjugates of the begin std::complex elements of a first container. Writes them at the end of a second container. This function is used to calculate the negative frequencies after a real fftw.
Definition: FFT.hpp:1582
Provides a class to manipulate 1d dynamic and generic arrays.
void split_radix_ifft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the split-radix ifft3d of a container with 3d iterators.
Definition: FFT.hpp:1251
void ifftw3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the ifftw3d of a container with 3d iterators.
Definition: FFT.hpp:2194
void ifft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the ifft3d of a container with 3d iterators.
Definition: FFT.hpp:2852
void complex_split_radix_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the complex split radix fft of a container (from FFTPACK.hpp).
Definition: FFT.hpp:828
void complex_split_radix_fft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the complex split-radix fft3d of a container with 3d iterators.
Definition: FFT.hpp:1303
void fftduplicate2d(ForwardIterator2D upper_left, ForwardIterator2D bottom_right)
Performs a duplicate of a 2D container, for use with the fft2d and ifft2d functions, in order to calculate the negative frequencies with the real fftw2d.
Definition: FFT.hpp:1613
void complex_fftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
complex fftw of a container (from fftw3.h).
Definition: FFT.hpp:1690
void dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the Discrete Cosinus Transform II (forward) of a container.
Definition: FFT.hpp:2993
unsigned int log2(unsigned int x)
Calculates the log2 of an integer.
Definition: FFT.hpp:135
void ifft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the inverse fft of a container.
Definition: FFT.hpp:2507
void idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the Discrete Cosinus Transform III (backward) of a container.
Definition: FFT.hpp:3029
void real_fftw3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the real fftw3d of a container with 3d iterators.
Definition: FFT.hpp:1982
void split_radix_ifft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the split radix fft backward of a container (from FFTPACK.hpp).
Definition: FFT.hpp:880
iterator end()
Returns a read/write iterator that points one past the last element in the Array. Iteration is done i...
Definition: Array.hpp:1084
Provides some statistics algorithms.
void radix2_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the radix2 (Cooley-Tukey) fft2d of a container with 2d iterators.
Definition: FFT.hpp:559
void real_fftw(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the real fftw of a container (from fftw3.h).
Definition: FFT.hpp:1641
void real_split_radix_fft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the real split-radix fft4d of a container with 4d iterators.
Definition: FFT.hpp:3137
void copy(_II first, _II last, _OI output_first)
Copy algorithm optimized for slip iterators.
Definition: copy_ext.hpp:177
Difference of Point2D class, specialization of DPoint<CoordType,DIM> with DIM = 2.
Definition: Array2d.hpp:129
void split_radix_idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the split radix Discrete Cosinus Transform III (i.e. idct) of a container (from FFTPACK...
Definition: FFT.hpp:1492
Numerical Vector class. This container statisfies the RandomAccessContainer concepts of the Standard ...
Definition: Vector.hpp:104
Imag functor. Return the imaginary part of x.
Definition: macros.hpp:112
unsigned int bitReverse(unsigned int x, const int log2n)
bitReverse function used in the Cooley-Tukey FFT algorithm
Definition: FFT.hpp:116
Provides some mathematical functors and constants.
iterator end()
Returns a read/write iterator that points one past the last element in the Array2d. Iteration is done in ordinary element order.
Definition: Array2d.hpp:2352
HyperVolume< T > sqrt(const HyperVolume< T > &M)
void complex_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the complex fft of a container.
Definition: FFT.hpp:2470
void resize(const size_type new_n, const T &val=T())
Resizes a Vector.
Definition: Vector.hpp:1465
void real_split_radix_fft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the real split-radix fft2d of a container with 2d iterators.
Definition: FFT.hpp:944
void fftw_idct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the Discrete Cosinus Transform III (backward) of a container (from fftw).
Definition: FFT.hpp:2362
void split_radix_dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the split radix Discrete Cosinus Transform II of a container (from FFTPACK.hpp).
Definition: FFT.hpp:1536
void real_fft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the real fft4d of a container with 4d iterators.
Definition: FFT.hpp:3823
void complex_fftw3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the complex fftw3d of a container with 3d iterators.
Definition: FFT.hpp:2037
void ifft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the ifft2d of a container with 2d iterators.
Definition: FFT.hpp:2629
void complex_fft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the complex fft3d of a container with 3d iterators.
Definition: FFT.hpp:2113
void fftshift3d(ForwardIterator3D front_upper_left, ForwardIterator3D back_bottom_right)
Performs a shift of a 2d container, for use with the fft2d and ifft2d functions, in order to move the...
Definition: FFT.hpp:230
void real_fftw2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the real fftw2d of a container with 2d iterators.
Definition: FFT.hpp:1778
void radix2_ifft1d2d(InputIterator1d begin1, OuputIterator1d begin2, std::size_t nb_lig, std::size_t nb_col)
Computes the inverse radix 2 fft2d of a container with 1d iterators (compatible with simple pointers)...
Definition: FFT.hpp:507
void real_fft(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the real fft of a container.
Definition: FFT.hpp:2430
Provides a class to manipulate Numerical Matrix.
void real_radix2_fft(InputIter_T a, Iter_T b, const int log2n)
Computes the Cooley-Tukey (radix2) real_fft of a container.
Definition: FFT.hpp:407
iterator begin()
Returns a read/write iterator that points to the first element in the Array2d. Iteration is done in o...
Definition: Array2d.hpp:2345
void radix2_ifft2d(InputBidirectionalIterator2d in_upper_left, InputBidirectionalIterator2d in_bottom_right, OutputBidirectionalIterator2d out_upper_left)
Computes the radix2 ifft2d of a container with 2d iterators.
Definition: FFT.hpp:620
void complex_fftw4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the complex fftw4d of a container with 4d iterators.
Definition: FFT.hpp:3563
This is a linear (one-dimensional) dynamic template container. This container statisfies the RandomAc...
Definition: Array.hpp:94
void split_radix_ifft4d(InputBidirectionalIterator4d in_first_front_upper_left, InputBidirectionalIterator4d in_last_back_bottom_right, OutputBidirectionalIterator4d out_first_front_upper_left)
Computes the split-radix ifft4d of a container with 4d iterators.
Definition: FFT.hpp:3293
void fftw_dct(InputIter in_begin, InputIter in_end, OutputIter out_begin)
Computes the Discrete Cosinus Transform II (forward) of a container (from fftw).
Definition: FFT.hpp:2327
void real_fft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the real fft3d of a container with 3d iterators.
Definition: FFT.hpp:2804
void radix2_ifft(InputIter_T a, Iter_T b, const int log2n)
Computes the inverse Cooley-Tukey (radix2) fft of a container.
Definition: FFT.hpp:351
void real_split_radix_fft3d(InputBidirectionalIterator3d in_front_upper_left, InputBidirectionalIterator3d in_back_bottom_right, OutputBidirectionalIterator3d out_front_upper_left)
Computes the real split-radix fft3d of a container with 3d iterators.
Definition: FFT.hpp:1180
Provides a class to manipulate numerical vectors.
This is a two-dimensional dynamic and generic container. This container statisfies the Bidirectionnal...
Definition: Array2d.hpp:135
iterator begin()
Returns a read/write iterator that points to the first element in the Vector. Iteration is done in or...
Definition: Vector.hpp:1474