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