1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Mark Borgerding mark a borgerding net
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.
namespace Eigen {
// FFTW uses non-const arguments
// so we must use ugly const_cast calls for all the args it uses
//
// This should be safe as long as
// 1. we use FFTW_ESTIMATE for all our planning
// see the FFTW docs section 4.3.2 "Planner Flags"
// 2. fftw_complex is compatible with std::complex
// This assumes std::complex<T> layout is array of size 2 with real,imag
template <typename T>
T * ei_fftw_cast(const T* p)
{
return const_cast<T*>( p);
}
fftw_complex * ei_fftw_cast( const std::complex<double> * p)
{
return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
}
fftwf_complex * ei_fftw_cast( const std::complex<float> * p)
{
return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
}
fftwl_complex * ei_fftw_cast( const std::complex<long double> * p)
{
return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
}
template <typename T>
struct ei_fftw_plan {};
template <>
struct ei_fftw_plan<float>
{
typedef float scalar_type;
typedef fftwf_complex complex_type;
fftwf_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
fftwf_execute_dft( m_plan, src,dst);
}
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
fftwf_execute_dft( m_plan, src,dst);
}
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
fftwf_execute_dft_r2c( m_plan,src,dst);
}
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
fftwf_execute_dft_c2r( m_plan, src,dst);
}
};
template <>
struct ei_fftw_plan<double>
{
typedef double scalar_type;
typedef fftw_complex complex_type;
fftw_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute_dft( m_plan, src,dst);
}
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
fftw_execute_dft( m_plan, src,dst);
}
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
fftw_execute_dft_r2c( m_plan,src,dst);
}
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
fftw_execute_dft_c2r( m_plan, src,dst);
}
};
template <>
struct ei_fftw_plan<long double>
{
typedef long double scalar_type;
typedef fftwl_complex complex_type;
fftwl_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
void fwd(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE);
fftwl_execute_dft( m_plan, src,dst);
}
void inv(complex_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE);
fftwl_execute_dft( m_plan, src,dst);
}
void fwd(complex_type * dst,scalar_type * src,int nfft) {
if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE);
fftwl_execute_dft_r2c( m_plan,src,dst);
}
void inv(scalar_type * dst,complex_type * src,int nfft) {
if (m_plan==NULL)
m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE);
fftwl_execute_dft_c2r( m_plan, src,dst);
}
};
template <typename _Scalar>
struct ei_fftw_impl
{
typedef _Scalar Scalar;
typedef std::complex<Scalar> Complex;
void clear()
{
m_plans.clear();
}
void fwd( Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
}
// real-to-complex forward FFT
void fwd( Complex * dst,const Scalar * src,int nfft)
{
get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft);
int nhbins=(nfft>>1)+1;
for (int k=nhbins;k < nfft; ++k )
dst[k] = conj(dst[nfft-k]);
}
// inverse complex-to-complex
void inv(Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
// scaling
Scalar s = 1./nfft;
for (int k=0;k<nfft;++k)
dst[k] *= s;
}
// half-complex to scalar
void inv( Scalar * dst,const Complex * src,int nfft)
{
get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
Scalar s = 1./nfft;
for (int k=0;k<nfft;++k)
dst[k] *= s;
}
private:
typedef ei_fftw_plan<Scalar> PlanData;
typedef std::map<int,PlanData> PlanMap;
PlanMap m_plans;
PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
{
bool inplace = (dst==src);
bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
int key = (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned;
return m_plans[key];
}
};
}
|
