Minors...Just ordering the libraries
[u/mrichter/AliRoot.git] / TPC / AliTPCMonitorFFT.cxx
CommitLineData
48265b32 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
ca7b8371 18Revision 1.1 2007/09/17 10:23:31 cvetan
19New TPC monitoring package from Stefan Kniege. The monitoring package can be started by running TPCMonitor.C macro located in macros folder.
20
48265b32 21*/
22
ca7b8371 23////////////////////////////////////////////////////////////////////////
24//
25// AliTPCMonitorFFT class
26//
27// Wrapper class to perform Fast Fourier Transformations.
28// The code is based on the Gnu Scientific Library.
29// See documentation of gsl for further details.
30//
31// Author: Stefan Kniege, IKF, Frankfurt
32//
33//
34/////////////////////////////////////////////////////////////////////////
35
36
48265b32 37#include "AliTPCMonitorFFT.h"
38
39
40ClassImp(AliTPCMonitorFFT)
41
42//__________________________________________________________________
43 AliTPCMonitorFFT::AliTPCMonitorFFT()
44{
45 // Constuctor
46}
47
48//__________________________________________________________________
49AliTPCMonitorFFT::~AliTPCMonitorFFT()
50{
51 // Destructor
52}
53
54
55//__________________________________________________________________
56Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
57{
58
59 // Wrapper function for forward fourier transformation
60 Int_t direction = 1;
61 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
62 return ret ;
63}
64
65//__________________________________________________________________
66Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
67{
68 // Wrapper function for backward fourier transformation
69 Int_t direction = -1;
70 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
71 return ret ;
72}
73
74//__________________________________________________________________
75Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
76{
77 // Wrapper function for inverse fourier transformation
78 Int_t direction = -1;
79 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
80
81 if (ret)
82 {
83 return ret;
84 }
85
86 const Double_t norm = 1.0 / n;
87 size_t i;
88 for (i = 0; i < n; i++)
89 {
90 REAL(data,stride,i) *= norm;
91 IMAG(data,stride,i) *= norm;
92 }
93
94 return ret;
95}
96
97//__________________________________________________________________
98Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data, Int_t stride, size_t n, Int_t sign)
99{
100 // Wrapper function for fourier transformation depending on sign ( forward 1;backward -1 )
101 Int_t result ;
102 size_t dual;
103 size_t bit;
104 size_t logn = 0;
105 Int_t status;
106
107 if (n == 1) /* identity operation */
108 {
109 return 0 ;
110 }
111 /* make sure that n is a power of 2 */
112
113 result = FFTBinaryLogn(n) ;
114
115 if (result == -1)
116 {
117 cout << " not a power of 2 " << endl;
118 }
119 else
120 {
121 logn = result ;
122 }
123
124 /* apply fft recursion */
125
126 dual = n / 2;
127
128 for (bit = 0; bit < logn; bit++)
129 {
130 Double_t w_real = 1.0;
131 Double_t w_imag = 0.0;
132
133 const Double_t theta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
134
135 const Double_t s = sin (theta);
136 const Double_t t = sin (theta / 2.0);
137 const Double_t s2 = 2.0 * t * t;
138
139 size_t a, b;
140
141 for (b = 0; b < dual; b++)
142 {
143 for (a = 0; a < n; a+= 2 * dual)
144 {
145 const size_t i = b + a;
146 const size_t j = b + a + dual;
147
148 const Double_t t1_real = REAL(data,stride,i) + REAL(data,stride,j);
149 const Double_t t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
150 const Double_t t2_real = REAL(data,stride,i) - REAL(data,stride,j);
151 const Double_t t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
152
153 REAL(data,stride,i) = t1_real;
154 IMAG(data,stride,i) = t1_imag;
155 REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
156 IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
157 }
158
159 /* trignometric recurrence for w-> exp(i theta) w */
160
161 {
162 const Double_t tmp_real = w_real - s * w_imag - s2 * w_real;
163 const Double_t tmp_imag = w_imag + s * w_real - s2 * w_imag;
164 w_real = tmp_real;
165 w_imag = tmp_imag;
166 }
167 }
168 dual /= 2;
169 }
170
171 /* bit reverse the ordering of output data for decimation in
172 frequency algorithm */
173
174 //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
175 status = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
176 return 0;
177}
178
179//__________________________________________________________________
180Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn)
181{
182 // Wrapper function from gnu scientific library
183 /* This is the Goldrader bit-reversal algorithm */
184
185 size_t i;
186 size_t j = 0;
187
188 logn = 0 ; /* not needed for this algorithm */
189
190 for (i = 0; i < n - 1; i++)
191 {
192 size_t k = n / 2 ;
193
194 if (i < j)
195 {
196 const Double_t tmp_real = REAL(data,stride,i);
197 const Double_t tmp_imag = IMAG(data,stride,i);
198 REAL(data,stride,i) = REAL(data,stride,j);
199 IMAG(data,stride,i) = IMAG(data,stride,j);
200 REAL(data,stride,j) = tmp_real;
201 IMAG(data,stride,j) = tmp_imag;
202 }
203
204 while (k <= j)
205 {
206 j = j - k ;
207 k = k / 2 ;
208 }
209
210 j += k ;
211 }
212
213 return 0;
214
215}
216
217//__________________________________________________________________
218Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n)
219{
220
221 // Return log on base 2
222 size_t ntest ;
223 size_t binary_logn = 0 ;
224 size_t k = 1;
225
226 while (k < n)
227 {
228 k *= 2;
229 binary_logn++;
230 }
231
232 ntest = (1 << binary_logn) ;
233
234 if (n != ntest )
235 {
236 return -1 ; /* n is not a power of 2 */
237 }
238
239 return binary_logn;
240}
241