Logo ROOT   6.10/09
Reference Guide
TFFTComplex.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 // TFFTComplex
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 // Computes complex input/output discrete Fourier transforms (DFT)
18 // in one or more dimensions. For the detailed information on the computed
19 // transforms please refer to the FFTW manual, chapter "What FFTW really computes".
20 //
21 // How to use it:
22 // 1) Create an instance of TFFTComplex - this will allocate input and output
23 // arrays (unless an in-place transform is specified)
24 // 2) Run the Init() function with the desired flags and settings
25 // 3) Set the data (via SetPoints(), SetPoint() or SetPointComplex() functions)
26 // 4) Run the Transform() function
27 // 5) Get the output (via GetPoints(), GetPoint() or GetPointComplex() functions)
28 // 6) Repeat steps 3)-5) as needed
29 //
30 // For a transform of the same size, but with different flags or sign, rerun the Init()
31 // function and continue with steps 3)-5)
32 // NOTE: 1) running Init() function will overwrite the input array! Don't set any data
33 // before running the Init() function
34 // 2) FFTW computes unnormalized transform, so doing a transform followed by
35 // its inverse will lead to the original array scaled by the transform size
36 //
37 //////////////////////////////////////////////////////////////////////////
38 
39 #include "TFFTComplex.h"
40 #include "fftw3.h"
41 #include "TComplex.h"
42 
43 
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 ///default
48 
50 {
51  fIn = 0;
52  fOut = 0;
53  fPlan = 0;
54  fN = 0;
55  fFlags = 0;
56  fNdim = 0;
57  fTotalSize = 0;
58  fSign = 1;
59 }
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 ///For 1d transforms
63 ///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array
64 
66 {
67  fIn = fftw_malloc(sizeof(fftw_complex) *n);
68  if (!inPlace)
69  fOut = fftw_malloc(sizeof(fftw_complex) * n);
70  else
71  fOut = 0;
72  fN = new Int_t[1];
73  fN[0] = n;
74  fTotalSize = n;
75  fNdim = 1;
76  fSign = 1;
77  fPlan = 0;
78  fFlags = 0;
79 }
80 
81 ////////////////////////////////////////////////////////////////////////////////
82 ///For multidim. transforms
83 ///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array
84 
86 {
87  fNdim = ndim;
88  fTotalSize = 1;
89  fN = new Int_t[fNdim];
90  for (Int_t i=0; i<fNdim; i++){
91  fN[i] = n[i];
92  fTotalSize*=n[i];
93  }
94  fIn = fftw_malloc(sizeof(fftw_complex)*fTotalSize);
95  if (!inPlace)
96  fOut = fftw_malloc(sizeof(fftw_complex) * fTotalSize);
97  else
98  fOut = 0;
99  fSign = 1;
100  fPlan = 0;
101  fFlags = 0;
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 ///Destroys the data arrays and the plan. However, some plan information stays around
106 ///until the root session is over, and is reused if other plans of the same size are
107 ///created
108 
110 {
111  fftw_destroy_plan((fftw_plan)fPlan);
112  fPlan = 0;
113  fftw_free((fftw_complex*)fIn);
114  if (fOut)
115  fftw_free((fftw_complex*)fOut);
116  if (fN)
117  delete [] fN;
118 }
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 ///Creates the fftw-plan
122 ///
123 ///NOTE: input and output arrays are overwritten during initialisation,
124 /// so don't set any points, before running this function!!!!!
125 ///
126 ///2nd parameter: +1
127 ///Argument kind is dummy and doesn't need to be specified
128 ///Possible flag_options:
129 ///"ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
130 /// performance
131 ///"M" (from "measure") - some time spend in finding the optimal way to do the transform
132 ///"P" (from "patient") - more time spend in finding the optimal way to do the transform
133 ///"EX" (from "exhaustive") - the most optimal way is found
134 ///This option should be chosen depending on how many transforms of the same size and
135 ///type are going to be done. Planning is only done once, for the first transform of this
136 ///size and type.
137 
138 void TFFTComplex::Init( Option_t *flags, Int_t sign,const Int_t* /*kind*/)
139 {
140  fSign = sign;
141  fFlags = flags;
142 
143  if (fPlan)
144  fftw_destroy_plan((fftw_plan)fPlan);
145  fPlan = 0;
146 
147  if (fOut)
148  fPlan = (void*)fftw_plan_dft(fNdim, fN, (fftw_complex*)fIn, (fftw_complex*)fOut, sign,MapFlag(flags));
149  else
150  fPlan = (void*)fftw_plan_dft(fNdim, fN, (fftw_complex*)fIn, (fftw_complex*)fIn, sign, MapFlag(flags));
151 }
152 
153 ////////////////////////////////////////////////////////////////////////////////
154 ///Computes the transform, specified in Init() function
155 
157 {
158  if (fPlan)
159  fftw_execute((fftw_plan)fPlan);
160  else {
161  Error("Transform", "transform not initialised");
162  return;
163  }
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 ///Copies the output(or input) into the argument array
168 
170 {
171  if (!fromInput){
172  for (Int_t i=0; i<2*fTotalSize; i+=2){
173  data[i] = ((fftw_complex*)fOut)[i/2][0];
174  data[i+1] = ((fftw_complex*)fOut)[i/2][1];
175  }
176  } else {
177  for (Int_t i=0; i<2*fTotalSize; i+=2){
178  data[i] = ((fftw_complex*)fIn)[i/2][0];
179  data[i+1] = ((fftw_complex*)fIn)[i/2][1];
180  }
181  }
182 }
183 
184 ////////////////////////////////////////////////////////////////////////////////
185 ///returns real and imaginary parts of the point #ipoint
186 
187 void TFFTComplex::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
188 {
189  if (fOut && !fromInput){
190  re = ((fftw_complex*)fOut)[ipoint][0];
191  im = ((fftw_complex*)fOut)[ipoint][1];
192  } else {
193  re = ((fftw_complex*)fIn)[ipoint][0];
194  im = ((fftw_complex*)fIn)[ipoint][1];
195  }
196 }
197 
198 ////////////////////////////////////////////////////////////////////////////////
199 ///For multidimensional transforms. Returns real and imaginary parts of the point #ipoint
200 
201 void TFFTComplex::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
202 {
203  Int_t ireal = ipoint[0];
204  for (Int_t i=0; i<fNdim-1; i++)
205  ireal=fN[i+1]*ireal + ipoint[i+1];
206 
207  if (fOut && !fromInput){
208  re = ((fftw_complex*)fOut)[ireal][0];
209  im = ((fftw_complex*)fOut)[ireal][1];
210  } else {
211  re = ((fftw_complex*)fIn)[ireal][0];
212  im = ((fftw_complex*)fIn)[ireal][1];
213  }
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 ///Copies real and imaginary parts of the output (input) into the argument arrays
218 
219 void TFFTComplex::GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput) const
220 {
221  if (fOut && !fromInput){
222  for (Int_t i=0; i<fTotalSize; i++){
223  re[i] = ((fftw_complex*)fOut)[i][0];
224  im[i] = ((fftw_complex*)fOut)[i][1];
225  }
226  } else {
227  for (Int_t i=0; i<fTotalSize; i++){
228  re[i] = ((fftw_complex*)fIn)[i][0];
229  im[i] = ((fftw_complex*)fIn)[i][1];
230  }
231  }
232 }
233 
234 ////////////////////////////////////////////////////////////////////////////////
235 ///Copies the output(input) into the argument array
236 
238 {
239  if (fOut && !fromInput){
240  for (Int_t i=0; i<fTotalSize; i+=2){
241  data[i] = ((fftw_complex*)fOut)[i/2][0];
242  data[i+1] = ((fftw_complex*)fOut)[i/2][1];
243  }
244  } else {
245  for (Int_t i=0; i<fTotalSize; i+=2){
246  data[i] = ((fftw_complex*)fIn)[i/2][0];
247  data[i+1] = ((fftw_complex*)fIn)[i/2][1];
248  }
249  }
250 }
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 ///sets real and imaginary parts of point # ipoint
254 
256 {
257  ((fftw_complex*)fIn)[ipoint][0]=re;
258  ((fftw_complex*)fIn)[ipoint][1]=im;
259 }
260 
261 ////////////////////////////////////////////////////////////////////////////////
262 ///For multidim. transforms. Sets real and imaginary parts of point # ipoint
263 
264 void TFFTComplex::SetPoint(const Int_t *ipoint, Double_t re, Double_t im)
265 {
266  Int_t ireal = ipoint[0];
267  for (Int_t i=0; i<fNdim-1; i++)
268  ireal=fN[i+1]*ireal + ipoint[i+1];
269 
270  ((fftw_complex*)fIn)[ireal][0]=re;
271  ((fftw_complex*)fIn)[ireal][1]=im;
272 }
273 
274 ////////////////////////////////////////////////////////////////////////////////
275 
277 {
278  ((fftw_complex*)fIn)[ipoint][0] = c.Re();
279  ((fftw_complex*)fIn)[ipoint][1] = c.Im();
280 }
281 
282 ////////////////////////////////////////////////////////////////////////////////
283 ///set all points. the values are copied. points should be ordered as follows:
284 ///[re_0, im_0, re_1, im_1, ..., re_n, im_n)
285 
287 {
288  for (Int_t i=0; i<2*fTotalSize-1; i+=2){
289  ((fftw_complex*)fIn)[i/2][0]=data[i];
290  ((fftw_complex*)fIn)[i/2][1]=data[i+1];
291  }
292 }
293 
294 ////////////////////////////////////////////////////////////////////////////////
295 ///set all points. the values are copied
296 
297 void TFFTComplex::SetPointsComplex(const Double_t *re_data, const Double_t *im_data)
298 {
299  if (!fIn){
300  Error("SetPointsComplex", "Size is not set yet");
301  return;
302  }
303  for (Int_t i=0; i<fTotalSize; i++){
304  ((fftw_complex*)fIn)[i][0]=re_data[i];
305  ((fftw_complex*)fIn)[i][1]=im_data[i];
306  }
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 ///allowed options:
311 ///"ES" - FFTW_ESTIMATE
312 ///"M" - FFTW_MEASURE
313 ///"P" - FFTW_PATIENT
314 ///"EX" - FFTW_EXHAUSTIVE
315 
317 {
318  TString opt = flag;
319  opt.ToUpper();
320  if (opt.Contains("ES"))
321  return FFTW_ESTIMATE;
322  if (opt.Contains("M"))
323  return FFTW_MEASURE;
324  if (opt.Contains("P"))
325  return FFTW_PATIENT;
326  if (opt.Contains("EX"))
327  return FFTW_EXHAUSTIVE;
328  return FFTW_ESTIMATE;
329 }
virtual void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const
Copies real and imaginary parts of the output (input) into the argument arrays.
void * fIn
Definition: TFFTComplex.h:48
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
sets real and imaginary parts of point # ipoint
virtual void Init(Option_t *flags, Int_t sign, const Int_t *)
Creates the fftw-plan.
const char Option_t
Definition: RtypesCore.h:62
virtual void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const
returns real and imaginary parts of the point #ipoint
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1112
Double_t Re() const
Definition: TComplex.h:43
Basic string class.
Definition: TString.h:129
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Int_t fTotalSize
Definition: TFFTComplex.h:52
Option_t * fFlags
Definition: TFFTComplex.h:55
virtual void SetPointComplex(Int_t ipoint, TComplex &c)
Int_t fNdim
Definition: TFFTComplex.h:51
Int_t fSign
Definition: TFFTComplex.h:54
virtual void Transform()
Computes the transform, specified in Init() function.
UInt_t MapFlag(Option_t *flag)
allowed options: "ES" - FFTW_ESTIMATE "M" - FFTW_MEASURE "P" - FFTW_PATIENT "EX" - FFTW_EXHAUSTIVE ...
TFFTComplex()
default
Definition: TFFTComplex.cxx:49
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
void * fOut
Definition: TFFTComplex.h:49
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const
Copies the output(or input) into the argument array.
virtual ~TFFTComplex()
Destroys the data arrays and the plan.
#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
virtual void SetPoints(const Double_t *data)
set all points.
Double_t Im() const
Definition: TComplex.h:44
Int_t * fN
Definition: TFFTComplex.h:53
void * fPlan
Definition: TFFTComplex.h:50
virtual void SetPointsComplex(const Double_t *re, const Double_t *im)
set all points. the values are copied
const Int_t n
Definition: legend1.C:16