Home | History | Annotate | Download | only in FFT 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2009 Mark Borgerding mark a borgerding net
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 namespace Eigen {
     11 
     12 namespace internal {
     13 
     14   // FFTW uses non-const arguments
     15   // so we must use ugly const_cast calls for all the args it uses
     16   //
     17   // This should be safe as long as
     18   // 1. we use FFTW_ESTIMATE for all our planning
     19   //       see the FFTW docs section 4.3.2 "Planner Flags"
     20   // 2. fftw_complex is compatible with std::complex
     21   //    This assumes std::complex<T> layout is array of size 2 with real,imag
     22   template <typename T>
     23   inline
     24   T * fftw_cast(const T* p)
     25   {
     26       return const_cast<T*>( p);
     27   }
     28 
     29   inline
     30   fftw_complex * fftw_cast( const std::complex<double> * p)
     31   {
     32       return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
     33   }
     34 
     35   inline
     36   fftwf_complex * fftw_cast( const std::complex<float> * p)
     37   {
     38       return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
     39   }
     40 
     41   inline
     42   fftwl_complex * fftw_cast( const std::complex<long double> * p)
     43   {
     44       return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
     45   }
     46 
     47   template <typename T>
     48   struct fftw_plan {};
     49 
     50   template <>
     51   struct fftw_plan<float>
     52   {
     53       typedef float scalar_type;
     54       typedef fftwf_complex complex_type;
     55       fftwf_plan m_plan;
     56       fftw_plan() :m_plan(NULL) {}
     57       ~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
     58 
     59       inline
     60       void fwd(complex_type * dst,complex_type * src,int nfft) {
     61           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     62           fftwf_execute_dft( m_plan, src,dst);
     63       }
     64       inline
     65       void inv(complex_type * dst,complex_type * src,int nfft) {
     66           if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     67           fftwf_execute_dft( m_plan, src,dst);
     68       }
     69       inline
     70       void fwd(complex_type * dst,scalar_type * src,int nfft) {
     71           if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     72           fftwf_execute_dft_r2c( m_plan,src,dst);
     73       }
     74       inline
     75       void inv(scalar_type * dst,complex_type * src,int nfft) {
     76           if (m_plan==NULL)
     77               m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     78           fftwf_execute_dft_c2r( m_plan, src,dst);
     79       }
     80 
     81       inline
     82       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
     83           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     84           fftwf_execute_dft( m_plan, src,dst);
     85       }
     86       inline
     87       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
     88           if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
     89           fftwf_execute_dft( m_plan, src,dst);
     90       }
     91 
     92   };
     93   template <>
     94   struct fftw_plan<double>
     95   {
     96       typedef double scalar_type;
     97       typedef fftw_complex complex_type;
     98       ::fftw_plan m_plan;
     99       fftw_plan() :m_plan(NULL) {}
    100       ~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
    101 
    102       inline
    103       void fwd(complex_type * dst,complex_type * src,int nfft) {
    104           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    105           fftw_execute_dft( m_plan, src,dst);
    106       }
    107       inline
    108       void inv(complex_type * dst,complex_type * src,int nfft) {
    109           if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    110           fftw_execute_dft( m_plan, src,dst);
    111       }
    112       inline
    113       void fwd(complex_type * dst,scalar_type * src,int nfft) {
    114           if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    115           fftw_execute_dft_r2c( m_plan,src,dst);
    116       }
    117       inline
    118       void inv(scalar_type * dst,complex_type * src,int nfft) {
    119           if (m_plan==NULL)
    120               m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    121           fftw_execute_dft_c2r( m_plan, src,dst);
    122       }
    123       inline
    124       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
    125           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    126           fftw_execute_dft( m_plan, src,dst);
    127       }
    128       inline
    129       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
    130           if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    131           fftw_execute_dft( m_plan, src,dst);
    132       }
    133   };
    134   template <>
    135   struct fftw_plan<long double>
    136   {
    137       typedef long double scalar_type;
    138       typedef fftwl_complex complex_type;
    139       fftwl_plan m_plan;
    140       fftw_plan() :m_plan(NULL) {}
    141       ~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
    142 
    143       inline
    144       void fwd(complex_type * dst,complex_type * src,int nfft) {
    145           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    146           fftwl_execute_dft( m_plan, src,dst);
    147       }
    148       inline
    149       void inv(complex_type * dst,complex_type * src,int nfft) {
    150           if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    151           fftwl_execute_dft( m_plan, src,dst);
    152       }
    153       inline
    154       void fwd(complex_type * dst,scalar_type * src,int nfft) {
    155           if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    156           fftwl_execute_dft_r2c( m_plan,src,dst);
    157       }
    158       inline
    159       void inv(scalar_type * dst,complex_type * src,int nfft) {
    160           if (m_plan==NULL)
    161               m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    162           fftwl_execute_dft_c2r( m_plan, src,dst);
    163       }
    164       inline
    165       void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
    166           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    167           fftwl_execute_dft( m_plan, src,dst);
    168       }
    169       inline
    170       void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
    171           if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
    172           fftwl_execute_dft( m_plan, src,dst);
    173       }
    174   };
    175 
    176   template <typename _Scalar>
    177   struct fftw_impl
    178   {
    179       typedef _Scalar Scalar;
    180       typedef std::complex<Scalar> Complex;
    181 
    182       inline
    183       void clear()
    184       {
    185         m_plans.clear();
    186       }
    187 
    188       // complex-to-complex forward FFT
    189       inline
    190       void fwd( Complex * dst,const Complex *src,int nfft)
    191       {
    192         get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
    193       }
    194 
    195       // real-to-complex forward FFT
    196       inline
    197       void fwd( Complex * dst,const Scalar * src,int nfft)
    198       {
    199           get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
    200       }
    201 
    202       // 2-d complex-to-complex
    203       inline
    204       void fwd2(Complex * dst, const Complex * src, int n0,int n1)
    205       {
    206           get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
    207       }
    208 
    209       // inverse complex-to-complex
    210       inline
    211       void inv(Complex * dst,const Complex  *src,int nfft)
    212       {
    213         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
    214       }
    215 
    216       // half-complex to scalar
    217       inline
    218       void inv( Scalar * dst,const Complex * src,int nfft)
    219       {
    220         get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
    221       }
    222 
    223       // 2-d complex-to-complex
    224       inline
    225       void inv2(Complex * dst, const Complex * src, int n0,int n1)
    226       {
    227         get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
    228       }
    229 
    230 
    231   protected:
    232       typedef fftw_plan<Scalar> PlanData;
    233 
    234       typedef std::map<int64_t,PlanData> PlanMap;
    235 
    236       PlanMap m_plans;
    237 
    238       inline
    239       PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
    240       {
    241           bool inplace = (dst==src);
    242           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
    243           int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
    244           return m_plans[key];
    245       }
    246 
    247       inline
    248       PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
    249       {
    250           bool inplace = (dst==src);
    251           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
    252           int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
    253           return m_plans[key];
    254       }
    255   };
    256 
    257 } // end namespace internal
    258 
    259 } // end namespace Eigen
    260 
    261 /* vim: set filetype=cpp et sw=2 ts=2 ai: */
    262