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