source: Sophya/trunk/SophyaExt/IFFTW/fftw3server.cc@ 3391

Last change on this file since 3391 was 3359, checked in by ansari, 18 years ago

Ajout methode version r_4 (float) a FFTWServer (sous condition de flag) - Reza 23/10/2007
File size: 13.3 KB
Line 
1#include "FFTW/fftw3.h"
2
3
4static void FillSizes4FFTW(const BaseArray & in, int * sz)
5{
6 int k1 = 0;
7 int k2 = 0;
8 for(k1=in.NbDimensions()-1; k1>=0; k1--) {
9 sz[k2] = in.Size(k1); k2++;
10 }
11}
12
13/* --Methode-- */
14//! Constructor - If preserve_input==true, input arrays will not be overwritten.
15FFTWServer::FFTWServer(bool preserve_input)
16#ifdef ALSO_FFTW_FLOAT_EXTSOP
17 : FFTServerInterface("FFTServer using FFTW3 package (r_4 and r_8)") ,
18#else
19 : FFTServerInterface("FFTServer using FFTW3 package (r_8 only)") ,
20#endif
21 ckR8("FFTWServer - ", true, false),
22 ckR4("FFTWServer - ", true, false),
23 _preserve_input(preserve_input)
24{
25}
26
27
28/* --Methode-- */
29FFTWServer::~FFTWServer()
30{
31}
32
33/* --Methode-- */
34FFTServerInterface * FFTWServer::Clone()
35{
36 return (new FFTWServer) ;
37}
38
39/* --Methode-- */
40void
41FFTWServer::FFTForward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
42{
43 int rank = ckR8.CheckResize(in, out);
44 if (rank == 1) { // One dimensional transform
45 fftw_plan p = fftw_plan_dft_1d(in.Size(), (fftw_complex *)in.Data(),
46 (fftw_complex *)out.Data(),
47 FFTW_FORWARD, FFTW_ESTIMATE);
48 if (p == NULL)
49 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) Error creating fftw_plan");
50 fftw_execute(p);
51 fftw_destroy_plan(p);
52 }
53 else { // Multi dimensional
54 if (in.NbDimensions() > MAXND_FFTW)
55 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
56 int sz[MAXND_FFTW];
57 FillSizes4FFTW(in, sz);
58 fftw_plan p = fftw_plan_dft(rank, sz,
59 (fftw_complex *)in.Data(), (fftw_complex *)out.Data(),
60 FFTW_FORWARD, FFTW_ESTIMATE);
61 if (p == NULL)
62 throw ParmError("FFTWServer::FFTForward( complex<r_8>, complex<r_8> ) Error creating fftw_plan");
63 fftw_execute(p);
64 fftw_destroy_plan(p);
65 }
66 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
67 return;
68}
69
70/* --Methode-- */
71void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
72{
73 int rank = ckR8.CheckResize(in, out);
74 if (rank == 1) { // One dimensional transform
75 fftw_plan p = fftw_plan_dft_1d(in.Size(), (fftw_complex *)in.Data(),
76 (fftw_complex *)out.Data(),
77 FFTW_BACKWARD, FFTW_ESTIMATE);
78 if (p == NULL)
79 throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) Error creating fftw_plan");
80 fftw_execute(p);
81 fftw_destroy_plan(p);
82 }
83 else { // Multi dimensional
84 if (in.NbDimensions() > MAXND_FFTW)
85 throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) rank > MAXND_FFTW !");
86 int sz[MAXND_FFTW];
87 FillSizes4FFTW(in, sz);
88 fftw_plan p = fftw_plan_dft(rank, sz,
89 (fftw_complex *)in.Data(), (fftw_complex *)out.Data(),
90 FFTW_BACKWARD, FFTW_ESTIMATE);
91 if (p == NULL)
92 throw ParmError("FFTWServer::FFTBackward( complex<r_8>, complex<r_8> ) Error creating fftw_plan");
93 fftw_execute(p);
94 fftw_destroy_plan(p);
95 }
96
97 return;
98}
99
100
101void FFTWServer::FFTForward(TArray< r_8 > & in, TArray< complex<r_8> > & out)
102{
103 int rank = ckR8.CheckResize(in, out);
104 if (rank == 1) { // One dimensional transform
105 fftw_plan p = fftw_plan_dft_r2c_1d(in.Size(), in.Data(),
106 (fftw_complex *)out.Data(),
107 FFTW_ESTIMATE);
108 if (p == NULL)
109 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan");
110 fftw_execute(p);
111 fftw_destroy_plan(p);
112
113
114 }
115 else { // Multi dimensional
116 if (in.NbDimensions() > MAXND_FFTW)
117 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) rank > MAXND_FFTW !");
118 int sz[MAXND_FFTW];
119 FillSizes4FFTW(in, sz);
120 fftw_plan p = fftw_plan_dft_r2c(rank, sz, in.Data(),
121 (fftw_complex *)out.Data(),
122 FFTW_ESTIMATE);
123 if (p == NULL)
124 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan");
125 fftw_execute(p);
126 fftw_destroy_plan(p);
127 }
128 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
129 return;
130
131}
132
133
134
135void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< r_8 > & out,
136 bool usoutsz)
137{
138 // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
139 // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
140 int rank = ckR8.CheckResize(in, out, usoutsz);
141 bool share = (_preserve_input) ? false : true;
142 TArray< complex<r_8> > inp(in, share);
143 if (rank == 1) { // One dimensional transform
144 fftw_plan p = fftw_plan_dft_c2r_1d(out.Size(), (fftw_complex *)inp.Data(),
145 out.Data(),
146 FFTW_ESTIMATE);
147 if (p == NULL)
148 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan");
149 fftw_execute(p);
150 fftw_destroy_plan(p);
151 }
152 else { // Multi dimensional
153 if (in.NbDimensions() > MAXND_FFTW)
154 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) rank > MAXND_FFTW !");
155 int sz[MAXND_FFTW];
156 FillSizes4FFTW(out, sz);
157 fftw_plan p = fftw_plan_dft_c2r(rank, sz, (fftw_complex *)inp.Data(),
158 out.Data(),
159 FFTW_ESTIMATE);
160 if (p == NULL)
161 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan");
162 fftw_execute(p);
163 fftw_destroy_plan(p);
164 }
165 return;
166}
167
168
169
170/* --Methode-- */
171void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa,
172 bool usz)
173{
174 TVector< complex<r_8> > in(ina);
175 TVector< r_8> out(outa);
176 int n = in.NElts();
177 r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
178 sa_size_t ncs;
179 if (usz) {
180 if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
181 throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
182 ncs = out.NElts();
183 }
184 else {
185 ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
186 ncs = 2*n-1 : ncs = 2*n-2;
187
188 if (out.NElts() != ncs) {
189 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
190 << " out.NElts()= " << out.NElts() << endl;
191 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
192 }
193 }
194 // if (ncs == 2*n-2) {
195 out(0) = in(0).real();
196 for(int k=1; k<n; k++) {
197 out(k) = in(k).real();
198 out(ncs-k) = in(k).imag();
199 }
200 if (ncs == 2*n-2)
201 out(n-1) = in(n-1).real();
202
203 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
204}
205
206
207/* --Methode-- */
208void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
209{
210 TVector< r_8> in(ina);
211 TVector< complex<r_8> > out(outa);
212
213 sa_size_t n = in.NElts();
214 sa_size_t ncs = n/2+1;
215 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
216 if (out.NElts() != ncs)
217 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
218
219 out(0) = complex<r_8> (in(0),0.);
220 for(int k=1; k<ncs; k++)
221 out(k) = complex<r_8>(in(k), in(n-k));
222 if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
223
224}
225
226
227#ifdef ALSO_FFTW_FLOAT_EXTSOP
228 /* ---------------------------------------------------------------------------
229 We define here single precision (float) version of FFTWServr methods
230 Needs the libfftw3f.a , in addition to libfftw3.a
231 --------------------------------------------------------------------------- */
232/* --Methode-- */
233void
234FFTWServer::FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
235{
236 int rank = ckR4.CheckResize(in, out);
237 if (rank == 1) { // One dimensional transform
238 fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
239 (fftwf_complex *)out.Data(),
240 FFTW_FORWARD, FFTW_ESTIMATE);
241 if (p == NULL)
242 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
243 fftwf_execute(p);
244 fftwf_destroy_plan(p);
245 }
246 else { // Multi dimensional
247 if (in.NbDimensions() > MAXND_FFTW)
248 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
249 int sz[MAXND_FFTW];
250 FillSizes4FFTW(in, sz);
251 fftwf_plan p = fftwf_plan_dft(rank, sz,
252 (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
253 FFTW_FORWARD, FFTW_ESTIMATE);
254 if (p == NULL)
255 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
256 fftwf_execute(p);
257 fftwf_destroy_plan(p);
258 }
259 if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
260 return;
261}
262
263/* --Methode-- */
264void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
265{
266 int rank = ckR4.CheckResize(in, out);
267 if (rank == 1) { // One dimensional transform
268 fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
269 (fftwf_complex *)out.Data(),
270 FFTW_BACKWARD, FFTW_ESTIMATE);
271 if (p == NULL)
272 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
273 fftwf_execute(p);
274 fftwf_destroy_plan(p);
275 }
276 else { // Multi dimensional
277 if (in.NbDimensions() > MAXND_FFTW)
278 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
279 int sz[MAXND_FFTW];
280 FillSizes4FFTW(in, sz);
281 fftwf_plan p = fftwf_plan_dft(rank, sz,
282 (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
283 FFTW_BACKWARD, FFTW_ESTIMATE);
284 if (p == NULL)
285 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
286 fftwf_execute(p);
287 fftwf_destroy_plan(p);
288 }
289
290 return;
291}
292
293
294void FFTWServer::FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out)
295{
296 int rank = ckR4.CheckResize(in, out);
297 if (rank == 1) { // One dimensional transform
298 fftwf_plan p = fftwf_plan_dft_r2c_1d(in.Size(), in.Data(),
299 (fftwf_complex *)out.Data(),
300 FFTW_ESTIMATE);
301 if (p == NULL)
302 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
303 fftwf_execute(p);
304 fftwf_destroy_plan(p);
305
306
307 }
308 else { // Multi dimensional
309 if (in.NbDimensions() > MAXND_FFTW)
310 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
311 int sz[MAXND_FFTW];
312 FillSizes4FFTW(in, sz);
313 fftwf_plan p = fftwf_plan_dft_r2c(rank, sz, in.Data(),
314 (fftwf_complex *)out.Data(),
315 FFTW_ESTIMATE);
316 if (p == NULL)
317 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
318 fftwf_execute(p);
319 fftwf_destroy_plan(p);
320 }
321 if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
322 return;
323
324}
325
326
327
328void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & out,
329 bool usoutsz)
330{
331 // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
332 // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
333 int rank = ckR4.CheckResize(in, out, usoutsz);
334 bool share = (_preserve_input) ? false : true;
335 TArray< complex<r_4> > inp(in, share);
336 if (rank == 1) { // One dimensional transform
337 fftwf_plan p = fftwf_plan_dft_c2r_1d(out.Size(), (fftwf_complex *)inp.Data(),
338 out.Data(),
339 FFTW_ESTIMATE);
340 if (p == NULL)
341 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
342 fftwf_execute(p);
343 fftwf_destroy_plan(p);
344 }
345 else { // Multi dimensional
346 if (in.NbDimensions() > MAXND_FFTW)
347 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
348 int sz[MAXND_FFTW];
349 FillSizes4FFTW(out, sz);
350 fftwf_plan p = fftwf_plan_dft_c2r(rank, sz, (fftwf_complex *)inp.Data(),
351 out.Data(),
352 FFTW_ESTIMATE);
353 if (p == NULL)
354 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
355 fftwf_execute(p);
356 fftwf_destroy_plan(p);
357 }
358 return;
359}
360
361
362
363/* --Methode-- */
364void FFTWServer::ReShapetoReal(TArray< complex<r_4> > const & ina, TArray< r_4 > & outa,
365 bool usz)
366{
367 TVector< complex<r_4> > in(ina);
368 TVector< r_4> out(outa);
369 int n = in.NElts();
370 r_4 thr = FFTArrayChecker<r_4>::ZeroThreshold();
371 sa_size_t ncs;
372 if (usz) {
373 if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
374 throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
375 ncs = out.NElts();
376 }
377 else {
378 ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
379 ncs = 2*n-1 : ncs = 2*n-2;
380
381 if (out.NElts() != ncs) {
382 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
383 << " out.NElts()= " << out.NElts() << endl;
384 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
385 }
386 }
387 // if (ncs == 2*n-2) {
388 out(0) = in(0).real();
389 for(int k=1; k<n; k++) {
390 out(k) = in(k).real();
391 out(ncs-k) = in(k).imag();
392 }
393 if (ncs == 2*n-2)
394 out(n-1) = in(n-1).real();
395
396 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
397}
398
399
400/* --Methode-- */
401void FFTWServer::ReShapetoCompl(TArray< r_4 > const & ina, TArray< complex<r_4> > & outa)
402{
403 TVector< r_4> in(ina);
404 TVector< complex<r_4> > out(outa);
405
406 sa_size_t n = in.NElts();
407 sa_size_t ncs = n/2+1;
408 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
409 if (out.NElts() != ncs)
410 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
411
412 out(0) = complex<r_4> (in(0),0.);
413 for(int k=1; k<ncs; k++)
414 out(k) = complex<r_4>(in(k), in(n-k));
415 if (n%2 == 0) out(ncs-1) = complex<r_4>(in(n/2), 0.);
416
417}
418
419#endif
Note: See TracBrowser for help on using the repository browser.