39 T * fftw_cast(
const T* p)
41 return const_cast<T*
>( p);
45 fftw_complex * fftw_cast(
const std::complex<double> * p)
47 return const_cast<fftw_complex*
>(
reinterpret_cast<const fftw_complex*
>(p) );
51 fftwf_complex * fftw_cast(
const std::complex<float> * p)
53 return const_cast<fftwf_complex*
>(
reinterpret_cast<const fftwf_complex*
>(p) );
57 fftwl_complex * fftw_cast(
const std::complex<long double> * p)
59 return const_cast<fftwl_complex*
>(
reinterpret_cast<const fftwl_complex*
>(p) );
66 struct fftw_plan<float>
68 typedef float scalar_type;
69 typedef fftwf_complex complex_type;
71 fftw_plan() :m_plan(NULL) {}
72 ~fftw_plan() {
if (m_plan) fftwf_destroy_plan(m_plan);}
75 void fwd(complex_type * dst,complex_type * src,
int nfft) {
76 if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
77 fftwf_execute_dft( m_plan, src,dst);
80 void inv(complex_type * dst,complex_type * src,
int nfft) {
81 if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
82 fftwf_execute_dft( m_plan, src,dst);
85 void fwd(complex_type * dst,scalar_type * src,
int nfft) {
86 if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
87 fftwf_execute_dft_r2c( m_plan,src,dst);
90 void inv(scalar_type * dst,complex_type * src,
int nfft) {
92 m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
93 fftwf_execute_dft_c2r( m_plan, src,dst);
97 void fwd2( complex_type * dst,complex_type * src,
int n0,
int n1) {
98 if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
99 fftwf_execute_dft( m_plan, src,dst);
102 void inv2( complex_type * dst,complex_type * src,
int n0,
int n1) {
103 if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
104 fftwf_execute_dft( m_plan, src,dst);
109 struct fftw_plan<double>
111 typedef double scalar_type;
112 typedef fftw_complex complex_type;
114 fftw_plan() :m_plan(NULL) {}
115 ~fftw_plan() {
if (m_plan) fftw_destroy_plan(m_plan);}
118 void fwd(complex_type * dst,complex_type * src,
int nfft) {
119 if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
120 fftw_execute_dft( m_plan, src,dst);
123 void inv(complex_type * dst,complex_type * src,
int nfft) {
124 if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
125 fftw_execute_dft( m_plan, src,dst);
128 void fwd(complex_type * dst,scalar_type * src,
int nfft) {
129 if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
130 fftw_execute_dft_r2c( m_plan,src,dst);
133 void inv(scalar_type * dst,complex_type * src,
int nfft) {
135 m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
136 fftw_execute_dft_c2r( m_plan, src,dst);
139 void fwd2( complex_type * dst,complex_type * src,
int n0,
int n1) {
140 if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
141 fftw_execute_dft( m_plan, src,dst);
144 void inv2( complex_type * dst,complex_type * src,
int n0,
int n1) {
145 if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
146 fftw_execute_dft( m_plan, src,dst);
150 struct fftw_plan<long double>
152 typedef long double scalar_type;
153 typedef fftwl_complex complex_type;
155 fftw_plan() :m_plan(NULL) {}
156 ~fftw_plan() {
if (m_plan) fftwl_destroy_plan(m_plan);}
159 void fwd(complex_type * dst,complex_type * src,
int nfft) {
160 if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
161 fftwl_execute_dft( m_plan, src,dst);
164 void inv(complex_type * dst,complex_type * src,
int nfft) {
165 if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
166 fftwl_execute_dft( m_plan, src,dst);
169 void fwd(complex_type * dst,scalar_type * src,
int nfft) {
170 if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
171 fftwl_execute_dft_r2c( m_plan,src,dst);
174 void inv(scalar_type * dst,complex_type * src,
int nfft) {
176 m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
177 fftwl_execute_dft_c2r( m_plan, src,dst);
180 void fwd2( complex_type * dst,complex_type * src,
int n0,
int n1) {
181 if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
182 fftwl_execute_dft( m_plan, src,dst);
185 void inv2( complex_type * dst,complex_type * src,
int n0,
int n1) {
186 if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
187 fftwl_execute_dft( m_plan, src,dst);
191 template <
typename _Scalar>
194 typedef _Scalar Scalar;
195 typedef std::complex<Scalar> Complex;
205 void fwd( Complex * dst,
const Complex *src,
int nfft)
207 get_plan(nfft,
false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
212 void fwd( Complex * dst,
const Scalar * src,
int nfft)
214 get_plan(nfft,
false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
219 void fwd2(Complex * dst,
const Complex * src,
int n0,
int n1)
221 get_plan(n0,n1,
false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
226 void inv(Complex * dst,
const Complex *src,
int nfft)
228 get_plan(nfft,
true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
233 void inv( Scalar * dst,
const Complex * src,
int nfft)
235 get_plan(nfft,
true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
240 void inv2(Complex * dst,
const Complex * src,
int n0,
int n1)
242 get_plan(n0,n1,
true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
247 typedef fftw_plan<Scalar> PlanData;
249 typedef std::map<int64_t,PlanData> PlanMap;
254 PlanData & get_plan(
int nfft,
bool inverse,
void * dst,
const void * src)
256 bool inplace = (dst==src);
257 bool aligned = ( (
reinterpret_cast<size_t>(src)&15) | (
reinterpret_cast<size_t>(dst)&15) ) == 0;
258 int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
263 PlanData & get_plan(
int n0,
int n1,
bool inverse,
void * dst,
const void * src)
265 bool inplace = (dst==src);
266 bool aligned = ( (
reinterpret_cast<size_t>(src)&15) | (
reinterpret_cast<size_t>(dst)&15) ) == 0;
267 int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;