Logo ROOT   6.18/05
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
15TVirtualFFT is an interface class for Fast Fourier Transforms.
16
17The default FFT library is FFTW. To use it, FFTW3 library should already
18be installed, and ROOT should be have fftw3 module enabled, with the directories
19of fftw3 include file and library specified (see installation instructions).
20Function SetDefaultFFT() allows to change the default library.
21
22## Available transform types:
23FFT:
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:
40Different types of transforms are specified by parameter kind of the SineCosine() static
41function. 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
52Formulas and detailed descriptions can be found in the chapter
53"What FFTW really computes" of the FFTW manual
54
55NOTE: 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:
61Call to the static function FFT returns a pointer to a fast Fourier transform
62with requested parameters. Call to the static function SineCosine returns a
63pointer 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~~~
80Different 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
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")) {
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}
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:59
#define gROOT
Definition: TROOT.h:414
#define R__LOCKGUARD(mutex)
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
const char * Data() const
Definition: TString.h:364
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
virtual ~TVirtualFFT()
destructor
Definition: TVirtualFFT.cxx:97
static void SetDefaultFFT(const char *name="")
static: set name of default fft
static void SetTransform(TVirtualFFT *fft)
static: set the current transfrom to parameter
virtual void Init(Option_t *flag, Int_t sign, const Int_t *kind)=0
static TVirtualFFT * FFT(Int_t ndim, Int_t *n, Option_t *option)
Returns a pointer to the FFT of requested size and type.
static TVirtualFFT * GetCurrentTransform()
static: return current fgFFT
virtual Int_t GetNdim() const =0
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.
virtual Option_t * GetTransformFlag() const =0
virtual Option_t * GetType() const =0
virtual Int_t GetSign() const =0
virtual Int_t * GetN() const =0
static TString fgDefault
Definition: TVirtualFFT.h:92
static TVirtualFFT * fgFFT
Definition: TVirtualFFT.h:91
static const char * GetDefaultFFT()
static: return the name of the default fft
const Int_t n
Definition: legend1.C:16