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 Eigen::numext::int64_t int64_t;
235 
236       typedef std::map<int64_t,PlanData> PlanMap;
237 
238       PlanMap m_plans;
239 
240       inline
241       PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
242       {
243           bool inplace = (dst==src);
244           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
245           int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
246           return m_plans[key];
247       }
248 
249       inline
250       PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
251       {
252           bool inplace = (dst==src);
253           bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
254           int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
255           return m_plans[key];
256       }
257   };
258 
259 } // end namespace internal
260 
261 } // end namespace Eigen
262 
263 /* vim: set filetype=cpp et sw=2 ts=2 ai: */
264