![]() |
Eigen-unsupported
3.3.3
|
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: */