source: Sophya/trunk/SophyaExt/IFFTW/fftwserver.cc@ 1528

Last change on this file since 1528 was 1424, checked in by ansari, 25 years ago

MAJ documentation - Reza 23/2/2001
File size: 11.6 KB
Line 
1#include "fftwserver.h"
2#include "FFTW/fftw.h"
3#include "FFTW/rfftw.h"
4
5/*!
6 \defgroup IFFTW IFFTW module
7 Module containing interface classes between Sophya objects
8 and FFTW Fourier transform package (see http://www.fftw.org )
9*/
10
11/*!
12 \class SOPHYA::FFTWServer
13 \ingroup IFFTW
14 An implementation of FFTServerInterface based on FFTW, double
15 precision arrays, using FFTW package, availabale from http://www.fftw.org.
16 Refer to FFTServerInterface for details about FFTServer operations.
17
18 \code
19 #include "fftwserver.h"
20 // ...
21 TMatrix<r_8> in(24,32);
22 TMatrix< complex<r_8> > out;
23 in = RandomSequence();
24 FFTWServer ffts;
25 ffts.setNormalize(true); // To have normalized transforms
26 cout << " FFTServer info string= " << ffts.getInfo() << endl;
27 cout << "in= " << in << endl;
28 cout << " Calling ffts.FFTForward(in, out) : " << endl;
29 ffts.FFTForward(in, out);
30 cout << "out= " << out << endl;
31 \endcode
32
33*/
34
35#define MAXND_FFTW 5
36
37namespace SOPHYA {
38
39class FFTWServerPlan {
40public:
41 FFTWServerPlan(int n, fftw_direction dir, bool fgreal=false);
42 FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal=false);
43 ~FFTWServerPlan();
44 void Recreate(int n);
45 void Recreate(int nd, int * nxyz);
46
47 static void FillSizes(const BaseArray & in, int * sz);
48
49 int _n; // Array dimension for 1-d arrays
50 int _nd; // Nb of dimensions for n-d arrays
51 int _nxyz[MAXND_FFTW]; // Array dimensions for n-d arrays
52 fftw_direction _dir;
53
54 fftw_plan p;
55 rfftw_plan rp;
56 fftwnd_plan pnd;
57 rfftwnd_plan rpnd;
58
59};
60
61} // Fin du namespace
62
63FFTWServerPlan::FFTWServerPlan(int n, fftw_direction dir, bool fgreal)
64{
65 if (n < 1)
66 throw ParmError("FFTWServerPlan: Array size <= 0 !");
67 p = NULL;
68 rp = NULL;
69 pnd = NULL;
70 rpnd = NULL;
71 _nd = -1;
72 for(int k=0; k<MAXND_FFTW; k++) _nxyz[k] = -10;
73 _n = n;
74 _dir = dir;
75 if (fgreal) {
76 rp = rfftw_create_plan(n, dir, FFTW_ESTIMATE);
77 if (rp == NULL)
78 throw AllocationError("FFTWServerPlan: failed to create plan (rp) !");
79 }
80 else {
81 p = fftw_create_plan(n, dir, FFTW_ESTIMATE);
82 if (p == NULL)
83 throw AllocationError("FFTWServerPlan: failed to create plan (p) !");
84 }
85
86}
87
88FFTWServerPlan::FFTWServerPlan(int nd, int * nxyz, fftw_direction dir, bool fgreal)
89{
90 int k;
91 if (nd > MAXND_FFTW)
92 throw ParmError("FFTWServerPlan: Array rank (nd) > MAXND_FFTW !");
93 p = NULL;
94 rp = NULL;
95 pnd = NULL;
96 rpnd = NULL;
97 _n = -10;
98 _nd = nd;
99 for(k=0; k<nd; k++) {
100 if (nxyz[k] < 1)
101 throw ParmError("FFTWServerPlan: One of the Array size <= 0 !");
102 _nxyz[k] = nxyz[k];
103 }
104 for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
105 _dir = dir;
106 if (fgreal) {
107 rpnd = rfftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
108 if (rpnd == NULL)
109 throw AllocationError("FFTWServerPlan: failed to create plan (rpnd) !");
110 }
111 else {
112 pnd = fftwnd_create_plan(_nd, _nxyz, dir, FFTW_ESTIMATE);
113 if (pnd == NULL)
114 throw AllocationError("FFTWServerPlan: failed to create plan (pnd) !");
115 }
116
117}
118
119FFTWServerPlan::~FFTWServerPlan()
120{
121 if (p) fftw_destroy_plan(p);
122 if (rp) rfftw_destroy_plan(rp);
123 if (pnd) fftwnd_destroy_plan(pnd);
124 if (rpnd) rfftwnd_destroy_plan(rpnd);
125}
126
127void
128FFTWServerPlan::Recreate(int n)
129{
130 if (n < 1)
131 throw ParmError("FFTWServerPlan::Recreate(n) n < 0 !");
132 if (_nd > 0)
133 throw ParmError("FFTWServerPlan::Recreate(n) Multi-dimensional plan ! > 0 !");
134 if (n == _n) return;
135 _n = n;
136 if (p) {
137 fftw_destroy_plan(p);
138 p = fftw_create_plan(n, _dir, FFTW_ESTIMATE);
139 if (p == NULL)
140 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (p) !");
141 }
142 else {
143 rfftw_destroy_plan(rp);
144 rp = rfftw_create_plan(n, _dir, FFTW_ESTIMATE);
145 if (rp == NULL)
146 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rp) !");
147 }
148
149}
150
151void
152FFTWServerPlan::Recreate(int nd, int * nxyz)
153{
154 if (_n > 0)
155 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) 1-dimensional plan !");
156 int k;
157 if (nd == _nd) {
158 bool samepl = true;
159 for (int k=0; k<nd; k++)
160 if (nxyz[k] != _nxyz[k]) samepl = false;
161 if (samepl) return;
162 }
163 if (nd > MAXND_FFTW)
164 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) Array rank (nd) > MAXND_FFTW !");
165 _nd = nd;
166 for(k=0; k<nd; k++) {
167 if (nxyz[k] < 1)
168 throw ParmError("FFTWServerPlan::Recreate(nd, nxyz) One of the Array size <= 0 !");
169 _nxyz[k] = nxyz[k];
170 }
171 for(k=nd; k<MAXND_FFTW; k++) _nxyz[k] = -10;
172 if (pnd) {
173 fftwnd_destroy_plan(pnd);
174 pnd = fftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
175 if (pnd == NULL)
176 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (pnd) !");
177 }
178 else {
179 rfftwnd_destroy_plan(rpnd);
180 rpnd = rfftwnd_create_plan(_nd, _nxyz, _dir, FFTW_ESTIMATE);
181 if (rpnd == NULL)
182 throw AllocationError("FFTWServerPlan::Recreate failed to create plan (rpnd) !");
183 }
184
185}
186
187void
188FFTWServerPlan::FillSizes(const BaseArray & in, int * sz)
189{
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}
196
197/* --Methode-- */
198FFTWServer::FFTWServer()
199 : FFTServerInterface("FFTServer using FFTW package")
200 , ckR8("FFTWServer - ", true, false)
201
202{
203 _p1df = NULL;
204 _p1db = NULL;
205 _pndf = NULL;
206 _pndb = NULL;
207
208 _p1drf = NULL;
209 _p1drb = NULL;
210 _pndrf = NULL;
211 _pndrb = NULL;
212}
213
214
215/* --Methode-- */
216FFTWServer::~FFTWServer()
217{
218 if (_p1df) delete _p1df ;
219 if (_p1db) delete _p1db ;
220 if (_pndf) delete _pndf ;
221 if (_pndb) delete _pndb ;
222
223 if (_p1drf) delete _p1drf ;
224 if (_p1drb) delete _p1drb ;
225 if (_pndrf) delete _pndrf ;
226 if (_pndrb) delete _pndrb ;
227}
228
229/* --Methode-- */
230FFTServerInterface * FFTWServer::Clone()
231{
232 return (new FFTWServer) ;
233}
234
235/* --Methode-- */
236void
237FFTWServer::FFTForward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
238{
239 int rank = ckR8.CheckResize(in, out);
240 if (rank == 1) { // One dimensional transform
241 if (_p1df) _p1df->Recreate(in.Size());
242 else _p1df = new FFTWServerPlan(in.Size(), FFTW_FORWARD, false);
243 fftw_one(_p1df->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
244 }
245 else { // Multi dimensional
246 if (in.NbDimensions() > MAXND_FFTW)
247 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
248 int sz[MAXND_FFTW];
249 FFTWServerPlan::FillSizes(in, sz);
250 // int k1 = 0;
251 // int k2 = 0;
252 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
253 // sz[k2] = in.Size(k1); k2++;
254 // }
255 if (_pndf) _pndf->Recreate(in.NbDimensions(), sz);
256 else _pndf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_FORWARD, false);
257 fftwnd_one(_pndf->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
258 }
259 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
260 return;
261}
262
263/* --Methode-- */
264void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< complex<r_8> > & out)
265{
266 int rank = ckR8.CheckResize(in, out);
267 if (rank == 1) { // One dimensional transform
268 if (_p1db) _p1db->Recreate(in.Size());
269 else _p1db = new FFTWServerPlan(in.Size(), FFTW_BACKWARD, false);
270 fftw_one(_p1db->p, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
271 }
272 else { // Multi dimensional
273 if (in.NbDimensions() > MAXND_FFTW)
274 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
275 int sz[MAXND_FFTW];
276 FFTWServerPlan::FillSizes(in, sz);
277 // int k1 = 0;
278 // int k2 = 0;
279 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
280 // sz[k2] = in.Size(k1); k2++;
281 // }
282 if (_pndb) _pndb->Recreate(in.NbDimensions(), sz);
283 else _pndb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_BACKWARD, false);
284 fftwnd_one(_pndb->pnd, (fftw_complex *)(in.Data()) , (fftw_complex *)(out.Data()) );
285 }
286
287 return;
288}
289
290
291void FFTWServer::FFTForward(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
292{
293 int rank = ckR8.CheckResize(in, out);
294 if (rank == 1) { // One dimensional transform
295 if (_p1drf) _p1drf->Recreate(in.Size());
296 else _p1drf = new FFTWServerPlan(in.Size(), FFTW_REAL_TO_COMPLEX, true);
297 TArray<r_8> outtemp;
298 outtemp.ReSize(in);
299 rfftw_one(_p1drf->rp, (fftw_real *)(in.Data()) , (fftw_real *)(outtemp.Data()));
300 ReShapetoCompl(outtemp, out);
301 }
302 else { // Multi dimensional
303 if (in.NbDimensions() > MAXND_FFTW)
304 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
305 int sz[MAXND_FFTW];
306 FFTWServerPlan::FillSizes(in, sz);
307 // int k1 = 0;
308 // int k2 = 0;
309 // for(k1=in.NbDimensions()-1; k1>=0; k1--) {
310 // sz[k2] = in.Size(k1); k2++;
311 // }
312 if (_pndrf) _pndrf->Recreate(in.NbDimensions(), sz);
313 else _pndrf = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_REAL_TO_COMPLEX, true);
314 rfftwnd_one_real_to_complex(_pndrf->rpnd, (fftw_real *)(in.Data()) ,
315 (fftw_complex *)(out.Data()) );
316 }
317 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
318 return;
319
320}
321
322
323
324void FFTWServer::FFTBackward(TArray< complex<r_8> > const & in, TArray< r_8 > & out,
325 bool usoutsz)
326{
327
328 int rank = ckR8.CheckResize(in, out, usoutsz);
329 if (rank == 1) { // One dimensional transform
330 TArray<r_8> intemp;
331 intemp.ReSize(out);
332 if (_p1drb) _p1drb->Recreate(out.Size());
333 else _p1drb = new FFTWServerPlan(out.Size(), FFTW_COMPLEX_TO_REAL, true);
334
335 ReShapetoReal(in, intemp);
336 // cerr << " DEBUG-FFTWServer::FFTBackward() in = \n" << in << endl;
337 // cerr << " DEBUG-FFTWServer::FFTBackward() intemp = \n" << intemp << endl;
338 rfftw_one(_p1drb->rp, (fftw_real *)(intemp.Data()) , (fftw_real *)(out.Data()));
339
340 }
341 else { // Multi dimensional
342 if (in.NbDimensions() > MAXND_FFTW)
343 throw ParmError("FFTWServer::FFTForward( complex<r_8>, r_8 ) rank > MAXND_FFTW !");
344 int sz[MAXND_FFTW];
345 FFTWServerPlan::FillSizes(out, sz);
346 // int k1 = 0;
347 // int k2 = 0;
348 // for(k1=out.NbDimensions()-1; k1>=0; k1--) {
349 // sz[k2] = out.Size(k1); k2++;
350 // }
351 if (_pndrb) _pndrb->Recreate(in.NbDimensions(), sz);
352 else _pndrb = new FFTWServerPlan(in.NbDimensions(), sz, FFTW_COMPLEX_TO_REAL, true);
353
354 rfftwnd_one_complex_to_real(_pndrb->rpnd, (fftw_complex *)(in.Data()) , (fftw_real *)(out.Data()) );
355 }
356 return;
357}
358
359
360
361/* --Methode-- */
362void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa)
363{
364 TVector< complex<r_8> > in(ina);
365 TVector< r_8> out(outa);
366 int n = in.NElts();
367 r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
368 sa_size_t ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
369 ncs = 2*n-1 : ncs = 2*n-2;
370
371 if (out.NElts() != ncs) {
372 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
373 << " out.NElts()= " << out.NElts() << endl;
374 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
375 }
376 // if (ncs == 2*n-2) {
377 out(0) = in(0).real();
378 for(int k=1; k<n; k++) {
379 out(k) = in(k).real();
380 out(ncs-k) = in(k).imag();
381 }
382 if (ncs == 2*n-2)
383 out(ncs-1) = in(n-1).real();
384
385 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
386}
387
388
389/* --Methode-- */
390void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
391{
392 TVector< r_8> in(ina);
393 TVector< complex<r_8> > out(outa);
394
395 sa_size_t n = in.NElts();
396 sa_size_t ncs = n/2+1;
397 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
398 if (out.NElts() != ncs)
399 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
400
401 out(0) = complex<r_8> (in(0),0.);
402 for(int k=1; k<ncs; k++)
403 out(k) = complex<r_8>(in(k), in(n-k));
404 if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
405
406}
407
Note: See TracBrowser for help on using the repository browser.