Logo ROOT  

Reference Guide
TSpectrumTransform.cxx
Go to the documentation of this file.
1// @(#)root/spectrum:$Id$
2// Author: Miroslav Morhac 25/09/06
3
4/** \class TSpectrumTransform
5 \ingroup Spectrum
6 \brief Advanced 1-dimensional orthogonal transform functions
7 \author Miroslav Morhac
8
9 Class to carry out transforms of 1D spectra, its filtering and
10 enhancement. It allows to calculate classic Fourier, Cosine, Sin,
11 Hartley, Walsh, Haar transforms as well as mixed transforms (Fourier-
12 Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sin-Walsh
13 and Sin-Haar). All the transforms are fast.
14
15 The algorithms in this class have been published in the following
16 references:
17
18 1. C.V. Hampton, B. Lian, Wm. C. McHarris: Fast-Fourier-transform
19 spectral enhancement techniques for gamma-ray spectroscopy.NIM A353(1994) 280-284.
20 2. Morhac M., Matousek V., New adaptive Cosine-Walsh transform and
21 its application to nuclear data compression, IEEE Transactions on
22 Signal Processing 48 (2000) 2693.
23 3. Morhac M., Matousek V., Data compression using new fast adaptive
24 Cosine-Haar transforms, Digital Signal Processing 8 (1998) 63.
25 4. Morhac M., Matousek V.: Multidimensional nuclear data compression
26 using fast adaptive Walsh-Haar transform. Acta Physica Slovaca 51 (2001) 307.
27 */
28
29#include "TSpectrumTransform.h"
30#include "TMath.h"
31
32ClassImp(TSpectrumTransform);
33
34////////////////////////////////////////////////////////////////////////////////
35///default constructor
36
37TSpectrumTransform::TSpectrumTransform()
38{
39 fSize=0;
40 fTransformType=kTransformCos;
41 fDegree=0;
42 fDirection=kTransformForward;
43 fXmin=0;
44 fXmax=0;
45 fFilterCoeff=0;
46 fEnhanceCoeff=0.5;
47}
48
49////////////////////////////////////////////////////////////////////////////////
50/// the constructor creates TSpectrumTransform object. Its size must be > than zero and must be power of 2.
51/// It sets default transform type to be Cosine transform. Transform parameters can be changed using setter functions.
52
53TSpectrumTransform::TSpectrumTransform(Int_t size):TNamed("SpectrumTransform", "Miroslav Morhac transformer")
54{
55 Int_t j,n;
56 if (size <= 0){
57 Error ("TSpectrumTransform","Invalid length, must be > than 0");
58 return;
59 }
60 j = 0;
61 n = 1;
62 for (; n < size;) {
63 j += 1;
64 n = n * 2;
65 }
66 if (n != size){
67 Error ("TSpectrumTransform","Invalid length, must be power of 2");
68 return;
69 }
70 fSize=size;
71 fTransformType=kTransformCos;
72 fDegree=0;
73 fDirection=kTransformForward;
74 fXmin=size/4;
75 fXmax=size-1;
76 fFilterCoeff=0;
77 fEnhanceCoeff=0.5;
78}
79
80
81////////////////////////////////////////////////////////////////////////////////
82/// Destructor
83
84TSpectrumTransform::~TSpectrumTransform()
85{
86}
87
88////////////////////////////////////////////////////////////////////////////////
89/// This function calculates Haar transform of a part of data
90/// Function parameters:
91/// - working_space-pointer to vector of transformed data
92/// - num-length of processed data
93/// - direction-forward or inverse transform
94
95void TSpectrumTransform::Haar(Double_t *working_space, int num, int direction)
96{
97 int i, ii, li, l2, l3, j, jj, jj1, lj, iter, m, jmin, jmax;
98 Double_t a, b, c, wlk;
99 Double_t val;
100 for (i = 0; i < num; i++)
101 working_space[i + num] = 0;
102 i = num;
103 iter = 0;
104 for (; i > 1;) {
105 iter += 1;
106 i = i / 2;
107 }
108 if (direction == kTransformForward) {
109 for (m = 1; m <= iter; m++) {
110 li = iter + 1 - m;
111 l2 = (Int_t) TMath::Power(2, li - 1);
112 for (i = 0; i < (2 * l2); i++) {
113 working_space[num + i] = working_space[i];
114 }
115 for (j = 0; j < l2; j++) {
116 l3 = l2 + j;
117 jj = 2 * j;
118 val = working_space[jj + num] + working_space[jj + 1 + num];
119 working_space[j] = val;
120 val = working_space[jj + num] - working_space[jj + 1 + num];
121 working_space[l3] = val;
122 }
123 }
124 }
125 val = working_space[0];
126 val = val / TMath::Sqrt(TMath::Power(2, iter));
127 working_space[0] = val;
128 val = working_space[1];
129 val = val / TMath::Sqrt(TMath::Power(2, iter));
130 working_space[1] = val;
131 for (ii = 2; ii <= iter; ii++) {
132 i = ii - 1;
133 wlk = 1 / TMath::Sqrt(TMath::Power(2, iter - i));
134 jmin = (Int_t) TMath::Power(2, i);
135 jmax = (Int_t) TMath::Power(2, ii) - 1;
136 for (j = jmin; j <= jmax; j++) {
137 val = working_space[j];
138 a = val;
139 a = a * wlk;
140 val = a;
141 working_space[j] = val;
142 }
143 }
144 if (direction == kTransformInverse) {
145 for (m = 1; m <= iter; m++) {
146 a = 2;
147 b = m - 1;
148 c = TMath::Power(a, b);
149 li = (Int_t) c;
150 for (i = 0; i < (2 * li); i++) {
151 working_space[i + num] = working_space[i];
152 }
153 for (j = 0; j < li; j++) {
154 lj = li + j;
155 jj = 2 * (j + 1) - 1;
156 jj1 = jj - 1;
157 val = working_space[j + num] - working_space[lj + num];
158 working_space[jj] = val;
159 val = working_space[j + num] + working_space[lj + num];
160 working_space[jj1] = val;
161 }
162 }
163 }
164 return;
165}
166
167////////////////////////////////////////////////////////////////////////////////
168/// This function calculates Walsh transform of a part of data
169/// Function parameters:
170/// - working_space-pointer to vector of transformed data
171/// - num-length of processed data
172
173void TSpectrumTransform::Walsh(Double_t *working_space, int num)
174{
175 int i, m, nump = 1, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter;
176 Double_t a;
177 Double_t val1, val2;
178 for (i = 0; i < num; i++)
179 working_space[i + num] = 0;
180 i = num;
181 iter = 0;
182 for (; i > 1;) {
183 iter += 1;
184 i = i / 2;
185 }
186 for (m = 1; m <= iter; m++) {
187 if (m == 1)
188 nump = 1;
189
190 else
191 nump = nump * 2;
192 mnum = num / nump;
193 mnum2 = mnum / 2;
194 for (mp = 0; mp < nump; mp++) {
195 ib = mp * mnum;
196 for (mp2 = 0; mp2 < mnum2; mp2++) {
197 mnum21 = mnum2 + mp2 + ib;
198 iba = ib + mp2;
199 val1 = working_space[iba];
200 val2 = working_space[mnum21];
201 working_space[iba + num] = val1 + val2;
202 working_space[mnum21 + num] = val1 - val2;
203 }
204 }
205 for (i = 0; i < num; i++) {
206 working_space[i] = working_space[i + num];
207 }
208 }
209 a = num;
210 a = TMath::Sqrt(a);
211 val2 = a;
212 for (i = 0; i < num; i++) {
213 val1 = working_space[i];
214 val1 = val1 / val2;
215 working_space[i] = val1;
216 }
217 return;
218}
219
220////////////////////////////////////////////////////////////////////////////////
221/// This function carries out bit-reverse reordering of data
222/// Function parameters:
223/// - working_space-pointer to vector of processed data
224/// - num-length of processed data
225
226void TSpectrumTransform::BitReverse(Double_t *working_space, int num)
227{
228 int ipower[26];
229 int i, ib, il, ibd, ip, ifac, i1;
230 for (i = 0; i < num; i++) {
231 working_space[i + num] = working_space[i];
232 }
233 for (i = 1; i <= num; i++) {
234 ib = i - 1;
235 il = 1;
236 lab9:ibd = ib / 2;
237 ipower[il - 1] = 1;
238 if (ib == (ibd * 2))
239 ipower[il - 1] = 0;
240 if (ibd == 0)
241 goto lab10;
242 ib = ibd;
243 il = il + 1;
244 goto lab9;
245 lab10:ip = 1;
246 ifac = num;
247 for (i1 = 1; i1 <= il; i1++) {
248 ifac = ifac / 2;
249 ip = ip + ifac * ipower[i1 - 1];
250 }
251 working_space[ip - 1] = working_space[i - 1 + num];
252 }
253 return;
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// This function calculates Fourier based transform of a part of data
258/// Function parameters:
259/// - working_space-pointer to vector of transformed data
260/// - num-length of processed data
261/// - hartley-1 if it is Hartley transform, 0 otherwise
262/// - direction-forward or inverse transform
263
264void TSpectrumTransform::Fourier(Double_t *working_space, int num, int hartley,
265 int direction, int zt_clear)
266{
267 int nxp2, nxp, i, j, k, m, iter, mxp, j1, j2, n1, n2, it;
268 Double_t a, b, c, d, sign, wpwr, arg, wr, wi, tr, ti, pi =
269 3.14159265358979323846;
270 Double_t val1, val2, val3, val4;
271 if (direction == kTransformForward && zt_clear == 0) {
272 for (i = 0; i < num; i++)
273 working_space[i + num] = 0;
274 }
275 i = num;
276 iter = 0;
277 for (; i > 1;) {
278 iter += 1;
279 i = i / 2;
280 }
281 sign = -1;
282 if (direction == kTransformInverse)
283 sign = 1;
284 nxp2 = num;
285 for (it = 1; it <= iter; it++) {
286 nxp = nxp2;
287 nxp2 = nxp / 2;
288 a = nxp2;
289 wpwr = pi / a;
290 for (m = 1; m <= nxp2; m++) {
291 a = m - 1;
292 arg = a * wpwr;
293 wr = TMath::Cos(arg);
294 wi = sign * TMath::Sin(arg);
295 for (mxp = nxp; mxp <= num; mxp += nxp) {
296 j1 = mxp - nxp + m;
297 j2 = j1 + nxp2;
298 val1 = working_space[j1 - 1];
299 val2 = working_space[j2 - 1];
300 val3 = working_space[j1 - 1 + num];
301 val4 = working_space[j2 - 1 + num];
302 a = val1;
303 b = val2;
304 c = val3;
305 d = val4;
306 tr = a - b;
307 ti = c - d;
308 a = a + b;
309 val1 = a;
310 working_space[j1 - 1] = val1;
311 c = c + d;
312 val1 = c;
313 working_space[j1 - 1 + num] = val1;
314 a = tr * wr - ti * wi;
315 val1 = a;
316 working_space[j2 - 1] = val1;
317 a = ti * wr + tr * wi;
318 val1 = a;
319 working_space[j2 - 1 + num] = val1;
320 }
321 }
322 }
323 n2 = num / 2;
324 n1 = num - 1;
325 j = 1;
326 for (i = 1; i <= n1; i++) {
327 if (i >= j)
328 goto lab55;
329 val1 = working_space[j - 1];
330 val2 = working_space[j - 1 + num];
331 val3 = working_space[i - 1];
332 working_space[j - 1] = val3;
333 working_space[j - 1 + num] = working_space[i - 1 + num];
334 working_space[i - 1] = val1;
335 working_space[i - 1 + num] = val2;
336 lab55: k = n2;
337 lab60: if (k >= j) goto lab65;
338 j = j - k;
339 k = k / 2;
340 goto lab60;
341 lab65: j = j + k;
342 }
343 a = num;
344 a = TMath::Sqrt(a);
345 for (i = 0; i < num; i++) {
346 if (hartley == 0) {
347 val1 = working_space[i];
348 b = val1;
349 b = b / a;
350 val1 = b;
351 working_space[i] = val1;
352 b = working_space[i + num];
353 b = b / a;
354 working_space[i + num] = b;
355 }
356
357 else {
358 b = working_space[i];
359 c = working_space[i + num];
360 b = (b + c) / a;
361 working_space[i] = b;
362 working_space[i + num] = 0;
363 }
364 }
365 if (hartley == 1 && direction == kTransformInverse) {
366 for (i = 1; i < num; i++)
367 working_space[num - i + num] = working_space[i];
368 working_space[0 + num] = working_space[0];
369 for (i = 0; i < num; i++) {
370 working_space[i] = working_space[i + num];
371 working_space[i + num] = 0;
372 }
373 }
374 return;
375}
376
377////////////////////////////////////////////////////////////////////////////////
378/// This function carries out bit-reverse reordering for Haar transform
379/// Function parameters:
380/// - working_space-pointer to vector of processed data
381/// - shift-shift of position of processing
382/// - start-initial position of processed data
383/// - num-length of processed data
384
385void TSpectrumTransform::BitReverseHaar(Double_t *working_space, int shift, int num,
386 int start)
387{
388 int ipower[26];
389 int i, ib, il, ibd, ip, ifac, i1;
390 for (i = 0; i < num; i++) {
391 working_space[i + shift + start] = working_space[i + start];
392 working_space[i + shift + start + 2 * shift] =
393 working_space[i + start + 2 * shift];
394 }
395 for (i = 1; i <= num; i++) {
396 ib = i - 1;
397 il = 1;
398 lab9: ibd = ib / 2;
399 ipower[il - 1] = 1;
400 if (ib == (ibd * 2))
401 ipower[il - 1] = 0;
402 if (ibd == 0)
403 goto lab10;
404 ib = ibd;
405 il = il + 1;
406 goto lab9;
407 lab10: ip = 1;
408 ifac = num;
409 for (i1 = 1; i1 <= il; i1++) {
410 ifac = ifac / 2;
411 ip = ip + ifac * ipower[i1 - 1];
412 }
413 working_space[ip - 1 + start] =
414 working_space[i - 1 + shift + start];
415 working_space[ip - 1 + start + 2 * shift] =
416 working_space[i - 1 + shift + start + 2 * shift];
417 }
418 return;
419}
420
421////////////////////////////////////////////////////////////////////////////////
422/// This function calculates generalized (mixed) transforms of different degrees
423/// Function parameters:
424/// - working_space-pointer to vector of transformed data
425/// - zt_clear-flag to clear imaginary data before staring
426/// - num-length of processed data
427/// - degree-degree of transform (see manual)
428/// - type-type of mixed transform (see manual)
429
430int TSpectrumTransform::GeneralExe(Double_t *working_space, int zt_clear, int num,
431 int degree, int type)
432{
433 int i, j, k, m, nump, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter,
434 mp2step, mppom, ring;
435 Double_t a, b, c, d, wpwr, arg, wr, wi, tr, ti, pi =
436 3.14159265358979323846;
437 Double_t val1, val2, val3, val4, a0oldr = 0, b0oldr = 0, a0r, b0r;
438 if (zt_clear == 0) {
439 for (i = 0; i < num; i++)
440 working_space[i + 2 * num] = 0;
441 }
442 i = num;
443 iter = 0;
444 for (; i > 1;) {
445 iter += 1;
446 i = i / 2;
447 }
448 a = num;
449 wpwr = 2.0 * pi / a;
450 nump = num;
451 mp2step = 1;
452 ring = num;
453 for (i = 0; i < iter - degree; i++)
454 ring = ring / 2;
455 for (m = 1; m <= iter; m++) {
456 nump = nump / 2;
457 mnum = num / nump;
458 mnum2 = mnum / 2;
459 if (m > degree
460 && (type == kTransformFourierHaar
461 || type == kTransformWalshHaar
462 || type == kTransformCosHaar
463 || type == kTransformSinHaar))
464 mp2step *= 2;
465 if (ring > 1)
466 ring = ring / 2;
467 for (mp = 0; mp < nump; mp++) {
468 if (type != kTransformWalshHaar) {
469 mppom = mp;
470 mppom = mppom % ring;
471 a = 0;
472 j = 1;
473 k = num / 4;
474 for (i = 0; i < (iter - 1); i++) {
475 if ((mppom & j) != 0)
476 a = a + k;
477 j = j * 2;
478 k = k / 2;
479 }
480 arg = a * wpwr;
481 wr = TMath::Cos(arg);
482 wi = TMath::Sin(arg);
483 }
484
485 else {
486 wr = 1;
487 wi = 0;
488 }
489 ib = mp * mnum;
490 for (mp2 = 0; mp2 < mnum2; mp2++) {
491 mnum21 = mnum2 + mp2 + ib;
492 iba = ib + mp2;
493 if (mp2 % mp2step == 0) {
494 a0r = a0oldr;
495 b0r = b0oldr;
496 a0r = 1 / TMath::Sqrt(2.0);
497 b0r = 1 / TMath::Sqrt(2.0);
498 }
499
500 else {
501 a0r = 1;
502 b0r = 0;
503 }
504 val1 = working_space[iba];
505 val2 = working_space[mnum21];
506 val3 = working_space[iba + 2 * num];
507 val4 = working_space[mnum21 + 2 * num];
508 a = val1;
509 b = val2;
510 c = val3;
511 d = val4;
512 tr = a * a0r + b * b0r;
513 val1 = tr;
514 working_space[num + iba] = val1;
515 ti = c * a0r + d * b0r;
516 val1 = ti;
517 working_space[num + iba + 2 * num] = val1;
518 tr =
519 a * b0r * wr - c * b0r * wi - b * a0r * wr + d * a0r * wi;
520 val1 = tr;
521 working_space[num + mnum21] = val1;
522 ti =
523 c * b0r * wr + a * b0r * wi - d * a0r * wr - b * a0r * wi;
524 val1 = ti;
525 working_space[num + mnum21 + 2 * num] = val1;
526 }
527 }
528 for (i = 0; i < num; i++) {
529 val1 = working_space[num + i];
530 working_space[i] = val1;
531 val1 = working_space[num + i + 2 * num];
532 working_space[i + 2 * num] = val1;
533 }
534 }
535 return (0);
536}
537
538////////////////////////////////////////////////////////////////////////////////
539/// This function calculates inverse generalized (mixed) transforms
540/// Function parameters:
541/// - working_space-pointer to vector of transformed data
542/// - num-length of processed data
543/// - degree-degree of transform (see manual)
544/// - type-type of mixed transform (see manual)
545
546int TSpectrumTransform::GeneralInv(Double_t *working_space, int num, int degree,
547 int type)
548{
549 int i, j, k, m, nump =
550 1, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter, mp2step, mppom,
551 ring;
552 Double_t a, b, c, d, wpwr, arg, wr, wi, tr, ti, pi =
553 3.14159265358979323846;
554 Double_t val1, val2, val3, val4, a0oldr = 0, b0oldr = 0, a0r, b0r;
555 i = num;
556 iter = 0;
557 for (; i > 1;) {
558 iter += 1;
559 i = i / 2;
560 }
561 a = num;
562 wpwr = 2.0 * pi / a;
563 mp2step = 1;
564 if (type == kTransformFourierHaar || type == kTransformWalshHaar
565 || type == kTransformCosHaar || type == kTransformSinHaar) {
566 for (i = 0; i < iter - degree; i++)
567 mp2step *= 2;
568 }
569 ring = 1;
570 for (m = 1; m <= iter; m++) {
571 if (m == 1)
572 nump = 1;
573
574 else
575 nump = nump * 2;
576 mnum = num / nump;
577 mnum2 = mnum / 2;
578 if (m > iter - degree + 1)
579 ring *= 2;
580 for (mp = nump - 1; mp >= 0; mp--) {
581 if (type != kTransformWalshHaar) {
582 mppom = mp;
583 mppom = mppom % ring;
584 a = 0;
585 j = 1;
586 k = num / 4;
587 for (i = 0; i < (iter - 1); i++) {
588 if ((mppom & j) != 0)
589 a = a + k;
590 j = j * 2;
591 k = k / 2;
592 }
593 arg = a * wpwr;
594 wr = TMath::Cos(arg);
595 wi = TMath::Sin(arg);
596 }
597
598 else {
599 wr = 1;
600 wi = 0;
601 }
602 ib = mp * mnum;
603 for (mp2 = 0; mp2 < mnum2; mp2++) {
604 mnum21 = mnum2 + mp2 + ib;
605 iba = ib + mp2;
606 if (mp2 % mp2step == 0) {
607 a0r = a0oldr;
608 b0r = b0oldr;
609 a0r = 1 / TMath::Sqrt(2.0);
610 b0r = 1 / TMath::Sqrt(2.0);
611 }
612
613 else {
614 a0r = 1;
615 b0r = 0;
616 }
617 val1 = working_space[iba];
618 val2 = working_space[mnum21];
619 val3 = working_space[iba + 2 * num];
620 val4 = working_space[mnum21 + 2 * num];
621 a = val1;
622 b = val2;
623 c = val3;
624 d = val4;
625 tr = a * a0r + b * wr * b0r + d * wi * b0r;
626 val1 = tr;
627 working_space[num + iba] = val1;
628 ti = c * a0r + d * wr * b0r - b * wi * b0r;
629 val1 = ti;
630 working_space[num + iba + 2 * num] = val1;
631 tr = a * b0r - b * wr * a0r - d * wi * a0r;
632 val1 = tr;
633 working_space[num + mnum21] = val1;
634 ti = c * b0r - d * wr * a0r + b * wi * a0r;
635 val1 = ti;
636 working_space[num + mnum21 + 2 * num] = val1;
637 }
638 }
639 if (m <= iter - degree
640 && (type == kTransformFourierHaar
641 || type == kTransformWalshHaar
642 || type == kTransformCosHaar
643 || type == kTransformSinHaar))
644 mp2step /= 2;
645 for (i = 0; i < num; i++) {
646 val1 = working_space[num + i];
647 working_space[i] = val1;
648 val1 = working_space[num + i + 2 * num];
649 working_space[i + 2 * num] = val1;
650 }
651 }
652 return (0);
653}
654
655////////////////////////////////////////////////////////////////////////////////
656/// This function transforms the source spectrum. The calling program
657/// should fill in input parameters.
658/// Transformed data are written into dest spectrum.
659///
660/// Function parameters:
661/// - source-pointer to the vector of source spectrum, its length should
662/// be size except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR
663/// transform. These need 2*size length to supply real and
664/// imaginary coefficients.
665/// - destVector-pointer to the vector of dest data, its length should be
666/// size except for direct FOURIER, FOUR-WALSH, FOUR-HAAR. These
667/// need 2*size length to store real and imaginary coefficients
668///
669/// ### Transform methods
670///
671/// Goal: to analyse experimental data using orthogonal transforms
672///
673/// - orthogonal transforms can be successfully used for the processing of
674/// nuclear spectra (not only)
675///
676/// - they can be used to remove high frequency noise, to increase
677/// signal-to-background ratio as well as to enhance low intensity components [1],
678/// to carry out e.g. Fourier analysis etc.
679///
680/// - we have implemented the function for the calculation of the commonly
681/// used orthogonal transforms as well as functions for the filtration and
682/// enhancement of experimental data
683///
684/// #### References:
685///
686/// [1] C.V. Hampton, B. Lian, Wm. C.
687/// McHarris: Fast-Fourier-transform spectral enhancement techniques for gamma-ray
688/// spectroscopy. NIM A353 (1994) 280-284.
689///
690/// [2] Morhac; M., Matouoek V.,
691/// New adaptive Cosine-Walsh transform and its application to nuclear data
692/// compression, IEEE Transactions on Signal Processing 48 (2000) 2693.
693///
694/// [3] Morhac; M., Matouoek V.,
695/// Data compression using new fast adaptive Cosine-Haar transforms, Digital Signal
696/// Processing 8 (1998) 63.
697///
698/// [4] Morhac; M., Matouoek V.:
699/// Multidimensional nuclear data compression using fast adaptive Walsh-Haar
700/// transform. Acta Physica Slovaca 51 (2001) 307.
701///
702/// ### Example - script Transform.c:
703///
704/// \image html spectrumtransform_transform_image002.jpg Fig. 1 Original gamma-ray spectrum
705/// \image html spectrumtransform_transform_image003.jpg Fig. 2 Transformed spectrum from Fig. 1 using Cosine transform
706///
707/// #### Script:
708/// Example to illustrate Transform function (class TSpectrumTransform).
709/// To execute this example, do:
710///
711/// `root > .x Transform.C`
712///
713/// ~~~ {.cpp}
714/// #include <TSpectrum>
715/// #include <TSpectrumTransform>
716/// void Transform() {
717/// Int_t i;
718/// Double_t nbins = 4096;
719/// Double_t xmin = 0;
720/// Double_t xmax = (Double_t)nbins;
721/// Double_t * source = new Double_t[nbins];
722/// Double_t * dest = new Double_t[nbins];
723/// TH1F *h = new TH1F("h","Transformed spectrum using Cosine transform",nbins,xmin,xmax);
724/// TFile *f = new TFile("spectra/TSpectrum.root");
725/// h=(TH1F*) f->Get("transform1;1");
726/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
727/// TCanvas *Transform1 = gROOT->GetListOfCanvases()->FindObject("Transform1");
728/// if (!Transform1) Transform1 = new TCanvas("Transform","Transform1",10,10,1000,700);
729/// TSpectrum *s = new TSpectrum();
730/// TSpectrumTransform *t = new TSpectrumTransform(4096);
731/// t->SetTransformType(t->kTransformCos,0);
732/// t->SetDirection(t->kTransformForward);
733/// t->Transform(source,dest);
734/// for (i = 0; i < nbins; i++) h->SetBinContent(i + 1,dest[i]);
735/// h->SetLineColor(kRed);
736/// h->Draw("L");
737/// }
738/// ~~~
739
740void TSpectrumTransform::Transform(const Double_t *source, Double_t *destVector)
741{
742 int i, j=0, k = 1, m, l;
743 Double_t val;
744 Double_t a, b, pi = 3.14159265358979323846;
745 Double_t *working_space = 0;
746 if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
747 if (fTransformType >= kTransformCosWalsh)
748 fDegree += 1;
749 k = (Int_t) TMath::Power(2, fDegree);
750 j = fSize / k;
751 }
752 switch (fTransformType) {
753 case kTransformHaar:
754 case kTransformWalsh:
755 working_space = new Double_t[2 * fSize];
756 break;
757 case kTransformCos:
758 case kTransformSin:
759 case kTransformFourier:
760 case kTransformHartley:
761 case kTransformFourierWalsh:
762 case kTransformFourierHaar:
763 case kTransformWalshHaar:
764 working_space = new Double_t[4 * fSize];
765 break;
766 case kTransformCosWalsh:
767 case kTransformCosHaar:
768 case kTransformSinWalsh:
769 case kTransformSinHaar:
770 working_space = new Double_t[8 * fSize];
771 break;
772 }
773 if (fDirection == kTransformForward) {
774 switch (fTransformType) {
775 case kTransformHaar:
776 for (i = 0; i < fSize; i++) {
777 working_space[i] = source[i];
778 }
779 Haar(working_space, fSize, fDirection);
780 for (i = 0; i < fSize; i++) {
781 destVector[i] = working_space[i];
782 }
783 break;
784 case kTransformWalsh:
785 for (i = 0; i < fSize; i++) {
786 working_space[i] = source[i];
787 }
788 Walsh(working_space, fSize);
789 BitReverse(working_space, fSize);
790 for (i = 0; i < fSize; i++) {
791 destVector[i] = working_space[i];
792 }
793 break;
794 case kTransformCos:
795 fSize = 2 * fSize;
796 for (i = 1; i <= (fSize / 2); i++) {
797 val = source[i - 1];
798 working_space[i - 1] = val;
799 working_space[fSize - i] = val;
800 }
801 Fourier(working_space, fSize, 0, kTransformForward, 0);
802 for (i = 0; i < fSize / 2; i++) {
803 a = pi * (Double_t) i / (Double_t) fSize;
804 a = TMath::Cos(a);
805 b = working_space[i];
806 a = b / a;
807 working_space[i] = a;
808 working_space[i + fSize] = 0;
809 } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
810 for (i = 0; i < fSize / 2; i++) {
811 destVector[i] = working_space[i];
812 }
813 break;
814 case kTransformSin:
815 fSize = 2 * fSize;
816 for (i = 1; i <= (fSize / 2); i++) {
817 val = source[i - 1];
818 working_space[i - 1] = val;
819 working_space[fSize - i] = -val;
820 }
821 Fourier(working_space, fSize, 0, kTransformForward, 0);
822 for (i = 0; i < fSize / 2; i++) {
823 a = pi * (Double_t) i / (Double_t) fSize;
824 a = TMath::Sin(a);
825 b = working_space[i];
826 if (a != 0)
827 a = b / a;
828 working_space[i - 1] = a;
829 working_space[i + fSize] = 0;
830 }
831 working_space[fSize / 2 - 1] =
832 working_space[fSize / 2] / TMath::Sqrt(2.0);
833 for (i = 0; i < fSize / 2; i++) {
834 destVector[i] = working_space[i];
835 }
836 break;
837 case kTransformFourier:
838 for (i = 0; i < fSize; i++) {
839 working_space[i] = source[i];
840 }
841 Fourier(working_space, fSize, 0, kTransformForward, 0);
842 for (i = 0; i < 2 * fSize; i++) {
843 destVector[i] = working_space[i];
844 }
845 break;
846 case kTransformHartley:
847 for (i = 0; i < fSize; i++) {
848 working_space[i] = source[i];
849 }
850 Fourier(working_space, fSize, 1, kTransformForward, 0);
851 for (i = 0; i < fSize; i++) {
852 destVector[i] = working_space[i];
853 }
854 break;
855 case kTransformFourierWalsh:
856 case kTransformFourierHaar:
857 case kTransformWalshHaar:
858 case kTransformCosWalsh:
859 case kTransformCosHaar:
860 case kTransformSinWalsh:
861 case kTransformSinHaar:
862 for (i = 0; i < fSize; i++) {
863 val = source[i];
864 if (fTransformType == kTransformCosWalsh
865 || fTransformType == kTransformCosHaar) {
866 j = (Int_t) TMath::Power(2, fDegree) / 2;
867 k = i / j;
868 k = 2 * k * j;
869 working_space[k + i % j] = val;
870 working_space[k + 2 * j - 1 - i % j] = val;
871 }
872
873 else if (fTransformType == kTransformSinWalsh
874 || fTransformType == kTransformSinHaar) {
875 j = (Int_t) TMath::Power(2, fDegree) / 2;
876 k = i / j;
877 k = 2 * k * j;
878 working_space[k + i % j] = val;
879 working_space[k + 2 * j - 1 - i % j] = -val;
880 }
881
882 else
883 working_space[i] = val;
884 }
885 if (fTransformType == kTransformFourierWalsh
886 || fTransformType == kTransformFourierHaar
887 || fTransformType == kTransformWalshHaar) {
888 for (i = 0; i < j; i++)
889 BitReverseHaar(working_space, fSize, k, i * k);
890 GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
891 }
892
893 else if (fTransformType == kTransformCosWalsh
894 || fTransformType == kTransformCosHaar) {
895 m = (Int_t) TMath::Power(2, fDegree);
896 l = 2 * fSize / m;
897 for (i = 0; i < l; i++)
898 BitReverseHaar(working_space, 2 * fSize, m, i * m);
899 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
900 for (i = 0; i < fSize; i++) {
901 k = i / j;
902 k = 2 * k * j;
903 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
904 a = TMath::Cos(a);
905 b = working_space[k + i % j];
906 if (i % j == 0)
907 a = b / TMath::Sqrt(2.0);
908
909 else
910 a = b / a;
911 working_space[i] = a;
912 working_space[i + 2 * fSize] = 0;
913 }
914 }
915
916 else if (fTransformType == kTransformSinWalsh
917 || fTransformType == kTransformSinHaar) {
918 m = (Int_t) TMath::Power(2, fDegree);
919 l = 2 * fSize / m;
920 for (i = 0; i < l; i++)
921 BitReverseHaar(working_space, 2 * fSize, m, i * m);
922 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
923 for (i = 0; i < fSize; i++) {
924 k = i / j;
925 k = 2 * k * j;
926 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
927 a = TMath::Cos(a);
928 b = working_space[j + k + i % j];
929 if (i % j == 0)
930 a = b / TMath::Sqrt(2.0);
931
932 else
933 a = b / a;
934 working_space[j + k / 2 - i % j - 1] = a;
935 working_space[i + 2 * fSize] = 0;
936 }
937 }
938 if (fTransformType > kTransformWalshHaar)
939 k = (Int_t) TMath::Power(2, fDegree - 1);
940
941 else
942 k = (Int_t) TMath::Power(2, fDegree);
943 j = fSize / k;
944 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
945 working_space[fSize + i] = working_space[l + i / j];
946 working_space[fSize + i + 2 * fSize] =
947 working_space[l + i / j + 2 * fSize];
948 }
949 for (i = 0; i < fSize; i++) {
950 working_space[i] = working_space[fSize + i];
951 working_space[i + 2 * fSize] =
952 working_space[fSize + i + 2 * fSize];
953 }
954 for (i = 0; i < fSize; i++) {
955 destVector[i] = working_space[i];
956 }
957 if (fTransformType == kTransformFourierWalsh
958 || fTransformType == kTransformFourierHaar) {
959 for (i = 0; i < fSize; i++) {
960 destVector[fSize + i] = working_space[i + 2 * fSize];
961 }
962 }
963 break;
964 }
965 }
966
967 else if (fDirection == kTransformInverse) {
968 switch (fTransformType) {
969 case kTransformHaar:
970 for (i = 0; i < fSize; i++) {
971 working_space[i] = source[i];
972 }
973 Haar(working_space, fSize, fDirection);
974 for (i = 0; i < fSize; i++) {
975 destVector[i] = working_space[i];
976 }
977 break;
978 case kTransformWalsh:
979 for (i = 0; i < fSize; i++) {
980 working_space[i] = source[i];
981 }
982 BitReverse(working_space, fSize);
983 Walsh(working_space, fSize);
984 for (i = 0; i < fSize; i++) {
985 destVector[i] = working_space[i];
986 }
987 break;
988 case kTransformCos:
989 for (i = 0; i < fSize; i++) {
990 working_space[i] = source[i];
991 }
992 fSize = 2 * fSize;
993 working_space[0] = working_space[0] * TMath::Sqrt(2.0);
994 for (i = 0; i < fSize / 2; i++) {
995 a = pi * (Double_t) i / (Double_t) fSize;
996 b = TMath::Sin(a);
997 a = TMath::Cos(a);
998 working_space[i + fSize] = (Double_t) working_space[i] * b;
999 working_space[i] = (Double_t) working_space[i] * a;
1000 } for (i = 2; i <= (fSize / 2); i++) {
1001 working_space[fSize - i + 1] = working_space[i - 1];
1002 working_space[fSize - i + 1 + fSize] =
1003 -working_space[i - 1 + fSize];
1004 }
1005 working_space[fSize / 2] = 0;
1006 working_space[fSize / 2 + fSize] = 0;
1007 Fourier(working_space, fSize, 0, kTransformInverse, 1);
1008 for (i = 0; i < fSize / 2; i++) {
1009 destVector[i] = working_space[i];
1010 }
1011 break;
1012 case kTransformSin:
1013 for (i = 0; i < fSize; i++) {
1014 working_space[i] = source[i];
1015 }
1016 fSize = 2 * fSize;
1017 working_space[fSize / 2] =
1018 working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
1019 for (i = fSize / 2 - 1; i > 0; i--) {
1020 a = pi * (Double_t) i / (Double_t) fSize;
1021 working_space[i + fSize] =
1022 -(Double_t) working_space[i - 1] * TMath::Cos(a);
1023 working_space[i] =
1024 (Double_t) working_space[i - 1] * TMath::Sin(a);
1025 } for (i = 2; i <= (fSize / 2); i++) {
1026 working_space[fSize - i + 1] = working_space[i - 1];
1027 working_space[fSize - i + 1 + fSize] =
1028 -working_space[i - 1 + fSize];
1029 }
1030 working_space[0] = 0;
1031 working_space[fSize] = 0;
1032 working_space[fSize / 2 + fSize] = 0;
1033 Fourier(working_space, fSize, 0, kTransformInverse, 0);
1034 for (i = 0; i < fSize / 2; i++) {
1035 destVector[i] = working_space[i];
1036 }
1037 break;
1038 case kTransformFourier:
1039 for (i = 0; i < 2 * fSize; i++) {
1040 working_space[i] = source[i];
1041 }
1042 Fourier(working_space, fSize, 0, kTransformInverse, 0);
1043 for (i = 0; i < fSize; i++) {
1044 destVector[i] = working_space[i];
1045 }
1046 break;
1047 case kTransformHartley:
1048 for (i = 0; i < fSize; i++) {
1049 working_space[i] = source[i];
1050 }
1051 Fourier(working_space, fSize, 1, kTransformInverse, 0);
1052 for (i = 0; i < fSize; i++) {
1053 destVector[i] = working_space[i];
1054 }
1055 break;
1056 case kTransformFourierWalsh:
1057 case kTransformFourierHaar:
1058 case kTransformWalshHaar:
1059 case kTransformCosWalsh:
1060 case kTransformCosHaar:
1061 case kTransformSinWalsh:
1062 case kTransformSinHaar:
1063 for (i = 0; i < fSize; i++) {
1064 working_space[i] = source[i];
1065 }
1066 if (fTransformType == kTransformFourierWalsh
1067 || fTransformType == kTransformFourierHaar) {
1068 for (i = 0; i < fSize; i++) {
1069 working_space[i + 2 * fSize] = source[fSize + i];
1070 }
1071 }
1072 if (fTransformType > kTransformWalshHaar)
1073 k = (Int_t) TMath::Power(2, fDegree - 1);
1074
1075 else
1076 k = (Int_t) TMath::Power(2, fDegree);
1077 j = fSize / k;
1078 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1079 working_space[fSize + l + i / j] = working_space[i];
1080 working_space[fSize + l + i / j + 2 * fSize] =
1081 working_space[i + 2 * fSize];
1082 }
1083 for (i = 0; i < fSize; i++) {
1084 working_space[i] = working_space[fSize + i];
1085 working_space[i + 2 * fSize] =
1086 working_space[fSize + i + 2 * fSize];
1087 }
1088 if (fTransformType == kTransformFourierWalsh
1089 || fTransformType == kTransformFourierHaar
1090 || fTransformType == kTransformWalshHaar) {
1091 GeneralInv(working_space, fSize, fDegree, fTransformType);
1092 for (i = 0; i < j; i++)
1093 BitReverseHaar(working_space, fSize, k, i * k);
1094 }
1095
1096 else if (fTransformType == kTransformCosWalsh
1097 || fTransformType == kTransformCosHaar) {
1098 j = (Int_t) TMath::Power(2, fDegree) / 2;
1099 m = (Int_t) TMath::Power(2, fDegree);
1100 l = 2 * fSize / m;
1101 for (i = 0; i < fSize; i++) {
1102 k = i / j;
1103 k = 2 * k * j;
1104 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1105 if (i % j == 0) {
1106 working_space[2 * fSize + k + i % j] =
1107 working_space[i] * TMath::Sqrt(2.0);
1108 working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
1109 }
1110
1111 else {
1112 b = TMath::Sin(a);
1113 a = TMath::Cos(a);
1114 working_space[4 * fSize + 2 * fSize + k + i % j] =
1115 -(Double_t) working_space[i] * b;
1116 working_space[2 * fSize + k + i % j] =
1117 (Double_t) working_space[i] * a;
1118 } } for (i = 0; i < fSize; i++) {
1119 k = i / j;
1120 k = 2 * k * j;
1121 if (i % j == 0) {
1122 working_space[2 * fSize + k + j] = 0;
1123 working_space[4 * fSize + 2 * fSize + k + j] = 0;
1124 }
1125
1126 else {
1127 working_space[2 * fSize + k + 2 * j - i % j] =
1128 working_space[2 * fSize + k + i % j];
1129 working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
1130 -working_space[4 * fSize + 2 * fSize + k + i % j];
1131 }
1132 }
1133 for (i = 0; i < 2 * fSize; i++) {
1134 working_space[i] = working_space[2 * fSize + i];
1135 working_space[4 * fSize + i] =
1136 working_space[4 * fSize + 2 * fSize + i];
1137 }
1138 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1139 m = (Int_t) TMath::Power(2, fDegree);
1140 l = 2 * fSize / m;
1141 for (i = 0; i < l; i++)
1142 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1143 }
1144
1145 else if (fTransformType == kTransformSinWalsh
1146 || fTransformType == kTransformSinHaar) {
1147 j = (Int_t) TMath::Power(2, fDegree) / 2;
1148 m = (Int_t) TMath::Power(2, fDegree);
1149 l = 2 * fSize / m;
1150 for (i = 0; i < fSize; i++) {
1151 k = i / j;
1152 k = 2 * k * j;
1153 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1154 if (i % j == 0) {
1155 working_space[2 * fSize + k + j + i % j] =
1156 working_space[j + k / 2 - i % j -
1157 1] * TMath::Sqrt(2.0);
1158 working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
1159 }
1160
1161 else {
1162 b = TMath::Sin(a);
1163 a = TMath::Cos(a);
1164 working_space[4 * fSize + 2 * fSize + k + j + i % j] =
1165 -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
1166 working_space[2 * fSize + k + j + i % j] =
1167 (Double_t) working_space[j + k / 2 - i % j - 1] * a;
1168 } } for (i = 0; i < fSize; i++) {
1169 k = i / j;
1170 k = 2 * k * j;
1171 if (i % j == 0) {
1172 working_space[2 * fSize + k] = 0;
1173 working_space[4 * fSize + 2 * fSize + k] = 0;
1174 }
1175
1176 else {
1177 working_space[2 * fSize + k + i % j] =
1178 working_space[2 * fSize + k + 2 * j - i % j];
1179 working_space[4 * fSize + 2 * fSize + k + i % j] =
1180 -working_space[4 * fSize + 2 * fSize + k + 2 * j -
1181 i % j];
1182 }
1183 }
1184 for (i = 0; i < 2 * fSize; i++) {
1185 working_space[i] = working_space[2 * fSize + i];
1186 working_space[4 * fSize + i] =
1187 working_space[4 * fSize + 2 * fSize + i];
1188 }
1189 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1190 for (i = 0; i < l; i++)
1191 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1192 }
1193 for (i = 0; i < fSize; i++) {
1194 if (fTransformType >= kTransformCosWalsh) {
1195 k = i / j;
1196 k = 2 * k * j;
1197 val = working_space[k + i % j];
1198 }
1199
1200 else
1201 val = working_space[i];
1202 destVector[i] = val;
1203 }
1204 break;
1205 }
1206 }
1207 delete[]working_space;
1208 return;
1209}
1210
1211////////////////////////////////////////////////////////////////////////////////
1212/// This function transforms the source spectrum. The calling program
1213/// should fill in input parameters. Then it sets transformed
1214/// coefficients in the given region (fXmin, fXmax) to the given
1215/// fFilterCoeff and transforms it back.
1216/// Filtered data are written into dest spectrum.
1217///
1218/// Function parameters:
1219/// - source-pointer to the vector of source spectrum, its length should
1220/// be size except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR
1221/// transform. These need 2*size length to supply real and
1222/// imaginary coefficients.
1223/// - destVector-pointer to the vector of dest data, its length should be
1224/// size except for direct FOURIER, FOUR-WALSH, FOUR-HAAR. These
1225/// need 2*size length to store real and imaginary coefficients
1226///
1227/// ### Example - script Filter.c:
1228///
1229/// \image html spectrumtransform_filter_image001.jpg Fig. 1 Original spectrum (black line) and filtered spectrum (red line) using Cosine transform and zonal filtration (channels 2048-4095 were set to 0)
1230///
1231/// #### Script:
1232///
1233/// Example to illustrate FilterZonal function (class TSpectrumTransform).
1234/// To execute this example, do:
1235///
1236/// `root > .x Filter.C`
1237///
1238/// ~~~ {.cpp}
1239/// #include <TSpectrum>
1240/// #include <TSpectrumTransform>
1241/// void Filter() {
1242/// Int_t i;
1243/// Double_t nbins = 4096;
1244/// Double_t xmin = 0;
1245/// Double_t xmax = (Double_t)nbins;
1246/// Double_t * source = new Double_t[nbins];
1247/// Double_t * dest = new Double_t[nbins];
1248/// TH1F *h = new TH1F("h","Zonal filtering using Cosine transform",nbins,xmin,xmax);
1249/// TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1250/// TFile *f = new TFile("spectra/TSpectrum.root");
1251/// h=(TH1F*) f->Get("transform1;1");
1252/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1253/// TCanvas *Transform1 = gROOT->GetListOfCanvases()->FindObject("Transform1");
1254/// if (!Transform1) Transform1 = new TCanvas("Transform","Transform1",10,10,1000,700);
1255/// h->SetAxisRange(700,1024);
1256/// h->Draw("L");
1257/// TSpectrum *s = new TSpectrum();
1258/// TSpectrumTransform *t = new TSpectrumTransform(4096);
1259/// t->SetTransformType(t->kTransformCos,0);
1260/// t->SetRegion(2048, 4095);
1261/// t->FilterZonal(source,dest);
1262/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,dest[i]);
1263/// d->SetLineColor(kRed);
1264/// d->Draw("SAME L");
1265/// }
1266/// ~~~
1267
1268void TSpectrumTransform::FilterZonal(const Double_t *source, Double_t *destVector)
1269{
1270 int i, j=0, k = 1, m, l;
1271 Double_t val;
1272 Double_t *working_space = 0;
1273 Double_t a, b, pi = 3.14159265358979323846, old_area, new_area;
1274 if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
1275 if (fTransformType >= kTransformCosWalsh)
1276 fDegree += 1;
1277 k = (Int_t) TMath::Power(2, fDegree);
1278 j = fSize / k;
1279 }
1280 switch (fTransformType) {
1281 case kTransformHaar:
1282 case kTransformWalsh:
1283 working_space = new Double_t[2 * fSize];
1284 break;
1285 case kTransformCos:
1286 case kTransformSin:
1287 case kTransformFourier:
1288 case kTransformHartley:
1289 case kTransformFourierWalsh:
1290 case kTransformFourierHaar:
1291 case kTransformWalshHaar:
1292 working_space = new Double_t[4 * fSize];
1293 break;
1294 case kTransformCosWalsh:
1295 case kTransformCosHaar:
1296 case kTransformSinWalsh:
1297 case kTransformSinHaar:
1298 working_space = new Double_t[8 * fSize];
1299 break;
1300 }
1301 switch (fTransformType) {
1302 case kTransformHaar:
1303 for (i = 0; i < fSize; i++) {
1304 working_space[i] = source[i];
1305 }
1306 Haar(working_space, fSize, kTransformForward);
1307 break;
1308 case kTransformWalsh:
1309 for (i = 0; i < fSize; i++) {
1310 working_space[i] = source[i];
1311 }
1312 Walsh(working_space, fSize);
1313 BitReverse(working_space, fSize);
1314 break;
1315 case kTransformCos:
1316 fSize = 2 * fSize;
1317 for (i = 1; i <= (fSize / 2); i++) {
1318 val = source[i - 1];
1319 working_space[i - 1] = val;
1320 working_space[fSize - i] = val;
1321 }
1322 Fourier(working_space, fSize, 0, kTransformForward, 0);
1323 for (i = 0; i < fSize / 2; i++) {
1324 a = pi * (Double_t) i / (Double_t) fSize;
1325 a = TMath::Cos(a);
1326 b = working_space[i];
1327 a = b / a;
1328 working_space[i] = a;
1329 working_space[i + fSize] = 0;
1330 } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
1331 fSize = fSize / 2;
1332 break;
1333 case kTransformSin:
1334 fSize = 2 * fSize;
1335 for (i = 1; i <= (fSize / 2); i++) {
1336 val = source[i - 1];
1337 working_space[i - 1] = val;
1338 working_space[fSize - i] = -val;
1339 }
1340 Fourier(working_space, fSize, 0, kTransformForward, 0);
1341 for (i = 0; i < fSize / 2; i++) {
1342 a = pi * (Double_t) i / (Double_t) fSize;
1343 a = TMath::Sin(a);
1344 b = working_space[i];
1345 if (a != 0)
1346 a = b / a;
1347 working_space[i - 1] = a;
1348 working_space[i + fSize] = 0;
1349 }
1350 working_space[fSize / 2 - 1] =
1351 working_space[fSize / 2] / TMath::Sqrt(2.0);
1352 fSize = fSize / 2;
1353 break;
1354 case kTransformFourier:
1355 for (i = 0; i < fSize; i++) {
1356 working_space[i] = source[i];
1357 }
1358 Fourier(working_space, fSize, 0, kTransformForward, 0);
1359 break;
1360 case kTransformHartley:
1361 for (i = 0; i < fSize; i++) {
1362 working_space[i] = source[i];
1363 }
1364 Fourier(working_space, fSize, 1, kTransformForward, 0);
1365 break;
1366 case kTransformFourierWalsh:
1367 case kTransformFourierHaar:
1368 case kTransformWalshHaar:
1369 case kTransformCosWalsh:
1370 case kTransformCosHaar:
1371 case kTransformSinWalsh:
1372 case kTransformSinHaar:
1373 for (i = 0; i < fSize; i++) {
1374 val = source[i];
1375 if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1376 j = (Int_t) TMath::Power(2, fDegree) / 2;
1377 k = i / j;
1378 k = 2 * k * j;
1379 working_space[k + i % j] = val;
1380 working_space[k + 2 * j - 1 - i % j] = val;
1381 }
1382
1383 else if (fTransformType == kTransformSinWalsh
1384 || fTransformType == kTransformSinHaar) {
1385 j = (Int_t) TMath::Power(2, fDegree) / 2;
1386 k = i / j;
1387 k = 2 * k * j;
1388 working_space[k + i % j] = val;
1389 working_space[k + 2 * j - 1 - i % j] = -val;
1390 }
1391
1392 else
1393 working_space[i] = val;
1394 }
1395 if (fTransformType == kTransformFourierWalsh
1396 || fTransformType == kTransformFourierHaar
1397 || fTransformType == kTransformWalshHaar) {
1398 for (i = 0; i < j; i++)
1399 BitReverseHaar(working_space, fSize, k, i * k);
1400 GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
1401 }
1402
1403 else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1404 m = (Int_t) TMath::Power(2, fDegree);
1405 l = 2 * fSize / m;
1406 for (i = 0; i < l; i++)
1407 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1408 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1409 for (i = 0; i < fSize; i++) {
1410 k = i / j;
1411 k = 2 * k * j;
1412 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1413 a = TMath::Cos(a);
1414 b = working_space[k + i % j];
1415 if (i % j == 0)
1416 a = b / TMath::Sqrt(2.0);
1417
1418 else
1419 a = b / a;
1420 working_space[i] = a;
1421 working_space[i + 2 * fSize] = 0;
1422 }
1423 }
1424
1425 else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
1426 m = (Int_t) TMath::Power(2, fDegree);
1427 l = 2 * fSize / m;
1428 for (i = 0; i < l; i++)
1429 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1430 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1431 for (i = 0; i < fSize; i++) {
1432 k = i / j;
1433 k = 2 * k * j;
1434 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1435 a = TMath::Cos(a);
1436 b = working_space[j + k + i % j];
1437 if (i % j == 0)
1438 a = b / TMath::Sqrt(2.0);
1439
1440 else
1441 a = b / a;
1442 working_space[j + k / 2 - i % j - 1] = a;
1443 working_space[i + 2 * fSize] = 0;
1444 }
1445 }
1446 if (fTransformType > kTransformWalshHaar)
1447 k = (Int_t) TMath::Power(2, fDegree - 1);
1448
1449 else
1450 k = (Int_t) TMath::Power(2, fDegree);
1451 j = fSize / k;
1452 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1453 working_space[fSize + i] = working_space[l + i / j];
1454 working_space[fSize + i + 2 * fSize] =
1455 working_space[l + i / j + 2 * fSize];
1456 }
1457 for (i = 0; i < fSize; i++) {
1458 working_space[i] = working_space[fSize + i];
1459 working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
1460 }
1461 break;
1462 }
1463 if (!working_space) return;
1464 for (i = 0, old_area = 0; i < fSize; i++) {
1465 old_area += working_space[i];
1466 }
1467 for (i = 0, new_area = 0; i < fSize; i++) {
1468 if (i >= fXmin && i <= fXmax)
1469 working_space[i] = fFilterCoeff;
1470 new_area += working_space[i];
1471 }
1472 if (new_area != 0) {
1473 a = old_area / new_area;
1474 for (i = 0; i < fSize; i++) {
1475 working_space[i] *= a;
1476 }
1477 }
1478 if (fTransformType == kTransformFourier) {
1479 for (i = 0, old_area = 0; i < fSize; i++) {
1480 old_area += working_space[i + fSize];
1481 }
1482 for (i = 0, new_area = 0; i < fSize; i++) {
1483 if (i >= fXmin && i <= fXmax)
1484 working_space[i + fSize] = fFilterCoeff;
1485 new_area += working_space[i + fSize];
1486 }
1487 if (new_area != 0) {
1488 a = old_area / new_area;
1489 for (i = 0; i < fSize; i++) {
1490 working_space[i + fSize] *= a;
1491 }
1492 }
1493 }
1494
1495 else if (fTransformType == kTransformFourierWalsh
1496 || fTransformType == kTransformFourierHaar) {
1497 for (i = 0, old_area = 0; i < fSize; i++) {
1498 old_area += working_space[i + 2 * fSize];
1499 }
1500 for (i = 0, new_area = 0; i < fSize; i++) {
1501 if (i >= fXmin && i <= fXmax)
1502 working_space[i + 2 * fSize] = fFilterCoeff;
1503 new_area += working_space[i + 2 * fSize];
1504 }
1505 if (new_area != 0) {
1506 a = old_area / new_area;
1507 for (i = 0; i < fSize; i++) {
1508 working_space[i + 2 * fSize] *= a;
1509 }
1510 }
1511 }
1512 switch (fTransformType) {
1513 case kTransformHaar:
1514 Haar(working_space, fSize, kTransformInverse);
1515 for (i = 0; i < fSize; i++) {
1516 destVector[i] = working_space[i];
1517 }
1518 break;
1519 case kTransformWalsh:
1520 BitReverse(working_space, fSize);
1521 Walsh(working_space, fSize);
1522 for (i = 0; i < fSize; i++) {
1523 destVector[i] = working_space[i];
1524 }
1525 break;
1526 case kTransformCos:
1527 fSize = 2 * fSize;
1528 working_space[0] = working_space[0] * TMath::Sqrt(2.0);
1529 for (i = 0; i < fSize / 2; i++) {
1530 a = pi * (Double_t) i / (Double_t) fSize;
1531 b = TMath::Sin(a);
1532 a = TMath::Cos(a);
1533 working_space[i + fSize] = (Double_t) working_space[i] * b;
1534 working_space[i] = (Double_t) working_space[i] * a;
1535 } for (i = 2; i <= (fSize / 2); i++) {
1536 working_space[fSize - i + 1] = working_space[i - 1];
1537 working_space[fSize - i + 1 + fSize] =
1538 -working_space[i - 1 + fSize];
1539 }
1540 working_space[fSize / 2] = 0;
1541 working_space[fSize / 2 + fSize] = 0;
1542 Fourier(working_space, fSize, 0, kTransformInverse, 1);
1543 for (i = 0; i < fSize / 2; i++) {
1544 destVector[i] = working_space[i];
1545 }
1546 break;
1547 case kTransformSin:
1548 fSize = 2 * fSize;
1549 working_space[fSize / 2] =
1550 working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
1551 for (i = fSize / 2 - 1; i > 0; i--) {
1552 a = pi * (Double_t) i / (Double_t) fSize;
1553 working_space[i + fSize] =
1554 -(Double_t) working_space[i - 1] * TMath::Cos(a);
1555 working_space[i] = (Double_t) working_space[i - 1] * TMath::Sin(a);
1556 } for (i = 2; i <= (fSize / 2); i++) {
1557 working_space[fSize - i + 1] = working_space[i - 1];
1558 working_space[fSize - i + 1 + fSize] =
1559 -working_space[i - 1 + fSize];
1560 }
1561 working_space[0] = 0;
1562 working_space[fSize] = 0;
1563 working_space[fSize / 2 + fSize] = 0;
1564 Fourier(working_space, fSize, 0, kTransformInverse, 0);
1565 for (i = 0; i < fSize / 2; i++) {
1566 destVector[i] = working_space[i];
1567 }
1568 break;
1569 case kTransformFourier:
1570 Fourier(working_space, fSize, 0, kTransformInverse, 0);
1571 for (i = 0; i < fSize; i++) {
1572 destVector[i] = working_space[i];
1573 }
1574 break;
1575 case kTransformHartley:
1576 Fourier(working_space, fSize, 1, kTransformInverse, 0);
1577 for (i = 0; i < fSize; i++) {
1578 destVector[i] = working_space[i];
1579 }
1580 break;
1581 case kTransformFourierWalsh:
1582 case kTransformFourierHaar:
1583 case kTransformWalshHaar:
1584 case kTransformCosWalsh:
1585 case kTransformCosHaar:
1586 case kTransformSinWalsh:
1587 case kTransformSinHaar:
1588 if (fTransformType > kTransformWalshHaar)
1589 k = (Int_t) TMath::Power(2, fDegree - 1);
1590
1591 else
1592 k = (Int_t) TMath::Power(2, fDegree);
1593 j = fSize / k;
1594 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1595 working_space[fSize + l + i / j] = working_space[i];
1596 working_space[fSize + l + i / j + 2 * fSize] =
1597 working_space[i + 2 * fSize];
1598 }
1599 for (i = 0; i < fSize; i++) {
1600 working_space[i] = working_space[fSize + i];
1601 working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
1602 }
1603 if (fTransformType == kTransformFourierWalsh
1604 || fTransformType == kTransformFourierHaar
1605 || fTransformType == kTransformWalshHaar) {
1606 GeneralInv(working_space, fSize, fDegree, fTransformType);
1607 for (i = 0; i < j; i++)
1608 BitReverseHaar(working_space, fSize, k, i * k);
1609 }
1610
1611 else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1612 j = (Int_t) TMath::Power(2, fDegree) / 2;
1613 m = (Int_t) TMath::Power(2, fDegree);
1614 l = 2 * fSize / m;
1615 for (i = 0; i < fSize; i++) {
1616 k = i / j;
1617 k = 2 * k * j;
1618 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1619 if (i % j == 0) {
1620 working_space[2 * fSize + k + i % j] =
1621 working_space[i] * TMath::Sqrt(2.0);
1622 working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
1623 }
1624
1625 else {
1626 b = TMath::Sin(a);
1627 a = TMath::Cos(a);
1628 working_space[4 * fSize + 2 * fSize + k + i % j] =
1629 -(Double_t) working_space[i] * b;
1630 working_space[2 * fSize + k + i % j] =
1631 (Double_t) working_space[i] * a;
1632 } } for (i = 0; i < fSize; i++) {
1633 k = i / j;
1634 k = 2 * k * j;
1635 if (i % j == 0) {
1636 working_space[2 * fSize + k + j] = 0;
1637 working_space[4 * fSize + 2 * fSize + k + j] = 0;
1638 }
1639
1640 else {
1641 working_space[2 * fSize + k + 2 * j - i % j] =
1642 working_space[2 * fSize + k + i % j];
1643 working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
1644 -working_space[4 * fSize + 2 * fSize + k + i % j];
1645 }
1646 }
1647 for (i = 0; i < 2 * fSize; i++) {
1648 working_space[i] = working_space[2 * fSize + i];
1649 working_space[4 * fSize + i] =
1650 working_space[4 * fSize + 2 * fSize + i];
1651 }
1652 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1653 m = (Int_t) TMath::Power(2, fDegree);
1654 l = 2 * fSize / m;
1655 for (i = 0; i < l; i++)
1656 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1657 }
1658
1659 else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
1660 j = (Int_t) TMath::Power(2, fDegree) / 2;
1661 m = (Int_t) TMath::Power(2, fDegree);
1662 l = 2 * fSize / m;
1663 for (i = 0; i < fSize; i++) {
1664 k = i / j;
1665 k = 2 * k * j;
1666 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1667 if (i % j == 0) {
1668 working_space[2 * fSize + k + j + i % j] =
1669 working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
1670 working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
1671 }
1672
1673 else {
1674 b = TMath::Sin(a);
1675 a = TMath::Cos(a);
1676 working_space[4 * fSize + 2 * fSize + k + j + i % j] =
1677 -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
1678 working_space[2 * fSize + k + j + i % j] =
1679 (Double_t) working_space[j + k / 2 - i % j - 1] * a;
1680 } } for (i = 0; i < fSize; i++) {
1681 k = i / j;
1682 k = 2 * k * j;
1683 if (i % j == 0) {
1684 working_space[2 * fSize + k] = 0;
1685 working_space[4 * fSize + 2 * fSize + k] = 0;
1686 }
1687
1688 else {
1689 working_space[2 * fSize + k + i % j] =
1690 working_space[2 * fSize + k + 2 * j - i % j];
1691 working_space[4 * fSize + 2 * fSize + k + i % j] =
1692 -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
1693 }
1694 }
1695 for (i = 0; i < 2 * fSize; i++) {
1696 working_space[i] = working_space[2 * fSize + i];
1697 working_space[4 * fSize + i] =
1698 working_space[4 * fSize + 2 * fSize + i];
1699 }
1700 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1701 for (i = 0; i < l; i++)
1702 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1703 }
1704 for (i = 0; i < fSize; i++) {
1705 if (fTransformType >= kTransformCosWalsh) {
1706 k = i / j;
1707 k = 2 * k * j;
1708 val = working_space[k + i % j];
1709 }
1710
1711 else
1712 val = working_space[i];
1713 destVector[i] = val;
1714 }
1715 break;
1716 }
1717 delete[]working_space;
1718 return;
1719}
1720
1721
1722/////////////////////////////////////////////////////////////////////////////////
1723/// This function transforms the source spectrum. The calling program
1724/// should fill in input parameters. Then it multiplies transformed
1725/// coefficients in the given region (fXmin, fXmax) by the given
1726/// fEnhanceCoeff and transforms it back
1727/// Processed data are written into dest spectrum.
1728///
1729/// Function parameters:
1730/// - source-pointer to the vector of source spectrum, its length should
1731/// be size except for inverse FOURIER, FOUR-WALSh, FOUR-HAAR
1732/// transform. These need 2*size length to supply real and
1733/// imaginary coefficients.
1734/// - destVector-pointer to the vector of dest data, its length should be
1735/// size except for direct FOURIER, FOUR-WALSh, FOUR-HAAR. These
1736/// need 2*size length to store real and imaginary coefficients
1737///
1738/// ### Example - script Enhance.c:
1739///
1740/// \image html spectrumtransform_enhance_image001.jpg Fig. 1 Original spectrum (black line) and enhanced spectrum (red line) using Cosine transform (channels 0-1024 were multiplied by 2)
1741///
1742/// #### Script:
1743///
1744/// Example to illustrate Enhance function (class TSpectrumTransform).
1745/// To execute this example, do:
1746///
1747/// `root > .x Enhance.C`
1748///
1749/// ~~~ {.cpp}
1750/// void Enhance() {
1751/// Int_t i;
1752/// Double_t nbins = 4096;
1753/// Double_t xmin = 0;
1754/// Double_t xmax = (Double_t)nbins;
1755/// Double_t * source = new Double_t[nbins];
1756/// Double_t * dest = new Double_t[nbins];
1757/// TH1F *h = new TH1F("h","Enhancement using Cosine transform",nbins,xmin,xmax);
1758/// TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1759/// TFile *f = new TFile("spectra/TSpectrum.root");
1760/// h=(TH1F*) f->Get("transform1;1");
1761/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1762/// TCanvas *Transform1 = gROOT->GetListOfCanvases()->FindObject("Transform1");
1763/// if (!Transform1) Transform1 = new TCanvas("Transform","Transform1",10,10,1000,700);
1764/// h->SetAxisRange(700,1024);
1765/// h->Draw("L");
1766/// TSpectrum *s = new TSpectrum();
1767/// TSpectrumTransform *t = new TSpectrumTransform(4096);
1768/// t->SetTransformType(t->kTransformCos,0);
1769/// t->SetRegion(0, 1024);
1770/// t->SetEnhanceCoeff(2);
1771/// t->Enhance(source,dest);
1772/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,dest[i]);
1773/// d->SetLineColor(kRed);
1774/// d->Draw("SAME L");
1775/// }
1776/// ~~~
1777
1778void TSpectrumTransform::Enhance(const Double_t *source, Double_t *destVector)
1779{
1780 int i, j=0, k = 1, m, l;
1781 Double_t val;
1782 Double_t *working_space = 0;
1783 Double_t a, b, pi = 3.14159265358979323846, old_area, new_area;
1784 if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
1785 if (fTransformType >= kTransformCosWalsh)
1786 fDegree += 1;
1787 k = (Int_t) TMath::Power(2, fDegree);
1788 j = fSize / k;
1789 }
1790 switch (fTransformType) {
1791 case kTransformHaar:
1792 case kTransformWalsh:
1793 working_space = new Double_t[2 * fSize];
1794 break;
1795 case kTransformCos:
1796 case kTransformSin:
1797 case kTransformFourier:
1798 case kTransformHartley:
1799 case kTransformFourierWalsh:
1800 case kTransformFourierHaar:
1801 case kTransformWalshHaar:
1802 working_space = new Double_t[4 * fSize];
1803 break;
1804 case kTransformCosWalsh:
1805 case kTransformCosHaar:
1806 case kTransformSinWalsh:
1807 case kTransformSinHaar:
1808 working_space = new Double_t[8 * fSize];
1809 break;
1810 }
1811 switch (fTransformType) {
1812 case kTransformHaar:
1813 for (i = 0; i < fSize; i++) {
1814 working_space[i] = source[i];
1815 }
1816 Haar(working_space, fSize, kTransformForward);
1817 break;
1818 case kTransformWalsh:
1819 for (i = 0; i < fSize; i++) {
1820 working_space[i] = source[i];
1821 }
1822 Walsh(working_space, fSize);
1823 BitReverse(working_space, fSize);
1824 break;
1825 case kTransformCos:
1826 fSize = 2 * fSize;
1827 for (i = 1; i <= (fSize / 2); i++) {
1828 val = source[i - 1];
1829 working_space[i - 1] = val;
1830 working_space[fSize - i] = val;
1831 }
1832 Fourier(working_space, fSize, 0, kTransformForward, 0);
1833 for (i = 0; i < fSize / 2; i++) {
1834 a = pi * (Double_t) i / (Double_t) fSize;
1835 a = TMath::Cos(a);
1836 b = working_space[i];
1837 a = b / a;
1838 working_space[i] = a;
1839 working_space[i + fSize] = 0;
1840 } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
1841 fSize = fSize / 2;
1842 break;
1843 case kTransformSin:
1844 fSize = 2 * fSize;
1845 for (i = 1; i <= (fSize / 2); i++) {
1846 val = source[i - 1];
1847 working_space[i - 1] = val;
1848 working_space[fSize - i] = -val;
1849 }
1850 Fourier(working_space, fSize, 0, kTransformForward, 0);
1851 for (i = 0; i < fSize / 2; i++) {
1852 a = pi * (Double_t) i / (Double_t) fSize;
1853 a = TMath::Sin(a);
1854 b = working_space[i];
1855 if (a != 0)
1856 a = b / a;
1857 working_space[i - 1] = a;
1858 working_space[i + fSize] = 0;
1859 }
1860 working_space[fSize / 2 - 1] =
1861 working_space[fSize / 2] / TMath::Sqrt(2.0);
1862 fSize = fSize / 2;
1863 break;
1864 case kTransformFourier:
1865 for (i = 0; i < fSize; i++) {
1866 working_space[i] = source[i];
1867 }
1868 Fourier(working_space, fSize, 0, kTransformForward, 0);
1869 break;
1870 case kTransformHartley:
1871 for (i = 0; i < fSize; i++) {
1872 working_space[i] = source[i];
1873 }
1874 Fourier(working_space, fSize, 1, kTransformForward, 0);
1875 break;
1876 case kTransformFourierWalsh:
1877 case kTransformFourierHaar:
1878 case kTransformWalshHaar:
1879 case kTransformCosWalsh:
1880 case kTransformCosHaar:
1881 case kTransformSinWalsh:
1882 case kTransformSinHaar:
1883 for (i = 0; i < fSize; i++) {
1884 val = source[i];
1885 if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1886 j = (Int_t) TMath::Power(2, fDegree) / 2;
1887 k = i / j;
1888 k = 2 * k * j;
1889 working_space[k + i % j] = val;
1890 working_space[k + 2 * j - 1 - i % j] = val;
1891 }
1892
1893 else if (fTransformType == kTransformSinWalsh
1894 || fTransformType == kTransformSinHaar) {
1895 j = (Int_t) TMath::Power(2, fDegree) / 2;
1896 k = i / j;
1897 k = 2 * k * j;
1898 working_space[k + i % j] = val;
1899 working_space[k + 2 * j - 1 - i % j] = -val;
1900 }
1901
1902 else
1903 working_space[i] = val;
1904 }
1905 if (fTransformType == kTransformFourierWalsh
1906 || fTransformType == kTransformFourierHaar
1907 || fTransformType == kTransformWalshHaar) {
1908 for (i = 0; i < j; i++)
1909 BitReverseHaar(working_space, fSize, k, i * k);
1910 GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
1911 }
1912
1913 else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1914 m = (Int_t) TMath::Power(2, fDegree);
1915 l = 2 * fSize / m;
1916 for (i = 0; i < l; i++)
1917 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1918 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1919 for (i = 0; i < fSize; i++) {
1920 k = i / j;
1921 k = 2 * k * j;
1922 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1923 a = TMath::Cos(a);
1924 b = working_space[k + i % j];
1925 if (i % j == 0)
1926 a = b / TMath::Sqrt(2.0);
1927
1928 else
1929 a = b / a;
1930 working_space[i] = a;
1931 working_space[i + 2 * fSize] = 0;
1932 }
1933 }
1934
1935 else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
1936 m = (Int_t) TMath::Power(2, fDegree);
1937 l = 2 * fSize / m;
1938 for (i = 0; i < l; i++)
1939 BitReverseHaar(working_space, 2 * fSize, m, i * m);
1940 GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1941 for (i = 0; i < fSize; i++) {
1942 k = i / j;
1943 k = 2 * k * j;
1944 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1945 a = TMath::Cos(a);
1946 b = working_space[j + k + i % j];
1947 if (i % j == 0)
1948 a = b / TMath::Sqrt(2.0);
1949
1950 else
1951 a = b / a;
1952 working_space[j + k / 2 - i % j - 1] = a;
1953 working_space[i + 2 * fSize] = 0;
1954 }
1955 }
1956 if (fTransformType > kTransformWalshHaar)
1957 k = (Int_t) TMath::Power(2, fDegree - 1);
1958
1959 else
1960 k = (Int_t) TMath::Power(2, fDegree);
1961 j = fSize / k;
1962 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1963 working_space[fSize + i] = working_space[l + i / j];
1964 working_space[fSize + i + 2 * fSize] =
1965 working_space[l + i / j + 2 * fSize];
1966 }
1967 for (i = 0; i < fSize; i++) {
1968 working_space[i] = working_space[fSize + i];
1969 working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
1970 }
1971 break;
1972 }
1973 if (!working_space) return;
1974 for (i = 0, old_area = 0; i < fSize; i++) {
1975 old_area += working_space[i];
1976 }
1977 for (i = 0, new_area = 0; i < fSize; i++) {
1978 if (i >= fXmin && i <= fXmax)
1979 working_space[i] *= fEnhanceCoeff;
1980 new_area += working_space[i];
1981 }
1982 if (new_area != 0) {
1983 a = old_area / new_area;
1984 for (i = 0; i < fSize; i++) {
1985 working_space[i] *= a;
1986 }
1987 }
1988 if (fTransformType == kTransformFourier) {
1989 for (i = 0, old_area = 0; i < fSize; i++) {
1990 old_area += working_space[i + fSize];
1991 }
1992 for (i = 0, new_area = 0; i < fSize; i++) {
1993 if (i >= fXmin && i <= fXmax)
1994 working_space[i + fSize] *= fEnhanceCoeff;
1995 new_area += working_space[i + fSize];
1996 }
1997 if (new_area != 0) {
1998 a = old_area / new_area;
1999 for (i = 0; i < fSize; i++) {
2000 working_space[i + fSize] *= a;
2001 }
2002 }
2003 }
2004
2005 else if (fTransformType == kTransformFourierWalsh
2006 || fTransformType == kTransformFourierHaar) {
2007 for (i = 0, old_area = 0; i < fSize; i++) {
2008 old_area += working_space[i + 2 * fSize];
2009 }
2010 for (i = 0, new_area = 0; i < fSize; i++) {
2011 if (i >= fXmin && i <= fXmax)
2012 working_space[i + 2 * fSize] *= fEnhanceCoeff;
2013 new_area += working_space[i + 2 * fSize];
2014 }
2015 if (new_area != 0) {
2016 a = old_area / new_area;
2017 for (i = 0; i < fSize; i++) {
2018 working_space[i + 2 * fSize] *= a;
2019 }
2020 }
2021 }
2022 switch (fTransformType) {
2023 case kTransformHaar:
2024 Haar(working_space, fSize, kTransformInverse);
2025 for (i = 0; i < fSize; i++) {
2026 destVector[i] = working_space[i];
2027 }
2028 break;
2029 case kTransformWalsh:
2030 BitReverse(working_space, fSize);
2031 Walsh(working_space, fSize);
2032 for (i = 0; i < fSize; i++) {
2033 destVector[i] = working_space[i];
2034 }
2035 break;
2036 case kTransformCos:
2037 fSize = 2 * fSize;
2038 working_space[0] = working_space[0] * TMath::Sqrt(2.0);
2039 for (i = 0; i < fSize / 2; i++) {
2040 a = pi * (Double_t) i / (Double_t) fSize;
2041 b = TMath::Sin(a);
2042 a = TMath::Cos(a);
2043 working_space[i + fSize] = (Double_t) working_space[i] * b;
2044 working_space[i] = (Double_t) working_space[i] * a;
2045 } for (i = 2; i <= (fSize / 2); i++) {
2046 working_space[fSize - i + 1] = working_space[i - 1];
2047 working_space[fSize - i + 1 + fSize] =
2048 -working_space[i - 1 + fSize];
2049 }
2050 working_space[fSize / 2] = 0;
2051 working_space[fSize / 2 + fSize] = 0;
2052 Fourier(working_space, fSize, 0, kTransformInverse, 1);
2053 for (i = 0; i < fSize / 2; i++) {
2054 destVector[i] = working_space[i];
2055 }
2056 break;
2057 case kTransformSin:
2058 fSize = 2 * fSize;
2059 working_space[fSize / 2] =
2060 working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
2061 for (i = fSize / 2 - 1; i > 0; i--) {
2062 a = pi * (Double_t) i / (Double_t) fSize;
2063 working_space[i + fSize] =
2064 -(Double_t) working_space[i - 1] * TMath::Cos(a);
2065 working_space[i] = (Double_t) working_space[i - 1] * TMath::Sin(a);
2066 } for (i = 2; i <= (fSize / 2); i++) {
2067 working_space[fSize - i + 1] = working_space[i - 1];
2068 working_space[fSize - i + 1 + fSize] =
2069 -working_space[i - 1 + fSize];
2070 }
2071 working_space[0] = 0;
2072 working_space[fSize] = 0;
2073 working_space[fSize / 2 + fSize] = 0;
2074 Fourier(working_space, fSize, 0, kTransformInverse, 0);
2075 for (i = 0; i < fSize / 2; i++) {
2076 destVector[i] = working_space[i];
2077 }
2078 break;
2079 case kTransformFourier:
2080 Fourier(working_space, fSize, 0, kTransformInverse, 0);
2081 for (i = 0; i < fSize; i++) {
2082 destVector[i] = working_space[i];
2083 }
2084 break;
2085 case kTransformHartley:
2086 Fourier(working_space, fSize, 1, kTransformInverse, 0);
2087 for (i = 0; i < fSize; i++) {
2088 destVector[i] = working_space[i];
2089 }
2090 break;
2091 case kTransformFourierWalsh:
2092 case kTransformFourierHaar:
2093 case kTransformWalshHaar:
2094 case kTransformCosWalsh:
2095 case kTransformCosHaar:
2096 case kTransformSinWalsh:
2097 case kTransformSinHaar:
2098 if (fTransformType > kTransformWalshHaar)
2099 k = (Int_t) TMath::Power(2, fDegree - 1);
2100
2101 else
2102 k = (Int_t) TMath::Power(2, fDegree);
2103 j = fSize / k;
2104 for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
2105 working_space[fSize + l + i / j] = working_space[i];
2106 working_space[fSize + l + i / j + 2 * fSize] =
2107 working_space[i + 2 * fSize];
2108 }
2109 for (i = 0; i < fSize; i++) {
2110 working_space[i] = working_space[fSize + i];
2111 working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
2112 }
2113 if (fTransformType == kTransformFourierWalsh
2114 || fTransformType == kTransformFourierHaar
2115 || fTransformType == kTransformWalshHaar) {
2116 GeneralInv(working_space, fSize, fDegree, fTransformType);
2117 for (i = 0; i < j; i++)
2118 BitReverseHaar(working_space, fSize, k, i * k);
2119 }
2120
2121 else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
2122 j = (Int_t) TMath::Power(2, fDegree) / 2;
2123 m = (Int_t) TMath::Power(2, fDegree);
2124 l = 2 * fSize / m;
2125 for (i = 0; i < fSize; i++) {
2126 k = i / j;
2127 k = 2 * k * j;
2128 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2129 if (i % j == 0) {
2130 working_space[2 * fSize + k + i % j] =
2131 working_space[i] * TMath::Sqrt(2.0);
2132 working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
2133 }
2134
2135 else {
2136 b = TMath::Sin(a);
2137 a = TMath::Cos(a);
2138 working_space[4 * fSize + 2 * fSize + k + i % j] =
2139 -(Double_t) working_space[i] * b;
2140 working_space[2 * fSize + k + i % j] =
2141 (Double_t) working_space[i] * a;
2142 } } for (i = 0; i < fSize; i++) {
2143 k = i / j;
2144 k = 2 * k * j;
2145 if (i % j == 0) {
2146 working_space[2 * fSize + k + j] = 0;
2147 working_space[4 * fSize + 2 * fSize + k + j] = 0;
2148 }
2149
2150 else {
2151 working_space[2 * fSize + k + 2 * j - i % j] =
2152 working_space[2 * fSize + k + i % j];
2153 working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
2154 -working_space[4 * fSize + 2 * fSize + k + i % j];
2155 }
2156 }
2157 for (i = 0; i < 2 * fSize; i++) {
2158 working_space[i] = working_space[2 * fSize + i];
2159 working_space[4 * fSize + i] =
2160 working_space[4 * fSize + 2 * fSize + i];
2161 }
2162 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2163 m = (Int_t) TMath::Power(2, fDegree);
2164 l = 2 * fSize / m;
2165 for (i = 0; i < l; i++)
2166 BitReverseHaar(working_space, 2 * fSize, m, i * m);
2167 }
2168
2169 else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
2170 j = (Int_t) TMath::Power(2, fDegree) / 2;
2171 m = (Int_t) TMath::Power(2, fDegree);
2172 l = 2 * fSize / m;
2173 for (i = 0; i < fSize; i++) {
2174 k = i / j;
2175 k = 2 * k * j;
2176 a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2177 if (i % j == 0) {
2178 working_space[2 * fSize + k + j + i % j] =
2179 working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
2180 working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
2181 }
2182
2183 else {
2184 b = TMath::Sin(a);
2185 a = TMath::Cos(a);
2186 working_space[4 * fSize + 2 * fSize + k + j + i % j] =
2187 -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
2188 working_space[2 * fSize + k + j + i % j] =
2189 (Double_t) working_space[j + k / 2 - i % j - 1] * a;
2190 } } for (i = 0; i < fSize; i++) {
2191 k = i / j;
2192 k = 2 * k * j;
2193 if (i % j == 0) {
2194 working_space[2 * fSize + k] = 0;
2195 working_space[4 * fSize + 2 * fSize + k] = 0;
2196 }
2197
2198 else {
2199 working_space[2 * fSize + k + i % j] =
2200 working_space[2 * fSize + k + 2 * j - i % j];
2201 working_space[4 * fSize + 2 * fSize + k + i % j] =
2202 -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
2203 }
2204 }
2205 for (i = 0; i < 2 * fSize; i++) {
2206 working_space[i] = working_space[2 * fSize + i];
2207 working_space[4 * fSize + i] =
2208 working_space[4 * fSize + 2 * fSize + i];
2209 }
2210 GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2211 for (i = 0; i < l; i++)
2212 BitReverseHaar(working_space, 2 * fSize, m, i * m);
2213 }
2214 for (i = 0; i < fSize; i++) {
2215 if (fTransformType >= kTransformCosWalsh) {
2216 k = i / j;
2217 k = 2 * k * j;
2218 val = working_space[k + i % j];
2219 }
2220
2221 else
2222 val = working_space[i];
2223 destVector[i] = val;
2224 }
2225 break;
2226 }
2227 delete[]working_space;
2228 return;
2229}
2230
2231////////////////////////////////////////////////////////////////////////////////
2232/// This function sets the following parameters for transform:
2233/// - transType - type of transform (Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar)
2234/// - degree - degree of mixed transform, applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar transforms
2235
2236void TSpectrumTransform::SetTransformType(Int_t transType, Int_t degree)
2237{
2238 Int_t j, n;
2239 j = 0;
2240 n = 1;
2241 for (; n < fSize;) {
2242 j += 1;
2243 n = n * 2;
2244 }
2245 if (transType < kTransformHaar || transType > kTransformSinHaar){
2246 Error ("TSpectrumTransform","Invalid type of transform");
2247 return;
2248 }
2249 if (transType >= kTransformFourierWalsh && transType <= kTransformSinHaar) {
2250 if (degree > j || degree < 1){
2251 Error ("TSpectrumTransform","Invalid degree of mixed transform");
2252 return;
2253 }
2254 }
2255 fTransformType=transType;
2256 fDegree=degree;
2257}
2258
2259////////////////////////////////////////////////////////////////////////////////
2260/// This function sets the filtering or enhancement region:
2261/// - xmin, xmax
2262
2263void TSpectrumTransform::SetRegion(Int_t xmin, Int_t xmax)
2264{
2265 if(xmin<0 || xmax < xmin || xmax >= fSize){
2266 Error("TSpectrumTransform", "Wrong range");
2267 return;
2268 }
2269 fXmin = xmin;
2270 fXmax = xmax;
2271}
2272
2273////////////////////////////////////////////////////////////////////////////////
2274/// This function sets the direction of the transform:
2275/// - direction (forward or inverse)
2276
2277void TSpectrumTransform::SetDirection(Int_t direction)
2278{
2279 if(direction != kTransformForward && direction != kTransformInverse){
2280 Error("TSpectrumTransform", "Wrong direction");
2281 return;
2282 }
2283 fDirection = direction;
2284}
2285
2286////////////////////////////////////////////////////////////////////////////////
2287/// This function sets the filter coefficient:
2288/// - filterCoeff - after the transform the filtered region (xmin, xmax) is replaced by this coefficient. Applies only for filtereng operation.
2289
2290void TSpectrumTransform::SetFilterCoeff(Double_t filterCoeff)
2291{
2292 fFilterCoeff = filterCoeff;
2293}
2294
2295////////////////////////////////////////////////////////////////////////////////
2296/// This function sets the enhancement coefficient:
2297/// - enhanceCoeff - after the transform the enhanced region (xmin, xmax) is multiplied by this coefficient. Applies only for enhancement operation.
2298
2299void TSpectrumTransform::SetEnhanceCoeff(Double_t enhanceCoeff)
2300{
2301 fEnhanceCoeff = enhanceCoeff;
2302}
d
#define d(i)
Definition: RSha256.hxx:102
b
#define b(i)
Definition: RSha256.hxx:100
c
#define c(i)
Definition: RSha256.hxx:101
Int_t
int Int_t
Definition: RtypesCore.h:43
Double_t
double Double_t
Definition: RtypesCore.h:57
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:361
type
int type
Definition: TGX11.cxx:120
xmin
float xmin
Definition: THbookFile.cxx:93
xmax
float xmax
Definition: THbookFile.cxx:93
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
TSpectrumTransform
Advanced 1-dimensional orthogonal transform functions.
Definition: TSpectrumTransform.h:18
TSpectrumTransform::Haar
void Haar(Double_t *working_space, Int_t num, Int_t direction)
This function calculates Haar transform of a part of data Function parameters:
Definition: TSpectrumTransform.cxx:95
TSpectrumTransform::~TSpectrumTransform
virtual ~TSpectrumTransform()
Destructor.
Definition: TSpectrumTransform.cxx:84
TSpectrumTransform::fTransformType
Int_t fTransformType
type of transformation (Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh,...
Definition: TSpectrumTransform.h:22
TSpectrumTransform::Transform
void Transform(const Double_t *source, Double_t *destVector)
This function transforms the source spectrum.
Definition: TSpectrumTransform.cxx:740
TSpectrumTransform::BitReverseHaar
void BitReverseHaar(Double_t *working_space, Int_t shift, Int_t num, Int_t start)
This function carries out bit-reverse reordering for Haar transform Function parameters:
Definition: TSpectrumTransform.cxx:385
TSpectrumTransform::SetRegion
void SetRegion(Int_t xmin, Int_t xmax)
This function sets the filtering or enhancement region:
Definition: TSpectrumTransform.cxx:2263
TSpectrumTransform::SetDirection
void SetDirection(Int_t direction)
This function sets the direction of the transform:
Definition: TSpectrumTransform.cxx:2277
TSpectrumTransform::Fourier
void Fourier(Double_t *working_space, Int_t num, Int_t hartley, Int_t direction, Int_t zt_clear)
This function calculates Fourier based transform of a part of data Function parameters:
Definition: TSpectrumTransform.cxx:264
TSpectrumTransform::Enhance
void Enhance(const Double_t *source, Double_t *destVector)
This function transforms the source spectrum.
Definition: TSpectrumTransform.cxx:1778
TSpectrumTransform::BitReverse
void BitReverse(Double_t *working_space, Int_t num)
This function carries out bit-reverse reordering of data Function parameters:
Definition: TSpectrumTransform.cxx:226
TSpectrumTransform::fSize
Int_t fSize
length of transformed data
Definition: TSpectrumTransform.h:21
TSpectrumTransform::fEnhanceCoeff
Double_t fEnhanceCoeff
multiplication coefficient applied in enhanced region;
Definition: TSpectrumTransform.h:28
TSpectrumTransform::SetEnhanceCoeff
void SetEnhanceCoeff(Double_t enhanceCoeff)
This function sets the enhancement coefficient:
Definition: TSpectrumTransform.cxx:2299
TSpectrumTransform::GeneralExe
Int_t GeneralExe(Double_t *working_space, Int_t zt_clear, Int_t num, Int_t degree, Int_t type)
This function calculates generalized (mixed) transforms of different degrees Function parameters:
Definition: TSpectrumTransform.cxx:430
TSpectrumTransform::Walsh
void Walsh(Double_t *working_space, Int_t num)
This function calculates Walsh transform of a part of data Function parameters:
Definition: TSpectrumTransform.cxx:173
TSpectrumTransform::GeneralInv
Int_t GeneralInv(Double_t *working_space, Int_t num, Int_t degree, Int_t type)
This function calculates inverse generalized (mixed) transforms Function parameters:
Definition: TSpectrumTransform.cxx:546
TSpectrumTransform::SetFilterCoeff
void SetFilterCoeff(Double_t filterCoeff)
This function sets the filter coefficient:
Definition: TSpectrumTransform.cxx:2290
TSpectrumTransform::fXmin
Int_t fXmin
first channel of filtered or enhanced region
Definition: TSpectrumTransform.h:25
TSpectrumTransform::fXmax
Int_t fXmax
last channel of filtered or enhanced region
Definition: TSpectrumTransform.h:26
TSpectrumTransform::fDirection
Int_t fDirection
forward or inverse transform
Definition: TSpectrumTransform.h:24
TSpectrumTransform::FilterZonal
void FilterZonal(const Double_t *source, Double_t *destVector)
This function transforms the source spectrum.
Definition: TSpectrumTransform.cxx:1268
TSpectrumTransform::fDegree
Int_t fDegree
degree of mixed transform, applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar,...
Definition: TSpectrumTransform.h:23
TSpectrumTransform::TSpectrumTransform
TSpectrumTransform()
default constructor
Definition: TSpectrumTransform.cxx:37
TSpectrumTransform::fFilterCoeff
Double_t fFilterCoeff
value set in the filtered region
Definition: TSpectrumTransform.h:27
TSpectrumTransform::SetTransformType
void SetTransformType(Int_t transType, Int_t degree)
This function sets the following parameters for transform:
Definition: TSpectrumTransform.cxx:2236
n
const Int_t n
Definition: legend1.C:16
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:67
TGeant4Unit::degree
static constexpr double degree
Definition: TGeant4SystemOfUnits.h:135
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
TMath::Cos
Double_t Cos(Double_t)
Definition: TMath.h:631
TMath::Sin
Double_t Sin(Double_t)
Definition: TMath.h:627
m
auto * m
Definition: textangle.C:8
l
auto * l
Definition: textangle.C:4
a
auto * a
Definition: textangle.C:12