src/FFT/ei_fftw_impl.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Mark Borgerding mark a borgerding net
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 namespace Eigen {
0011
0012 namespace internal {
0013
0014   // FFTW uses non-const arguments
0015   // so we must use ugly const_cast calls for all the args it uses
0016   //
0017   // This should be safe as long as
0018   // 1. we use FFTW_ESTIMATE for all our planning
0019   //       see the FFTW docs section 4.3.2 "Planner Flags"
0020   // 2. fftw_complex is compatible with std::complex
0021   //    This assumes std::complex<T> layout is array of size 2 with real,imag
0022   template <typename T>
0023   inline
0024   T * fftw_cast(const T* p)
0025   {
0026       return const_cast<T*>( p);
0027   }
0028
0029   inline
0030   fftw_complex * fftw_cast( const std::complex<double> * p)
0031   {
0032       return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
0033   }
0034
0035   inline
0036   fftwf_complex * fftw_cast( const std::complex<float> * p)
0037   {
0038       return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
0039   }
0040
0041   inline
0042   fftwl_complex * fftw_cast( const std::complex<long double> * p)
0043   {
0044       return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
0045   }
0046
0047   template <typename T>
0048   struct fftw_plan {};
0049
0050   template <>
0051   struct fftw_plan<float>
0052   {
0053       typedef float scalar_type;
0054       typedef fftwf_complex complex_type;
0055       fftwf_plan m_plan;
0056       fftw_plan() :m_plan(NULL) {}
0057       ~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
0058
0059       inline
0060       void fwd(complex_type * dst,complex_type * src,int nfft) {
0061           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0062           fftwf_execute_dft( m_plan, src,dst);
0063       }
0064       inline
0065       void inv(complex_type * dst,complex_type * src,int nfft) {
0066           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0067           fftwf_execute_dft( m_plan, src,dst);
0068       }
0069       inline
0070       void fwd(complex_type * dst,scalar_type * src,int nfft) {
0071           if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0072           fftwf_execute_dft_r2c( m_plan,src,dst);
0073       }
0074       inline
0075       void inv(scalar_type * dst,complex_type * src,int nfft) {
0076           if (m_plan==NULL)
0077               m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0078           fftwf_execute_dft_c2r( m_plan, src,dst);
0079       }
0080
0081       inline
0082       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
0083           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0084           fftwf_execute_dft( m_plan, src,dst);
0085       }
0086       inline
0087       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
0088           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0089           fftwf_execute_dft( m_plan, src,dst);
0090       }
0091
0092   };
0093   template <>
0094   struct fftw_plan<double>
0095   {
0096       typedef double scalar_type;
0097       typedef fftw_complex complex_type;
0098       ::fftw_plan m_plan;
0099       fftw_plan() :m_plan(NULL) {}
0100       ~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
0101
0102       inline
0103       void fwd(complex_type * dst,complex_type * src,int nfft) {
0104           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0105           fftw_execute_dft( m_plan, src,dst);
0106       }
0107       inline
0108       void inv(complex_type * dst,complex_type * src,int nfft) {
0109           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0110           fftw_execute_dft( m_plan, src,dst);
0111       }
0112       inline
0113       void fwd(complex_type * dst,scalar_type * src,int nfft) {
0114           if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0115           fftw_execute_dft_r2c( m_plan,src,dst);
0116       }
0117       inline
0118       void inv(scalar_type * dst,complex_type * src,int nfft) {
0119           if (m_plan==NULL)
0120               m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0121           fftw_execute_dft_c2r( m_plan, src,dst);
0122       }
0123       inline
0124       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
0125           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0126           fftw_execute_dft( m_plan, src,dst);
0127       }
0128       inline
0129       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
0130           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0131           fftw_execute_dft( m_plan, src,dst);
0132       }
0133   };
0134   template <>
0135   struct fftw_plan<long double>
0136   {
0137       typedef long double scalar_type;
0138       typedef fftwl_complex complex_type;
0139       fftwl_plan m_plan;
0140       fftw_plan() :m_plan(NULL) {}
0141       ~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
0142
0143       inline
0144       void fwd(complex_type * dst,complex_type * src,int nfft) {
0145           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0146           fftwl_execute_dft( m_plan, src,dst);
0147       }
0148       inline
0149       void inv(complex_type * dst,complex_type * src,int nfft) {
0150           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0151           fftwl_execute_dft( m_plan, src,dst);
0152       }
0153       inline
0154       void fwd(complex_type * dst,scalar_type * src,int nfft) {
0155           if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0156           fftwl_execute_dft_r2c( m_plan,src,dst);
0157       }
0158       inline
0159       void inv(scalar_type * dst,complex_type * src,int nfft) {
0160           if (m_plan==NULL)
0161               m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0162           fftwl_execute_dft_c2r( m_plan, src,dst);
0163       }
0164       inline
0165       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
0166           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0167           fftwl_execute_dft( m_plan, src,dst);
0168       }
0169       inline
0170       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
0171           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
0172           fftwl_execute_dft( m_plan, src,dst);
0173       }
0174   };
0175
0176   template <typename _Scalar>
0177   struct fftw_impl
0178   {
0179       typedef _Scalar Scalar;
0180       typedef std::complex<Scalar> Complex;
0181
0182       inline
0183       void clear()
0184       {
0185         m_plans.clear();
0186       }
0187
0188       // complex-to-complex forward FFT
0189       inline
0190       void fwd( Complex * dst,const Complex *src,int nfft)
0191       {
0192         get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
0193       }
0194
0195       // real-to-complex forward FFT
0196       inline
0197       void fwd( Complex * dst,const Scalar * src,int nfft)
0198       {
0199           get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
0200       }
0201
0202       // 2-d complex-to-complex
0203       inline
0204       void fwd2(Complex * dst, const Complex * src, int n0,int n1)
0205       {
0206           get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
0207       }
0208
0209       // inverse complex-to-complex
0210       inline
0211       void inv(Complex * dst,const Complex  *src,int nfft)
0212       {
0213         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
0214       }
0215
0216       // half-complex to scalar
0217       inline
0218       void inv( Scalar * dst,const Complex * src,int nfft)
0219       {
0220         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
0221       }
0222
0223       // 2-d complex-to-complex
0224       inline
0225       void inv2(Complex * dst, const Complex * src, int n0,int n1)
0226       {
0227         get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
0228       }
0229
0230
0231   protected:
0232       typedef fftw_plan<Scalar> PlanData;
0233
0234       typedef Eigen::numext::int64_t int64_t;
0235
0236       typedef std::map<int64_t,PlanData> PlanMap;
0237
0238       PlanMap m_plans;
0239
0240       inline
0241       PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
0242       {
0243           bool inplace = (dst==src);
0244           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
0245           int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
0246           return m_plans[key];
0247       }
0248
0249       inline
0250       PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
0251       {
0252           bool inplace = (dst==src);
0253           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
0254           int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
0255           return m_plans[key];
0256       }
0257   };
0258
0259 } // end namespace internal
0260
0261 } // end namespace Eigen