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