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 fftw_cast(const T * p)24 T * fftw_cast(const T* p) 25 { 26 return const_cast<T*>( p); 27 } 28 29 inline fftw_cast(const std::complex<double> * p)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 fftw_cast(const std::complex<float> * p)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 fftw_cast(const std::complex<long double> * p)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