Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFFTComplex.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 TFFTComplex
14///
15/// One of the interface classes to the FFTW package, can be used directly
16/// or via the TVirtualFFT class. Only the basic interface of FFTW is implemented.
17/// Computes complex input/output discrete Fourier transforms (DFT)
18/// in one or more dimensions. For the detailed information on the computed
19/// transforms please refer to the FFTW manual, chapter "What FFTW really computes".
20///
21/// How to use it:
22///
23/// 1. Create an instance of TFFTComplex - this will allocate input and output
24/// arrays (unless an in-place transform is specified)
25/// 2. Run the Init() function with the desired flags and settings
26/// 3. Set the data (via SetPoints(), SetPoint() or SetPointComplex() functions)
27/// 4. Run the Transform() function
28/// 5. Get the output (via GetPoints(), GetPoint() or GetPointComplex() functions)
29/// 6. Repeat steps 3)-5) as needed
30///
31/// For a transform of the same size, but with different flags or sign, rerun the Init()
32/// function and continue with steps 3)-5)
33///
34/// NOTE:
35/// 1. running Init() function will overwrite the input array! Don't set any data
36/// before running the Init() function
37/// 2. FFTW computes unnormalized transform, so doing a transform followed by
38/// its inverse will lead to the original array scaled by the transform size
39///
40////////////////////////////////////////////////////////////////////////////////
41
42#include "TFFTComplex.h"
43#include "fftw3.h"
44#include "TComplex.h"
45
46
48
49////////////////////////////////////////////////////////////////////////////////
50///default
51
53{
54 fIn = nullptr;
55 fOut = nullptr;
56 fPlan = nullptr;
57 fN = nullptr;
58 fNdim = 0;
59 fTotalSize = 0;
60 fSign = 1;
61}
62
63////////////////////////////////////////////////////////////////////////////////
64///For 1d transforms
65///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array
66
68{
69 fIn = fftw_malloc(sizeof(fftw_complex) *n);
70 if (!inPlace)
71 fOut = fftw_malloc(sizeof(fftw_complex) * n);
72 else
73 fOut = nullptr;
74 fN = new Int_t[1];
75 fN[0] = n;
76 fTotalSize = n;
77 fNdim = 1;
78 fSign = 1;
79 fPlan = nullptr;
80}
81
82////////////////////////////////////////////////////////////////////////////////
83///For multidim. transforms
84///Allocates memory for the input array, and, if inPlace = kFALSE, for the output array
85
87{
88 fNdim = ndim;
89 fTotalSize = 1;
90 fN = new Int_t[fNdim];
91 for (Int_t i=0; i<fNdim; i++){
92 fN[i] = n[i];
93 fTotalSize*=n[i];
94 }
95 fIn = fftw_malloc(sizeof(fftw_complex)*fTotalSize);
96 if (!inPlace)
97 fOut = fftw_malloc(sizeof(fftw_complex) * fTotalSize);
98 else
99 fOut = nullptr;
100 fSign = 1;
101 fPlan = nullptr;
102}
103
104////////////////////////////////////////////////////////////////////////////////
105///Destroys the data arrays and the plan. However, some plan information stays around
106///until the root session is over, and is reused if other plans of the same size are
107///created
108
110{
111 fftw_destroy_plan((fftw_plan)fPlan);
112 fPlan = nullptr;
113 fftw_free((fftw_complex*)fIn);
114 if (fOut)
115 fftw_free((fftw_complex*)fOut);
116 if (fN)
117 delete [] fN;
118}
119
120////////////////////////////////////////////////////////////////////////////////
121///Creates the fftw-plan
122///
123///NOTE: input and output arrays are overwritten during initialisation,
124/// so don't set any points, before running this function!!!!!
125///
126///2nd parameter: +1
127///
128///Argument kind is dummy and doesn't need to be specified
129///
130///Possible flag_options:
131/// - "ES" (from "estimate") - no time in preparing the transform, but probably sub-optimal
132/// performance
133/// - "M" (from "measure") - some time spend in finding the optimal way to do the transform
134/// - "P" (from "patient") - more time spend in finding the optimal way to do the transform
135/// - "EX" (from "exhaustive") - the most optimal way is found
136///This option should be chosen depending on how many transforms of the same size and
137///type are going to be done. Planning is only done once, for the first transform of this
138///size and type.
139
140void TFFTComplex::Init( Option_t *flags, Int_t sign,const Int_t* /*kind*/)
141{
142 fSign = sign;
143 fFlags = flags;
144
145 if (fPlan)
146 fftw_destroy_plan((fftw_plan)fPlan);
147 fPlan = nullptr;
148
149 if (fOut)
150 fPlan = (void*)fftw_plan_dft(fNdim, fN, (fftw_complex*)fIn, (fftw_complex*)fOut, sign,MapFlag(flags));
151 else
152 fPlan = (void*)fftw_plan_dft(fNdim, fN, (fftw_complex*)fIn, (fftw_complex*)fIn, sign, MapFlag(flags));
153}
154
155////////////////////////////////////////////////////////////////////////////////
156///Computes the transform, specified in Init() function
157
159{
160 if (fPlan)
161 fftw_execute((fftw_plan)fPlan);
162 else {
163 Error("Transform", "transform not initialised");
164 return;
165 }
166}
167
168////////////////////////////////////////////////////////////////////////////////
169///Copies the output(or input) into the argument array
170
172{
173 if (!fromInput){
174 for (Int_t i=0; i<2*fTotalSize; i+=2){
175 data[i] = ((fftw_complex*)fOut)[i/2][0];
176 data[i+1] = ((fftw_complex*)fOut)[i/2][1];
177 }
178 } else {
179 for (Int_t i=0; i<2*fTotalSize; i+=2){
180 data[i] = ((fftw_complex*)fIn)[i/2][0];
181 data[i+1] = ((fftw_complex*)fIn)[i/2][1];
182 }
183 }
184}
185
186////////////////////////////////////////////////////////////////////////////////
187///returns real and imaginary parts of the point `#ipoint`
188
189void TFFTComplex::GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
190{
191 if (fOut && !fromInput){
192 re = ((fftw_complex*)fOut)[ipoint][0];
193 im = ((fftw_complex*)fOut)[ipoint][1];
194 } else {
195 re = ((fftw_complex*)fIn)[ipoint][0];
196 im = ((fftw_complex*)fIn)[ipoint][1];
197 }
198}
199
200////////////////////////////////////////////////////////////////////////////////
201///For multidimensional transforms. Returns real and imaginary parts of the point `#ipoint`
202
203void TFFTComplex::GetPointComplex(const Int_t *ipoint, Double_t &re, Double_t &im, Bool_t fromInput) const
204{
205 Int_t ireal = ipoint[0];
206 for (Int_t i=0; i<fNdim-1; i++)
207 ireal=fN[i+1]*ireal + ipoint[i+1];
208
209 if (fOut && !fromInput){
210 re = ((fftw_complex*)fOut)[ireal][0];
211 im = ((fftw_complex*)fOut)[ireal][1];
212 } else {
213 re = ((fftw_complex*)fIn)[ireal][0];
214 im = ((fftw_complex*)fIn)[ireal][1];
215 }
216}
217
218////////////////////////////////////////////////////////////////////////////////
219///Copies real and imaginary parts of the output (input) into the argument arrays
220
222{
223 if (fOut && !fromInput){
224 for (Int_t i=0; i<fTotalSize; i++){
225 re[i] = ((fftw_complex*)fOut)[i][0];
226 im[i] = ((fftw_complex*)fOut)[i][1];
227 }
228 } else {
229 for (Int_t i=0; i<fTotalSize; i++){
230 re[i] = ((fftw_complex*)fIn)[i][0];
231 im[i] = ((fftw_complex*)fIn)[i][1];
232 }
233 }
234}
235
236////////////////////////////////////////////////////////////////////////////////
237///Copies the output(input) into the argument array
238
240{
241 if (fOut && !fromInput){
242 for (Int_t i=0; i<fTotalSize; i+=2){
243 data[i] = ((fftw_complex*)fOut)[i/2][0];
244 data[i+1] = ((fftw_complex*)fOut)[i/2][1];
245 }
246 } else {
247 for (Int_t i=0; i<fTotalSize; i+=2){
248 data[i] = ((fftw_complex*)fIn)[i/2][0];
249 data[i+1] = ((fftw_complex*)fIn)[i/2][1];
250 }
251 }
252}
253
254////////////////////////////////////////////////////////////////////////////////
255///sets real and imaginary parts of point # ipoint
256
258{
259 ((fftw_complex*)fIn)[ipoint][0]=re;
260 ((fftw_complex*)fIn)[ipoint][1]=im;
261}
262
263////////////////////////////////////////////////////////////////////////////////
264///For multidim. transforms. Sets real and imaginary parts of point # ipoint
265
266void TFFTComplex::SetPoint(const Int_t *ipoint, Double_t re, Double_t im)
267{
268 Int_t ireal = ipoint[0];
269 for (Int_t i=0; i<fNdim-1; i++)
270 ireal=fN[i+1]*ireal + ipoint[i+1];
271
272 ((fftw_complex*)fIn)[ireal][0]=re;
273 ((fftw_complex*)fIn)[ireal][1]=im;
274}
275
276////////////////////////////////////////////////////////////////////////////////
277
279{
280 ((fftw_complex*)fIn)[ipoint][0] = c.Re();
281 ((fftw_complex*)fIn)[ipoint][1] = c.Im();
282}
283
284////////////////////////////////////////////////////////////////////////////////
285///set all points. the values are copied. points should be ordered as follows:
286///[re_0, im_0, re_1, im_1, ..., re_n, im_n)
287
289{
290 for (Int_t i=0; i<2*fTotalSize-1; i+=2){
291 ((fftw_complex*)fIn)[i/2][0]=data[i];
292 ((fftw_complex*)fIn)[i/2][1]=data[i+1];
293 }
294}
295
296////////////////////////////////////////////////////////////////////////////////
297///set all points. the values are copied
298
299void TFFTComplex::SetPointsComplex(const Double_t *re_data, const Double_t *im_data)
300{
301 if (!fIn){
302 Error("SetPointsComplex", "Size is not set yet");
303 return;
304 }
305 for (Int_t i=0; i<fTotalSize; i++){
306 ((fftw_complex*)fIn)[i][0]=re_data[i];
307 ((fftw_complex*)fIn)[i][1]=im_data[i];
308 }
309}
310
311////////////////////////////////////////////////////////////////////////////////
312///allowed options:
313/// - "ES" - FFTW_ESTIMATE
314/// - "M" - FFTW_MEASURE
315/// - "P" - FFTW_PATIENT
316/// - "EX" - FFTW_EXHAUSTIVE
317
319{
320 TString opt = flag;
321 opt.ToUpper();
322 if (opt.Contains("ES"))
323 return FFTW_ESTIMATE;
324 if (opt.Contains("M"))
325 return FFTW_MEASURE;
326 if (opt.Contains("P"))
327 return FFTW_PATIENT;
328 if (opt.Contains("EX"))
329 return FFTW_EXHAUSTIVE;
330 return FFTW_ESTIMATE;
331}
#define c(i)
Definition RSha256.hxx:101
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 TFFTComplex.h:20
Int_t fTotalSize
Definition TFFTComplex.h:26
void * fPlan
Definition TFFTComplex.h:24
UInt_t MapFlag(Option_t *flag)
allowed options:
void SetPoints(const Double_t *data) override
set all points.
void SetPointComplex(Int_t ipoint, TComplex &c) override
void GetPointsComplex(Double_t *re, Double_t *im, Bool_t fromInput=kFALSE) const override
Copies real and imaginary parts of the output (input) into the argument arrays.
void GetPoints(Double_t *data, Bool_t fromInput=kFALSE) const override
Copies the output(or input) into the argument array.
void SetPointsComplex(const Double_t *re, const Double_t *im) override
set all points. the values are copied
void Transform() override
Computes the transform, specified in Init() function.
TString fFlags
Definition TFFTComplex.h:29
TFFTComplex()
default
void Init(Option_t *flags, Int_t sign, const Int_t *) override
Creates the fftw-plan.
void * fOut
Definition TFFTComplex.h:23
Int_t * fN
Definition TFFTComplex.h:27
~TFFTComplex() override
Destroys the data arrays and the plan.
void GetPointComplex(Int_t ipoint, Double_t &re, Double_t &im, Bool_t fromInput=kFALSE) const override
returns real and imaginary parts of the point #ipoint
void SetPoint(Int_t ipoint, Double_t re, Double_t im=0) override
sets real and imaginary parts of point # ipoint
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