Logo ROOT  
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 "TError.h"
87 
90 
92 
93 ////////////////////////////////////////////////////////////////////////////////
94 ///destructor
95 
97 {
98  if (this==fgFFT)
99  fgFFT = 0;
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 ///Returns a pointer to the FFT of requested size and type.
104 ///
105 /// \param[in] ndim number of transform dimensions
106 /// \param[in] n sizes of each dimension (an array at least ndim long)
107 /// \param [in] option consists of 3 parts - flag option and an option to create a new TVirtualFFT
108 /// 1. transform type option:
109 /// Available transform types are:
110 /// C2CForward, C2CBackward, C2R, R2C, R2HC, HC2R, DHT
111 /// see class description for details
112 /// 2. flag option: choosing how much time should be spent in planning the transform:
113 /// Possible options:
114 /// - "ES" (from "estimate") - no time in preparing the transform,
115 /// but probably sub-optimal performance
116 /// - "M" (from "measure") - some time spend in finding the optimal way
117 /// to do the transform
118 /// - "P" (from "patient") - more time spend in finding the optimal way
119 /// to do the transform
120 /// - "EX" (from "exhaustive") - the most optimal way is found
121 /// This option should be chosen depending on how many transforms of the
122 /// same size and type are going to be done.
123 /// Planning is only done once, for the first transform of this size and type.
124 /// 3. option allowing to choose between the global fgFFT and a new TVirtualFFT object
125 /// "" - default, changes and returns the global fgFFT variable
126 /// "K" (from "keep")- without touching the global fgFFT,
127 /// creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
128 ///
129 /// Examples of valid options: "R2C ES K", "C2CF M", "DHT P K", etc.
130 
132 {
133 
134  Int_t inputtype=0, currenttype=0;
135  TString opt = option;
136  opt.ToUpper();
137  //find the tranform flag
138  Option_t *flag;
139  flag = "ES";
140  if (opt.Contains("ES")) flag = "ES";
141  if (opt.Contains("M")) flag = "M";
142  if (opt.Contains("P")) flag = "P";
143  if (opt.Contains("EX")) flag = "EX";
144 
145  Int_t ndiff = 0;
146 
147  if (!opt.Contains("K")) {
148  if (fgFFT){
149  //if the global transform exists, check if it should be changed
150  if (fgFFT->GetNdim()!=ndim)
151  ndiff++;
152  else {
153  Int_t *ncurrent = fgFFT->GetN();
154  for (Int_t i=0; i<ndim; i++){
155  if (n[i]!=ncurrent[i])
156  ndiff++;
157  }
158  }
159  Option_t *t = fgFFT->GetType();
160  if (!opt.Contains(t)) {
161  if (opt.Contains("HC") || opt.Contains("DHT"))
162  inputtype = 1;
163  if (strcmp(t,"R2HC")==0 || strcmp(t,"HC2R")==0 || strcmp(t,"DHT")==0)
164  currenttype=1;
165 
166  if (!(inputtype==1 && currenttype==1))
167  ndiff++;
168  }
169  if (ndiff>0){
170  delete fgFFT;
171  fgFFT = 0;
172  }
173  }
174  }
175 
176  Int_t sign = 0;
177  if (opt.Contains("C2CB") || opt.Contains("C2R"))
178  sign = 1;
179  if (opt.Contains("C2CF") || opt.Contains("R2C"))
180  sign = -1;
181 
182  TVirtualFFT *fft = 0;
183  if (opt.Contains("K") || !fgFFT) {
184 
186 
187  TPluginHandler *h;
188  TString pluginname;
189  if (fgDefault.Length()==0) fgDefault="fftw";
190  if (strcmp(fgDefault.Data(),"fftw")==0) {
191  if (opt.Contains("C2C")) pluginname = "fftwc2c";
192  if (opt.Contains("C2R")) pluginname = "fftwc2r";
193  if (opt.Contains("R2C")) pluginname = "fftwr2c";
194  if (opt.Contains("HC") || opt.Contains("DHT")) pluginname = "fftwr2r";
195  if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
196  if (h->LoadPlugin()==-1) {
197  ::Error("TVirtualFFT::FFT", "handler not found");
198  return 0;
199  }
200  fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
201  if (!fft) {
202  ::Error("TVirtualFFT::FFT", "plugin failed to create TVirtualFFT object");
203  return 0;
204  }
205  Int_t *kind = new Int_t[1];
206  if (pluginname=="fftwr2r") {
207  if (opt.Contains("R2HC")) kind[0] = 10;
208  if (opt.Contains("HC2R")) kind[0] = 11;
209  if (opt.Contains("DHT")) kind[0] = 12;
210  }
211  fft->Init(flag, sign, kind);
212  if (!opt.Contains("K")) {
213  fgFFT = fft;
214  }
215  delete [] kind;
216  return fft;
217  }
218  else {
219  ::Error("TVirtualFFT::FFT", "plugin not found");
220  return 0;
221  }
222  }
223  } else {
224 
226 
227  //if the global transform already exists and just needs to be reinitialised
228  //with different parameters
229  if (fgFFT->GetSign()!=sign || !opt.Contains(fgFFT->GetTransformFlag()) || !opt.Contains(fgFFT->GetType())) {
230  Int_t *kind = new Int_t[1];
231  if (inputtype==1) {
232  if (opt.Contains("R2HC")) kind[0] = 10;
233  if (opt.Contains("HC2R")) kind[0] = 11;
234  if (opt.Contains("DHT")) kind[0] = 12;
235  }
236  fgFFT->Init(flag, sign, kind);
237  delete [] kind;
238  }
239  }
240  return fgFFT;
241 }
242 
243 ////////////////////////////////////////////////////////////////////////////////
244 ///Returns a pointer to a sine or cosine transform of requested size and kind
245 ///
246 ///Parameters:
247 /// \param [in] ndim number of transform dimensions
248 /// \param [in] n sizes of each dimension (an array at least ndim long)
249 /// \param [in] r2rkind transform kind for each dimension
250 /// 4 different kinds of sine and cosine transforms are available
251 /// - DCT-I - kind=0
252 /// - DCT-II - kind=1
253 /// - DCT-III - kind=2
254 /// - DCT-IV - kind=3
255 /// - DST-I - kind=4
256 /// - DST-II - kind=5
257 /// - DST-III - kind=6
258 /// - DST-IV - kind=7
259 /// \param [in] option : consists of 2 parts
260 /// - flag option and an option to create a new TVirtualFFT
261 /// - flag option: choosing how much time should be spent in planning the transform:
262 /// Possible options:
263 /// - "ES" (from "estimate") - no time in preparing the transform,
264 /// but probably sub-optimal performance
265 /// - "M" (from "measure") - some time spend in finding the optimal way
266 /// to do the transform
267 /// - "P" (from "patient") - more time spend in finding the optimal way
268 /// to do the transform
269 /// - "EX" (from "exhaustive") - the most optimal way is found
270 /// This option should be chosen depending on how many transforms of the
271 /// same size and type are going to be done.
272 /// Planning is only done once, for the first transform of this size and type.
273 /// - option allowing to choose between the global fgFFT and a new TVirtualFFT object
274 /// - "" - default, changes and returns the global fgFFT variable
275 /// - "K" (from "keep")- without touching the global fgFFT,
276 /// creates and returns a new TVirtualFFT*. User is then responsible for deleting it.
277 /// Examples of valid options: "ES K", "EX", etc
278 
280 {
281  TString opt = option;
282  //find the tranform flag
283  Option_t *flag;
284  flag = "ES";
285  if (opt.Contains("ES")) flag = "ES";
286  if (opt.Contains("M")) flag = "M";
287  if (opt.Contains("P")) flag = "P";
288  if (opt.Contains("EX")) flag = "EX";
289 
290  if (!opt.Contains("K")) {
291  if (fgFFT){
292  Int_t ndiff = 0;
293  if (fgFFT->GetNdim()!=ndim || strcmp(fgFFT->GetType(),"R2R")!=0)
294  ndiff++;
295  else {
296  Int_t *ncurrent = fgFFT->GetN();
297  for (Int_t i=0; i<ndim; i++) {
298  if (n[i] != ncurrent[i])
299  ndiff++;
300  }
301 
302  }
303  if (ndiff>0) {
304  delete fgFFT;
305  fgFFT = 0;
306  }
307  }
308  }
309  TVirtualFFT *fft = 0;
310 
312 
313  if (!fgFFT || opt.Contains("K")) {
314  TPluginHandler *h;
315  TString pluginname;
316  if (fgDefault.Length()==0) fgDefault="fftw";
317  if (strcmp(fgDefault.Data(),"fftw")==0) {
318  pluginname = "fftwr2r";
319  if ((h=gROOT->GetPluginManager()->FindHandler("TVirtualFFT", pluginname))) {
320  if (h->LoadPlugin()==-1){
321  ::Error("TVirtualFFT::SineCosine", "handler not found");
322  return 0;
323  }
324  fft = (TVirtualFFT*)h->ExecPlugin(3, ndim, n, kFALSE);
325  if (!fft) {
326  ::Error("TVirtualFFT::SineCosine", "plugin failed to create TVirtualFFT object");
327  return 0;
328  }
329  fft->Init(flag, 0, r2rkind);
330  if (!opt.Contains("K"))
331  fgFFT = fft;
332  return fft;
333  } else {
334  ::Error("TVirtualFFT::SineCosine", "handler not found");
335  return 0;
336  }
337  }
338  }
339 
340  //if (fgFFT->GetTransformFlag()!=flag)
341  fgFFT->Init(flag,0, r2rkind);
342  return fgFFT;
343 }
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// static: return current fgFFT
347 
349 {
350  if (fgFFT)
351  return fgFFT;
352  else{
353  ::Warning("TVirtualFFT::GetCurrentTransform", "fgFFT is not defined yet");
354  return 0;
355  }
356 }
357 
358 ////////////////////////////////////////////////////////////////////////////////
359 /// static: set the current transfrom to parameter
360 
362 {
363  fgFFT = fft;
364 }
365 
366 ////////////////////////////////////////////////////////////////////////////////
367 /// static: return the name of the default fft
368 
370 {
371  return fgDefault.Data();
372 }
373 
374 ////////////////////////////////////////////////////////////////////////////////
375 /// static: set name of default fft
376 
378 {
379  if (fgDefault == name) return;
380  delete fgFFT;
381  fgFFT = 0;
382  fgDefault = name;
383 }
n
const Int_t n
Definition: legend1.C:16
TVirtualFFT::SetTransform
static void SetTransform(TVirtualFFT *fft)
static: set the current transfrom to parameter
Definition: TVirtualFFT.cxx:361
Option_t
const char Option_t
Definition: RtypesCore.h:66
TVirtualFFT::GetCurrentTransform
static TVirtualFFT * GetCurrentTransform()
static: return current fgFFT
Definition: TVirtualFFT.cxx:348
TString::Data
const char * Data() const
Definition: TString.h:369
TVirtualFFT::SetDefaultFFT
static void SetDefaultFFT(const char *name="")
static: set name of default fft
Definition: TVirtualFFT.cxx:377
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TVirtualFFT::fgFFT
static TVirtualFFT * fgFFT
Definition: TVirtualFFT.h:91
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:890
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
TVirtualFFT::~TVirtualFFT
virtual ~TVirtualFFT()
destructor
Definition: TVirtualFFT.cxx:96
TString::Length
Ssiz_t Length() const
Definition: TString.h:410
TVirtualFFT.h
TString
Basic string class.
Definition: TString.h:136
TROOT.h
TString::ToUpper
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
TPluginHandler
Definition: TPluginManager.h:101
TVirtualFFT::fgDefault
static TString fgDefault
Definition: TVirtualFFT.h:92
TVirtualFFT::GetN
virtual Int_t * GetN() const =0
gROOTMutex
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:61
h
#define h(i)
Definition: RSha256.hxx:106
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TVirtualFFT::GetTransformFlag
virtual Option_t * GetTransformFlag() const =0
R__LOCKGUARD
#define R__LOCKGUARD(mutex)
Definition: TVirtualMutex.h:104
TVirtualFFT::GetSign
virtual Int_t GetSign() const =0
TVirtualFFT::Init
virtual void Init(Option_t *flag, Int_t sign, const Int_t *kind)=0
TVirtualFFT::SineCosine
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.
Definition: TVirtualFFT.cxx:279
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:876
TVirtualFFT::FFT
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
Definition: TVirtualFFT.cxx:131
TPluginManager.h
TVirtualFFT::GetNdim
virtual Int_t GetNdim() const =0
name
char name[80]
Definition: TGX11.cxx:110
TVirtualFFT::GetDefaultFFT
static const char * GetDefaultFFT()
static: return the name of the default fft
Definition: TVirtualFFT.cxx:369
TVirtualFFT::GetType
virtual Option_t * GetType() const =0
TVirtualFFT
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
gROOT
#define gROOT
Definition: TROOT.h:406
int
TError.h