ei_fftw_impl.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. 
00003 //
00004 // Copyright (C) 2009 Mark Borgerding mark a borgerding net
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 namespace Eigen { 
00011 
00012 namespace internal {
00013 
00014   // FFTW uses non-const arguments
00015   // so we must use ugly const_cast calls for all the args it uses
00016   //
00017   // This should be safe as long as 
00018   // 1. we use FFTW_ESTIMATE for all our planning
00019   //       see the FFTW docs section 4.3.2 "Planner Flags"
00020   // 2. fftw_complex is compatible with std::complex
00021   //    This assumes std::complex<T> layout is array of size 2 with real,imag
00022   template <typename T> 
00023   inline 
00024   T * fftw_cast(const T* p)
00025   { 
00026       return const_cast<T*>( p); 
00027   }
00028 
00029   inline 
00030   fftw_complex * fftw_cast( const std::complex<double> * p)
00031   {
00032       return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) ); 
00033   }
00034 
00035   inline 
00036   fftwf_complex * fftw_cast( const std::complex<float> * p)
00037   { 
00038       return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) ); 
00039   }
00040 
00041   inline 
00042   fftwl_complex * fftw_cast( const std::complex<long double> * p)
00043   { 
00044       return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) ); 
00045   }
00046 
00047   template <typename T> 
00048   struct fftw_plan {};
00049 
00050   template <> 
00051   struct fftw_plan<float>
00052   {
00053       typedef float scalar_type;
00054       typedef fftwf_complex complex_type;
00055       fftwf_plan m_plan;
00056       fftw_plan() :m_plan(NULL) {}
00057       ~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
00058 
00059       inline
00060       void fwd(complex_type * dst,complex_type * src,int nfft) {
00061           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00062           fftwf_execute_dft( m_plan, src,dst);
00063       }
00064       inline
00065       void inv(complex_type * dst,complex_type * src,int nfft) {
00066           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00067           fftwf_execute_dft( m_plan, src,dst);
00068       }
00069       inline
00070       void fwd(complex_type * dst,scalar_type * src,int nfft) {
00071           if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00072           fftwf_execute_dft_r2c( m_plan,src,dst);
00073       }
00074       inline
00075       void inv(scalar_type * dst,complex_type * src,int nfft) {
00076           if (m_plan==NULL)
00077               m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00078           fftwf_execute_dft_c2r( m_plan, src,dst);
00079       }
00080 
00081       inline 
00082       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
00083           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00084           fftwf_execute_dft( m_plan, src,dst);
00085       }
00086       inline 
00087       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
00088           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00089           fftwf_execute_dft( m_plan, src,dst);
00090       }
00091 
00092   };
00093   template <> 
00094   struct fftw_plan<double>
00095   {
00096       typedef double scalar_type;
00097       typedef fftw_complex complex_type;
00098       ::fftw_plan m_plan;
00099       fftw_plan() :m_plan(NULL) {}
00100       ~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
00101 
00102       inline
00103       void fwd(complex_type * dst,complex_type * src,int nfft) {
00104           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00105           fftw_execute_dft( m_plan, src,dst);
00106       }
00107       inline
00108       void inv(complex_type * dst,complex_type * src,int nfft) {
00109           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00110           fftw_execute_dft( m_plan, src,dst);
00111       }
00112       inline
00113       void fwd(complex_type * dst,scalar_type * src,int nfft) {
00114           if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00115           fftw_execute_dft_r2c( m_plan,src,dst);
00116       }
00117       inline
00118       void inv(scalar_type * dst,complex_type * src,int nfft) {
00119           if (m_plan==NULL)
00120               m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00121           fftw_execute_dft_c2r( m_plan, src,dst);
00122       }
00123       inline 
00124       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
00125           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00126           fftw_execute_dft( m_plan, src,dst);
00127       }
00128       inline 
00129       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
00130           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00131           fftw_execute_dft( m_plan, src,dst);
00132       }
00133   };
00134   template <> 
00135   struct fftw_plan<long double>
00136   {
00137       typedef long double scalar_type;
00138       typedef fftwl_complex complex_type;
00139       fftwl_plan m_plan;
00140       fftw_plan() :m_plan(NULL) {}
00141       ~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
00142 
00143       inline
00144       void fwd(complex_type * dst,complex_type * src,int nfft) {
00145           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00146           fftwl_execute_dft( m_plan, src,dst);
00147       }
00148       inline
00149       void inv(complex_type * dst,complex_type * src,int nfft) {
00150           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00151           fftwl_execute_dft( m_plan, src,dst);
00152       }
00153       inline
00154       void fwd(complex_type * dst,scalar_type * src,int nfft) {
00155           if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00156           fftwl_execute_dft_r2c( m_plan,src,dst);
00157       }
00158       inline
00159       void inv(scalar_type * dst,complex_type * src,int nfft) {
00160           if (m_plan==NULL)
00161               m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00162           fftwl_execute_dft_c2r( m_plan, src,dst);
00163       }
00164       inline 
00165       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
00166           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00167           fftwl_execute_dft( m_plan, src,dst);
00168       }
00169       inline 
00170       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
00171           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
00172           fftwl_execute_dft( m_plan, src,dst);
00173       }
00174   };
00175 
00176   template <typename _Scalar>
00177   struct fftw_impl
00178   {
00179       typedef _Scalar Scalar;
00180       typedef std::complex<Scalar> Complex;
00181 
00182       inline
00183       void clear() 
00184       {
00185         m_plans.clear();
00186       }
00187 
00188       // complex-to-complex forward FFT
00189       inline
00190       void fwd( Complex * dst,const Complex *src,int nfft)
00191       {
00192         get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
00193       }
00194 
00195       // real-to-complex forward FFT
00196       inline
00197       void fwd( Complex * dst,const Scalar * src,int nfft) 
00198       {
00199           get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
00200       }
00201 
00202       // 2-d complex-to-complex
00203       inline
00204       void fwd2(Complex * dst, const Complex * src, int n0,int n1)
00205       {
00206           get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
00207       }
00208 
00209       // inverse complex-to-complex
00210       inline
00211       void inv(Complex * dst,const Complex  *src,int nfft)
00212       {
00213         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
00214       }
00215 
00216       // half-complex to scalar
00217       inline
00218       void inv( Scalar * dst,const Complex * src,int nfft) 
00219       {
00220         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
00221       }
00222 
00223       // 2-d complex-to-complex
00224       inline
00225       void inv2(Complex * dst, const Complex * src, int n0,int n1)
00226       {
00227         get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
00228       }
00229 
00230 
00231   protected:
00232       typedef fftw_plan<Scalar> PlanData;
00233 
00234       typedef std::map<int64_t,PlanData> PlanMap;
00235 
00236       PlanMap m_plans;
00237 
00238       inline
00239       PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
00240       {
00241           bool inplace = (dst==src);
00242           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
00243           int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
00244           return m_plans[key];
00245       }
00246 
00247       inline
00248       PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
00249       {
00250           bool inplace = (dst==src);
00251           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
00252           int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
00253           return m_plans[key];
00254       }
00255   };
00256 
00257 } // end namespace internal
00258 
00259 } // end namespace Eigen
00260 
00261 /* vim: set filetype=cpp et sw=2 ts=2 ai: */


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:07