source: pacpussensors/trunk/Vislab/lib3dv/eigen/unsupported/Eigen/src/FFT/ei_fftw_impl.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 9.0 KB
Line 
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
10namespace Eigen {
11
12namespace 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: */
Note: See TracBrowser for help on using the repository browser.