Logo ROOT  
Reference Guide
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 = 0;
72 fOut = 0;
73 fPlan = 0;
74 fN = 0;
75 fKind = 0;
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 = 0;
88 else fOut = fftw_malloc(sizeof(Double_t)*n);
89
90 fPlan = 0;
91 fNdim = 1;
92 fN = new Int_t[1];
93 fN[0] = n;
94 fKind = 0;
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 = 0;
109 fPlan = 0;
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 = 0;
119}
120
121////////////////////////////////////////////////////////////////////////////////
122///clean-up
123
125{
126 fftw_destroy_plan((fftw_plan)fPlan);
127 fPlan = 0;
128 fftw_free(fIn);
129 fIn = 0;
130 if (fOut){
131 fftw_free(fOut);
132 }
133 fOut = 0;
134 if (fN)
135 delete [] fN;
136 fN = 0;
137 if (fKind)
138 fftw_free((fftw_r2r_kind*)fKind);
139 fKind = 0;
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 = 0;
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
226void TFFTReal::GetPoints(Double_t *data, Bool_t fromInput) const
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 tranforms. 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 0;
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}
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
One of the interface classes to the FFTW package, can be used directly or via the TVirtualFFT class.
Definition: TFFTReal.h:20
virtual void Transform()
Computes the transform, specified in Init() function.
Definition: TFFTReal.cxx:197
Int_t fNdim
Definition: TFFTReal.h:25
virtual void Init(Option_t *flags, Int_t sign, const Int_t *kind)
Creates the fftw-plan.
Definition: TFFTReal.cxx:176
virtual Option_t * GetType() const
Returns the type of the transform.
Definition: TFFTReal.cxx:210
virtual void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const
Copies the output (or input) points into the provided array, that should be big enough.
Definition: TFFTReal.cxx:226
void * fPlan
Definition: TFFTReal.h:24
UInt_t MapFlag(Option_t *flag)
allowed options:
Definition: TFFTReal.cxx:403
Int_t fTotalSize
Definition: TFFTReal.h:26
virtual Double_t * GetPointsReal(Bool_t fromInput=kFALSE) const
Returns the output (or input) array we have 4 different cases:
Definition: TFFTReal.cxx:296
virtual void GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const
Only for input of HC2R and output of R2HC and for 1d.
Definition: TFFTReal.cxx:282
Int_t * fN
Definition: TFFTReal.h:27
TString fFlags
Definition: TFFTReal.h:29
void * fOut
Definition: TFFTReal.h:23
virtual void SetPoints(const Double_t *data)
Sets all points.
Definition: TFFTReal.cxx:348
TFFTReal()
default
Definition: TFFTReal.cxx:69
void * fIn
Definition: TFFTReal.h:22
virtual ~TFFTReal()
clean-up
Definition: TFFTReal.cxx:124
Int_t MapOptions(const Int_t *kind)
transfers the r2r_kind parameters to fftw type
Definition: TFFTReal.cxx:357
void * fKind
Definition: TFFTReal.h:28
virtual Double_t GetPointReal(Int_t ipoint, Bool_t fromInput=kFALSE) const
For 1d tranforms. Returns point #ipoint.
Definition: TFFTReal.cxx:236
virtual void SetPoint(Int_t ipoint, Double_t re, Double_t im=0)
Definition: TFFTReal.cxx:311
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
Basic string class.
Definition: TString.h:131
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
const Int_t n
Definition: legend1.C:16