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

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

Correction remplacement rank par NbDimension() lors de FFT tableaux multi-dim avec certaines tailles= 1 - Reza(+cmv) 29/11/2007
File size: 13.4 KB
RevLine 
[3000]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)
[3359]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)
[3000]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);
[3412]58 fftw_plan p = fftw_plan_dft(in.NbDimensions(), sz,
[3000]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);
[3412]88 fftw_plan p = fftw_plan_dft(in.NbDimensions(), sz,
[3000]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(),
[3412]106 (fftw_complex *)out.Data(), FFTW_ESTIMATE);
[3000]107 if (p == NULL)
108 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan");
109 fftw_execute(p);
110 fftw_destroy_plan(p);
111
112
113 }
114 else { // Multi dimensional
115 if (in.NbDimensions() > MAXND_FFTW)
116 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) rank > MAXND_FFTW !");
117 int sz[MAXND_FFTW];
118 FillSizes4FFTW(in, sz);
[3412]119 fftw_plan p = fftw_plan_dft_r2c(in.NbDimensions(), sz, in.Data(),
120 (fftw_complex *)out.Data(), FFTW_ESTIMATE);
[3000]121 if (p == NULL)
122 throw ParmError("FFTWServer::FFTForward(r_8, complex<r_8> ) Error creating fftw_plan");
123 fftw_execute(p);
124 fftw_destroy_plan(p);
125 }
126 if(this->getNormalize()) out=out/complex<r_8>((double)in.Size(),0.);
127 return;
128
129}
130
131
132
133void FFTWServer::FFTBackward(TArray< complex<r_8> > & in, TArray< r_8 > & out,
134 bool usoutsz)
135{
136 // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
137 // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
138 int rank = ckR8.CheckResize(in, out, usoutsz);
139 bool share = (_preserve_input) ? false : true;
140 TArray< complex<r_8> > inp(in, share);
141 if (rank == 1) { // One dimensional transform
142 fftw_plan p = fftw_plan_dft_c2r_1d(out.Size(), (fftw_complex *)inp.Data(),
[3412]143 out.Data(), FFTW_ESTIMATE);
[3000]144 if (p == NULL)
145 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan");
146 fftw_execute(p);
147 fftw_destroy_plan(p);
148 }
149 else { // Multi dimensional
150 if (in.NbDimensions() > MAXND_FFTW)
151 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) rank > MAXND_FFTW !");
152 int sz[MAXND_FFTW];
153 FillSizes4FFTW(out, sz);
[3412]154 fftw_plan p = fftw_plan_dft_c2r(out.NbDimensions(), sz, (fftw_complex *)inp.Data(),
155 out.Data(), FFTW_ESTIMATE);
[3000]156 if (p == NULL)
157 throw ParmError("FFTWServer::FFTBackward(r_8, complex<r_8> ) Error creating fftw_plan");
158 fftw_execute(p);
159 fftw_destroy_plan(p);
160 }
161 return;
162}
163
164
165
166/* --Methode-- */
167void FFTWServer::ReShapetoReal(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa,
168 bool usz)
169{
170 TVector< complex<r_8> > in(ina);
171 TVector< r_8> out(outa);
172 int n = in.NElts();
173 r_8 thr = FFTArrayChecker<r_8>::ZeroThreshold();
174 sa_size_t ncs;
175 if (usz) {
176 if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
177 throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
178 ncs = out.NElts();
179 }
180 else {
181 ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
182 ncs = 2*n-1 : ncs = 2*n-2;
183
184 if (out.NElts() != ncs) {
185 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
186 << " out.NElts()= " << out.NElts() << endl;
187 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
188 }
189 }
190 // if (ncs == 2*n-2) {
191 out(0) = in(0).real();
192 for(int k=1; k<n; k++) {
193 out(k) = in(k).real();
194 out(ncs-k) = in(k).imag();
195 }
196 if (ncs == 2*n-2)
197 out(n-1) = in(n-1).real();
198
199 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
200}
201
202
203/* --Methode-- */
204void FFTWServer::ReShapetoCompl(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
205{
206 TVector< r_8> in(ina);
207 TVector< complex<r_8> > out(outa);
208
209 sa_size_t n = in.NElts();
210 sa_size_t ncs = n/2+1;
211 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
212 if (out.NElts() != ncs)
213 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
214
215 out(0) = complex<r_8> (in(0),0.);
216 for(int k=1; k<ncs; k++)
217 out(k) = complex<r_8>(in(k), in(n-k));
218 if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n/2), 0.);
219
220}
221
[3359]222
223#ifdef ALSO_FFTW_FLOAT_EXTSOP
224 /* ---------------------------------------------------------------------------
225 We define here single precision (float) version of FFTWServr methods
226 Needs the libfftw3f.a , in addition to libfftw3.a
227 --------------------------------------------------------------------------- */
228/* --Methode-- */
229void
230FFTWServer::FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
231{
232 int rank = ckR4.CheckResize(in, out);
233 if (rank == 1) { // One dimensional transform
234 fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
235 (fftwf_complex *)out.Data(),
236 FFTW_FORWARD, FFTW_ESTIMATE);
237 if (p == NULL)
238 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
239 fftwf_execute(p);
240 fftwf_destroy_plan(p);
241 }
242 else { // Multi dimensional
243 if (in.NbDimensions() > MAXND_FFTW)
244 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
245 int sz[MAXND_FFTW];
246 FillSizes4FFTW(in, sz);
[3412]247 fftwf_plan p = fftwf_plan_dft(in.NbDimensions(), sz,
[3359]248 (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
249 FFTW_FORWARD, FFTW_ESTIMATE);
250 if (p == NULL)
251 throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
252 fftwf_execute(p);
253 fftwf_destroy_plan(p);
254 }
255 if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
256 return;
257}
258
259/* --Methode-- */
260void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
261{
262 int rank = ckR4.CheckResize(in, out);
263 if (rank == 1) { // One dimensional transform
264 fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
265 (fftwf_complex *)out.Data(),
266 FFTW_BACKWARD, FFTW_ESTIMATE);
267 if (p == NULL)
268 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
269 fftwf_execute(p);
270 fftwf_destroy_plan(p);
271 }
272 else { // Multi dimensional
273 if (in.NbDimensions() > MAXND_FFTW)
274 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
275 int sz[MAXND_FFTW];
276 FillSizes4FFTW(in, sz);
[3412]277 fftwf_plan p = fftwf_plan_dft(in.NbDimensions(), sz,
[3359]278 (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
279 FFTW_BACKWARD, FFTW_ESTIMATE);
280 if (p == NULL)
281 throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
282 fftwf_execute(p);
283 fftwf_destroy_plan(p);
284 }
285
286 return;
287}
288
289
290void FFTWServer::FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out)
291{
292 int rank = ckR4.CheckResize(in, out);
293 if (rank == 1) { // One dimensional transform
294 fftwf_plan p = fftwf_plan_dft_r2c_1d(in.Size(), in.Data(),
[3412]295 (fftwf_complex *)out.Data(), FFTW_ESTIMATE);
[3359]296 if (p == NULL)
297 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
298 fftwf_execute(p);
299 fftwf_destroy_plan(p);
300
301
302 }
303 else { // Multi dimensional
304 if (in.NbDimensions() > MAXND_FFTW)
305 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
306 int sz[MAXND_FFTW];
307 FillSizes4FFTW(in, sz);
[3412]308 fftwf_plan p = fftwf_plan_dft_r2c(in.NbDimensions(), sz, in.Data(),
309 (fftwf_complex *)out.Data(), FFTW_ESTIMATE);
[3359]310 if (p == NULL)
311 throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
312 fftwf_execute(p);
313 fftwf_destroy_plan(p);
314 }
315 if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
316 return;
317
318}
319
320
321
322void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & out,
323 bool usoutsz)
324{
325 // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
326 // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
327 int rank = ckR4.CheckResize(in, out, usoutsz);
328 bool share = (_preserve_input) ? false : true;
329 TArray< complex<r_4> > inp(in, share);
330 if (rank == 1) { // One dimensional transform
331 fftwf_plan p = fftwf_plan_dft_c2r_1d(out.Size(), (fftwf_complex *)inp.Data(),
[3412]332 out.Data(), FFTW_ESTIMATE);
[3359]333 if (p == NULL)
334 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
335 fftwf_execute(p);
336 fftwf_destroy_plan(p);
337 }
338 else { // Multi dimensional
339 if (in.NbDimensions() > MAXND_FFTW)
340 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
341 int sz[MAXND_FFTW];
342 FillSizes4FFTW(out, sz);
[3412]343 fftwf_plan p = fftwf_plan_dft_c2r(out.NbDimensions(), sz, (fftwf_complex *)inp.Data(),
344 out.Data(), FFTW_ESTIMATE);
[3359]345 if (p == NULL)
346 throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
347 fftwf_execute(p);
348 fftwf_destroy_plan(p);
349 }
350 return;
351}
352
353
354
355/* --Methode-- */
356void FFTWServer::ReShapetoReal(TArray< complex<r_4> > const & ina, TArray< r_4 > & outa,
357 bool usz)
358{
359 TVector< complex<r_4> > in(ina);
360 TVector< r_4> out(outa);
361 int n = in.NElts();
362 r_4 thr = FFTArrayChecker<r_4>::ZeroThreshold();
363 sa_size_t ncs;
364 if (usz) {
365 if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
366 throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
367 ncs = out.NElts();
368 }
369 else {
370 ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
371 ncs = 2*n-1 : ncs = 2*n-2;
372
373 if (out.NElts() != ncs) {
374 cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
375 << " out.NElts()= " << out.NElts() << endl;
376 throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
377 }
378 }
379 // if (ncs == 2*n-2) {
380 out(0) = in(0).real();
381 for(int k=1; k<n; k++) {
382 out(k) = in(k).real();
383 out(ncs-k) = in(k).imag();
384 }
385 if (ncs == 2*n-2)
386 out(n-1) = in(n-1).real();
387
388 // for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
389}
390
391
392/* --Methode-- */
393void FFTWServer::ReShapetoCompl(TArray< r_4 > const & ina, TArray< complex<r_4> > & outa)
394{
395 TVector< r_4> in(ina);
396 TVector< complex<r_4> > out(outa);
397
398 sa_size_t n = in.NElts();
399 sa_size_t ncs = n/2+1;
400 sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
401 if (out.NElts() != ncs)
402 throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
403
404 out(0) = complex<r_4> (in(0),0.);
405 for(int k=1; k<ncs; k++)
406 out(k) = complex<r_4>(in(k), in(n-k));
407 if (n%2 == 0) out(ncs-1) = complex<r_4>(in(n/2), 0.);
408
409}
410
411#endif
Note: See TracBrowser for help on using the repository browser.