ROOT   6.10/09 Reference Guide
TFFTReal.cxx
Go to the documentation of this file.
1 // @(#)root/fft:$Id$
2 // Author: Anna Kreshuk 07/4/2006
3
4 /*************************************************************************
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11
12 //////////////////////////////////////////////////////////////////////////
13 //
14 // TFFTReal
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 transforms called r2r in FFTW manual:
19 // - transforms of real input and output in "halfcomplex" format i.e.
20 // real and imaginary parts for a transform of size n stored as
21 // (r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1)
22 // - discrete Hartley transform
23 // - sine and cosine transforms (DCT-I,II,III,IV and DST-I,II,III,IV)
24 // For the detailed information on the computed
25 // transforms please refer to the FFTW manual, chapter "What FFTW really computes".
26 //
27 // How to use it:
28 // 1) Create an instance of TFFTReal - this will allocate input and output
29 // arrays (unless an in-place transform is specified)
30 // 2) Run the Init() function with the desired flags and settings (see function
31 // comments for possible kind parameters)
32 // 3) Set the data (via SetPoints()or SetPoint() functions)
33 // 4) Run the Transform() function
34 // 5) Get the output (via GetPoints() or GetPoint() functions)
35 // 6) Repeat steps 3)-5) as needed
36 // For a transform of the same size, but of different kind (or with different flags),
37 // rerun the Init() function and continue with steps 3)-5)
38 //
39 // NOTE: 1) running Init() function will overwrite the input array! Don't set any data
40 // before running the Init() function!
41 // 2) FFTW computes unnormalized transform, so doing a transform followed by
42 // its inverse will lead to the original array scaled BY:
43 // - transform size (N) for R2HC, HC2R, DHT transforms
44 // - 2*(N-1) for DCT-I (REDFT00)
45 // - 2*(N+1) for DST-I (RODFT00)
46 // - 2*N for the remaining transforms
47 // Transform inverses:
48 // R2HC<-->HC2R
49 // DHT<-->DHT
50 // DCT-I<-->DCT-I
51 // DCT-II<-->DCT-III
52 // DCT-IV<-->DCT-IV
53 // DST-I<-->DST-I
54 // DST-II<-->DST-III
55 // DST-IV<-->DST-IV
56 //
57 //////////////////////////////////////////////////////////////////////////
58
59 #include "TFFTReal.h"
60 #include "fftw3.h"
61
63
64 ////////////////////////////////////////////////////////////////////////////////
65 ///default
66
68 {
69  fIn = 0;
70  fOut = 0;
71  fPlan = 0;
72  fN = 0;
73  fKind = 0;
74  fFlags = 0;
75  fNdim = 0;
76  fTotalSize = 0;
77 }
78
79 ////////////////////////////////////////////////////////////////////////////////
80 ///For 1d transforms
81 ///n here is the physical size of the transform (see FFTW manual for more details)
82
84 {
85  fIn = fftw_malloc(sizeof(Double_t)*n);
86  if (inPlace) fOut = 0;
87  else fOut = fftw_malloc(sizeof(Double_t)*n);
88
89  fPlan = 0;
90  fNdim = 1;
91  fN = new Int_t[1];
92  fN[0] = n;
93  fKind = 0;
94  fTotalSize = n;
95  fFlags = 0;
96 }
97
98 ////////////////////////////////////////////////////////////////////////////////
99 ///For multidimensional transforms
100 ///1st parameter is the # of dimensions,
101 ///2nd is the sizes (physical) of the transform in each dimension
102
104 {
105  fTotalSize = 1;
106  fNdim = ndim;
107  fN = new Int_t[ndim];
108  fKind = 0;
109  fPlan = 0;
110  fFlags = 0;
111  for (Int_t i=0; i<ndim; i++){
112  fTotalSize*=n[i];
113  fN[i] = n[i];
114  }
115  fIn = fftw_malloc(sizeof(Double_t)*fTotalSize);
116  if (!inPlace)
117  fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
118  else
119  fOut = 0;
120 }
121
122 ////////////////////////////////////////////////////////////////////////////////
123 ///clean-up
124
126 {
127  fftw_destroy_plan((fftw_plan)fPlan);
128  fPlan = 0;
129  fftw_free(fIn);
130  fIn = 0;
131  if (fOut){
132  fftw_free(fOut);
133  }
134  fOut = 0;
135  if (fN)
136  delete [] fN;
137  fN = 0;
138  if (fKind)
139  fftw_free((fftw_r2r_kind*)fKind);
140  fKind = 0;
141 }
142
143 ////////////////////////////////////////////////////////////////////////////////
144 ///Creates the fftw-plan
145 ///
146 ///NOTE: input and output arrays are overwritten during initialisation,
147 /// so don't set any points, before running this function!!!!!
148 ///
149 ///1st parameter:
150 /// Possible flag_options:
151 /// "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
152 /// performance
153 /// "M" (from "measure") - some time spend in finding the optimal way to do the transform
154 /// "P" (from "patient") - more time spend in finding the optimal way to do the transform
155 /// "EX" (from "exhaustive") - the most optimal way is found
156 /// This option should be chosen depending on how many transforms of the same size and
157 /// type are going to be done. Planning is only done once, for the first transform of this
158 /// size and type.
159 ///2nd parameter is dummy and doesn't need to be specified
160 ///3rd parameter- transform kind for each dimension
161 /// 4 different kinds of sine and cosine transforms are available
162 /// DCT-I - kind=0
163 /// DCT-II - kind=1
164 /// DCT-III - kind=2
165 /// DCT-IV - kind=3
166 /// DST-I - kind=4
167 /// DST-II - kind=5
168 /// DSTIII - kind=6
169 /// DSTIV - kind=7
170
171 void TFFTReal::Init( Option_t* flags,Int_t /*sign*/, const Int_t *kind)
172 {
173  if (fPlan)
174  fftw_destroy_plan((fftw_plan)fPlan);
175  fPlan = 0;
176
177  if (!fKind)
178  fKind = (fftw_r2r_kind*)fftw_malloc(sizeof(fftw_r2r_kind)*fNdim);
179
180  if (MapOptions(kind)){
181  if (fOut)
182  fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fOut, (fftw_r2r_kind*)fKind, MapFlag(flags));
183  else
184  fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fIn, (fftw_r2r_kind*)fKind, MapFlag(flags));
185  fFlags = flags;
186  }
187 }
188
189 ////////////////////////////////////////////////////////////////////////////////
190 ///Computes the transform, specified in Init() function
191
193 {
194  if (fPlan)
195  fftw_execute((fftw_plan)fPlan);
196  else {
197  Error("Transform", "transform hasn't been initialised");
198  return;
199  }
200 }
201
202 ////////////////////////////////////////////////////////////////////////////////
203 ///Returns the type of the transform
204
206 {
207  if (!fKind) {
208  Error("GetType", "Type not defined yet (kind not set)");
209  return "";
210  }
211  if (((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC) return "R2HC";
212  if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R) return "HC2R";
213  if (((fftw_r2r_kind*)fKind)[0]==FFTW_DHT) return "DHT";
214  else return "R2R";
215 }
216
217 ////////////////////////////////////////////////////////////////////////////////
218 ///Copies the output (or input) points into the provided array, that should
219 ///be big enough
220
221 void TFFTReal::GetPoints(Double_t *data, Bool_t fromInput) const
222 {
223  const Double_t * array = GetPointsReal(fromInput);
224  if (!array) return;
225  std::copy(array, array+fTotalSize, data);
226 }
227
228 ////////////////////////////////////////////////////////////////////////////////
229 ///For 1d tranforms. Returns point #ipoint
230
231 Double_t TFFTReal::GetPointReal(Int_t ipoint, Bool_t fromInput) const
232 {
233  if (ipoint<0 || ipoint>fTotalSize){
234  Error("GetPointReal", "No such point");
235  return 0;
236  }
237  const Double_t * array = GetPointsReal(fromInput);
238  return ( array ) ? array[ipoint] : 0;
239 }
240
241 ////////////////////////////////////////////////////////////////////////////////
242 ///For multidim.transforms. Returns point #ipoint
243
244 Double_t TFFTReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
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  const Double_t * array = GetPointsReal(fromInput);
251  return ( array ) ? array[ireal] : 0;
252 }
253
254 ////////////////////////////////////////////////////////////////////////////////
255 ///Only for input of HC2R and output of R2HC
256
257 void TFFTReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
258 {
259  const Double_t * array = GetPointsReal(fromInput);
260  if (!array) return;
261  if ( ( ((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC && !fromInput ) ||
262  ( ((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R && fromInput ) )
263  {
264  if (ipoint<fN[0]/2+1){
265  re = array[ipoint];
266  im = array[fN[0]-ipoint];
267  } else {
268  re = array[fN[0]-ipoint];
269  im = -array[ipoint];
270  }
271  if ((fN[0]%2)==0 && ipoint==fN[0]/2) im = 0;
272  }
273 }
274 ////////////////////////////////////////////////////////////////////////////////
275 ///Only for input of HC2R and output of R2HC and for 1d
276
277 void TFFTReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
278 {
279  GetPointComplex(ipoint[0], re, im, fromInput);
280 }
281
282 ////////////////////////////////////////////////////////////////////////////////
283 ///Returns the output (or input) array
284
286 {
287  // we have 4 different cases
288  // fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
289  // fromInput = false; fOut = NULL (transformed is in place) : return fIn
290  // fromInput = true; fOut = !NULL : return fIn
291  // fromInput = true; fOut = NULL return an error since input array is overwritten
292
293  if (!fromInput && fOut)
294  return (Double_t*)fOut;
295  else if (fromInput && !fOut) {
296  Error("GetPointsReal","Input array was destroyed");
297  return 0;
298  }
299  //R__ASSERT(fIn);
300  return (Double_t*)fIn;
301 }
302
303 ////////////////////////////////////////////////////////////////////////////////
304
306 {
307  if (ipoint<0 || ipoint>fTotalSize){
308  Error("SetPoint", "illegal point index");
309  return;
310  }
311  if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R){
312  if ((fN[0]%2)==0 && ipoint==fN[0]/2)
313  ((Double_t*)fIn)[ipoint] = re;
314  else {
315  ((Double_t*)fIn)[ipoint] = re;
316  ((Double_t*)fIn)[fN[0]-ipoint]=im;
317  }
318  }
319  else
320  ((Double_t*)fIn)[ipoint]=re;
321 }
322
323 ////////////////////////////////////////////////////////////////////////////////
324 ///Since multidimensional R2HC and HC2R transforms are not supported,
325 ///third parameter is dummy
326
327 void TFFTReal::SetPoint(const Int_t *ipoint, Double_t re, Double_t /*im*/)
328 {
329  Int_t ireal = ipoint[0];
330  for (Int_t i=0; i<fNdim-1; i++)
331  ireal=fN[i+1]*ireal + ipoint[i+1];
332  if (ireal < 0 || ireal >fTotalSize){
333  Error("SetPoint", "illegal point index");
334  return;
335  }
336  ((Double_t*)fIn)[ireal]=re;
337 }
338
339 ////////////////////////////////////////////////////////////////////////////////
340 ///Sets all points
341
343 {
344  for (Int_t i=0; i<fTotalSize; i++)
345  ((Double_t*)fIn)[i] = data[i];
346 }
347
348 ////////////////////////////////////////////////////////////////////////////////
349 ///transfers the r2r_kind parameters to fftw type
350
352 {
353  if (kind[0] == 10){
354  if (fNdim>1){
355  Error("Init", "Multidimensional R2HC transforms are not supported, use R2C interface instead");
356  return 0;
357  }
358  ((fftw_r2r_kind*)fKind)[0] = FFTW_R2HC;
359  }
360  else if (kind[0] == 11) {
361  if (fNdim>1){
362  Error("Init", "Multidimensional HC2R transforms are not supported, use C2R interface instead");
363  return 0;
364  }
365  ((fftw_r2r_kind*)fKind)[0] = FFTW_HC2R;
366  }
367  else if (kind[0] == 12) {
368  for (Int_t i=0; i<fNdim; i++)
369  ((fftw_r2r_kind*)fKind)[i] = FFTW_DHT;
370  }
371  else {
372  for (Int_t i=0; i<fNdim; i++){
373  switch (kind[i]) {
374  case 0: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT00; break;
375  case 1: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT01; break;
376  case 2: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT10; break;
377  case 3: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT11; break;
378  case 4: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT00; break;
379  case 5: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT01; break;
380  case 6: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT10; break;
381  case 7: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT11; break;
382  default:
383  ((fftw_r2r_kind*)fKind)[i] = FFTW_R2HC; break;
384  }
385  }
386  }
387  return 1;
388 }
389
390 ////////////////////////////////////////////////////////////////////////////////
391 ///allowed options:
392 ///"ES" - FFTW_ESTIMATE
393 ///"M" - FFTW_MEASURE
394 ///"P" - FFTW_PATIENT
395 ///"EX" - FFTW_EXHAUSTIVE
396
398 {
399  TString opt = flag;
400  opt.ToUpper();
401  if (opt.Contains("ES"))
402  return FFTW_ESTIMATE;
403  if (opt.Contains("M"))
404  return FFTW_MEASURE;
405  if (opt.Contains("P"))
406  return FFTW_PATIENT;
407  if (opt.Contains("EX"))
408  return FFTW_EXHAUSTIVE;
409  return FFTW_ESTIMATE;
410 }
virtual void SetPoints(const Double_t *data)
Sets all points.
Definition: TFFTReal.cxx:342
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const
For 1d tranforms. Returns point #ipoint.
Definition: TFFTReal.cxx:231
const char Option_t
Definition: RtypesCore.h:62
Int_t * fN
Definition: TFFTReal.h:73
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1112
void * fOut
Definition: TFFTReal.h:69
virtual Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const
Returns the output (or input) array.
Definition: TFFTReal.cxx:285
Basic string class.
Definition: TString.h:129
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual Option_t * GetType() const
Returns the type of the transform.
Definition: TFFTReal.cxx:205
virtual void GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const
Only for input of HC2R and output of R2HC and for 1d.
Definition: TFFTReal.cxx:277
void * fIn
Definition: TFFTReal.h:68
Int_t MapOptions(const Int_t *kind)
transfers the r2r_kind parameters to fftw type
Definition: TFFTReal.cxx:351
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE ...
Definition: TFFTReal.cxx:397
virtual void Init(Option_t *flags, Int_t sign, const Int_t *kind)
Creates the fftw-plan.
Definition: TFFTReal.cxx:171
void * fKind
Definition: TFFTReal.h:74
unsigned int UInt_t
Definition: RtypesCore.h:42
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
Definition: TFFTReal.cxx:305
TFFTReal()
default
Definition: TFFTReal.cxx:67
Int_t fNdim
Definition: TFFTReal.h:71
virtual ~TFFTReal()
clean-up
Definition: TFFTReal.cxx:125
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:572
Int_t fTotalSize
Definition: TFFTReal.h:72
Option_t * fFlags
Definition: TFFTReal.h:75
virtual void Transform()
Computes the transform, specified in Init() function.
Definition: TFFTReal.cxx:192
const Int_t n
Definition: legend1.C:16
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const
Copies the output (or input) points into the provided array, that should be big enough.
Definition: TFFTReal.cxx:221
void * fPlan
Definition: TFFTReal.h:70