Logo ROOT   6.13/01
Reference Guide
TFFTReal.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 // 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 
62 ClassImp(TFFTReal);
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 ///default
66 
68 {
69  fIn = 0;
70  fOut = 0;
71  fPlan = 0;
72  fN = 0;
73  fKind = 0;
74  fNdim = 0;
75  fTotalSize = 0;
76 }
77 
78 ////////////////////////////////////////////////////////////////////////////////
79 ///For 1d transforms
80 ///n here is the physical size of the transform (see FFTW manual for more details)
81 
82 TFFTReal::TFFTReal(Int_t n, Bool_t inPlace)
83 {
84  fIn = fftw_malloc(sizeof(Double_t)*n);
85  if (inPlace) fOut = 0;
86  else fOut = fftw_malloc(sizeof(Double_t)*n);
87 
88  fPlan = 0;
89  fNdim = 1;
90  fN = new Int_t[1];
91  fN[0] = n;
92  fKind = 0;
93  fTotalSize = n;
94 }
95 
96 ////////////////////////////////////////////////////////////////////////////////
97 ///For multidimensional transforms
98 ///1st parameter is the # of dimensions,
99 ///2nd is the sizes (physical) of the transform in each dimension
100 
101 TFFTReal::TFFTReal(Int_t ndim, Int_t *n, Bool_t inPlace)
102 {
103  fTotalSize = 1;
104  fNdim = ndim;
105  fN = new Int_t[ndim];
106  fKind = 0;
107  fPlan = 0;
108  for (Int_t i=0; i<ndim; i++){
109  fTotalSize*=n[i];
110  fN[i] = n[i];
111  }
112  fIn = fftw_malloc(sizeof(Double_t)*fTotalSize);
113  if (!inPlace)
114  fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
115  else
116  fOut = 0;
117 }
118 
119 ////////////////////////////////////////////////////////////////////////////////
120 ///clean-up
121 
123 {
124  fftw_destroy_plan((fftw_plan)fPlan);
125  fPlan = 0;
126  fftw_free(fIn);
127  fIn = 0;
128  if (fOut){
129  fftw_free(fOut);
130  }
131  fOut = 0;
132  if (fN)
133  delete [] fN;
134  fN = 0;
135  if (fKind)
136  fftw_free((fftw_r2r_kind*)fKind);
137  fKind = 0;
138 }
139 
140 ////////////////////////////////////////////////////////////////////////////////
141 ///Creates the fftw-plan
142 ///
143 ///NOTE: input and output arrays are overwritten during initialisation,
144 /// so don't set any points, before running this function!!!!!
145 ///
146 ///1st parameter:
147 /// Possible flag_options:
148 /// "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
149 /// performance
150 /// "M" (from "measure") - some time spend in finding the optimal way to do the transform
151 /// "P" (from "patient") - more time spend in finding the optimal way to do the transform
152 /// "EX" (from "exhaustive") - the most optimal way is found
153 /// This option should be chosen depending on how many transforms of the same size and
154 /// type are going to be done. Planning is only done once, for the first transform of this
155 /// size and type.
156 ///2nd parameter is dummy and doesn't need to be specified
157 ///3rd parameter- transform kind for each dimension
158 /// 4 different kinds of sine and cosine transforms are available
159 /// DCT-I - kind=0
160 /// DCT-II - kind=1
161 /// DCT-III - kind=2
162 /// DCT-IV - kind=3
163 /// DST-I - kind=4
164 /// DST-II - kind=5
165 /// DSTIII - kind=6
166 /// DSTIV - kind=7
167 
168 void TFFTReal::Init( Option_t* flags,Int_t /*sign*/, const Int_t *kind)
169 {
170  if (fPlan)
171  fftw_destroy_plan((fftw_plan)fPlan);
172  fPlan = 0;
173 
174  if (!fKind)
175  fKind = (fftw_r2r_kind*)fftw_malloc(sizeof(fftw_r2r_kind)*fNdim);
176 
177  if (MapOptions(kind)){
178  if (fOut)
179  fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fOut, (fftw_r2r_kind*)fKind, MapFlag(flags));
180  else
181  fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fIn, (fftw_r2r_kind*)fKind, MapFlag(flags));
182  fFlags = flags;
183  }
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 ///Computes the transform, specified in Init() function
188 
190 {
191  if (fPlan)
192  fftw_execute((fftw_plan)fPlan);
193  else {
194  Error("Transform", "transform hasn't been initialised");
195  return;
196  }
197 }
198 
199 ////////////////////////////////////////////////////////////////////////////////
200 ///Returns the type of the transform
201 
202 Option_t *TFFTReal::GetType() const
203 {
204  if (!fKind) {
205  Error("GetType", "Type not defined yet (kind not set)");
206  return "";
207  }
208  if (((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC) return "R2HC";
209  if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R) return "HC2R";
210  if (((fftw_r2r_kind*)fKind)[0]==FFTW_DHT) return "DHT";
211  else return "R2R";
212 }
213 
214 ////////////////////////////////////////////////////////////////////////////////
215 ///Copies the output (or input) points into the provided array, that should
216 ///be big enough
217 
218 void TFFTReal::GetPoints(Double_t *data, Bool_t fromInput) const
219 {
220  const Double_t * array = GetPointsReal(fromInput);
221  if (!array) return;
222  std::copy(array, array+fTotalSize, data);
223 }
224 
225 ////////////////////////////////////////////////////////////////////////////////
226 ///For 1d tranforms. Returns point #ipoint
227 
228 Double_t TFFTReal::GetPointReal(Int_t ipoint, Bool_t fromInput) const
229 {
230  if (ipoint<0 || ipoint>fTotalSize){
231  Error("GetPointReal", "No such point");
232  return 0;
233  }
234  const Double_t * array = GetPointsReal(fromInput);
235  return ( array ) ? array[ipoint] : 0;
236 }
237 
238 ////////////////////////////////////////////////////////////////////////////////
239 ///For multidim.transforms. Returns point #ipoint
240 
241 Double_t TFFTReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
242 {
243  Int_t ireal = ipoint[0];
244  for (Int_t i=0; i<fNdim-1; i++)
245  ireal=fN[i+1]*ireal + ipoint[i+1];
246 
247  const Double_t * array = GetPointsReal(fromInput);
248  return ( array ) ? array[ireal] : 0;
249 }
250 
251 ////////////////////////////////////////////////////////////////////////////////
252 ///Only for input of HC2R and output of R2HC
253 
254 void TFFTReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
255 {
256  const Double_t * array = GetPointsReal(fromInput);
257  if (!array) return;
258  if ( ( ((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC && !fromInput ) ||
259  ( ((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R && fromInput ) )
260  {
261  if (ipoint<fN[0]/2+1){
262  re = array[ipoint];
263  im = array[fN[0]-ipoint];
264  } else {
265  re = array[fN[0]-ipoint];
266  im = -array[ipoint];
267  }
268  if ((fN[0]%2)==0 && ipoint==fN[0]/2) im = 0;
269  }
270 }
271 ////////////////////////////////////////////////////////////////////////////////
272 ///Only for input of HC2R and output of R2HC and for 1d
273 
274 void TFFTReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
275 {
276  GetPointComplex(ipoint[0], re, im, fromInput);
277 }
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 ///Returns the output (or input) array
281 
282 Double_t* TFFTReal::GetPointsReal(Bool_t fromInput) const
283 {
284  // we have 4 different cases
285  // fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
286  // fromInput = false; fOut = NULL (transformed is in place) : return fIn
287  // fromInput = true; fOut = !NULL : return fIn
288  // fromInput = true; fOut = NULL return an error since input array is overwritten
289 
290  if (!fromInput && fOut)
291  return (Double_t*)fOut;
292  else if (fromInput && !fOut) {
293  Error("GetPointsReal","Input array was destroyed");
294  return 0;
295  }
296  //R__ASSERT(fIn);
297  return (Double_t*)fIn;
298 }
299 
300 ////////////////////////////////////////////////////////////////////////////////
301 
302 void TFFTReal::SetPoint(Int_t ipoint, Double_t re, Double_t im)
303 {
304  if (ipoint<0 || ipoint>fTotalSize){
305  Error("SetPoint", "illegal point index");
306  return;
307  }
308  if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R){
309  if ((fN[0]%2)==0 && ipoint==fN[0]/2)
310  ((Double_t*)fIn)[ipoint] = re;
311  else {
312  ((Double_t*)fIn)[ipoint] = re;
313  ((Double_t*)fIn)[fN[0]-ipoint]=im;
314  }
315  }
316  else
317  ((Double_t*)fIn)[ipoint]=re;
318 }
319 
320 ////////////////////////////////////////////////////////////////////////////////
321 ///Since multidimensional R2HC and HC2R transforms are not supported,
322 ///third parameter is dummy
323 
324 void TFFTReal::SetPoint(const Int_t *ipoint, Double_t re, Double_t /*im*/)
325 {
326  Int_t ireal = ipoint[0];
327  for (Int_t i=0; i<fNdim-1; i++)
328  ireal=fN[i+1]*ireal + ipoint[i+1];
329  if (ireal < 0 || ireal >fTotalSize){
330  Error("SetPoint", "illegal point index");
331  return;
332  }
333  ((Double_t*)fIn)[ireal]=re;
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 ///Sets all points
338 
339 void TFFTReal::SetPoints(const Double_t *data)
340 {
341  for (Int_t i=0; i<fTotalSize; i++)
342  ((Double_t*)fIn)[i] = data[i];
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 ///transfers the r2r_kind parameters to fftw type
347 
348 Int_t TFFTReal::MapOptions(const Int_t *kind)
349 {
350  if (kind[0] == 10){
351  if (fNdim>1){
352  Error("Init", "Multidimensional R2HC transforms are not supported, use R2C interface instead");
353  return 0;
354  }
355  ((fftw_r2r_kind*)fKind)[0] = FFTW_R2HC;
356  }
357  else if (kind[0] == 11) {
358  if (fNdim>1){
359  Error("Init", "Multidimensional HC2R transforms are not supported, use C2R interface instead");
360  return 0;
361  }
362  ((fftw_r2r_kind*)fKind)[0] = FFTW_HC2R;
363  }
364  else if (kind[0] == 12) {
365  for (Int_t i=0; i<fNdim; i++)
366  ((fftw_r2r_kind*)fKind)[i] = FFTW_DHT;
367  }
368  else {
369  for (Int_t i=0; i<fNdim; i++){
370  switch (kind[i]) {
371  case 0: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT00; break;
372  case 1: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT01; break;
373  case 2: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT10; break;
374  case 3: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT11; break;
375  case 4: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT00; break;
376  case 5: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT01; break;
377  case 6: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT10; break;
378  case 7: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT11; break;
379  default:
380  ((fftw_r2r_kind*)fKind)[i] = FFTW_R2HC; break;
381  }
382  }
383  }
384  return 1;
385 }
386 
387 ////////////////////////////////////////////////////////////////////////////////
388 ///allowed options:
389 ///"ES" - FFTW_ESTIMATE
390 ///"M" - FFTW_MEASURE
391 ///"P" - FFTW_PATIENT
392 ///"EX" - FFTW_EXHAUSTIVE
393 
394 UInt_t TFFTReal::MapFlag(Option_t *flag)
395 {
396  TString opt = flag;
397  opt.ToUpper();
398  if (opt.Contains("ES"))
399  return FFTW_ESTIMATE;
400  if (opt.Contains("M"))
401  return FFTW_MEASURE;
402  if (opt.Contains("P"))
403  return FFTW_PATIENT;
404  if (opt.Contains("EX"))
405  return FFTW_EXHAUSTIVE;
406  return FFTW_ESTIMATE;
407 }
virtual void SetPoints(const Double_t *data)
Sets all points.
Definition: TFFTReal.cxx:339
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const
For 1d tranforms. Returns point #ipoint.
Definition: TFFTReal.cxx:228
Int_t * fN
Definition: TFFTReal.h:74
void * fOut
Definition: TFFTReal.h:70
virtual Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const
Returns the output (or input) array.
Definition: TFFTReal.cxx:282
virtual Option_t * GetType() const
Returns the type of the transform.
Definition: TFFTReal.cxx:202
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:274
void * fIn
Definition: TFFTReal.h:69
Int_t MapOptions(const Int_t *kind)
transfers the r2r_kind parameters to fftw type
Definition: TFFTReal.cxx:348
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE ...
Definition: TFFTReal.cxx:394
virtual void Init(Option_t *flags, Int_t sign, const Int_t *kind)
Creates the fftw-plan.
Definition: TFFTReal.cxx:168
void * fKind
Definition: TFFTReal.h:75
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
Definition: TFFTReal.cxx:302
TFFTReal()
default
Definition: TFFTReal.cxx:67
Int_t fNdim
Definition: TFFTReal.h:72
virtual ~TFFTReal()
clean-up
Definition: TFFTReal.cxx:122
Int_t fTotalSize
Definition: TFFTReal.h:73
TString fFlags
Definition: TFFTReal.h:76
virtual void Transform()
Computes the transform, specified in Init() function.
Definition: TFFTReal.cxx:189
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:218
void * fPlan
Definition: TFFTReal.h:71