Logo ROOT  
Reference Guide
TFFTRealComplex.cxx
Go to the documentation of this file.
1// @(#)root/fft:$Id$
2// Author: Anna Kreshuk 07/4/2006
3
4/*************************************************************************
5 * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12////////////////////////////////////////////////////////////////////////////////
13///
14/// \class TFFTRealComplex
15///
16/// One of the interface classes to the FFTW package, can be used directly
17/// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented.
18///
19/// Computes a real input/complex output discrete Fourier transform in 1 or more
20/// dimensions. However, only out-of-place transforms are now supported for transforms
21/// in more than 1 dimension. For detailed information about the computed transforms,
22/// please refer to the FFTW manual
23///
24/// How to use it:
25/// 1. Create an instance of TFFTRealComplex - this will allocate input and output
26/// arrays (unless an in-place transform is specified)
27/// 2. Run the Init() function with the desired flags and settings (see function
28/// comments for possible kind parameters)
29/// 3. Set the data (via SetPoints()or SetPoint() functions)
30/// 4. Run the Transform() function
31/// 5. Get the output (via GetPoints() or GetPoint() functions)
32/// 6. Repeat steps 3)-5) as needed
33/// For a transform of the same size, but with different flags,
34/// rerun the Init() function and continue with steps 3)-5)
35///
36/// NOTE:
37/// 1. running Init() function will overwrite the input array! Don't set any data
38/// before running the Init() function
39/// 2. FFTW computes unnormalized transform, so doing a transform followed by
40/// its inverse will lead to the original array scaled by the transform size
41///
42///
43/////////////////////////////////////////////////////////////////////////////////
44
45#include "TFFTRealComplex.h"
46#include "fftw3.h"
47#include "TComplex.h"
48
49
51
52////////////////////////////////////////////////////////////////////////////////
53///default
54
56{
57 fIn = 0;
58 fOut = 0;
59 fPlan = 0;
60 fN = 0;
61 fNdim = 0;
62 fTotalSize = 0;
63}
64
65////////////////////////////////////////////////////////////////////////////////
66///For 1d transforms
67///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array
68
70{
71
72 if (!inPlace){
73 fIn = fftw_malloc(sizeof(Double_t)*n);
74 fOut = fftw_malloc(sizeof(fftw_complex)*(n/2+1));
75 } else {
76 fIn = fftw_malloc(sizeof(Double_t)*(2*(n/2+1)));
77 fOut = 0;
78 }
79 fN = new Int_t[1];
80 fN[0] = n;
81 fTotalSize = n;
82 fNdim = 1;
83 fPlan = 0;
84}
85
86////////////////////////////////////////////////////////////////////////////////
87///For ndim-dimensional transforms
88///Second argument contains sizes of the transform in each dimension
89
91{
92 if (ndim>1 && inPlace==kTRUE){
93 Error("TFFTRealComplex", "multidimensional in-place r2c transforms are not implemented yet");
94 return;
95 }
96 fNdim = ndim;
97 fTotalSize = 1;
98 fN = new Int_t[fNdim];
99 for (Int_t i=0; i<fNdim; i++){
100 fN[i] = n[i];
101 fTotalSize*=n[i];
102 }
103 Int_t sizeout = Int_t(Double_t(fTotalSize)*(n[ndim-1]/2+1)/n[ndim-1]);
104 if (!inPlace){
105 fIn = fftw_malloc(sizeof(Double_t)*fTotalSize);
106 fOut = fftw_malloc(sizeof(fftw_complex)*sizeout);
107 } else {
108 fIn = fftw_malloc(sizeof(Double_t)*(2*sizeout));
109 fOut = 0;
110 }
111 fPlan = 0;
112}
113
114////////////////////////////////////////////////////////////////////////////////
115///Destroys the data arrays and the plan. However, some plan information stays around
116///until the root session is over, and is reused if other plans of the same size are
117///created
118
120{
121 fftw_destroy_plan((fftw_plan)fPlan);
122 fPlan = 0;
123 fftw_free(fIn);
124 fIn = 0;
125 if (fOut)
126 fftw_free((fftw_complex*)fOut);
127 fOut = 0;
128 if (fN)
129 delete [] fN;
130 fN = 0;
131}
132
133////////////////////////////////////////////////////////////////////////////////
134///Creates the fftw-plan
135///
136///NOTE: input and output arrays are overwritten during initialisation,
137/// so don't set any points, before running this function!!!!!
138///
139///Arguments sign and kind are dummy and not need to be specified
140///Possible flag_options:
141///
142/// - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
143/// performance
144/// - "M" (from "measure") - some time spend in finding the optimal way to do the transform
145/// - "P" (from "patient") - more time spend in finding the optimal way to do the transform
146/// - "EX" (from "exhaustive") - the most optimal way is found
147///
148///This option should be chosen depending on how many transforms of the same size and
149///type are going to be done. Planning is only done once, for the first transform of this
150///size and type.
151
152void TFFTRealComplex::Init(Option_t *flags,Int_t /*sign*/, const Int_t* /*kind*/)
153{
154 fFlags = flags;
155
156 if (fPlan)
157 fftw_destroy_plan((fftw_plan)fPlan);
158 fPlan = 0;
159
160 if (fOut)
161 fPlan = (void*)fftw_plan_dft_r2c(fNdim, fN, (Double_t*)fIn, (fftw_complex*)fOut,MapFlag(flags));
162 else
163 fPlan = (void*)fftw_plan_dft_r2c(fNdim, fN, (Double_t*)fIn, (fftw_complex*)fIn, MapFlag(flags));
164}
165
166////////////////////////////////////////////////////////////////////////////////
167///Computes the transform, specified in Init() function
168
170{
171
172 if (fPlan){
173 fftw_execute((fftw_plan)fPlan);
174 }
175 else {
176 Error("Transform", "transform hasn't been initialised");
177 return;
178 }
179}
180
181////////////////////////////////////////////////////////////////////////////////
182///Fills the array data with the computed transform.
183///Only (roughly) a half of the transform is copied (exactly the output of FFTW),
184///the rest being Hermitian symmetric with the first half
185
186void TFFTRealComplex::GetPoints(Double_t *data, Bool_t fromInput) const
187{
188 if (fromInput){
189 for (Int_t i=0; i<fTotalSize; i++)
190 data[i] = ((Double_t*)fIn)[i];
191 } else {
192 Int_t realN = 2*Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
193 if (fOut){
194 for (Int_t i=0; i<realN; i+=2){
195 data[i] = ((fftw_complex*)fOut)[i/2][0];
196 data[i+1] = ((fftw_complex*)fOut)[i/2][1];
197 }
198 }
199 else {
200 for (Int_t i=0; i<realN; i++)
201 data[i] = ((Double_t*)fIn)[i];
202 }
203 }
204}
205
206////////////////////////////////////////////////////////////////////////////////
207///Returns the real part of the point #ipoint from the output or the point #ipoint
208///from the input
209
211{
212 if (fOut && !fromInput){
213 Warning("GetPointReal", "Output is complex. Only real part returned");
214 return ((fftw_complex*)fOut)[ipoint][0];
215 }
216 else
217 return ((Double_t*)fIn)[ipoint];
218}
219
220////////////////////////////////////////////////////////////////////////////////
221///Returns the real part of the point #ipoint from the output or the point #ipoint
222///from the input
223
224Double_t TFFTRealComplex::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
225{
226 Int_t ireal = ipoint[0];
227 for (Int_t i=0; i<fNdim-1; i++)
228 ireal=fN[i+1]*ireal + ipoint[i+1];
229
230 if (fOut && !fromInput){
231 Warning("GetPointReal", "Output is complex. Only real part returned");
232 return ((fftw_complex*)fOut)[ireal][0];
233 }
234 else
235 return ((Double_t*)fIn)[ireal];
236}
237
238
239////////////////////////////////////////////////////////////////////////////////
240///Returns the point #ipoint.
241///For 1d, if ipoint > fN/2+1 (the point is in the Hermitian symmetric part), it is still
242///returned. For >1d, only the first (roughly)half of points can be returned
243///For 2d, see function GetPointComplex(Int_t *ipoint,...)
244
245void TFFTRealComplex::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
246{
247 if (fromInput){
248 re = ((Double_t*)fIn)[ipoint];
249 } else {
250 if (fNdim==1){
251 if (fOut){
252 if (ipoint<fN[0]/2+1){
253 re = ((fftw_complex*)fOut)[ipoint][0];
254 im = ((fftw_complex*)fOut)[ipoint][1];
255 } else {
256 re = ((fftw_complex*)fOut)[fN[0]-ipoint][0];
257 im = -((fftw_complex*)fOut)[fN[0]-ipoint][1];
258 }
259 } else {
260 if (ipoint<fN[0]/2+1){
261 re = ((Double_t*)fIn)[2*ipoint];
262 im = ((Double_t*)fIn)[2*ipoint+1];
263 } else {
264 re = ((Double_t*)fIn)[2*(fN[0]-ipoint)];
265 im = ((Double_t*)fIn)[2*(fN[0]-ipoint)+1];
266 }
267 }
268 }
269 else {
270 Int_t realN = 2*Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
271 if (ipoint>realN/2){
272 Error("GetPointComplex", "Illegal index value");
273 return;
274 }
275 if (fOut){
276 re = ((fftw_complex*)fOut)[ipoint][0];
277 im = ((fftw_complex*)fOut)[ipoint][1];
278 } else {
279 re = ((Double_t*)fIn)[2*ipoint];
280 im = ((Double_t*)fIn)[2*ipoint+1];
281 }
282 }
283 }
284}
285////////////////////////////////////////////////////////////////////////////////
286///For multidimensional transforms. Returns the point #ipoint.
287///In case of transforms of more than 2 dimensions,
288///only points from the first (roughly)half are returned, the rest being Hermitian symmetric
289
290void TFFTRealComplex::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
291{
292
293 Int_t ireal = ipoint[0];
294 for (Int_t i=0; i<fNdim-2; i++)
295 ireal=fN[i+1]*ireal + ipoint[i+1];
296 //special treatment of the last dimension
297 ireal = (fN[fNdim-1]/2+1)*ireal + ipoint[fNdim-1];
298
299 if (fromInput){
300 re = ((Double_t*)fIn)[ireal];
301 return;
302 }
303 if (fNdim==1){
304 if (fOut){
305 if (ipoint[0] <fN[0]/2+1){
306 re = ((fftw_complex*)fOut)[ipoint[0]][0];
307 im = ((fftw_complex*)fOut)[ipoint[0]][1];
308 } else {
309 re = ((fftw_complex*)fOut)[fN[0]-ipoint[0]][0];
310 im = -((fftw_complex*)fOut)[fN[0]-ipoint[0]][1];
311 }
312 } else {
313 if (ireal <fN[0]/2+1){
314 re = ((Double_t*)fIn)[2*ipoint[0]];
315 im = ((Double_t*)fIn)[2*ipoint[0]+1];
316 } else {
317 re = ((Double_t*)fIn)[2*(fN[0]-ipoint[0])];
318 im = ((Double_t*)fIn)[2*(fN[0]-ipoint[0])+1];
319 }
320 }
321 }
322 else if (fNdim==2){
323 if (fOut){
324 if (ipoint[1]<fN[1]/2+1){
325 re = ((fftw_complex*)fOut)[ipoint[1]+(fN[1]/2+1)*ipoint[0]][0];
326 im = ((fftw_complex*)fOut)[ipoint[1]+(fN[1]/2+1)*ipoint[0]][1];
327 } else {
328 if (ipoint[0]==0){
329 re = ((fftw_complex*)fOut)[fN[1]-ipoint[1]][0];
330 im = -((fftw_complex*)fOut)[fN[1]-ipoint[1]][1];
331 } else {
332 re = ((fftw_complex*)fOut)[(fN[1]-ipoint[1])+(fN[1]/2+1)*(fN[0]-ipoint[0])][0];
333 im = -((fftw_complex*)fOut)[(fN[1]-ipoint[1])+(fN[1]/2+1)*(fN[0]-ipoint[0])][1];
334 }
335 }
336 } else {
337 if (ipoint[1]<fN[1]/2+1){
338 re = ((Double_t*)fIn)[2*(ipoint[1]+(fN[1]/2+1)*ipoint[0])];
339 im = ((Double_t*)fIn)[2*(ipoint[1]+(fN[1]/2+1)*ipoint[0])+1];
340 } else {
341 if (ipoint[0]==0){
342 re = ((Double_t*)fIn)[2*(fN[1]-ipoint[1])];
343 im = -((Double_t*)fIn)[2*(fN[1]-ipoint[1])+1];
344 } else {
345 re = ((Double_t*)fIn)[2*((fN[1]-ipoint[1])+(fN[1]/2+1)*(fN[0]-ipoint[0]))];
346 im = -((Double_t*)fIn)[2*((fN[1]-ipoint[1])+(fN[1]/2+1)*(fN[0]-ipoint[0]))+1];
347 }
348 }
349 }
350 }
351 else {
352
353 if (fOut){
354 re = ((fftw_complex*)fOut)[ireal][0];
355 im = ((fftw_complex*)fOut)[ireal][1];
356 } else {
357 re = ((Double_t*)fIn)[2*ireal];
358 im = ((Double_t*)fIn)[2*ireal+1];
359 }
360 }
361}
362
363
364////////////////////////////////////////////////////////////////////////////////
365///Returns the input array// One of the interface classes to the FFTW package, can be used directly
366/// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented.
367
369{
370 if (!fromInput){
371 Error("GetPointsReal", "Output array is complex");
372 return 0;
373 }
374 return (Double_t*)fIn;
375}
376
377////////////////////////////////////////////////////////////////////////////////
378///Fills the argument arrays with the real and imaginary parts of the computed transform.
379///Only (roughly) a half of the transform is copied, the rest being Hermitian
380///symmetric with the first half
381
383{
384 Int_t nreal;
385 if (fOut && !fromInput){
386 nreal = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
387 for (Int_t i=0; i<nreal; i++){
388 re[i] = ((fftw_complex*)fOut)[i][0];
389 im[i] = ((fftw_complex*)fOut)[i][1];
390 }
391 } else if (fromInput) {
392 for (Int_t i=0; i<fTotalSize; i++){
393 re[i] = ((Double_t *)fIn)[i];
394 im[i] = 0;
395 }
396 }
397 else {
398 nreal = 2*Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
399 for (Int_t i=0; i<nreal; i+=2){
400 re[i/2] = ((Double_t*)fIn)[i];
401 im[i/2] = ((Double_t*)fIn)[i+1];
402 }
403 }
404}
405
406////////////////////////////////////////////////////////////////////////////////
407///Fills the argument arrays with the real and imaginary parts of the computed transform.
408///Only (roughly) a half of the transform is copied, the rest being Hermitian
409///symmetric with the first half
410
412{
413 Int_t nreal;
414
415 if (fOut && !fromInput){
416 nreal = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
417
418 for (Int_t i=0; i<nreal; i+=2){
419 data[i] = ((fftw_complex*)fOut)[i/2][0];
420 data[i+1] = ((fftw_complex*)fOut)[i/2][1];
421 }
422 } else if (fromInput){
423 for (Int_t i=0; i<fTotalSize; i+=2){
424 data[i] = ((Double_t*)fIn)[i/2];
425 data[i+1] = 0;
426 }
427 }
428 else {
429
430 nreal = 2*Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
431 for (Int_t i=0; i<nreal; i++)
432 data[i] = ((Double_t*)fIn)[i];
433 }
434}
435
436////////////////////////////////////////////////////////////////////////////////
437///Set the point #ipoint
438
440{
441 ((Double_t *)fIn)[ipoint] = re;
442}
443
444////////////////////////////////////////////////////////////////////////////////
445///For multidimensional transforms. Set the point #ipoint
446
447void TFFTRealComplex::SetPoint(const Int_t *ipoint, Double_t re, Double_t /*im*/)
448{
449 Int_t ireal = ipoint[0];
450 for (Int_t i=0; i<fNdim-1; i++)
451 ireal=fN[i+1]*ireal + ipoint[i+1];
452 ((Double_t*)fIn)[ireal]=re;
453}
454
455////////////////////////////////////////////////////////////////////////////////
456///Set all input points
457
459{
460 for (Int_t i=0; i<fTotalSize; i++){
461 ((Double_t*)fIn)[i]=data[i];
462 }
463}
464
465////////////////////////////////////////////////////////////////////////////////
466///Sets the point #ipoint (only the real part of the argument is taken)
467
469{
470 ((Double_t *)fIn)[ipoint]=c.Re();
471}
472
473////////////////////////////////////////////////////////////////////////////////
474///Set all points. Only the real array is used
475
477{
478 for (Int_t i=0; i<fTotalSize; i++)
479 ((Double_t *)fIn)[i] = re[i];
480}
481
482////////////////////////////////////////////////////////////////////////////////
483///allowed options:
484///"ES"
485///"M"
486///"P"
487///"EX"
488
490{
491
492 TString opt = flag;
493 opt.ToUpper();
494 if (opt.Contains("ES"))
495 return FFTW_ESTIMATE;
496 if (opt.Contains("M"))
497 return FFTW_MEASURE;
498 if (opt.Contains("P"))
499 return FFTW_PATIENT;
500 if (opt.Contains("EX"))
501 return FFTW_EXHAUSTIVE;
502 return FFTW_ESTIMATE;
503}
#define c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
One of the interface classes to the FFTW package, can be used directly or via the TVirtualFFT class.
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" "M" "P" "EX"
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
Set the point #ipoint.
virtual void SetPoints(const Double_t *data)
Set all input points.
virtual void Init(Option_t *flags, Int_t, const Int_t *)
Creates the fftw-plan.
virtual void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const
Returns the point #ipoint.
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const
Returns the real part of the point #ipoint from the output or the point #ipoint from the input.
virtual void Transform()
Computes the transform, specified in Init() function.
virtual ~TFFTRealComplex()
Destroys the data arrays and the plan.
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const
Fills the array data with the computed transform.
virtual void SetPointComplex(Int_t ipoint, TComplex &c)
Sets the point #ipoint (only the real part of the argument is taken)
virtual void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const
Fills the argument arrays with the real and imaginary parts of the computed transform.
TFFTRealComplex()
default
virtual Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const
Returns the input array// One of the interface classes to the FFTW package, can be used directly or v...
virtual void SetPointsComplex(const Double_t *re, const Double_t *im)
Set all points. Only the real array is used.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
Basic string class.
Definition: TString.h:131
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
const Int_t n
Definition: legend1.C:16