Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFFTReal.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/// \class TFFTReal
14/// One of the interface classes to the FFTW package, can be used directly
15/// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented.
16///
17/// Computes transforms called r2r in FFTW manual:
18/// - transforms of real input and output in "halfcomplex" format i.e.
19/// real and imaginary parts for a transform of size n stored as
20/// (r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1)
21/// - discrete Hartley transform
22/// - sine and cosine transforms (DCT-I,II,III,IV and DST-I,II,III,IV)
23/// For the detailed information on the computed
24/// transforms please refer to the FFTW manual, chapter "What FFTW really computes".
25///
26/// How to use it:
27/// 1. Create an instance of TFFTReal - this will allocate input and output
28/// arrays (unless an in-place transform is specified)
29/// 2. Run the Init() function with the desired flags and settings (see function
30/// comments for possible kind parameters)
31/// 3. Set the data (via SetPoints()or SetPoint() functions)
32/// 4. Run the Transform() function
33/// 5. Get the output (via GetPoints() or GetPoint() functions)
34/// 6. Repeat steps 3)-5) as needed
35///
36/// For a transform of the same size, but of different kind (or with different flags),
37/// rerun the Init() function and continue with steps 3)-5)
38///
39/// NOTE:
40/// 1. running Init() function will overwrite the input array! Don't set any data
41/// before running the Init() function!
42/// 2. FFTW computes unnormalized transform, so doing a transform followed by
43/// its inverse will lead to the original array scaled BY:
44/// - transform size (N) for R2HC, HC2R, DHT transforms
45/// - 2*(N-1) for DCT-I (REDFT00)
46/// - 2*(N+1) for DST-I (RODFT00)
47/// - 2*N for the remaining transforms
48///
49/// Transform inverses:
50/// - R2HC<-->HC2R
51/// - DHT<-->DHT
52/// - DCT-I<-->DCT-I
53/// - DCT-II<-->DCT-III
54/// - DCT-IV<-->DCT-IV
55/// - DST-I<-->DST-I
56/// - DST-II<-->DST-III
57/// - DST-IV<-->DST-IV
58///
59////////////////////////////////////////////////////////////////////////////////
60
61#include "TFFTReal.h"
62#include "fftw3.h"
63
65
66////////////////////////////////////////////////////////////////////////////////
67///default
68
70{
71 fIn = nullptr;
72 fOut = nullptr;
73 fPlan = nullptr;
74 fN = nullptr;
75 fKind = nullptr;
76 fNdim = 0;
77 fTotalSize = 0;
78}
79
80////////////////////////////////////////////////////////////////////////////////
81///For 1d transforms
82///n here is the physical size of the transform (see FFTW manual for more details)
83
85{
86 fIn = fftw_malloc(sizeof(Double_t)*n);
87 if (inPlace) fOut = nullptr;
88 else fOut = fftw_malloc(sizeof(Double_t)*n);
89
90 fPlan = nullptr;
91 fNdim = 1;
92 fN = new Int_t[1];
93 fN[0] = n;
94 fKind = nullptr;
95 fTotalSize = n;
96}
97
98////////////////////////////////////////////////////////////////////////////////
99///For multidimensional transforms
100///1st parameter is the # of dimensions,
101///2nd is the sizes (physical) of the transform in each dimension
102
104{
105 fTotalSize = 1;
106 fNdim = ndim;
107 fN = new Int_t[ndim];
108 fKind = nullptr;
109 fPlan = nullptr;
110 for (Int_t i=0; i<ndim; i++){
111 fTotalSize*=n[i];
112 fN[i] = n[i];
113 }
114 fIn = fftw_malloc(sizeof(Double_t)*fTotalSize);
115 if (!inPlace)
116 fOut = fftw_malloc(sizeof(Double_t)*fTotalSize);
117 else
118 fOut = nullptr;
119}
120
121////////////////////////////////////////////////////////////////////////////////
122///clean-up
123
125{
126 fftw_destroy_plan((fftw_plan)fPlan);
127 fPlan = nullptr;
128 fftw_free(fIn);
129 fIn = nullptr;
130 if (fOut){
131 fftw_free(fOut);
132 }
133 fOut = nullptr;
134 if (fN)
135 delete [] fN;
136 fN = nullptr;
137 if (fKind)
138 fftw_free((fftw_r2r_kind*)fKind);
139 fKind = nullptr;
140}
141
142////////////////////////////////////////////////////////////////////////////////
143///Creates the fftw-plan
144///
145///NOTE: input and output arrays are overwritten during initialisation,
146/// so don't set any points, before running this function!!!!!
147///
148/// #### 1st parameter:
149/// Possible flag_options:
150///
151/// - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
152/// performance
153/// - "M" (from "measure") - some time spend in finding the optimal way to do the transform
154/// - "P" (from "patient") - more time spend in finding the optimal way to do the transform
155/// - "EX" (from "exhaustive") - the most optimal way is found
156///
157/// This option should be chosen depending on how many transforms of the same size and
158/// type are going to be done. Planning is only done once, for the first transform of this
159/// size and type.
160///
161/// #### 2nd parameter:
162/// is dummy and doesn't need to be specified
163///
164/// #### 3rd parameter:
165/// transform kind for each dimension
166/// 4 different kinds of sine and cosine transforms are available
167/// - DCT-I - kind=0
168/// - DCT-II - kind=1
169/// - DCT-III - kind=2
170/// - DCT-IV - kind=3
171/// - DST-I - kind=4
172/// - DST-II - kind=5
173/// - DSTIII - kind=6
174/// - DSTIV - kind=7
175
176void TFFTReal::Init( Option_t* flags,Int_t /*sign*/, const Int_t *kind)
177{
178 if (fPlan)
179 fftw_destroy_plan((fftw_plan)fPlan);
180 fPlan = nullptr;
181
182 if (!fKind)
183 fKind = (fftw_r2r_kind*)fftw_malloc(sizeof(fftw_r2r_kind)*fNdim);
184
185 if (MapOptions(kind)){
186 if (fOut)
187 fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fOut, (fftw_r2r_kind*)fKind, MapFlag(flags));
188 else
189 fPlan = (void*)fftw_plan_r2r(fNdim, fN, (Double_t*)fIn, (Double_t*)fIn, (fftw_r2r_kind*)fKind, MapFlag(flags));
190 fFlags = flags;
191 }
192}
193
194////////////////////////////////////////////////////////////////////////////////
195///Computes the transform, specified in Init() function
196
198{
199 if (fPlan)
200 fftw_execute((fftw_plan)fPlan);
201 else {
202 Error("Transform", "transform hasn't been initialised");
203 return;
204 }
205}
206
207////////////////////////////////////////////////////////////////////////////////
208///Returns the type of the transform
209
211{
212 if (!fKind) {
213 Error("GetType", "Type not defined yet (kind not set)");
214 return "";
215 }
216 if (((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC) return "R2HC";
217 if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R) return "HC2R";
218 if (((fftw_r2r_kind*)fKind)[0]==FFTW_DHT) return "DHT";
219 else return "R2R";
220}
221
222////////////////////////////////////////////////////////////////////////////////
223///Copies the output (or input) points into the provided array, that should
224///be big enough
225
227{
228 const Double_t * array = GetPointsReal(fromInput);
229 if (!array) return;
230 std::copy(array, array+fTotalSize, data);
231}
232
233////////////////////////////////////////////////////////////////////////////////
234///For 1d transformations. Returns point `#ipoint`
235
237{
238 if (ipoint<0 || ipoint>fTotalSize){
239 Error("GetPointReal", "No such point");
240 return 0;
241 }
242 const Double_t * array = GetPointsReal(fromInput);
243 return ( array ) ? array[ipoint] : 0;
244}
245
246////////////////////////////////////////////////////////////////////////////////
247///For multidim.transforms. Returns point `#ipoint`
248
249Double_t TFFTReal::GetPointReal(const Int_t *ipoint, Bool_t fromInput) const
250{
251 Int_t ireal = ipoint[0];
252 for (Int_t i=0; i<fNdim-1; i++)
253 ireal=fN[i+1]*ireal + ipoint[i+1];
254
255 const Double_t * array = GetPointsReal(fromInput);
256 return ( array ) ? array[ireal] : 0;
257}
258
259////////////////////////////////////////////////////////////////////////////////
260///Only for input of HC2R and output of R2HC
261
262void TFFTReal::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
263{
264 const Double_t * array = GetPointsReal(fromInput);
265 if (!array) return;
266 if ( ( ((fftw_r2r_kind*)fKind)[0]==FFTW_R2HC && !fromInput ) ||
267 ( ((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R && fromInput ) )
268 {
269 if (ipoint<fN[0]/2+1){
270 re = array[ipoint];
271 im = array[fN[0]-ipoint];
272 } else {
273 re = array[fN[0]-ipoint];
274 im = -array[ipoint];
275 }
276 if ((fN[0]%2)==0 && ipoint==fN[0]/2) im = 0;
277 }
278}
279////////////////////////////////////////////////////////////////////////////////
280///Only for input of HC2R and output of R2HC and for 1d
281
282void TFFTReal::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
283{
284 GetPointComplex(ipoint[0], re, im, fromInput);
285}
286
287////////////////////////////////////////////////////////////////////////////////
288///Returns the output (or input) array
289/// we have 4 different cases:
290///
291/// - fromInput = false; fOut = !NULL (transformed is not in place) : return fOut
292/// - fromInput = false; fOut = NULL (transformed is in place) : return fIn
293/// - fromInput = true; fOut = !NULL : return fIn
294/// - fromInput = true; fOut = NULL return an error since input array is overwritten
295
297{
298
299 if (!fromInput && fOut)
300 return (Double_t*)fOut;
301 else if (fromInput && !fOut) {
302 Error("GetPointsReal","Input array was destroyed");
303 return nullptr;
304 }
305 //R__ASSERT(fIn);
306 return (Double_t*)fIn;
307}
308
309////////////////////////////////////////////////////////////////////////////////
310
312{
313 if (ipoint<0 || ipoint>fTotalSize){
314 Error("SetPoint", "illegal point index");
315 return;
316 }
317 if (((fftw_r2r_kind*)fKind)[0]==FFTW_HC2R){
318 if ((fN[0]%2)==0 && ipoint==fN[0]/2)
319 ((Double_t*)fIn)[ipoint] = re;
320 else {
321 ((Double_t*)fIn)[ipoint] = re;
322 ((Double_t*)fIn)[fN[0]-ipoint]=im;
323 }
324 }
325 else
326 ((Double_t*)fIn)[ipoint]=re;
327}
328
329////////////////////////////////////////////////////////////////////////////////
330///Since multidimensional R2HC and HC2R transforms are not supported,
331///third parameter is dummy
332
333void TFFTReal::SetPoint(const Int_t *ipoint, Double_t re, Double_t /*im*/)
334{
335 Int_t ireal = ipoint[0];
336 for (Int_t i=0; i<fNdim-1; i++)
337 ireal=fN[i+1]*ireal + ipoint[i+1];
338 if (ireal < 0 || ireal >fTotalSize){
339 Error("SetPoint", "illegal point index");
340 return;
341 }
342 ((Double_t*)fIn)[ireal]=re;
343}
344
345////////////////////////////////////////////////////////////////////////////////
346///Sets all points
347
349{
350 for (Int_t i=0; i<fTotalSize; i++)
351 ((Double_t*)fIn)[i] = data[i];
352}
353
354////////////////////////////////////////////////////////////////////////////////
355///transfers the r2r_kind parameters to fftw type
356
358{
359 if (kind[0] == 10){
360 if (fNdim>1){
361 Error("Init", "Multidimensional R2HC transforms are not supported, use R2C interface instead");
362 return 0;
363 }
364 ((fftw_r2r_kind*)fKind)[0] = FFTW_R2HC;
365 }
366 else if (kind[0] == 11) {
367 if (fNdim>1){
368 Error("Init", "Multidimensional HC2R transforms are not supported, use C2R interface instead");
369 return 0;
370 }
371 ((fftw_r2r_kind*)fKind)[0] = FFTW_HC2R;
372 }
373 else if (kind[0] == 12) {
374 for (Int_t i=0; i<fNdim; i++)
375 ((fftw_r2r_kind*)fKind)[i] = FFTW_DHT;
376 }
377 else {
378 for (Int_t i=0; i<fNdim; i++){
379 switch (kind[i]) {
380 case 0: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT00; break;
381 case 1: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT01; break;
382 case 2: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT10; break;
383 case 3: ((fftw_r2r_kind*)fKind)[i] = FFTW_REDFT11; break;
384 case 4: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT00; break;
385 case 5: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT01; break;
386 case 6: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT10; break;
387 case 7: ((fftw_r2r_kind*)fKind)[i] = FFTW_RODFT11; break;
388 default:
389 ((fftw_r2r_kind*)fKind)[i] = FFTW_R2HC; break;
390 }
391 }
392 }
393 return 1;
394}
395
396////////////////////////////////////////////////////////////////////////////////
397///allowed options:
398/// - "ES" - FFTW_ESTIMATE
399/// - "M" - FFTW_MEASURE
400/// - "P" - FFTW_PATIENT
401/// - "EX" - FFTW_EXHAUSTIVE
402
404{
405 TString opt = flag;
406 opt.ToUpper();
407 if (opt.Contains("ES"))
408 return FFTW_ESTIMATE;
409 if (opt.Contains("M"))
410 return FFTW_MEASURE;
411 if (opt.Contains("P"))
412 return FFTW_PATIENT;
413 if (opt.Contains("EX"))
414 return FFTW_EXHAUSTIVE;
415 return FFTW_ESTIMATE;
416}
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
One of the interface classes to the FFTW package, can be used directly or via the TVirtualFFT class.
Definition TFFTReal.h:20
Int_t fNdim
Definition TFFTReal.h:25
void Transform() override
Computes the transform, specified in Init() function.
Definition TFFTReal.cxx:197
void SetPoint(Int_t ipoint, Double_t re, Double_t im=0) override
Definition TFFTReal.cxx:311
void * fPlan
Definition TFFTReal.h:24
UInt_t MapFlag(Option_t *flag)
allowed options:
Definition TFFTReal.cxx:403
Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const override
Returns the output (or input) array we have 4 different cases:
Definition TFFTReal.cxx:296
Int_t fTotalSize
Definition TFFTReal.h:26
void SetPoints(const Double_t *data) override
Sets all points.
Definition TFFTReal.cxx:348
void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const override
Copies the output (or input) points into the provided array, that should be big enough.
Definition TFFTReal.cxx:226
Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const override
For 1d transformations. Returns point #ipoint
Definition TFFTReal.cxx:236
Int_t * fN
Definition TFFTReal.h:27
Option_t * GetType() const override
Returns the type of the transform.
Definition TFFTReal.cxx:210
TString fFlags
Definition TFFTReal.h:29
void * fOut
Definition TFFTReal.h:23
TFFTReal()
default
Definition TFFTReal.cxx:69
void * fIn
Definition TFFTReal.h:22
Int_t MapOptions(const Int_t *kind)
transfers the r2r_kind parameters to fftw type
Definition TFFTReal.cxx:357
void Init(Option_t *flags, Int_t sign, const Int_t *kind) override
Creates the fftw-plan.
Definition TFFTReal.cxx:176
void * fKind
Definition TFFTReal.h:28
void GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const override
Only for input of HC2R and output of R2HC and for 1d.
Definition TFFTReal.cxx:282
~TFFTReal() override
clean-up
Definition TFFTReal.cxx:124
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:976
Basic string class.
Definition TString.h:139
void ToUpper()
Change string to upper case.
Definition TString.cxx:1195
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
const Int_t n
Definition legend1.C:16