removed TPC CA tracker
[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$
18*/
19
20#include "AliTPCMonitorFFT.h"
21
22
23ClassImp(AliTPCMonitorFFT)
24
25//__________________________________________________________________
26 AliTPCMonitorFFT::AliTPCMonitorFFT()
27{
28 // Constuctor
29}
30
31//__________________________________________________________________
32AliTPCMonitorFFT::~AliTPCMonitorFFT()
33{
34 // Destructor
35}
36
37
38//__________________________________________________________________
39Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
40{
41
42 // Wrapper function for forward fourier transformation
43 Int_t direction = 1;
44 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
45 return ret ;
46}
47
48//__________________________________________________________________
49Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
50{
51 // Wrapper function for backward fourier transformation
52 Int_t direction = -1;
53 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
54 return ret ;
55}
56
57//__________________________________________________________________
58Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
59{
60 // Wrapper function for inverse fourier transformation
61 Int_t direction = -1;
62 Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
63
64 if (ret)
65 {
66 return ret;
67 }
68
69 const Double_t norm = 1.0 / n;
70 size_t i;
71 for (i = 0; i < n; i++)
72 {
73 REAL(data,stride,i) *= norm;
74 IMAG(data,stride,i) *= norm;
75 }
76
77 return ret;
78}
79
80//__________________________________________________________________
81Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data, Int_t stride, size_t n, Int_t sign)
82{
83 // Wrapper function for fourier transformation depending on sign ( forward 1;backward -1 )
84 Int_t result ;
85 size_t dual;
86 size_t bit;
87 size_t logn = 0;
88 Int_t status;
89
90 if (n == 1) /* identity operation */
91 {
92 return 0 ;
93 }
94 /* make sure that n is a power of 2 */
95
96 result = FFTBinaryLogn(n) ;
97
98 if (result == -1)
99 {
100 cout << " not a power of 2 " << endl;
101 }
102 else
103 {
104 logn = result ;
105 }
106
107 /* apply fft recursion */
108
109 dual = n / 2;
110
111 for (bit = 0; bit < logn; bit++)
112 {
113 Double_t w_real = 1.0;
114 Double_t w_imag = 0.0;
115
116 const Double_t theta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
117
118 const Double_t s = sin (theta);
119 const Double_t t = sin (theta / 2.0);
120 const Double_t s2 = 2.0 * t * t;
121
122 size_t a, b;
123
124 for (b = 0; b < dual; b++)
125 {
126 for (a = 0; a < n; a+= 2 * dual)
127 {
128 const size_t i = b + a;
129 const size_t j = b + a + dual;
130
131 const Double_t t1_real = REAL(data,stride,i) + REAL(data,stride,j);
132 const Double_t t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
133 const Double_t t2_real = REAL(data,stride,i) - REAL(data,stride,j);
134 const Double_t t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
135
136 REAL(data,stride,i) = t1_real;
137 IMAG(data,stride,i) = t1_imag;
138 REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
139 IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
140 }
141
142 /* trignometric recurrence for w-> exp(i theta) w */
143
144 {
145 const Double_t tmp_real = w_real - s * w_imag - s2 * w_real;
146 const Double_t tmp_imag = w_imag + s * w_real - s2 * w_imag;
147 w_real = tmp_real;
148 w_imag = tmp_imag;
149 }
150 }
151 dual /= 2;
152 }
153
154 /* bit reverse the ordering of output data for decimation in
155 frequency algorithm */
156
157 //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
158 status = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
159 return 0;
160}
161
162//__________________________________________________________________
163Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn)
164{
165 // Wrapper function from gnu scientific library
166 /* This is the Goldrader bit-reversal algorithm */
167
168 size_t i;
169 size_t j = 0;
170
171 logn = 0 ; /* not needed for this algorithm */
172
173 for (i = 0; i < n - 1; i++)
174 {
175 size_t k = n / 2 ;
176
177 if (i < j)
178 {
179 const Double_t tmp_real = REAL(data,stride,i);
180 const Double_t tmp_imag = IMAG(data,stride,i);
181 REAL(data,stride,i) = REAL(data,stride,j);
182 IMAG(data,stride,i) = IMAG(data,stride,j);
183 REAL(data,stride,j) = tmp_real;
184 IMAG(data,stride,j) = tmp_imag;
185 }
186
187 while (k <= j)
188 {
189 j = j - k ;
190 k = k / 2 ;
191 }
192
193 j += k ;
194 }
195
196 return 0;
197
198}
199
200//__________________________________________________________________
201Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n)
202{
203
204 // Return log on base 2
205 size_t ntest ;
206 size_t binary_logn = 0 ;
207 size_t k = 1;
208
209 while (k < n)
210 {
211 k *= 2;
212 binary_logn++;
213 }
214
215 ntest = (1 << binary_logn) ;
216
217 if (n != ntest )
218 {
219 return -1 ; /* n is not a power of 2 */
220 }
221
222 return binary_logn;
223}
224