00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 namespace Eigen {
00026
00027 namespace internal {
00028
00029
00030
00031
00032
00033
00034
00035
00036
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
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
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
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
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
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
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 }
00273
00274 }
00275
00276