ROOT  6.06/09
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  fFlags = 0;
62  fNdim = 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  if (!inPlace){
72  fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1));
73  fOut = fftw_malloc(sizeof(Double_t)*n);
74 
75  } else {
76  fIn = fftw_malloc(sizeof(fftw_complex)*(n/2+1));
77  fOut = 0;
78  }
79  fNdim = 1;
80  fN = new Int_t[1];
81  fN[0] = n;
82  fPlan = 0;
83  fTotalSize = n;
84  fFlags = 0;
85 }
86 
87 ////////////////////////////////////////////////////////////////////////////////
88 ///For ndim-dimensional transforms
89 ///Second argurment contains sizes of the transform in each dimension
90 
92 {
93  fNdim = ndim;
94  fTotalSize = 1;
95  fN = new Int_t[fNdim];
96  for (Int_t i=0; i<fNdim; i++){
97  fN[i] = n[i];
98  fTotalSize*=n[i];
99  }
100  Int_t sizein = Int_t(Double_t(fTotalSize)*(n[ndim-1]/2+1)/n[ndim-1]);
101  if (!inPlace){
102  fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
103  fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
104  } else {
105  fIn = fftw_malloc(sizeof(fftw_complex)*sizein);
106  fOut = 0;
107  }
108  fPlan = 0;
109  fFlags = 0;
110 }
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((fftw_complex*)fIn);
123  if (fOut)
124  fftw_free(fOut);
125  fIn = 0;
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 
149 void TFFTComplexReal::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_c2r(fNdim, fN,(fftw_complex*)fIn,(Double_t*)fOut, MapFlag(flags));
159  else
160  fPlan = (void*)fftw_plan_dft_c2r(fNdim, fN, (fftw_complex*)fIn, (Double_t*)fIn, MapFlag(flags));
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 ///Computes the transform, specified in Init() function
165 
167 {
168  if (fPlan)
169  fftw_execute((fftw_plan)fPlan);
170  else {
171  Error("Transform", "transform was not initialized");
172  return;
173  }
174 }
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 ///Fills the argument array with the computed transform
178 /// Works only for output (input array is destroyed in a C2R transform)
179 
180 void TFFTComplexReal::GetPoints(Double_t *data, Bool_t fromInput) const
181 {
182  if (fromInput){
183  Error("GetPoints", "Input array has been destroyed");
184  return;
185  }
186  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
187  std::copy(array,array+fTotalSize, data);
188 }
189 
190 ////////////////////////////////////////////////////////////////////////////////
191 ///Returns the point #ipoint
192 /// Works only for output (input array is destroyed in a C2R transform)
193 
195 {
196  if (fromInput){
197  Error("GetPointReal", "Input array has been destroyed");
198  return 0;
199  }
200  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
201  return array[ipoint];
202 }
203 
204 ////////////////////////////////////////////////////////////////////////////////
205 ///For multidimensional transforms. Returns the point #ipoint
206 /// Works only for output (input array is destroyed in a C2R transform)
207 
208 Double_t TFFTComplexReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
209 {
210  if (fromInput){
211  Error("GetPointReal", "Input array has been destroyed");
212  return 0;
213  }
214  Int_t ireal = ipoint[0];
215  for (Int_t i=0; i<fNdim-1; i++)
216  ireal=fN[i+1]*ireal + ipoint[i+1];
217 
218  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
219  return array[ireal];
220 }
221 
222 
223 ////////////////////////////////////////////////////////////////////////////////
224 /// Works only for output (input array is destroyed in a C2R transform)
225 
226 void TFFTComplexReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
227 {
228  if (fromInput){
229  Error("GetPointComplex", "Input array has been destroyed");
230  return;
231  }
232  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
233  re = array[ipoint];
234  im = 0;
235 }
236 
237 ////////////////////////////////////////////////////////////////////////////////
238 ///For multidimensional transforms. Returns the point #ipoint
239 /// Works only for output (input array is destroyed in a C2R transform)
240 
241 void TFFTComplexReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
242 {
243  if (fromInput){
244  Error("GetPointComplex", "Input array has been destroyed");
245  return;
246  }
247  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
248 
249  Int_t ireal = ipoint[0];
250  for (Int_t i=0; i<fNdim-1; i++)
251  ireal=fN[i+1]*ireal + ipoint[i+1];
252 
253  re = array[ireal];
254  im = 0;
255 }
256 ////////////////////////////////////////////////////////////////////////////////
257 ///Returns the array of computed transform
258 /// Works only for output (input array is destroyed in a C2R transform)
259 
261 {
262  // we have 2 different cases
263  // fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
264  // fromInput = false; fOut = NULL (transformed is in place) : return fIn
265 
266  if (fromInput) {
267  Error("GetPointsReal","Input array was destroyed");
268  return 0;
269  }
270  return (Double_t*) ( (fOut) ? fOut : fIn );
271 }
272 
273 ////////////////////////////////////////////////////////////////////////////////
274 ///Fills the argument array with the computed transform
275 /// Works only for output (input array is destroyed in a C2R transform)
276 
278 {
279  if (fromInput){
280  Error("GetPointsComplex", "Input array has been destroyed");
281  return;
282  }
283  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
284  for (Int_t i=0; i<fTotalSize; i++){
285  re[i] = array[i];
286  im[i] = 0;
287  }
288 }
289 
290 ////////////////////////////////////////////////////////////////////////////////
291 ///Fills the argument array with the computed transform.
292 /// Works only for output (input array is destroyed in a C2R transform)
293 
295 {
296  if (fromInput){
297  Error("GetPointsComplex", "Input array has been destroyed");
298  return;
299  }
300  const Double_t * array = (const Double_t*) ( (fOut) ? fOut : fIn );
301  for (Int_t i=0; i<fTotalSize; i+=2){
302  data[i] = array[i/2];
303  data[i+1]=0;
304  }
305 }
306 
307 ////////////////////////////////////////////////////////////////////////////////
308 ///since the input must be complex-Hermitian, if the ipoint > n/2, the according
309 ///point before n/2 is set to (re, -im)
310 
312 {
313  if (ipoint <= fN[0]/2){
314  ((fftw_complex*)fIn)[ipoint][0] = re;
315  ((fftw_complex*)fIn)[ipoint][1] = im;
316  } else {
317  ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][0] = re;
318  ((fftw_complex*)fOut)[2*(fN[0]/2)-ipoint][1] = -im;
319  }
320 }
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 ///Set the point #ipoint. Since the input is Hermitian, only the first (roughly)half of
324 ///the points have to be set.
325 
327 {
328  Int_t ireal = ipoint[0];
329  for (Int_t i=0; i<fNdim-2; i++){
330  ireal=fN[i+1]*ireal + ipoint[i+1];
331  }
332  ireal = (fN[fNdim-1]/2+1)*ireal+ipoint[fNdim-1];
333  Int_t realN = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
334 
335  if (ireal > realN){
336  Error("SetPoint", "Illegal index value");
337  return;
338  }
339  ((fftw_complex*)fIn)[ireal][0] = re;
340  ((fftw_complex*)fIn)[ireal][1] = im;
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 ///since the input must be complex-Hermitian, if the ipoint > n/2, the according
345 ///point before n/2 is set to (re, -im)
346 
348 {
349  if (ipoint <= fN[0]/2){
350  ((fftw_complex*)fIn)[ipoint][0] = c.Re();
351  ((fftw_complex*)fIn)[ipoint][1] = c.Im();
352  } else {
353  ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][0] = c.Re();
354  ((fftw_complex*)fIn)[2*(fN[0]/2)-ipoint][1] = -c.Im();
355  }
356 }
357 
358 ////////////////////////////////////////////////////////////////////////////////
359 ///set all points. the values are copied. points should be ordered as follows:
360 ///[re_0, im_0, re_1, im_1, ..., re_n, im_n)
361 
363 {
364  Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
365 
366  for (Int_t i=0; i<2*(sizein); i+=2){
367  ((fftw_complex*)fIn)[i/2][0]=data[i];
368  ((fftw_complex*)fIn)[i/2][1]=data[i+1];
369  }
370 }
371 
372 ////////////////////////////////////////////////////////////////////////////////
373 ///Set all points. The values are copied.
374 
376 {
377  Int_t sizein = Int_t(Double_t(fTotalSize)*(fN[fNdim-1]/2+1)/fN[fNdim-1]);
378  for (Int_t i=0; i<sizein; i++){
379  ((fftw_complex*)fIn)[i][0]=re[i];
380  ((fftw_complex*)fIn)[i][1]=im[i];
381  }
382 }
383 
384 ////////////////////////////////////////////////////////////////////////////////
385 ///allowed options:
386 ///"ES" - FFTW_ESTIMATE
387 ///"M" - FFTW_MEASURE
388 ///"P" - FFTW_PATIENT
389 ///"EX" - FFTW_EXHAUSTIVE
390 
392 {
393  TString opt = flag;
394  opt.ToUpper();
395  if (opt.Contains("ES"))
396  return FFTW_ESTIMATE;
397  if (opt.Contains("M"))
398  return FFTW_MEASURE;
399  if (opt.Contains("P"))
400  return FFTW_PATIENT;
401  if (opt.Contains("EX"))
402  return FFTW_EXHAUSTIVE;
403  return FFTW_ESTIMATE;
404 }
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 ...
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...
const char Option_t
Definition: RtypesCore.h:62
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) ...
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1101
Double_t Re() const
Definition: TComplex.h:46
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual void Init(Option_t *flags, Int_t, const Int_t *)
Creates the fftw-plan.
virtual ~TFFTComplexReal()
Destroys the data arrays and the plan.
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 SetPoints(const Double_t *data)
set all points.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
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 void Transform()
Computes the transform, specified in Init() function.
Option_t * fFlags
unsigned int UInt_t
Definition: RtypesCore.h:42
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE ...
virtual void SetPointsComplex(const Double_t *re, const Double_t *im)
Set all points. The values are copied.
double Double_t
Definition: RtypesCore.h:55
Double_t Im() const
Definition: TComplex.h:47
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...
ClassImp(TFFTComplexReal) TFFTComplexReal
default
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
const Int_t n
Definition: legend1.C:16
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 ...