Logo ROOT  

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