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