1 | #include "fftwserver.h"
|
---|
2 | #include "FFTW/fftw.h"
|
---|
3 | #include "FFTW/rfftw.h"
|
---|
4 |
|
---|
5 |
|
---|
6 | #define MAXND_FFTW 5
|
---|
7 |
|
---|
8 | class FFTWServerPlan {
|
---|
9 | public:
|
---|
10 | FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false);
|
---|
11 | FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false);
|
---|
12 | ~FFTWServerPlan();
|
---|
13 | void Recreate(int n);
|
---|
14 | void Recreate(int nd, int * nxyz);
|
---|
15 |
|
---|
16 | static void FillSizes(const BaseArray & in, int * sz);
|
---|
17 |
|
---|
18 | int _n; // Array dimension for 1-d arrays
|
---|
19 | int _nd; // Nb of dimensions for n-d arrays
|
---|
20 | int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays
|
---|
21 | fftw_direction _dir;
|
---|
22 |
|
---|
23 | fftw_plan p;
|
---|
24 | rfftw_plan rp;
|
---|
25 | fftwnd_plan pnd;
|
---|
26 | rfftwnd_plan rpnd;
|
---|
27 |
|
---|
28 | };
|
---|
29 |
|
---|
30 | FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal)
|
---|
31 | {
|
---|
32 | if (n < 1)
|
---|
33 | throw ParmError("FFTWServerPlan: Array size <= 0 !");
|
---|
34 | p = NULL;
|
---|
35 | rp = NULL;
|
---|
36 | pnd = NULL;
|
---|
37 | rpnd = NULL;
|
---|
38 | _nd = -1;
|
---|
39 | for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
---|
40 | _n = n;
|
---|
41 | _dir = dir;
|
---|
42 | if (fgreal) rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE);
|
---|
43 | else p = fftw_create_plan(n, dir, FFTW_ESTIMATE);
|
---|
44 | }
|
---|
45 |
|
---|
46 | FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal)
|
---|
47 | {
|
---|
48 | int k;
|
---|
49 | if (nd > MAXND_FFTW)
|
---|
50 | throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !");
|
---|
51 | p = NULL;
|
---|
52 | rp = NULL;
|
---|
53 | pnd = NULL;
|
---|
54 | rpnd = NULL;
|
---|
55 | _n = -10;
|
---|
56 | _nd = nd;
|
---|
57 | for(k=0; k<nd; k++) {
|
---|
58 | if (nxyz[k] < 1)
|
---|
59 | throw ParmError("FFTWServerPlan: One of the Array size <= 0 !");
|
---|
60 | _nxyz[k] = nxyz[k];
|
---|
61 | }
|
---|
62 | for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
---|
63 | _dir = dir;
|
---|
64 | if (fgreal) rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
|
---|
65 | else pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
|
---|
66 | }
|
---|
67 |
|
---|
68 | FFTWServerPlan::~FFTWServerPlan()
|
---|
69 | {
|
---|
70 | if (p) fftw_destroy_plan(p);
|
---|
71 | if (rp) rfftw_destroy_plan(rp);
|
---|
72 | if (pnd) fftwnd_destroy_plan(pnd);
|
---|
73 | if (rpnd) rfftwnd_destroy_plan(rpnd);
|
---|
74 | }
|
---|
75 |
|
---|
76 | void
|
---|
77 | FFTWServerPlan::Recreate(int n)
|
---|
78 | {
|
---|
79 | if (n < 1)
|
---|
80 | throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !");
|
---|
81 | if (_nd > 0)
|
---|
82 | throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !");
|
---|
83 | if (n == _n) return;
|
---|
84 | _n = n;
|
---|
85 | if (p) {
|
---|
86 | fftw_destroy_plan(p);
|
---|
87 | p = fftw_create_plan(n, _dir, FFTW_ESTIMATE);
|
---|
88 | }
|
---|
89 | else {
|
---|
90 | rfftw_destroy_plan(rp);
|
---|
91 | rp = rfftw_create_plan(n, _dir, FFTW_ESTIMATE);
|
---|
92 | }
|
---|
93 | }
|
---|
94 |
|
---|
95 | void
|
---|
96 | FFTWServerPlan::Recreate(int nd, int * nxyz)
|
---|
97 | {
|
---|
98 | if (_n > 0)
|
---|
99 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !");
|
---|
100 | int k;
|
---|
101 | if (nd == _nd) {
|
---|
102 | bool samepl = true;
|
---|
103 | for (int k=0; k<nd; k++)
|
---|
104 | if (nxyz[k] != _nxyz[k]) samepl = false;
|
---|
105 | if (samepl) return;
|
---|
106 | }
|
---|
107 | if (nd > MAXND_FFTW)
|
---|
108 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !");
|
---|
109 | _nd = nd;
|
---|
110 | for(k=0; k<nd; k++) {
|
---|
111 | if (nxyz[k] < 1)
|
---|
112 | throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !");
|
---|
113 | _nxyz[k] = nxyz[k];
|
---|
114 | }
|
---|
115 | for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
|
---|
116 | if (pnd) {
|
---|
117 | fftwnd_destroy_plan(pnd);
|
---|
118 | pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
|
---|
119 | }
|
---|
120 | else {
|
---|
121 | rfftwnd_destroy_plan(rpnd);
|
---|
122 | rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
|
---|
123 | }
|
---|
124 |
|
---|
125 | }
|
---|
126 |
|
---|
127 | void
|
---|
128 | FFTWServerPlan::FillSizes(const BaseArray & in, int * sz)
|
---|
129 | {
|
---|
130 | int k1 = 0;
|
---|
131 | int k2 = 0;
|
---|
132 | for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
---|
133 | sz[k2] = in.Size(k1); k2++;
|
---|
134 | }
|
---|
135 | }
|
---|
136 |
|
---|
137 | /* --Methode-- */
|
---|
138 | FFTWServer::FFTWServer()
|
---|
139 | : FFTServerInterface("FFTServer using FFTW package")
|
---|
140 | , ckR8("FFTWServer - ", true, false)
|
---|
141 |
|
---|
142 | {
|
---|
143 | _p1df = NULL;
|
---|
144 | _p1db = NULL;
|
---|
145 | _pndf = NULL;
|
---|
146 | _pndb = NULL;
|
---|
147 |
|
---|
148 | _p1drf = NULL;
|
---|
149 | _p1drb = NULL;
|
---|
150 | _pndrf = NULL;
|
---|
151 | _pndrb = NULL;
|
---|
152 | }
|
---|
153 |
|
---|
154 |
|
---|
155 | /* --Methode-- */
|
---|
156 | FFTWServer::~FFTWServer()
|
---|
157 | {
|
---|
158 | if (_p1df) delete _p1df ;
|
---|
159 | if (_p1db) delete _p1db ;
|
---|
160 | if (_pndf) delete _pndf ;
|
---|
161 | if (_pndb) delete _pndb ;
|
---|
162 |
|
---|
163 | if (_p1drf) delete _p1drf ;
|
---|
164 | if (_p1drb) delete _p1drb ;
|
---|
165 | if (_pndrf) delete _pndrf ;
|
---|
166 | if (_pndrb) delete _pndrb ;
|
---|
167 | }
|
---|
168 |
|
---|
169 | /* --Methode-- */
|
---|
170 | FFTServerInterface * FFTWServer::Clone()
|
---|
171 | {
|
---|
172 | return (new FFTWServer) ;
|
---|
173 | }
|
---|
174 |
|
---|
175 | /* --Methode-- */
|
---|
176 | void
|
---|
177 | FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
|
---|
178 | {
|
---|
179 | int rank = ckR8.CheckResize(in, out);
|
---|
180 | if (rank == 1) { // One dimensional transform
|
---|
181 | if (_p1df) _p1df->Recreate(in.Size());
|
---|
182 | else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false);
|
---|
183 | fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
---|
184 | }
|
---|
185 | else { // Multi dimensional
|
---|
186 | if (in.NbDimensions() > MAXND_FFTW)
|
---|
187 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
---|
188 | int sz[MAXND_FFTW];
|
---|
189 | FFTWServerPlan::FillSizes(in, sz);
|
---|
190 | // int k1 = 0;
|
---|
191 | // int k2 = 0;
|
---|
192 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
---|
193 | // sz[k2] = in.Size(k1); k2++;
|
---|
194 | // }
|
---|
195 | if (_pndf) _pndf->Recreate(in.NbDimensions(), sz);
|
---|
196 | else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false);
|
---|
197 | fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
---|
198 | }
|
---|
199 | if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
|
---|
200 | return;
|
---|
201 | }
|
---|
202 |
|
---|
203 | /* --Methode-- */
|
---|
204 | void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
|
---|
205 | {
|
---|
206 | int rank = ckR8.CheckResize(in, out);
|
---|
207 | if (rank == 1) { // One dimensional transform
|
---|
208 | if (_p1db) _p1db->Recreate(in.Size());
|
---|
209 | else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false);
|
---|
210 | fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
---|
211 | }
|
---|
212 | else { // Multi dimensional
|
---|
213 | if (in.NbDimensions() > MAXND_FFTW)
|
---|
214 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
---|
215 | int sz[MAXND_FFTW];
|
---|
216 | FFTWServerPlan::FillSizes(in, sz);
|
---|
217 | // int k1 = 0;
|
---|
218 | // int k2 = 0;
|
---|
219 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
---|
220 | // sz[k2] = in.Size(k1); k2++;
|
---|
221 | // }
|
---|
222 | if (_pndb) _pndb->Recreate(in.NbDimensions(), sz);
|
---|
223 | else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false);
|
---|
224 | fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
|
---|
225 | }
|
---|
226 |
|
---|
227 | return;
|
---|
228 | }
|
---|
229 |
|
---|
230 |
|
---|
231 | void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
|
---|
232 | {
|
---|
233 | int rank = ckR8.CheckResize(in, out);
|
---|
234 | if (rank == 1) { // One dimensional transform
|
---|
235 | if (_p1drf) _p1drf->Recreate(in.Size());
|
---|
236 | else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true);
|
---|
237 | TArray<r_8> outtemp;
|
---|
238 | outtemp.ReSize(in);
|
---|
239 | rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data()));
|
---|
240 | ReShapetoCompl(outtemp, out);
|
---|
241 | }
|
---|
242 | else { // Multi dimensional
|
---|
243 | if (in.NbDimensions() > MAXND_FFTW)
|
---|
244 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
|
---|
245 | int sz[MAXND_FFTW];
|
---|
246 | FFTWServerPlan::FillSizes(in, sz);
|
---|
247 | // int k1 = 0;
|
---|
248 | // int k2 = 0;
|
---|
249 | // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
|
---|
250 | // sz[k2] = in.Size(k1); k2++;
|
---|
251 | // }
|
---|
252 | if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz);
|
---|
253 | else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true);
|
---|
254 | rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) ,
|
---|
255 | (fftw_complex *)(out.Data()) );
|
---|
256 | }
|
---|
257 | if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
|
---|
258 | return;
|
---|
259 |
|
---|
260 | }
|
---|
261 |
|
---|
262 |
|
---|
263 |
|
---|
264 | void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out,
|
---|
265 | bool usoutsz)
|
---|
266 | {
|
---|
267 |
|
---|
268 | int rank = ckR8.CheckResize(in, out, usoutsz);
|
---|
269 | if (rank == 1) { // One dimensional transform
|
---|
270 | TArray<r_8> intemp;
|
---|
271 | intemp.ReSize(out);
|
---|
272 | if (_p1drb) _p1drb->Recreate(out.Size());
|
---|
273 | else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
|
---|
274 |
|
---|
275 | ReShapetoReal(in, intemp);
|
---|
276 | // cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
|
---|
277 | // cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
|
---|
278 | rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(out.Data()));
|
---|
279 |
|
---|
280 | }
|
---|
281 | else { // Multi dimensional
|
---|
282 | if (in.NbDimensions() > MAXND_FFTW)
|
---|
283 | throw ParmError("FFTWServer::FFTForward( complex<r_8>, r_8 ) rank > MAXND_FFTW !");
|
---|
284 | int sz[MAXND_FFTW];
|
---|
285 | FFTWServerPlan::FillSizes(out, sz);
|
---|
286 | // int k1 = 0;
|
---|
287 | // int k2 = 0;
|
---|
288 | // for(k1=out.NbDimensions()-1; k1>=0; k1--) {
|
---|
289 | // sz[k2] = out.Size(k1); k2++;
|
---|
290 | // }
|
---|
291 | if (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
|
---|
292 | else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
|
---|
293 |
|
---|
294 | rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
|
---|
295 | }
|
---|
296 | return;
|
---|
297 | }
|
---|
298 |
|
---|
299 |
|
---|
300 |
|
---|
301 | /* --Methode-- */
|
---|
302 | void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa)
|
---|
303 | {
|
---|
304 | TVector< complex<r_8> > in(ina);
|
---|
305 | TVector< r_8> out(outa);
|
---|
306 | int n = in.NElts();
|
---|
307 | r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
|
---|
308 | sa_size_t ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
|
---|
309 | ncs = 2*n-1 : ncs = 2*n-2;
|
---|
310 |
|
---|
311 | if (out.NElts() != ncs) {
|
---|
312 | cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
|
---|
313 | << " out.NElts()= " << out.NElts() << endl;
|
---|
314 | throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
|
---|
315 | }
|
---|
316 | // if (ncs == 2*n-2) {
|
---|
317 | out(0) = in(0).real();
|
---|
318 | for(int k=1; k<n; k++) {
|
---|
319 | out(k) = in(k).real();
|
---|
320 | out(ncs-k) = in(k).imag();
|
---|
321 | }
|
---|
322 | if (ncs == 2*n-2)
|
---|
323 | out(ncs-1) = in(n-1).real();
|
---|
324 |
|
---|
325 | // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
|
---|
326 | }
|
---|
327 |
|
---|
328 |
|
---|
329 | /* --Methode-- */
|
---|
330 | void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
|
---|
331 | {
|
---|
332 | TVector< r_8> in(ina);
|
---|
333 | TVector< complex<r_8> > out(outa);
|
---|
334 |
|
---|
335 | sa_size_t n = in.NElts();
|
---|
336 | sa_size_t ncs = n/2+1;
|
---|
337 | sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
|
---|
338 | if (out.NElts() != ncs)
|
---|
339 | throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
|
---|
340 |
|
---|
341 | out(0) = complex<r_8> (in(0),0.);
|
---|
342 | for(int k=1; k<ncs; k++)
|
---|
343 | out(k) = complex<r_8>(in(k), in(n-k));
|
---|
344 | if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
|
---|
345 |
|
---|
346 | }
|
---|
347 |
|
---|