Logo ROOT   6.16/01
Reference Guide
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//
14// TFFTComplexReal
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 the inverse of the real-to-complex transforms (class TFFTRealComplex)
20// taking complex input (storing the non-redundant half of a logically Hermitian array)
21// to real output (see FFTW manual for more details)
22//
23// How to use it:
24// 1) Create an instance of TFFTComplexReal - this will allocate input and output
25// arrays (unless an in-place transform is specified)
26// 2) Run the Init() function with the desired flags and settings
27// 3) Set the data (via SetPoints(), SetPoint() or SetPointComplex() functions)
28// 4) Run the Transform() function
29// 5) Get the output (via GetPoints(), GetPoint() or GetPointReal() functions)
30// 6) Repeat steps 3)-5) as needed
31//
32// For a transform of the same size, but with different flags, rerun the Init()
33// function and continue with steps 3)-5)
34// NOTE: 1) running Init() function will overwrite the input array! Don't set any data
35// before running the Init() function
36// 2) FFTW computes unnormalized transform, so doing a transform followed by
37// its inverse will lead to the original array scaled by the transform size
38//
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 = 0;
57 fOut = 0;
58 fPlan = 0;
59 fN = 0;
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 = 0;
77 }
78 fNdim = 1;
79 fN = new Int_t[1];
80 fN[0] = n;
81 fPlan = 0;
82 fTotalSize = n;
83}
84
85////////////////////////////////////////////////////////////////////////////////
86///For ndim-dimensional transforms
87///Second argurment 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 = 0;
105 }
106 fPlan = 0;
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 = 0;
119 fftw_free((fftw_complex*)fIn);
120 if (fOut)
121 fftw_free(fOut);
122 fIn = 0;
123 fOut = 0;
124 if (fN)
125 delete [] fN;
126 fN = 0;
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///"ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
138/// performanc
139///"M" (from "measure") - some time spend in finding the optimal way to do the transform
140///"P" (from "patient") - more time spend in finding the optimal way to do the transform
141///"EX" (from "exhaustive") - the most optimal way is found
142///This option should be chosen depending on how many transforms of the same size and
143///type are going to be done. Planning is only done once, for the first transform of this
144///size and type.
145
146void TFFTComplexReal::Init( Option_t *flags, Int_t /*sign*/,const Int_t* /*kind*/)
147{
148 fFlags = flags;
149
150 if (fPlan)
151 fftw_destroy_plan((fftw_plan)fPlan);
152 fPlan = 0;
153
154 if (fOut)
155 fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN,(fftw_complex*)fIn,(Double_t*)fOut, MapFlag(flags));
156 else
157 fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN, (fftw_complex*)fIn, (Double_t*)fIn, MapFlag(flags));
158}
159
160////////////////////////////////////////////////////////////////////////////////
161///Computes the transform, specified in Init() function
162
164{
165 if (fPlan)
166 fftw_execute((fftw_plan)fPlan);
167 else {
168 Error("Transform", "transform was not initialized");
169 return;
170 }
171}
172
173////////////////////////////////////////////////////////////////////////////////
174///Fills the argument array with the computed transform
175/// Works only for output (input array is destroyed in a C2R transform)
176
178{
179 if (fromInput){
180 Error("GetPoints", "Input array has been destroyed");
181 return;
182 }
183 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
184 std::copy(array,array+fTotalSize, data);
185}
186
187////////////////////////////////////////////////////////////////////////////////
188///Returns the point #ipoint
189/// Works only for output (input array is destroyed in a C2R transform)
190
192{
193 if (fromInput){
194 Error("GetPointReal", "Input array has been destroyed");
195 return 0;
196 }
197 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
198 return array[ipoint];
199}
200
201////////////////////////////////////////////////////////////////////////////////
202///For multidimensional transforms. Returns the point #ipoint
203/// Works only for output (input array is destroyed in a C2R transform)
204
205Double_t TFFTComplexReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
206{
207 if (fromInput){
208 Error("GetPointReal", "Input array has been destroyed");
209 return 0;
210 }
211 Int_t ireal = ipoint[0];
212 for (Int_t i=0; i<fNdim-1; i++)
213 ireal=fN[i+1]*ireal + ipoint[i+1];
214
215 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
216 return array[ireal];
217}
218
219
220////////////////////////////////////////////////////////////////////////////////
221/// Works only for output (input array is destroyed in a C2R transform)
222
223void TFFTComplexReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
224{
225 if (fromInput){
226 Error("GetPointComplex", "Input array has been destroyed");
227 return;
228 }
229 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
230 re = array[ipoint];
231 im = 0;
232}
233
234////////////////////////////////////////////////////////////////////////////////
235///For multidimensional transforms. Returns the point #ipoint
236/// Works only for output (input array is destroyed in a C2R transform)
237
238void TFFTComplexReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
239{
240 if (fromInput){
241 Error("GetPointComplex", "Input array has been destroyed");
242 return;
243 }
244 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
245
246 Int_t ireal = ipoint[0];
247 for (Int_t i=0; i<fNdim-1; i++)
248 ireal=fN[i+1]*ireal + ipoint[i+1];
249
250 re = array[ireal];
251 im = 0;
252}
253////////////////////////////////////////////////////////////////////////////////
254///Returns the array of computed transform
255/// Works only for output (input array is destroyed in a C2R transform)
256
258{
259 // we have 2 different cases
260 // fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
261 // fromInput = false; fOut = NULL (transformed is in place) : return fIn
262
263 if (fromInput) {
264 Error("GetPointsReal","Input array was destroyed");
265 return 0;
266 }
267 return (Double_t*) ( (fOut) ? fOut : fIn );
268}
269
270////////////////////////////////////////////////////////////////////////////////
271///Fills the argument array with the computed transform
272/// Works only for output (input array is destroyed in a C2R transform)
273
275{
276 if (fromInput){
277 Error("GetPointsComplex", "Input array has been destroyed");
278 return;
279 }
280 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
281 for (Int_t i=0; i<fTotalSize; i++){
282 re[i] = array[i];
283 im[i] = 0;
284 }
285}
286
287////////////////////////////////////////////////////////////////////////////////
288///Fills the argument array with the computed transform.
289/// Works only for output (input array is destroyed in a C2R transform)
290
292{
293 if (fromInput){
294 Error("GetPointsComplex", "Input array has been destroyed");
295 return;
296 }
297 const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
298 for (Int_t i=0; i<fTotalSize; i+=2){
299 data[i] = array[i/2];
300 data[i+1]=0;
301 }
302}
303
304////////////////////////////////////////////////////////////////////////////////
305///since the input must be complex-Hermitian, if the ipoint > n/2, the according
306///point before n/2 is set to (re, -im)
307
309{
310 if (ipoint <= fN[0]/2){
311 ((fftw_complex*)fIn)[ipoint][0] = re;
312 ((fftw_complex*)fIn)[ipoint][1] = im;
313 } else {
314 ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][0] = re;
315 ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][1] = -im;
316 }
317}
318
319////////////////////////////////////////////////////////////////////////////////
320///Set the point #ipoint. Since the input is Hermitian, only the first (roughly)half of
321///the points have to be set.
322
324{
325 Int_t ireal = ipoint[0];
326 for (Int_t i=0; i<fNdim-2; i++){
327 ireal=fN[i+1]*ireal + ipoint[i+1];
328 }
329 ireal = (fN[fNdim-1]/2+1)*ireal+ipoint[fNdim-1];
330 Int_t realN = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
331
332 if (ireal > realN){
333 Error("SetPoint", "Illegal index value");
334 return;
335 }
336 ((fftw_complex*)fIn)[ireal][0] = re;
337 ((fftw_complex*)fIn)[ireal][1] = im;
338}
339
340////////////////////////////////////////////////////////////////////////////////
341///since the input must be complex-Hermitian, if the ipoint > n/2, the according
342///point before n/2 is set to (re, -im)
343
345{
346 if (ipoint <= fN[0]/2){
347 ((fftw_complex*)fIn)[ipoint][0] = c.Re();
348 ((fftw_complex*)fIn)[ipoint][1] = c.Im();
349 } else {
350 ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][0] = c.Re();
351 ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][1] = -c.Im();
352 }
353}
354
355////////////////////////////////////////////////////////////////////////////////
356///set all points. the values are copied. points should be ordered as follows:
357///[re_0, im_0, re_1, im_1, ..., re_n, im_n)
358
360{
361 Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
362
363 for (Int_t i=0; i<2*(sizein); i+=2){
364 ((fftw_complex*)fIn)[i/2][0]=data[i];
365 ((fftw_complex*)fIn)[i/2][1]=data[i+1];
366 }
367}
368
369////////////////////////////////////////////////////////////////////////////////
370///Set all points. The values are copied.
371
373{
374 Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
375 for (Int_t i=0; i<sizein; i++){
376 ((fftw_complex*)fIn)[i][0]=re[i];
377 ((fftw_complex*)fIn)[i][1]=im[i];
378 }
379}
380
381////////////////////////////////////////////////////////////////////////////////
382///allowed options:
383///"ES" - FFTW_ESTIMATE
384///"M" - FFTW_MEASURE
385///"P" - FFTW_PATIENT
386///"EX" - FFTW_EXHAUSTIVE
387
389{
390 TString opt = flag;
391 opt.ToUpper();
392 if (opt.Contains("ES"))
393 return FFTW_ESTIMATE;
394 if (opt.Contains("M"))
395 return FFTW_MEASURE;
396 if (opt.Contains("P"))
397 return FFTW_PATIENT;
398 if (opt.Contains("EX"))
399 return FFTW_EXHAUSTIVE;
400 return FFTW_ESTIMATE;
401}
#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 char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:363
virtual Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const
Returns the array of computed transform Works only for output (input array is destroyed in a C2R tran...
virtual void SetPointComplex(Int_t ipoint, TComplex &c)
since the input must be complex-Hermitian, if the ipoint > n/2, the according point before n/2 is set...
virtual void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const
Fills the argument array with the computed transform Works only for output (input array is destroyed ...
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const
Returns the point #ipoint Works only for output (input array is destroyed in a C2R transform)
virtual void SetPoints(const Double_t *data)
set all points.
virtual void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const
Works only for output (input array is destroyed in a C2R transform)
virtual void Transform()
Computes the transform, specified in Init() function.
virtual void Init(Option_t *flags, Int_t, const Int_t *)
Creates the fftw-plan.
TFFTComplexReal()
default
virtual void SetPointsComplex(const Double_t *re, const Double_t *im)
Set all points. The values are copied.
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
since the input must be complex-Hermitian, if the ipoint > n/2, the according point before n/2 is set...
virtual ~TFFTComplexReal()
Destroys the data arrays and the plan.
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const
Fills the argument array with the computed transform Works only for output (input array is destroyed ...
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