Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 = nullptr;
58 fOut = nullptr;
59 fPlan = nullptr;
60 fN = nullptr;
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 = nullptr;
78 }
79 fN = new Int_t[1];
80 fN[0] = n;
81 fTotalSize = n;
82 fNdim = 1;
83 fPlan = nullptr;
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 = nullptr;
110 }
111 fPlan = nullptr;
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 = nullptr;
123 fftw_free(fIn);
124 fIn = nullptr;
125 if (fOut)
126 fftw_free((fftw_complex*)fOut);
127 fOut = nullptr;
128 if (fN)
129 delete [] fN;
130 fN = nullptr;
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 = nullptr;
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
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
222/// `#ipoint` 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 nullptr;
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:45
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
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"
void Init(Option_t *flags, Int_t, const Int_t *) override
Creates the fftw-plan.
void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const override
Returns the point #ipoint.
void SetPoints(const Double_t *data) override
Set all input points.
Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const override
Returns the input array// One of the interface classes to the FFTW package, can be used directly or v...
Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const override
Returns the real part of the point #ipoint from the output or the point #ipoint from the input.
void SetPointComplex(Int_t ipoint, TComplex &c) override
Sets the point #ipoint (only the real part of the argument is taken)
void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const override
Fills the argument arrays with the real and imaginary parts of the computed transform.
void SetPointsComplex(const Double_t *re, const Double_t *im) override
Set all points. Only the real array is used.
void SetPoint(Int_t ipoint, Double_t re, Double_t im=0) override
Set the point #ipoint
void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const override
Fills the array data with the computed transform.
void Transform() override
Computes the transform, specified in Init() function.
~TFFTRealComplex() override
Destroys the data arrays and the plan.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:962
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:976
Basic string class.
Definition TString.h:139
void ToUpper()
Change string to upper case.
Definition TString.cxx:1195
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
const Int_t n
Definition legend1.C:16