Logo ROOT  

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