Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFFTComplexReal.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/// \class TFFTComplexReal
14///
15/// One of the interface classes to the FFTW package, can be used directly
16/// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented.
17///
18/// Computes the inverse of the real-to-complex transforms (class TFFTRealComplex)
19/// taking complex input (storing the non-redundant half of a logically Hermitian array)
20/// to real output (see FFTW manual for more details)
21///
22/// How to use it:
23/// 1. Create an instance of TFFTComplexReal - this will allocate input and output
24/// arrays (unless an in-place transform is specified)
25/// 2. Run the Init() function with the desired flags and settings
26/// 3. Set the data (via SetPoints(), SetPoint() or SetPointComplex() functions)
27/// 4. Run the Transform() function
28/// 5. Get the output (via GetPoints(), GetPoint() or GetPointReal() functions)
29/// 6. Repeat steps 3)-5) as needed
30///
31/// For a transform of the same size, but with different flags, rerun the Init()
32/// function and continue with steps 3)-5)
33///
34/// NOTE:
35/// 1. running Init() function will overwrite the input array! Don't set any data
36/// before running the Init() function
37/// 2. FFTW computes unnormalized transform, so doing a transform followed by
38/// its inverse will lead to the original array scaled by the transform size
39/// 3. In Complex to Real transform the input array is destroyed. It cannot then
40/// be retrieved when using the Get's methods.
41///
42////////////////////////////////////////////////////////////////////////////////
43
44#include "TFFTComplexReal.h"
45#include "fftw3.h"
46#include "TComplex.h"
47
48
50
51////////////////////////////////////////////////////////////////////////////////
52///default
53
55{
56 fIn = nullptr;
57 fOut = nullptr;
58 fPlan = nullptr;
59 fN = nullptr;
60 fTotalSize = 0;
61 fNdim = 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 if (!inPlace){
71 fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1));
72 fOut = fftw_malloc(sizeof(Double_t)*n);
73
74 } else {
75 fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1));
76 fOut = nullptr;
77 }
78 fNdim = 1;
79 fN = new Int_t[1];
80 fN[0] = n;
81 fPlan = nullptr;
82 fTotalSize = n;
83}
84
85////////////////////////////////////////////////////////////////////////////////
86///For ndim-dimensional transforms
87///Second argument contains sizes of the transform in each dimension
88
90{
91 fNdim = ndim;
92 fTotalSize = 1;
93 fN = new Int_t[fNdim];
94 for (Int_t i=0; i<fNdim; i++){
95 fN[i] = n[i];
96 fTotalSize*=n[i];
97 }
98 Int_t sizein = Int_t(Double_t(fTotalSize)*(n[ndim-1]/2+1)/n[ndim-1]);
99 if (!inPlace){
100 fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
101 fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
102 } else {
103 fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
104 fOut = nullptr;
105 }
106 fPlan = nullptr;
107}
108
109
110////////////////////////////////////////////////////////////////////////////////
111///Destroys the data arrays and the plan. However, some plan information stays around
112///until the root session is over, and is reused if other plans of the same size are
113///created
114
116{
117 fftw_destroy_plan((fftw_plan)fPlan);
118 fPlan = nullptr;
119 fftw_free((fftw_complex*)fIn);
120 if (fOut)
121 fftw_free(fOut);
122 fIn = nullptr;
123 fOut = nullptr;
124 if (fN)
125 delete [] fN;
126 fN = nullptr;
127}
128
129////////////////////////////////////////////////////////////////////////////////
130///Creates the fftw-plan
131///
132///NOTE: input and output arrays are overwritten during initialisation,
133/// so don't set any points, before running this function!!!!!
134///
135///Arguments sign and kind are dummy and not need to be specified
136///Possible flag_options:
137///
138/// - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
139/// performance
140/// - "M" (from "measure") - some time spend in finding the optimal way to do the transform
141/// - "P" (from "patient") - more time spend in finding the optimal way to do the transform
142/// - "EX" (from "exhaustive") - the most optimal way is found
143///
144///This option should be chosen depending on how many transforms of the same size and
145///type are going to be done. Planning is only done once, for the first transform of this
146///size and type.
147
148void TFFTComplexReal::Init( Option_t *flags, Int_t /*sign*/,const Int_t* /*kind*/)
149{
150 fFlags = flags;
151
152 if (fPlan)
153 fftw_destroy_plan((fftw_plan)fPlan);
154 fPlan = nullptr;
155
156 if (fOut)
157 fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN,(fftw_complex*)fIn,(Double_t*)fOut, MapFlag(flags));
158 else
159 fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN, (fftw_complex*)fIn, (Double_t*)fIn, MapFlag(flags));
160}
161
162////////////////////////////////////////////////////////////////////////////////
163///Computes the transform, specified in Init() function
164
166{
167 if (fPlan)
168 fftw_execute((fftw_plan)fPlan);
169 else {
170 Error("Transform", "transform was not initialized");
171 return;
172 }
173}
174
175////////////////////////////////////////////////////////////////////////////////
176///Fills the argument array with the computed transform
177/// Works only for output (input array is destroyed in a C2R transform)
178
180{
181 if (fromInput){
182 Error("GetPoints", "Input array has been destroyed");
183 return;
184 }
185 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
186 std::copy(array,array+fTotalSize, data);
187}
188
189////////////////////////////////////////////////////////////////////////////////
190///Returns the point `#ipoint`
191/// Works only for output (input array is destroyed in a C2R transform)
192
194{
195 if (fromInput){
196 Error("GetPointReal", "Input array has been destroyed");
197 return 0;
198 }
199 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
200 return array[ipoint];
201}
202
203////////////////////////////////////////////////////////////////////////////////
204///For multidimensional transforms. Returns the point `#ipoint`
205/// Works only for output (input array is destroyed in a C2R transform)
206
207Double_t TFFTComplexReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
208{
209 if (fromInput){
210 Error("GetPointReal", "Input array has been destroyed");
211 return 0;
212 }
213 Int_t ireal = ipoint[0];
214 for (Int_t i=0; i<fNdim-1; i++)
215 ireal=fN[i+1]*ireal + ipoint[i+1];
216
217 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
218 return array[ireal];
219}
220
221
222////////////////////////////////////////////////////////////////////////////////
223/// Works only for output (input array is destroyed in a C2R transform)
224
225void TFFTComplexReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
226{
227 if (fromInput){
228 Error("GetPointComplex", "Input array has been destroyed");
229 return;
230 }
231 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
232 re = array[ipoint];
233 im = 0;
234}
235
236////////////////////////////////////////////////////////////////////////////////
237///For multidimensional transforms. Returns the point `#ipoint`
238/// Works only for output (input array is destroyed in a C2R transform)
239
240void TFFTComplexReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
241{
242 if (fromInput){
243 Error("GetPointComplex", "Input array has been destroyed");
244 return;
245 }
246 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
247
248 Int_t ireal = ipoint[0];
249 for (Int_t i=0; i<fNdim-1; i++)
250 ireal=fN[i+1]*ireal + ipoint[i+1];
251
252 re = array[ireal];
253 im = 0;
254}
255////////////////////////////////////////////////////////////////////////////////
256///Returns the array of computed transform
257/// Works only for output (input array is destroyed in a C2R transform)
258
260{
261 // we have 2 different cases
262 // fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
263 // fromInput = false; fOut = NULL (transformed is in place) : return fIn
264
265 if (fromInput) {
266 Error("GetPointsReal","Input array was destroyed");
267 return nullptr;
268 }
269 return (Double_t*) ( (fOut) ? fOut : fIn );
270}
271
272////////////////////////////////////////////////////////////////////////////////
273///Fills the argument array with the computed transform
274/// Works only for output (input array is destroyed in a C2R transform)
275
277{
278 if (fromInput){
279 Error("GetPointsComplex", "Input array has been destroyed");
280 return;
281 }
282 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
283 for (Int_t i=0; i<fTotalSize; i++){
284 re[i] = array[i];
285 im[i] = 0;
286 }
287}
288
289////////////////////////////////////////////////////////////////////////////////
290///Fills the argument array with the computed transform.
291/// Works only for output (input array is destroyed in a C2R transform)
292
294{
295 if (fromInput){
296 Error("GetPointsComplex", "Input array has been destroyed");
297 return;
298 }
299 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
300 for (Int_t i=0; i<fTotalSize; i+=2){
301 data[i] = array[i/2];
302 data[i+1]=0;
303 }
304}
305
306////////////////////////////////////////////////////////////////////////////////
307///since the input must be complex-Hermitian, if the ipoint > n/2, the according
308///point before n/2 is set to (re, -im)
309
311{
312 if (ipoint <= fN[0]/2){
313 ((fftw_complex*)fIn)[ipoint][0] = re;
314 ((fftw_complex*)fIn)[ipoint][1] = im;
315 } else {
316 ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][0] = re;
317 ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][1] = -im;
318 }
319}
320
321////////////////////////////////////////////////////////////////////////////////
322///Set the point `#ipoint`. Since the input is Hermitian, only the first (roughly) half of
323///the points have to be set.
324
326{
327 Int_t ireal = ipoint[0];
328 for (Int_t i=0; i<fNdim-2; i++){
329 ireal=fN[i+1]*ireal + ipoint[i+1];
330 }
331 ireal = (fN[fNdim-1]/2+1)*ireal+ipoint[fNdim-1];
332 Int_t realN = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
333
334 if (ireal > realN){
335 Error("SetPoint", "Illegal index value");
336 return;
337 }
338 ((fftw_complex*)fIn)[ireal][0] = re;
339 ((fftw_complex*)fIn)[ireal][1] = im;
340}
341
342////////////////////////////////////////////////////////////////////////////////
343///since the input must be complex-Hermitian, if the ipoint > n/2, the according
344///point before n/2 is set to (re, -im)
345
347{
348 if (ipoint <= fN[0]/2){
349 ((fftw_complex*)fIn)[ipoint][0] = c.Re();
350 ((fftw_complex*)fIn)[ipoint][1] = c.Im();
351 } else {
352 ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][0] = c.Re();
353 ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][1] = -c.Im();
354 }
355}
356
357////////////////////////////////////////////////////////////////////////////////
358///set all points. the values are copied. points should be ordered as follows:
359///[re_0, im_0, re_1, im_1, ..., re_n, im_n)
360
362{
363 Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
364
365 for (Int_t i=0; i<2*(sizein); i+=2){
366 ((fftw_complex*)fIn)[i/2][0]=data[i];
367 ((fftw_complex*)fIn)[i/2][1]=data[i+1];
368 }
369}
370
371////////////////////////////////////////////////////////////////////////////////
372///Set all points. The values are copied.
373
375{
376 Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
377 for (Int_t i=0; i<sizein; i++){
378 ((fftw_complex*)fIn)[i][0]=re[i];
379 ((fftw_complex*)fIn)[i][1]=im[i];
380 }
381}
382
383////////////////////////////////////////////////////////////////////////////////
384///allowed options:
385///"ES" - FFTW_ESTIMATE
386///"M" - FFTW_MEASURE
387///"P" - FFTW_PATIENT
388///"EX" - FFTW_EXHAUSTIVE
389
391{
392 TString opt = flag;
393 opt.ToUpper();
394 if (opt.Contains("ES"))
395 return FFTW_ESTIMATE;
396 if (opt.Contains("M"))
397 return FFTW_MEASURE;
398 if (opt.Contains("P"))
399 return FFTW_PATIENT;
400 if (opt.Contains("EX"))
401 return FFTW_EXHAUSTIVE;
402 return FFTW_ESTIMATE;
403}
#define c(i)
Definition RSha256.hxx:101
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
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.
void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const override
Fills the argument array with the computed transform Works only for output (input array is destroyed ...
void SetPointComplex(Int_t ipoint, TComplex &c) override
since the input must be complex-Hermitian, if the ipoint > n/2, the according point before n/2 is set...
void SetPointsComplex(const Double_t *re, const Double_t *im) override
Set all points. The values are copied.
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE
void Transform() override
Computes the transform, specified in Init() function.
void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const override
Fills the argument array with the computed transform Works only for output (input array is destroyed ...
~TFFTComplexReal() override
Destroys the data arrays and the plan.
Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const override
Returns the array of computed transform Works only for output (input array is destroyed in a C2R tran...
void SetPoints(const Double_t *data) override
set all points.
void SetPoint(Int_t ipoint, Double_t re, Double_t im=0) override
since the input must be complex-Hermitian, if the ipoint > n/2, the according point before n/2 is set...
void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const override
Works only for output (input array is destroyed in a C2R transform)
Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const override
Returns the point #ipoint Works only for output (input array is destroyed in a C2R transform)
void Init(Option_t *flags, Int_t, const Int_t *) override
Creates the fftw-plan.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
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