Logo ROOT   6.12/07
Reference Guide
TVirtualFFT.cxx
Go to the documentation of this file.
1 // @(#)root/base:$Id$
2 // Author: Anna Kreshuk 10/04/2006
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, 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 /** \class TVirtualFFT
13 \ingroup Base
14 
15 TVirtualFFT is an interface class for Fast Fourier Transforms.
16 
17 The default FFT library is FFTW. To use it, FFTW3 library should already
18 be installed, and ROOT should be have fftw3 module enabled, with the directories
19 of fftw3 include file and library specified (see installation instructions).
20 Function SetDefaultFFT() allows to change the default library.
21 
22 ## Available transform types:
23 FFT:
24  - "C2CFORWARD" - a complex input/output discrete Fourier transform (DFT)
25  in one or more dimensions, -1 in the exponent
26  - "C2CBACKWARD"- a complex input/output discrete Fourier transform (DFT)
27  in one or more dimensions, +1 in the exponent
28  - "R2C" - a real-input/complex-output discrete Fourier transform (DFT)
29  in one or more dimensions,
30  - "C2R" - inverse transforms to "R2C", taking complex input
31  (storing the non-redundant half of a logically Hermitian array)
32  to real output
33  - "R2HC" - a real-input DFT with output in ¡Èhalfcomplex¡É format,
34  i.e. real and imaginary parts for a transform of size n stored as
35  r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1
36  - "HC2R" - computes the reverse of FFTW_R2HC, above
37  - "DHT" - computes a discrete Hartley transform
38 
39 ## Sine/cosine transforms:
40 Different types of transforms are specified by parameter kind of the SineCosine() static
41 function. 4 different kinds of sine and cosine transforms are available
42 
43  - DCT-I (REDFT00 in FFTW3 notation)- kind=0
44  - DCT-II (REDFT01 in FFTW3 notation)- kind=1
45  - DCT-III(REDFT10 in FFTW3 notation)- kind=2
46  - DCT-IV (REDFT11 in FFTW3 notation)- kind=3
47  - DST-I (RODFT00 in FFTW3 notation)- kind=4
48  - DST-II (RODFT01 in FFTW3 notation)- kind=5
49  - DST-III(RODFT10 in FFTW3 notation)- kind=6
50  - DST-IV (RODFT11 in FFTW3 notation)- kind=7
51 
52 Formulas and detailed descriptions can be found in the chapter
53 "What FFTW really computes" of the FFTW manual
54 
55 NOTE: FFTW computes unnormalized transforms, so doing a transform, followed by its
56  inverse will give the original array, multiplied by normalization constant
57  (transform size(N) for FFT, 2*(N-1) for DCT-I, 2*(N+1) for DST-I, 2*N for
58  other sine/cosine transforms)
59 
60 ## How to use it:
61 Call to the static function FFT returns a pointer to a fast Fourier transform
62 with requested parameters. Call to the static function SineCosine returns a
63 pointer to a sine or cosine transform with requested parameters. Example:
64 ~~~ {.cpp}
65 {
66  Int_t N = 10; Double_t *in = new Double_t[N];
67  TVirtualFFT *fftr2c = TVirtualFFT::FFT(1, &N, "R2C");
68  fftr2c->SetPoints(in);
69  fftr2c->Transform();
70  Double_t re, im;
71  for (Int_t i=0; i<N; i++)
72  fftr2c->GetPointComplex(i, re, im);
73  ...
74  fftr2c->SetPoints(in2);
75  ...
76  fftr2c->SetPoints(in3);
77  ...
78 }
79 ~~~
80 Different options are explained in the function comments
81 */
82 
83 #include "TROOT.h"
84 #include "TVirtualFFT.h"
85 #include "TPluginManager.h"
86 #include "TEnv.h"
87 #include "TError.h"
88 
91 
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 ///destructor
96 
98 {
99  if (this==fgFFT)
100  fgFFT = 0;
101 }
102 
103 ////////////////////////////////////////////////////////////////////////////////
104 ///Returns a pointer to the FFT of requested size and type.
105 ///
106 /// \param[in] ndim number of transform dimensions
107 /// \param[in] n sizes of each dimension (an array at least ndim long)
108 /// \param [in] option consists of 3 parts - flag option and an option to create a new TVirtualFFT
109 /// 1. transform type option:
110 /// Available transform types are:
111 /// C2CForward, C2CBackward, C2R, R2C, R2HC, HC2R, DHT
112 /// see class description for details
113 /// 2. flag option: choosing how much time should be spent in planning the transform:
114 /// Possible options:
115 /// - "ES" (from "estimate") - no time in preparing the transform,
116 /// but probably sub-optimal performance
117 /// - "M" (from "measure") - some time spend in finding the optimal way
118 /// to do the transform
119 /// - "P" (from "patient") - more time spend in finding the optimal way
120 /// to do the transform
121 /// - "EX" (from "exhaustive") - the most optimal way is found
122 /// This option should be chosen depending on how many transforms of the
123 /// same size and type are going to be done.
124 /// Planning is only done once, for the first transform of this size and type.
125 /// 3. option allowing to choose between the global fgFFT and a new TVirtualFFT object
126 /// "" - default, changes and returns the global fgFFT variable
127 /// "K" (from "keep")- without touching the global fgFFT,
128 /// creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
129 ///
130 /// Examples of valid options: "R2C ES K", "C2CF M", "DHT P K", etc.
131 
133 {
134 
135  Int_t inputtype=0, currenttype=0;
136  TString opt = option;
137  opt.ToUpper();
138  //find the tranform flag
139  Option_t *flag;
140  flag = "ES";
141  if (opt.Contains("ES")) flag = "ES";
142  if (opt.Contains("M")) flag = "M";
143  if (opt.Contains("P")) flag = "P";
144  if (opt.Contains("EX")) flag = "EX";
145 
146  Int_t ndiff = 0;
147 
148  if (!opt.Contains("K")) {
149  if (fgFFT){
150  //if the global transform exists, check if it should be changed
151  if (fgFFT->GetNdim()!=ndim)
152  ndiff++;
153  else {
154  Int_t *ncurrent = fgFFT->GetN();
155  for (Int_t i=0; i<ndim; i++){
156  if (n[i]!=ncurrent[i])
157  ndiff++;
158  }
159  }
160  Option_t *t = fgFFT->GetType();
161  if (!opt.Contains(t)) {
162  if (opt.Contains("HC") || opt.Contains("DHT"))
163  inputtype = 1;
164  if (strcmp(t,"R2HC")==0 || strcmp(t,"HC2R")==0 || strcmp(t,"DHT")==0)
165  currenttype=1;
166 
167  if (!(inputtype==1 && currenttype==1))
168  ndiff++;
169  }
170  if (ndiff>0){
171  delete fgFFT;
172  fgFFT = 0;
173  }
174  }
175  }
176 
177  Int_t sign = 0;
178  if (opt.Contains("C2CB") || opt.Contains("C2R"))
179  sign = 1;
180  if (opt.Contains("C2CF") || opt.Contains("R2C"))
181  sign = -1;
182 
183  TVirtualFFT *fft = 0;
184  if (opt.Contains("K") || !fgFFT) {
185 
187 
188  TPluginHandler *h;
189  TString pluginname;
190  if (fgDefault.Length()==0) fgDefault="fftw";
191  if (strcmp(fgDefault.Data(),"fftw")==0) {
192  if (opt.Contains("C2C")) pluginname = "fftwc2c";
193  if (opt.Contains("C2R")) pluginname = "fftwc2r";
194  if (opt.Contains("R2C")) pluginname = "fftwr2c";
195  if (opt.Contains("HC") || opt.Contains("DHT")) pluginname = "fftwr2r";
196  if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
197  if (h->LoadPlugin()==-1) {
198  ::Error("TVirtualFFT::FFT", "handler not found");
199  return 0;
200  }
201  fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
202  if (!fft) {
203  ::Error("TVirtualFFT::FFT", "plugin failed to create TVirtualFFT object");
204  return 0;
205  }
206  Int_t *kind = new Int_t[1];
207  if (pluginname=="fftwr2r") {
208  if (opt.Contains("R2HC")) kind[0] = 10;
209  if (opt.Contains("HC2R")) kind[0] = 11;
210  if (opt.Contains("DHT")) kind[0] = 12;
211  }
212  fft->Init(flag, sign, kind);
213  if (!opt.Contains("K")) {
214  fgFFT = fft;
215  }
216  delete [] kind;
217  return fft;
218  }
219  else {
220  ::Error("TVirtualFFT::FFT", "plugin not found");
221  return 0;
222  }
223  }
224  } else {
225 
227 
228  //if the global transform already exists and just needs to be reinitialised
229  //with different parameters
230  if (fgFFT->GetSign()!=sign || !opt.Contains(fgFFT->GetTransformFlag()) || !opt.Contains(fgFFT->GetType())) {
231  Int_t *kind = new Int_t[1];
232  if (inputtype==1) {
233  if (opt.Contains("R2HC")) kind[0] = 10;
234  if (opt.Contains("HC2R")) kind[0] = 11;
235  if (opt.Contains("DHT")) kind[0] = 12;
236  }
237  fgFFT->Init(flag, sign, kind);
238  delete [] kind;
239  }
240  }
241  return fgFFT;
242 }
243 
244 ////////////////////////////////////////////////////////////////////////////////
245 ///Returns a pointer to a sine or cosine transform of requested size and kind
246 ///
247 ///Parameters:
248 /// \param [in] ndim number of transform dimensions
249 /// \param [in] n sizes of each dimension (an array at least ndim long)
250 /// \param [in] r2rkind transform kind for each dimension
251 /// 4 different kinds of sine and cosine transforms are available
252 /// - DCT-I - kind=0
253 /// - DCT-II - kind=1
254 /// - DCT-III - kind=2
255 /// - DCT-IV - kind=3
256 /// - DST-I - kind=4
257 /// - DST-II - kind=5
258 /// - DST-III - kind=6
259 /// - DST-IV - kind=7
260 /// \param [in] option : consists of 2 parts
261 /// - flag option and an option to create a new TVirtualFFT
262 /// - flag option: choosing how much time should be spent in planning the transform:
263 /// Possible options:
264 /// - "ES" (from "estimate") - no time in preparing the transform,
265 /// but probably sub-optimal performance
266 /// - "M" (from "measure") - some time spend in finding the optimal way
267 /// to do the transform
268 /// - "P" (from "patient") - more time spend in finding the optimal way
269 /// to do the transform
270 /// - "EX" (from "exhaustive") - the most optimal way is found
271 /// This option should be chosen depending on how many transforms of the
272 /// same size and type are going to be done.
273 /// Planning is only done once, for the first transform of this size and type.
274 /// - option allowing to choose between the global fgFFT and a new TVirtualFFT object
275 /// - "" - default, changes and returns the global fgFFT variable
276 /// - "K" (from "keep")- without touching the global fgFFT,
277 /// creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
278 /// Examples of valid options: "ES K", "EX", etc
279 
281 {
282  TString opt = option;
283  //find the tranform flag
284  Option_t *flag;
285  flag = "ES";
286  if (opt.Contains("ES")) flag = "ES";
287  if (opt.Contains("M")) flag = "M";
288  if (opt.Contains("P")) flag = "P";
289  if (opt.Contains("EX")) flag = "EX";
290 
291  if (!opt.Contains("K")) {
292  if (fgFFT){
293  Int_t ndiff = 0;
294  if (fgFFT->GetNdim()!=ndim || strcmp(fgFFT->GetType(),"R2R")!=0)
295  ndiff++;
296  else {
297  Int_t *ncurrent = fgFFT->GetN();
298  for (Int_t i=0; i<ndim; i++) {
299  if (n[i] != ncurrent[i])
300  ndiff++;
301  }
302 
303  }
304  if (ndiff>0) {
305  delete fgFFT;
306  fgFFT = 0;
307  }
308  }
309  }
310  TVirtualFFT *fft = 0;
311 
313 
314  if (!fgFFT || opt.Contains("K")) {
315  TPluginHandler *h;
316  TString pluginname;
317  if (fgDefault.Length()==0) fgDefault="fftw";
318  if (strcmp(fgDefault.Data(),"fftw")==0) {
319  pluginname = "fftwr2r";
320  if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
321  if (h->LoadPlugin()==-1){
322  ::Error("TVirtualFFT::SineCosine", "handler not found");
323  return 0;
324  }
325  fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
326  if (!fft) {
327  ::Error("TVirtualFFT::SineCosine", "plugin failed to create TVirtualFFT object");
328  return 0;
329  }
330  fft->Init(flag, 0, r2rkind);
331  if (!opt.Contains("K"))
332  fgFFT = fft;
333  return fft;
334  } else {
335  ::Error("TVirtualFFT::SineCosine", "handler not found");
336  return 0;
337  }
338  }
339  }
340 
341  //if (fgFFT->GetTransformFlag()!=flag)
342  fgFFT->Init(flag,0, r2rkind);
343  return fgFFT;
344 }
345 
346 ////////////////////////////////////////////////////////////////////////////////
347 /// static: return current fgFFT
348 
350 {
351  if (fgFFT)
352  return fgFFT;
353  else{
354  ::Warning("TVirtualFFT::GetCurrentTransform", "fgFFT is not defined yet");
355  return 0;
356  }
357 }
358 
359 ////////////////////////////////////////////////////////////////////////////////
360 /// static: set the current transfrom to parameter
361 
363 {
364  fgFFT = fft;
365 }
366 
367 ////////////////////////////////////////////////////////////////////////////////
368 /// static: return the name of the default fft
369 
371 {
372  return fgDefault.Data();
373 }
374 
375 ////////////////////////////////////////////////////////////////////////////////
376 /// static: set name of default fft
377 
379 {
380  if (fgDefault == name) return;
381  delete fgFFT;
382  fgFFT = 0;
383  fgDefault = name;
384 }
static void SetTransform(TVirtualFFT *fft)
static: set the current transfrom to parameter
static const char * GetDefaultFFT()
static: return the name of the default fft
virtual ~TVirtualFFT()
destructor
Definition: TVirtualFFT.cxx:97
const char Option_t
Definition: RtypesCore.h:62
TH1 * h
Definition: legend2.C:5
static TVirtualFFT * GetCurrentTransform()
static: return current fgFFT
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1112
#define gROOT
Definition: TROOT.h:402
Int_t LoadPlugin()
Load the plugin library for this handler.
Basic string class.
Definition: TString.h:125
int Int_t
Definition: RtypesCore.h:41
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:57
static TString fgDefault
Definition: TVirtualFFT.h:92
static void SetDefaultFFT(const char *name="")
static: set name of default fft
static TVirtualFFT * fgFFT
Definition: TVirtualFFT.h:91
virtual Int_t * GetN() const =0
virtual void Init(Option_t *flag, Int_t sign, const Int_t *kind)=0
Long_t ExecPlugin(int nargs, const T &... params)
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
Ssiz_t Length() const
Definition: TString.h:386
virtual Option_t * GetTransformFlag() const =0
const Bool_t kFALSE
Definition: RtypesCore.h:88
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
static TVirtualFFT * SineCosine(Int_t ndim, Int_t *n, Int_t *r2rkind, Option_t *option)
Returns a pointer to a sine or cosine transform of requested size and kind.
#define ClassImp(name)
Definition: Rtypes.h:359
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:570
virtual Int_t GetSign() const =0
#define R__LOCKGUARD(mutex)
virtual Int_t GetNdim() const =0
const Int_t n
Definition: legend1.C:16
char name[80]
Definition: TGX11.cxx:109
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual Option_t * GetType() const =0
const char * Data() const
Definition: TString.h:345