]>
Commit | Line | Data |
---|---|---|
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$ | |
fb3305d1 | 18 | Revision 1.2 2007/10/12 13:36:27 cvetan |
19 | Coding convention fixes from Stefan | |
20 | ||
ca7b8371 | 21 | Revision 1.1 2007/09/17 10:23:31 cvetan |
22 | New TPC monitoring package from Stefan Kniege. The monitoring package can be started by running TPCMonitor.C macro located in macros folder. | |
23 | ||
48265b32 | 24 | */ |
25 | ||
ca7b8371 | 26 | //////////////////////////////////////////////////////////////////////// |
fb3305d1 | 27 | //// |
28 | //// AliTPCMonitorFFT class | |
29 | //// | |
30 | //// Wrapper class to perform Fast Fourier Transformations. | |
31 | //// The code is based on the Gnu Scientific Library. | |
32 | //// See documentation of gsl for further details. | |
33 | //// | |
34 | //// Author: Stefan Kniege, IKF, Frankfurt | |
35 | //// | |
36 | //// | |
ca7b8371 | 37 | ///////////////////////////////////////////////////////////////////////// |
38 | ||
39 | ||
48265b32 | 40 | #include "AliTPCMonitorFFT.h" |
41 | ||
42 | ||
43 | ClassImp(AliTPCMonitorFFT) | |
44 | ||
45 | //__________________________________________________________________ | |
46 | AliTPCMonitorFFT::AliTPCMonitorFFT() | |
47 | { | |
48 | // Constuctor | |
49 | } | |
50 | ||
51 | //__________________________________________________________________ | |
52 | AliTPCMonitorFFT::~AliTPCMonitorFFT() | |
53 | { | |
54 | // Destructor | |
55 | } | |
56 | ||
57 | ||
58 | //__________________________________________________________________ | |
59 | Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n) | |
60 | { | |
61 | ||
62 | // Wrapper function for forward fourier transformation | |
63 | Int_t direction = 1; | |
64 | Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction); | |
65 | return ret ; | |
66 | } | |
67 | ||
68 | //__________________________________________________________________ | |
69 | Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n) | |
70 | { | |
71 | // Wrapper function for backward fourier transformation | |
72 | Int_t direction = -1; | |
73 | Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction); | |
74 | return ret ; | |
75 | } | |
76 | ||
77 | //__________________________________________________________________ | |
78 | Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n) | |
79 | { | |
80 | // Wrapper function for inverse fourier transformation | |
81 | Int_t direction = -1; | |
82 | Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction); | |
83 | ||
84 | if (ret) | |
85 | { | |
86 | return ret; | |
87 | } | |
88 | ||
fb3305d1 | 89 | const Double_t knorm = 1.0 / n; |
48265b32 | 90 | size_t i; |
91 | for (i = 0; i < n; i++) | |
92 | { | |
fb3305d1 | 93 | REAL(data,stride,i) *= knorm; |
94 | IMAG(data,stride,i) *= knorm; | |
48265b32 | 95 | } |
96 | ||
97 | return ret; | |
98 | } | |
99 | ||
100 | //__________________________________________________________________ | |
101 | Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data, Int_t stride, size_t n, Int_t sign) | |
102 | { | |
103 | // Wrapper function for fourier transformation depending on sign ( forward 1;backward -1 ) | |
104 | Int_t result ; | |
105 | size_t dual; | |
106 | size_t bit; | |
107 | size_t logn = 0; | |
108 | Int_t status; | |
109 | ||
110 | if (n == 1) /* identity operation */ | |
111 | { | |
112 | return 0 ; | |
113 | } | |
114 | /* make sure that n is a power of 2 */ | |
115 | ||
116 | result = FFTBinaryLogn(n) ; | |
117 | ||
118 | if (result == -1) | |
119 | { | |
120 | cout << " not a power of 2 " << endl; | |
121 | } | |
122 | else | |
123 | { | |
124 | logn = result ; | |
125 | } | |
126 | ||
127 | /* apply fft recursion */ | |
128 | ||
129 | dual = n / 2; | |
130 | ||
131 | for (bit = 0; bit < logn; bit++) | |
132 | { | |
fb3305d1 | 133 | Double_t wreal = 1.0; |
134 | Double_t wimag = 0.0; | |
48265b32 | 135 | |
fb3305d1 | 136 | const Double_t ktheta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual)); |
48265b32 | 137 | |
fb3305d1 | 138 | const Double_t ks = sin (ktheta); |
139 | const Double_t kt = sin (ktheta / 2.0); | |
140 | const Double_t ks2 = 2.0 * kt * kt; | |
48265b32 | 141 | |
142 | size_t a, b; | |
143 | ||
144 | for (b = 0; b < dual; b++) | |
145 | { | |
146 | for (a = 0; a < n; a+= 2 * dual) | |
147 | { | |
fb3305d1 | 148 | const size_t ki = b + a; |
149 | const size_t kj = b + a + dual; | |
48265b32 | 150 | |
fb3305d1 | 151 | const Double_t kt1real = REAL(data,stride,ki) + REAL(data,stride,kj); |
152 | const Double_t kt1imag = IMAG(data,stride,ki) + IMAG(data,stride,kj); | |
153 | const Double_t kt2real = REAL(data,stride,ki) - REAL(data,stride,kj); | |
154 | const Double_t kt2imag = IMAG(data,stride,ki) - IMAG(data,stride,kj); | |
155 | ||
156 | REAL(data,stride,ki) = kt1real; | |
157 | IMAG(data,stride,ki) = kt1imag; | |
158 | REAL(data,stride,kj) = wreal*kt2real - wimag * kt2imag; | |
159 | IMAG(data,stride,kj) = wreal*kt2imag + wimag * kt2real; | |
48265b32 | 160 | } |
161 | ||
fb3305d1 | 162 | /* trignometric recurrence for w-> exp(i ktheta) w */ |
48265b32 | 163 | |
164 | { | |
fb3305d1 | 165 | const Double_t ktmpreal = wreal - ks * wimag - ks2 * wreal; |
166 | const Double_t ktmpimag = wimag + ks * wreal - ks2 * wimag; | |
167 | wreal = ktmpreal; | |
168 | wimag = ktmpimag; | |
48265b32 | 169 | } |
170 | } | |
171 | dual /= 2; | |
172 | } | |
173 | ||
174 | /* bit reverse the ordering of output data for decimation in | |
175 | frequency algorithm */ | |
176 | ||
177 | //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ; | |
178 | status = ComplexBitReverseOrderWrap(data, stride, n, logn) ; | |
179 | return 0; | |
180 | } | |
181 | ||
182 | //__________________________________________________________________ | |
fb3305d1 | 183 | Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn) const |
48265b32 | 184 | { |
185 | // Wrapper function from gnu scientific library | |
186 | /* This is the Goldrader bit-reversal algorithm */ | |
187 | ||
188 | size_t i; | |
189 | size_t j = 0; | |
190 | ||
191 | logn = 0 ; /* not needed for this algorithm */ | |
192 | ||
193 | for (i = 0; i < n - 1; i++) | |
194 | { | |
195 | size_t k = n / 2 ; | |
196 | ||
197 | if (i < j) | |
198 | { | |
fb3305d1 | 199 | const Double_t ktmpreal = REAL(data,stride,i); |
200 | const Double_t ktmpimag = IMAG(data,stride,i); | |
48265b32 | 201 | REAL(data,stride,i) = REAL(data,stride,j); |
202 | IMAG(data,stride,i) = IMAG(data,stride,j); | |
fb3305d1 | 203 | REAL(data,stride,j) = ktmpreal; |
204 | IMAG(data,stride,j) = ktmpimag; | |
48265b32 | 205 | } |
206 | ||
207 | while (k <= j) | |
208 | { | |
209 | j = j - k ; | |
210 | k = k / 2 ; | |
211 | } | |
212 | ||
213 | j += k ; | |
214 | } | |
215 | ||
216 | return 0; | |
217 | ||
218 | } | |
219 | ||
220 | //__________________________________________________________________ | |
fb3305d1 | 221 | Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n) const |
48265b32 | 222 | { |
223 | ||
224 | // Return log on base 2 | |
225 | size_t ntest ; | |
fb3305d1 | 226 | size_t binarylogn = 0 ; |
48265b32 | 227 | size_t k = 1; |
228 | ||
229 | while (k < n) | |
230 | { | |
231 | k *= 2; | |
fb3305d1 | 232 | binarylogn++; |
48265b32 | 233 | } |
234 | ||
fb3305d1 | 235 | ntest = (1 << binarylogn) ; |
48265b32 | 236 | |
237 | if (n != ntest ) | |
238 | { | |
239 | return -1 ; /* n is not a power of 2 */ | |
240 | } | |
241 | ||
fb3305d1 | 242 | return binarylogn; |
48265b32 | 243 | } |
244 |