]>
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$ | |
18 | */ | |
19 | ||
20 | #include "AliTPCMonitorFFT.h" | |
21 | ||
22 | ||
23 | ClassImp(AliTPCMonitorFFT) | |
24 | ||
25 | //__________________________________________________________________ | |
26 | AliTPCMonitorFFT::AliTPCMonitorFFT() | |
27 | { | |
28 | // Constuctor | |
29 | } | |
30 | ||
31 | //__________________________________________________________________ | |
32 | AliTPCMonitorFFT::~AliTPCMonitorFFT() | |
33 | { | |
34 | // Destructor | |
35 | } | |
36 | ||
37 | ||
38 | //__________________________________________________________________ | |
39 | Int_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 | //__________________________________________________________________ | |
49 | Int_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 | //__________________________________________________________________ | |
58 | Int_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 | //__________________________________________________________________ | |
81 | Int_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 | //__________________________________________________________________ | |
163 | Int_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 | //__________________________________________________________________ | |
201 | Int_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 |