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
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 "TError.h"
87
90
92
93////////////////////////////////////////////////////////////////////////////////
94///destructor
95
97{
98 if (this==fgFFT)
99 fgFFT = nullptr;
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 = nullptr;
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 = nullptr;
183 if (opt.Contains("K") || !fgFFT) {
184
186
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 nullptr;
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 nullptr;
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 nullptr;
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 = nullptr;
306 }
307 }
308 }
309 TVirtualFFT *fft = nullptr;
310
312
313 if (!fgFFT || opt.Contains("K")) {
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 nullptr;
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 nullptr;
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 nullptr;
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 nullptr;
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 = nullptr;
382 fgDefault = name;
383}
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:61
#define gROOT
Definition: TROOT.h:404
#define R__LOCKGUARD(mutex)
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
Basic string class.
Definition: TString.h:136
Ssiz_t Length() const
Definition: TString.h:410
const char * Data() const
Definition: TString.h:369
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1163
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
TVirtualFFT is an interface class for Fast Fourier Transforms.
Definition: TVirtualFFT.h:88
virtual ~TVirtualFFT()
destructor
Definition: TVirtualFFT.cxx:96
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