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