Logo ROOT  

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