Logo ROOT  

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