Logo ROOT  

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