ROOT  6.06/08
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 
743 void TSpectrumTransform::Transform(const Double_t *source, Double_t *destVector)
744 {
745 /*
746 <div class=Section1>
747 
748 <p class=MsoNormal><b><span style='font-size:14.0pt'>Transform methods</span></b></p>
749 
750 <p class=MsoNormal style='text-align:justify'><i>&nbsp;</i></p>
751 
752 <p class=MsoNormal style='text-align:justify'><i>Goal: to analyze experimental
753 data using orthogonal transforms</i></p>
754 
755 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
756 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
757 </span>orthogonal transforms can be successfully used for the processing of
758 nuclear spectra (not only) </p>
759 
760 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
761 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
762 </span>they can be used to remove high frequency noise, to increase
763 signal-to-background ratio as well as to enhance low intensity components [1],
764 to carry out e.g. Fourier analysis etc. </p>
765 
766 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
767 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
768 </span>we have implemented the function for the calculation of the commonly
769 used orthogonal transforms as well as functions for the filtration and
770 enhancement of experimental data</p>
771 
772 <p class=MsoNormal><i>&nbsp;</i></p>
773 
774 <p class=MsoNormal><i>Function:</i></p>
775 
776 <p class=MsoNormal><b>void TSpectrumTransform::Transform(const <a
777 href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> *fSource,
778 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a>
779 *fDest)</b></p>
780 
781 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
782 
783 <p class=MsoNormal style='text-align:justify'>This function transforms the
784 source spectrum according to the given input parameters. Transformed data are
785 written into dest spectrum. Before the Transform function is called the class
786 must be created by constructor and the type of the transform as well as some
787 other parameters must be set using a set of setter functions.</p>
788 
789 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
790 
791 <p class=MsoNormal><i><span style='color:red'>Member variables of
792 TSpectrumTransform class:</span></i></p>
793 
794 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'> <b>fSource</b>-pointer
795 to the vector of source spectrum. Its length should be equal to the “fSize”
796 parameter except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR transforms. These
797 need “2*fSize” length to supply real and imaginary coefficients.                   </p>
798 
799 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'><b>fDest</b>-pointer
800 to the vector of destination spectrum. Its length should be equal to the
801 “fSize” parameter except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR transforms.
802 These need “2*fSize” length to store real and imaginary coefficients. </p>
803 
804 <p class=MsoNormal style='text-align:justify'>        <b>fSize</b>-basic length
805 of the source and dest spectrum. <span style='color:fuchsia'>It should be power
806 of 2.</span></p>
807 
808 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify;text-indent:
809 -2.85pt'><b>fType</b>-type of transform</p>
810 
811 <p class=MsoNormal style='text-align:justify'>            Classic transforms:</p>
812 
813 <p class=MsoNormal style='text-align:justify'>                        kTransformHaar
814 </p>
815 
816 <p class=MsoNormal style='text-align:justify'>                        kTransformWalsh
817 </p>
818 
819 <p class=MsoNormal style='text-align:justify'>                        kTransformCos
820 </p>
821 
822 <p class=MsoNormal style='text-align:justify'>                        kTransformSin
823 </p>
824 
825 <p class=MsoNormal style='text-align:justify'>                        kTransformFourier
826 </p>
827 
828 <p class=MsoNormal style='text-align:justify'>                        kTransformHartey
829 </p>
830 
831 <p class=MsoNormal style='text-align:justify'>            Mixed transforms:</p>
832 
833 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierWalsh
834 </p>
835 
836 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierHaar
837 </p>
838 
839 <p class=MsoNormal style='text-align:justify'>                        kTransformWalshHaar
840 </p>
841 
842 <p class=MsoNormal style='text-align:justify'>                        kTransformCosWalsh
843 </p>
844 
845 <p class=MsoNormal style='text-align:justify'>                        kTransformCosHaar
846 </p>
847 
848 <p class=MsoNormal style='text-align:justify'>                        kTransformSinWalsh
849 </p>
850 
851 <p class=MsoNormal style='text-align:justify'>                        kTransformSinHaar
852 </p>
853 
854 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDirection</b>-direction-transform
855 direction (forward, inverse)</p>
856 
857 <p class=MsoNormal style='text-align:justify'>                        kTransformForward
858 </p>
859 
860 <p class=MsoNormal style='text-align:justify'>                        kTransformInverse
861 </p>
862 
863 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDegree</b>-applies
864 only for mixed transforms [2], [3], [4]. </p>
865 
866 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'>                
867 <span style='color:fuchsia'> Allowed range  <sub><img border=0 width=100
868 height=27 src="gif/spectrumtransform_transform_image001.gif"></sub>. </span></p>
869 
870 <p class=MsoNormal style='text-align:justify'><b><i>References:</i></b></p>
871 
872 <p class=MsoNormal style='text-align:justify'>[1] C.V. Hampton, B. Lian, Wm. C.
873 McHarris: Fast-Fourier-transform spectral enhancement techniques for gamma-ray
874 spectroscopy. NIM A353 (1994) 280-284. </p>
875 
876 <p class=MsoNormal style='text-align:justify'>[2] Morhá&#269; M., Matoušek V.,
877 New adaptive Cosine-Walsh  transform and its application to nuclear data
878 compression, IEEE Transactions on Signal Processing 48 (2000) 2693.  </p>
879 
880 <p class=MsoNormal style='text-align:justify'>[3] Morhá&#269; M., Matoušek V.,
881 Data compression using new fast adaptive Cosine-Haar transforms, Digital Signal
882 Processing 8 (1998) 63. </p>
883 
884 <p class=MsoNormal style='text-align:justify'>[4] Morhá&#269; M., Matoušek V.:
885 Multidimensional nuclear data compression using fast adaptive Walsh-Haar
886 transform. Acta Physica Slovaca 51 (2001) 307. </p>
887 
888 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
889 
890 <p class=MsoNormal style='text-align:justify'><i>Example  – script Transform.c:</i></p>
891 
892 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'><img
893 width=600 height=324 src="gif/spectrumtransform_transform_image002.jpg"></span></p>
894 
895 <p class=MsoNormal><b>Fig. 1 Original gamma-ray spectrum</b></p>
896 
897 <p class=MsoNormal><b><span style='font-size:14.0pt'><img border=0 width=601
898 height=402 src="gif/spectrumtransform_transform_image003.jpg"></span></b></p>
899 
900 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'>&nbsp;</span></p>
901 
902 <p class=MsoNormal><b>Fig. 2 Transformed spectrum from Fig. 1 using Cosine
903 transform</b></p>
904 
905 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>&nbsp;</span></b></p>
906 
907 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
908 
909 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
910 Transform function (class TSpectrumTransform).</span></p>
911 
912 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
913 do</span></p>
914 
915 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Transform.C</span></p>
916 
917 <p class=MsoNormal><span style='font-size:10.0pt'>   </span></p>
918 
919 <p class=MsoNormal><span style='font-size:10.0pt'>#include &lt;TSpectrum&gt;</span></p>
920 
921 <p class=MsoNormal><span style='font-size:10.0pt'>#include
922 &lt;TSpectrumTransform&gt;</span></p>
923 
924 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
925 
926 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Transform() {</span></p>
927 
928 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
929 
930 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
931 4096;</span></p>
932 
933 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
934 0;</span></p>
935 
936 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
937 (Double_t)nbins;</span></p>
938 
939 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
940 style='font-size:10.0pt'>Double_t * source = new Double_t[nbins];</span></p>
941 
942 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t * dest = new
943 Double_t[nbins];   </span></p>
944 
945 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new TH1F(&quot;h&quot;,&quot;Transformed
946 spectrum using Cosine transform&quot;,nbins,xmin,xmax);</span></p>
947 
948 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
949 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
950 
951 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
952 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
953 
954 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
955 i++) source[i]=h-&gt;GetBinContent(i + 1);         </span></p>
956 
957 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 =
958 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
959 
960 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
961 Transform1 = new
962 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
963 
964 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
965 TSpectrum();</span></p>
966 
967 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
968 new TSpectrumTransform(4096);</span></p>
969 
970 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
971 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
972 
973 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
974 t-&gt;SetDirection(t-&gt;kTransformForward);</span></p>
975 
976 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
977 style='font-size:10.0pt'>t-&gt;Transform(source,dest);</span></p>
978 
979 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
980 i++) h-&gt;SetBinContent(i + 1,dest[i]);   </span></p>
981 
982 <p class=MsoNormal><span style='font-size:10.0pt'>  
983 h-&gt;SetLineColor(kRed);      </span></p>
984 
985 <p class=MsoNormal><span style='font-size:10.0pt'>   h-&gt;Draw(&quot;L&quot;);</span></p>
986 
987 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
988 
989 </div>
990 
991  */
992  int i, j=0, k = 1, m, l;
993  Double_t val;
994  Double_t a, b, pi = 3.14159265358979323846;
995  Double_t *working_space = 0;
996  if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
997  if (fTransformType >= kTransformCosWalsh)
998  fDegree += 1;
999  k = (Int_t) TMath::Power(2, fDegree);
1000  j = fSize / k;
1001  }
1002  switch (fTransformType) {
1003  case kTransformHaar:
1004  case kTransformWalsh:
1005  working_space = new Double_t[2 * fSize];
1006  break;
1007  case kTransformCos:
1008  case kTransformSin:
1009  case kTransformFourier:
1010  case kTransformHartley:
1011  case kTransformFourierWalsh:
1012  case kTransformFourierHaar:
1013  case kTransformWalshHaar:
1014  working_space = new Double_t[4 * fSize];
1015  break;
1016  case kTransformCosWalsh:
1017  case kTransformCosHaar:
1018  case kTransformSinWalsh:
1019  case kTransformSinHaar:
1020  working_space = new Double_t[8 * fSize];
1021  break;
1022  }
1023  if (fDirection == kTransformForward) {
1024  switch (fTransformType) {
1025  case kTransformHaar:
1026  for (i = 0; i < fSize; i++) {
1027  working_space[i] = source[i];
1028  }
1029  Haar(working_space, fSize, fDirection);
1030  for (i = 0; i < fSize; i++) {
1031  destVector[i] = working_space[i];
1032  }
1033  break;
1034  case kTransformWalsh:
1035  for (i = 0; i < fSize; i++) {
1036  working_space[i] = source[i];
1037  }
1038  Walsh(working_space, fSize);
1039  BitReverse(working_space, fSize);
1040  for (i = 0; i < fSize; i++) {
1041  destVector[i] = working_space[i];
1042  }
1043  break;
1044  case kTransformCos:
1045  fSize = 2 * fSize;
1046  for (i = 1; i <= (fSize / 2); i++) {
1047  val = source[i - 1];
1048  working_space[i - 1] = val;
1049  working_space[fSize - i] = val;
1050  }
1051  Fourier(working_space, fSize, 0, kTransformForward, 0);
1052  for (i = 0; i < fSize / 2; i++) {
1053  a = pi * (Double_t) i / (Double_t) fSize;
1054  a = TMath::Cos(a);
1055  b = working_space[i];
1056  a = b / a;
1057  working_space[i] = a;
1058  working_space[i + fSize] = 0;
1059  } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
1060  for (i = 0; i < fSize / 2; i++) {
1061  destVector[i] = working_space[i];
1062  }
1063  break;
1064  case kTransformSin:
1065  fSize = 2 * fSize;
1066  for (i = 1; i <= (fSize / 2); i++) {
1067  val = source[i - 1];
1068  working_space[i - 1] = val;
1069  working_space[fSize - i] = -val;
1070  }
1071  Fourier(working_space, fSize, 0, kTransformForward, 0);
1072  for (i = 0; i < fSize / 2; i++) {
1073  a = pi * (Double_t) i / (Double_t) fSize;
1074  a = TMath::Sin(a);
1075  b = working_space[i];
1076  if (a != 0)
1077  a = b / a;
1078  working_space[i - 1] = a;
1079  working_space[i + fSize] = 0;
1080  }
1081  working_space[fSize / 2 - 1] =
1082  working_space[fSize / 2] / TMath::Sqrt(2.0);
1083  for (i = 0; i < fSize / 2; i++) {
1084  destVector[i] = working_space[i];
1085  }
1086  break;
1087  case kTransformFourier:
1088  for (i = 0; i < fSize; i++) {
1089  working_space[i] = source[i];
1090  }
1091  Fourier(working_space, fSize, 0, kTransformForward, 0);
1092  for (i = 0; i < 2 * fSize; i++) {
1093  destVector[i] = working_space[i];
1094  }
1095  break;
1096  case kTransformHartley:
1097  for (i = 0; i < fSize; i++) {
1098  working_space[i] = source[i];
1099  }
1100  Fourier(working_space, fSize, 1, kTransformForward, 0);
1101  for (i = 0; i < fSize; i++) {
1102  destVector[i] = working_space[i];
1103  }
1104  break;
1105  case kTransformFourierWalsh:
1106  case kTransformFourierHaar:
1107  case kTransformWalshHaar:
1108  case kTransformCosWalsh:
1109  case kTransformCosHaar:
1110  case kTransformSinWalsh:
1111  case kTransformSinHaar:
1112  for (i = 0; i < fSize; i++) {
1113  val = source[i];
1114  if (fTransformType == kTransformCosWalsh
1115  || fTransformType == kTransformCosHaar) {
1116  j = (Int_t) TMath::Power(2, fDegree) / 2;
1117  k = i / j;
1118  k = 2 * k * j;
1119  working_space[k + i % j] = val;
1120  working_space[k + 2 * j - 1 - i % j] = val;
1121  }
1122 
1123  else if (fTransformType == kTransformSinWalsh
1124  || fTransformType == kTransformSinHaar) {
1125  j = (Int_t) TMath::Power(2, fDegree) / 2;
1126  k = i / j;
1127  k = 2 * k * j;
1128  working_space[k + i % j] = val;
1129  working_space[k + 2 * j - 1 - i % j] = -val;
1130  }
1131 
1132  else
1133  working_space[i] = val;
1134  }
1135  if (fTransformType == kTransformFourierWalsh
1136  || fTransformType == kTransformFourierHaar
1137  || fTransformType == kTransformWalshHaar) {
1138  for (i = 0; i < j; i++)
1139  BitReverseHaar(working_space, fSize, k, i * k);
1140  GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
1141  }
1142 
1143  else if (fTransformType == kTransformCosWalsh
1144  || fTransformType == kTransformCosHaar) {
1145  m = (Int_t) TMath::Power(2, fDegree);
1146  l = 2 * fSize / m;
1147  for (i = 0; i < l; i++)
1148  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1149  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1150  for (i = 0; i < fSize; i++) {
1151  k = i / j;
1152  k = 2 * k * j;
1153  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1154  a = TMath::Cos(a);
1155  b = working_space[k + i % j];
1156  if (i % j == 0)
1157  a = b / TMath::Sqrt(2.0);
1158 
1159  else
1160  a = b / a;
1161  working_space[i] = a;
1162  working_space[i + 2 * fSize] = 0;
1163  }
1164  }
1165 
1166  else if (fTransformType == kTransformSinWalsh
1167  || fTransformType == kTransformSinHaar) {
1168  m = (Int_t) TMath::Power(2, fDegree);
1169  l = 2 * fSize / m;
1170  for (i = 0; i < l; i++)
1171  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1172  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1173  for (i = 0; i < fSize; i++) {
1174  k = i / j;
1175  k = 2 * k * j;
1176  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1177  a = TMath::Cos(a);
1178  b = working_space[j + k + i % j];
1179  if (i % j == 0)
1180  a = b / TMath::Sqrt(2.0);
1181 
1182  else
1183  a = b / a;
1184  working_space[j + k / 2 - i % j - 1] = a;
1185  working_space[i + 2 * fSize] = 0;
1186  }
1187  }
1188  if (fTransformType > kTransformWalshHaar)
1189  k = (Int_t) TMath::Power(2, fDegree - 1);
1190 
1191  else
1192  k = (Int_t) TMath::Power(2, fDegree);
1193  j = fSize / k;
1194  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1195  working_space[fSize + i] = working_space[l + i / j];
1196  working_space[fSize + i + 2 * fSize] =
1197  working_space[l + i / j + 2 * fSize];
1198  }
1199  for (i = 0; i < fSize; i++) {
1200  working_space[i] = working_space[fSize + i];
1201  working_space[i + 2 * fSize] =
1202  working_space[fSize + i + 2 * fSize];
1203  }
1204  for (i = 0; i < fSize; i++) {
1205  destVector[i] = working_space[i];
1206  }
1207  if (fTransformType == kTransformFourierWalsh
1208  || fTransformType == kTransformFourierHaar) {
1209  for (i = 0; i < fSize; i++) {
1210  destVector[fSize + i] = working_space[i + 2 * fSize];
1211  }
1212  }
1213  break;
1214  }
1215  }
1216 
1217  else if (fDirection == kTransformInverse) {
1218  switch (fTransformType) {
1219  case kTransformHaar:
1220  for (i = 0; i < fSize; i++) {
1221  working_space[i] = source[i];
1222  }
1223  Haar(working_space, fSize, fDirection);
1224  for (i = 0; i < fSize; i++) {
1225  destVector[i] = working_space[i];
1226  }
1227  break;
1228  case kTransformWalsh:
1229  for (i = 0; i < fSize; i++) {
1230  working_space[i] = source[i];
1231  }
1232  BitReverse(working_space, fSize);
1233  Walsh(working_space, fSize);
1234  for (i = 0; i < fSize; i++) {
1235  destVector[i] = working_space[i];
1236  }
1237  break;
1238  case kTransformCos:
1239  for (i = 0; i < fSize; i++) {
1240  working_space[i] = source[i];
1241  }
1242  fSize = 2 * fSize;
1243  working_space[0] = working_space[0] * TMath::Sqrt(2.0);
1244  for (i = 0; i < fSize / 2; i++) {
1245  a = pi * (Double_t) i / (Double_t) fSize;
1246  b = TMath::Sin(a);
1247  a = TMath::Cos(a);
1248  working_space[i + fSize] = (Double_t) working_space[i] * b;
1249  working_space[i] = (Double_t) working_space[i] * a;
1250  } for (i = 2; i <= (fSize / 2); i++) {
1251  working_space[fSize - i + 1] = working_space[i - 1];
1252  working_space[fSize - i + 1 + fSize] =
1253  -working_space[i - 1 + fSize];
1254  }
1255  working_space[fSize / 2] = 0;
1256  working_space[fSize / 2 + fSize] = 0;
1257  Fourier(working_space, fSize, 0, kTransformInverse, 1);
1258  for (i = 0; i < fSize / 2; i++) {
1259  destVector[i] = working_space[i];
1260  }
1261  break;
1262  case kTransformSin:
1263  for (i = 0; i < fSize; i++) {
1264  working_space[i] = source[i];
1265  }
1266  fSize = 2 * fSize;
1267  working_space[fSize / 2] =
1268  working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
1269  for (i = fSize / 2 - 1; i > 0; i--) {
1270  a = pi * (Double_t) i / (Double_t) fSize;
1271  working_space[i + fSize] =
1272  -(Double_t) working_space[i - 1] * TMath::Cos(a);
1273  working_space[i] =
1274  (Double_t) working_space[i - 1] * TMath::Sin(a);
1275  } for (i = 2; i <= (fSize / 2); i++) {
1276  working_space[fSize - i + 1] = working_space[i - 1];
1277  working_space[fSize - i + 1 + fSize] =
1278  -working_space[i - 1 + fSize];
1279  }
1280  working_space[0] = 0;
1281  working_space[fSize] = 0;
1282  working_space[fSize / 2 + fSize] = 0;
1283  Fourier(working_space, fSize, 0, kTransformInverse, 0);
1284  for (i = 0; i < fSize / 2; i++) {
1285  destVector[i] = working_space[i];
1286  }
1287  break;
1288  case kTransformFourier:
1289  for (i = 0; i < 2 * fSize; i++) {
1290  working_space[i] = source[i];
1291  }
1292  Fourier(working_space, fSize, 0, kTransformInverse, 0);
1293  for (i = 0; i < fSize; i++) {
1294  destVector[i] = working_space[i];
1295  }
1296  break;
1297  case kTransformHartley:
1298  for (i = 0; i < fSize; i++) {
1299  working_space[i] = source[i];
1300  }
1301  Fourier(working_space, fSize, 1, kTransformInverse, 0);
1302  for (i = 0; i < fSize; i++) {
1303  destVector[i] = working_space[i];
1304  }
1305  break;
1306  case kTransformFourierWalsh:
1307  case kTransformFourierHaar:
1308  case kTransformWalshHaar:
1309  case kTransformCosWalsh:
1310  case kTransformCosHaar:
1311  case kTransformSinWalsh:
1312  case kTransformSinHaar:
1313  for (i = 0; i < fSize; i++) {
1314  working_space[i] = source[i];
1315  }
1316  if (fTransformType == kTransformFourierWalsh
1317  || fTransformType == kTransformFourierHaar) {
1318  for (i = 0; i < fSize; i++) {
1319  working_space[i + 2 * fSize] = source[fSize + i];
1320  }
1321  }
1322  if (fTransformType > kTransformWalshHaar)
1323  k = (Int_t) TMath::Power(2, fDegree - 1);
1324 
1325  else
1326  k = (Int_t) TMath::Power(2, fDegree);
1327  j = fSize / k;
1328  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1329  working_space[fSize + l + i / j] = working_space[i];
1330  working_space[fSize + l + i / j + 2 * fSize] =
1331  working_space[i + 2 * fSize];
1332  }
1333  for (i = 0; i < fSize; i++) {
1334  working_space[i] = working_space[fSize + i];
1335  working_space[i + 2 * fSize] =
1336  working_space[fSize + i + 2 * fSize];
1337  }
1338  if (fTransformType == kTransformFourierWalsh
1339  || fTransformType == kTransformFourierHaar
1340  || fTransformType == kTransformWalshHaar) {
1341  GeneralInv(working_space, fSize, fDegree, fTransformType);
1342  for (i = 0; i < j; i++)
1343  BitReverseHaar(working_space, fSize, k, i * k);
1344  }
1345 
1346  else if (fTransformType == kTransformCosWalsh
1347  || fTransformType == kTransformCosHaar) {
1348  j = (Int_t) TMath::Power(2, fDegree) / 2;
1349  m = (Int_t) TMath::Power(2, fDegree);
1350  l = 2 * fSize / m;
1351  for (i = 0; i < fSize; i++) {
1352  k = i / j;
1353  k = 2 * k * j;
1354  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1355  if (i % j == 0) {
1356  working_space[2 * fSize + k + i % j] =
1357  working_space[i] * TMath::Sqrt(2.0);
1358  working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
1359  }
1360 
1361  else {
1362  b = TMath::Sin(a);
1363  a = TMath::Cos(a);
1364  working_space[4 * fSize + 2 * fSize + k + i % j] =
1365  -(Double_t) working_space[i] * b;
1366  working_space[2 * fSize + k + i % j] =
1367  (Double_t) working_space[i] * a;
1368  } } for (i = 0; i < fSize; i++) {
1369  k = i / j;
1370  k = 2 * k * j;
1371  if (i % j == 0) {
1372  working_space[2 * fSize + k + j] = 0;
1373  working_space[4 * fSize + 2 * fSize + k + j] = 0;
1374  }
1375 
1376  else {
1377  working_space[2 * fSize + k + 2 * j - i % j] =
1378  working_space[2 * fSize + k + i % j];
1379  working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
1380  -working_space[4 * fSize + 2 * fSize + k + i % j];
1381  }
1382  }
1383  for (i = 0; i < 2 * fSize; i++) {
1384  working_space[i] = working_space[2 * fSize + i];
1385  working_space[4 * fSize + i] =
1386  working_space[4 * fSize + 2 * fSize + i];
1387  }
1388  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1389  m = (Int_t) TMath::Power(2, fDegree);
1390  l = 2 * fSize / m;
1391  for (i = 0; i < l; i++)
1392  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1393  }
1394 
1395  else if (fTransformType == kTransformSinWalsh
1396  || fTransformType == kTransformSinHaar) {
1397  j = (Int_t) TMath::Power(2, fDegree) / 2;
1398  m = (Int_t) TMath::Power(2, fDegree);
1399  l = 2 * fSize / m;
1400  for (i = 0; i < fSize; i++) {
1401  k = i / j;
1402  k = 2 * k * j;
1403  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1404  if (i % j == 0) {
1405  working_space[2 * fSize + k + j + i % j] =
1406  working_space[j + k / 2 - i % j -
1407  1] * TMath::Sqrt(2.0);
1408  working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
1409  }
1410 
1411  else {
1412  b = TMath::Sin(a);
1413  a = TMath::Cos(a);
1414  working_space[4 * fSize + 2 * fSize + k + j + i % j] =
1415  -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
1416  working_space[2 * fSize + k + j + i % j] =
1417  (Double_t) working_space[j + k / 2 - i % j - 1] * a;
1418  } } for (i = 0; i < fSize; i++) {
1419  k = i / j;
1420  k = 2 * k * j;
1421  if (i % j == 0) {
1422  working_space[2 * fSize + k] = 0;
1423  working_space[4 * fSize + 2 * fSize + k] = 0;
1424  }
1425 
1426  else {
1427  working_space[2 * fSize + k + i % j] =
1428  working_space[2 * fSize + k + 2 * j - i % j];
1429  working_space[4 * fSize + 2 * fSize + k + i % j] =
1430  -working_space[4 * fSize + 2 * fSize + k + 2 * j -
1431  i % j];
1432  }
1433  }
1434  for (i = 0; i < 2 * fSize; i++) {
1435  working_space[i] = working_space[2 * fSize + i];
1436  working_space[4 * fSize + i] =
1437  working_space[4 * fSize + 2 * fSize + i];
1438  }
1439  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
1440  for (i = 0; i < l; i++)
1441  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1442  }
1443  for (i = 0; i < fSize; i++) {
1444  if (fTransformType >= kTransformCosWalsh) {
1445  k = i / j;
1446  k = 2 * k * j;
1447  val = working_space[k + i % j];
1448  }
1449 
1450  else
1451  val = working_space[i];
1452  destVector[i] = val;
1453  }
1454  break;
1455  }
1456  }
1457  delete[]working_space;
1458  return;
1459 }
1460 
1461 //////////FilterZonal FUNCTION - CALCULATES DIFFERENT 1-D ORTHOGONAL TRANSFORMS, SETS GIVEN REGION TO FILTER COEFFICIENT AND TRANSFORMS IT BACK//////
1462 
1463 ////////////////////////////////////////////////////////////////////////////////
1464 /////////////////////////////////////////////////////////////////////////////////
1465 /// ONE-DIMENSIONAL FILTER ZONAL FUNCTION
1466 /// This function transforms the source spectrum. The calling program
1467 /// should fill in input parameters. Then it sets transformed
1468 /// coefficients in the given region (fXmin, fXmax) to the given
1469 /// fFilterCoeff and transforms it back.
1470 /// Filtered data are written into dest spectrum.
1471 ///
1472 /// Function parameters:
1473 /// source-pointer to the vector of source spectrum, its length should
1474 /// be size except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR
1475 /// transform. These need 2*size length to supply real and
1476 /// imaginary coefficients.
1477 /// destVector-pointer to the vector of dest data, its length should be
1478 /// size except for direct FOURIER, FOUR-WALSH, FOUR-HAAR. These
1479 /// need 2*size length to store real and imaginary coefficients
1480 ///
1481 /////////////////////////////////////////////////////////////////////////////////
1482 ///
1483 
1484 void TSpectrumTransform::FilterZonal(const Double_t *source, Double_t *destVector)
1485 {
1486 /*
1487 <div class=Section2>
1488 
1489 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of zonal filtering</span></b></p>
1490 
1491 <p class=MsoNormal><i>&nbsp;</i></p>
1492 
1493 <p class=MsoNormal><i>Function:</i></p>
1494 
1495 <p class=MsoNormal><b>void TSpectrumTransform::FilterZonal(const <a
1496 href="http://root.cern.ch/root/html/ListOfTypes.html#Double_t">Double_t</a> *fSource,
1497 <a href="http://root.cern.ch/root/html/ListOfTypes.html#Double_t">Double_t</a> *fDest)</b></p>
1498 
1499 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
1500 
1501 <p class=MsoNormal style='text-align:justify'>This function transforms the
1502 source spectrum (for details see Transform function). Before the FilterZonal
1503 function is called the class must be created by constructor and the type of the
1504 transform as well as some other parameters must be set using a set of setter
1505 functions. The FilterZonal function sets transformed coefficients in the given
1506 region (fXmin, fXmax) to the given fFilterCoeff and transforms it back.
1507 Filtered data are written into dest spectrum. </p>
1508 
1509 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'>&nbsp;</span></i></p>
1510 
1511 <p class=MsoNormal style='text-align:justify'><i>Example – script Filter.c:</i></p>
1512 
1513 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
1514 border=0 width=601 height=402 src="gif/spectrumtransform_filter_image001.jpg"></span></i></p>
1515 
1516 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum
1517 (black line) and filtered spectrum (red line) using Cosine transform and zonal
1518 filtration (channels 2048-4095 were set to 0) </b></p>
1519 
1520 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
1521 
1522 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
1523 
1524 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
1525 FilterZonal function (class TSpectrumTransform).</span></p>
1526 
1527 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
1528 do</span></p>
1529 
1530 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Filter.C</span></p>
1531 
1532 <p class=MsoNormal><span style='font-size:10.0pt'>   </span></p>
1533 
1534 <p class=MsoNormal><span style='font-size:10.0pt'>#include &lt;TSpectrum&gt;</span></p>
1535 
1536 <p class=MsoNormal><span style='font-size:10.0pt'>#include
1537 &lt;TSpectrumTransform&gt;</span></p>
1538 
1539 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
1540 
1541 <p class=MsoNormal><span style='font-size:10.0pt'>void Filter() {</span></p>
1542 
1543 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
1544 
1545 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
1546 4096;</span></p>
1547 
1548 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
1549 0;</span></p>
1550 
1551 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
1552 (Double_t)nbins;</span></p>
1553 
1554 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
1555 style='font-size:10.0pt'>Double_t * source = new Double_t[nbins];</span></p>
1556 
1557 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t * dest = new
1558 Double_t[nbins];   </span></p>
1559 
1560 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new
1561 TH1F(&quot;h&quot;,&quot;Zonal filtering using Cosine
1562 transform&quot;,nbins,xmin,xmax);</span></p>
1563 
1564 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
1565 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);         </span></p>
1566 
1567 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
1568 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
1569 
1570 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
1571 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
1572 
1573 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
1574 i++) source[i]=h-&gt;GetBinContent(i + 1);     </span></p>
1575 
1576 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 =
1577 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
1578 
1579 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
1580 Transform1 = new
1581 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
1582 
1583 <p class=MsoNormal><span style='font-size:10.0pt'>  
1584 h-&gt;SetAxisRange(700,1024);   </span></p>
1585 
1586 <p class=MsoNormal><span style='font-size:10.0pt'>  
1587 h-&gt;Draw(&quot;L&quot;);   </span></p>
1588 
1589 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
1590 TSpectrum();</span></p>
1591 
1592 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
1593 new TSpectrumTransform(4096);</span></p>
1594 
1595 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
1596 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
1597 
1598 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1599 t-&gt;SetRegion(2048, 4095);</span></p>
1600 
1601 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1602 t-&gt;FilterZonal(source,dest);     </span></p>
1603 
1604 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
1605 style='font-size:10.0pt'>for (i = 0; i &lt; nbins; i++) d-&gt;SetBinContent(i +
1606 1,dest[i]);</span></p>
1607 
1608 <p class=MsoNormal><span style='font-size:10.0pt'>  
1609 d-&gt;SetLineColor(kRed);   </span></p>
1610 
1611 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
1612 L&quot;);</span></p>
1613 
1614 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
1615 
1616 </div>
1617 
1618  */
1619  int i, j=0, k = 1, m, l;
1620  Double_t val;
1621  Double_t *working_space = 0;
1622  Double_t a, b, pi = 3.14159265358979323846, old_area, new_area;
1623  if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
1624  if (fTransformType >= kTransformCosWalsh)
1625  fDegree += 1;
1626  k = (Int_t) TMath::Power(2, fDegree);
1627  j = fSize / k;
1628  }
1629  switch (fTransformType) {
1630  case kTransformHaar:
1631  case kTransformWalsh:
1632  working_space = new Double_t[2 * fSize];
1633  break;
1634  case kTransformCos:
1635  case kTransformSin:
1636  case kTransformFourier:
1637  case kTransformHartley:
1638  case kTransformFourierWalsh:
1639  case kTransformFourierHaar:
1640  case kTransformWalshHaar:
1641  working_space = new Double_t[4 * fSize];
1642  break;
1643  case kTransformCosWalsh:
1644  case kTransformCosHaar:
1645  case kTransformSinWalsh:
1646  case kTransformSinHaar:
1647  working_space = new Double_t[8 * fSize];
1648  break;
1649  }
1650  switch (fTransformType) {
1651  case kTransformHaar:
1652  for (i = 0; i < fSize; i++) {
1653  working_space[i] = source[i];
1654  }
1655  Haar(working_space, fSize, kTransformForward);
1656  break;
1657  case kTransformWalsh:
1658  for (i = 0; i < fSize; i++) {
1659  working_space[i] = source[i];
1660  }
1661  Walsh(working_space, fSize);
1662  BitReverse(working_space, fSize);
1663  break;
1664  case kTransformCos:
1665  fSize = 2 * fSize;
1666  for (i = 1; i <= (fSize / 2); i++) {
1667  val = source[i - 1];
1668  working_space[i - 1] = val;
1669  working_space[fSize - i] = val;
1670  }
1671  Fourier(working_space, fSize, 0, kTransformForward, 0);
1672  for (i = 0; i < fSize / 2; i++) {
1673  a = pi * (Double_t) i / (Double_t) fSize;
1674  a = TMath::Cos(a);
1675  b = working_space[i];
1676  a = b / a;
1677  working_space[i] = a;
1678  working_space[i + fSize] = 0;
1679  } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
1680  fSize = fSize / 2;
1681  break;
1682  case kTransformSin:
1683  fSize = 2 * fSize;
1684  for (i = 1; i <= (fSize / 2); i++) {
1685  val = source[i - 1];
1686  working_space[i - 1] = val;
1687  working_space[fSize - i] = -val;
1688  }
1689  Fourier(working_space, fSize, 0, kTransformForward, 0);
1690  for (i = 0; i < fSize / 2; i++) {
1691  a = pi * (Double_t) i / (Double_t) fSize;
1692  a = TMath::Sin(a);
1693  b = working_space[i];
1694  if (a != 0)
1695  a = b / a;
1696  working_space[i - 1] = a;
1697  working_space[i + fSize] = 0;
1698  }
1699  working_space[fSize / 2 - 1] =
1700  working_space[fSize / 2] / TMath::Sqrt(2.0);
1701  fSize = fSize / 2;
1702  break;
1703  case kTransformFourier:
1704  for (i = 0; i < fSize; i++) {
1705  working_space[i] = source[i];
1706  }
1707  Fourier(working_space, fSize, 0, kTransformForward, 0);
1708  break;
1709  case kTransformHartley:
1710  for (i = 0; i < fSize; i++) {
1711  working_space[i] = source[i];
1712  }
1713  Fourier(working_space, fSize, 1, kTransformForward, 0);
1714  break;
1715  case kTransformFourierWalsh:
1716  case kTransformFourierHaar:
1717  case kTransformWalshHaar:
1718  case kTransformCosWalsh:
1719  case kTransformCosHaar:
1720  case kTransformSinWalsh:
1721  case kTransformSinHaar:
1722  for (i = 0; i < fSize; i++) {
1723  val = source[i];
1724  if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1725  j = (Int_t) TMath::Power(2, fDegree) / 2;
1726  k = i / j;
1727  k = 2 * k * j;
1728  working_space[k + i % j] = val;
1729  working_space[k + 2 * j - 1 - i % j] = val;
1730  }
1731 
1732  else if (fTransformType == kTransformSinWalsh
1733  || fTransformType == kTransformSinHaar) {
1734  j = (Int_t) TMath::Power(2, fDegree) / 2;
1735  k = i / j;
1736  k = 2 * k * j;
1737  working_space[k + i % j] = val;
1738  working_space[k + 2 * j - 1 - i % j] = -val;
1739  }
1740 
1741  else
1742  working_space[i] = val;
1743  }
1744  if (fTransformType == kTransformFourierWalsh
1745  || fTransformType == kTransformFourierHaar
1746  || fTransformType == kTransformWalshHaar) {
1747  for (i = 0; i < j; i++)
1748  BitReverseHaar(working_space, fSize, k, i * k);
1749  GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
1750  }
1751 
1752  else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1753  m = (Int_t) TMath::Power(2, fDegree);
1754  l = 2 * fSize / m;
1755  for (i = 0; i < l; i++)
1756  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1757  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1758  for (i = 0; i < fSize; i++) {
1759  k = i / j;
1760  k = 2 * k * j;
1761  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1762  a = TMath::Cos(a);
1763  b = working_space[k + i % j];
1764  if (i % j == 0)
1765  a = b / TMath::Sqrt(2.0);
1766 
1767  else
1768  a = b / a;
1769  working_space[i] = a;
1770  working_space[i + 2 * fSize] = 0;
1771  }
1772  }
1773 
1774  else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
1775  m = (Int_t) TMath::Power(2, fDegree);
1776  l = 2 * fSize / m;
1777  for (i = 0; i < l; i++)
1778  BitReverseHaar(working_space, 2 * fSize, m, i * m);
1779  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
1780  for (i = 0; i < fSize; i++) {
1781  k = i / j;
1782  k = 2 * k * j;
1783  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1784  a = TMath::Cos(a);
1785  b = working_space[j + k + i % j];
1786  if (i % j == 0)
1787  a = b / TMath::Sqrt(2.0);
1788 
1789  else
1790  a = b / a;
1791  working_space[j + k / 2 - i % j - 1] = a;
1792  working_space[i + 2 * fSize] = 0;
1793  }
1794  }
1795  if (fTransformType > kTransformWalshHaar)
1796  k = (Int_t) TMath::Power(2, fDegree - 1);
1797 
1798  else
1799  k = (Int_t) TMath::Power(2, fDegree);
1800  j = fSize / k;
1801  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1802  working_space[fSize + i] = working_space[l + i / j];
1803  working_space[fSize + i + 2 * fSize] =
1804  working_space[l + i / j + 2 * fSize];
1805  }
1806  for (i = 0; i < fSize; i++) {
1807  working_space[i] = working_space[fSize + i];
1808  working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
1809  }
1810  break;
1811  }
1812  if (!working_space) return;
1813  for (i = 0, old_area = 0; i < fSize; i++) {
1814  old_area += working_space[i];
1815  }
1816  for (i = 0, new_area = 0; i < fSize; i++) {
1817  if (i >= fXmin && i <= fXmax)
1818  working_space[i] = fFilterCoeff;
1819  new_area += working_space[i];
1820  }
1821  if (new_area != 0) {
1822  a = old_area / new_area;
1823  for (i = 0; i < fSize; i++) {
1824  working_space[i] *= a;
1825  }
1826  }
1827  if (fTransformType == kTransformFourier) {
1828  for (i = 0, old_area = 0; i < fSize; i++) {
1829  old_area += working_space[i + fSize];
1830  }
1831  for (i = 0, new_area = 0; i < fSize; i++) {
1832  if (i >= fXmin && i <= fXmax)
1833  working_space[i + fSize] = fFilterCoeff;
1834  new_area += working_space[i + fSize];
1835  }
1836  if (new_area != 0) {
1837  a = old_area / new_area;
1838  for (i = 0; i < fSize; i++) {
1839  working_space[i + fSize] *= a;
1840  }
1841  }
1842  }
1843 
1844  else if (fTransformType == kTransformFourierWalsh
1845  || fTransformType == kTransformFourierHaar) {
1846  for (i = 0, old_area = 0; i < fSize; i++) {
1847  old_area += working_space[i + 2 * fSize];
1848  }
1849  for (i = 0, new_area = 0; i < fSize; i++) {
1850  if (i >= fXmin && i <= fXmax)
1851  working_space[i + 2 * fSize] = fFilterCoeff;
1852  new_area += working_space[i + 2 * fSize];
1853  }
1854  if (new_area != 0) {
1855  a = old_area / new_area;
1856  for (i = 0; i < fSize; i++) {
1857  working_space[i + 2 * fSize] *= a;
1858  }
1859  }
1860  }
1861  switch (fTransformType) {
1862  case kTransformHaar:
1863  Haar(working_space, fSize, kTransformInverse);
1864  for (i = 0; i < fSize; i++) {
1865  destVector[i] = working_space[i];
1866  }
1867  break;
1868  case kTransformWalsh:
1869  BitReverse(working_space, fSize);
1870  Walsh(working_space, fSize);
1871  for (i = 0; i < fSize; i++) {
1872  destVector[i] = working_space[i];
1873  }
1874  break;
1875  case kTransformCos:
1876  fSize = 2 * fSize;
1877  working_space[0] = working_space[0] * TMath::Sqrt(2.0);
1878  for (i = 0; i < fSize / 2; i++) {
1879  a = pi * (Double_t) i / (Double_t) fSize;
1880  b = TMath::Sin(a);
1881  a = TMath::Cos(a);
1882  working_space[i + fSize] = (Double_t) working_space[i] * b;
1883  working_space[i] = (Double_t) working_space[i] * a;
1884  } for (i = 2; i <= (fSize / 2); i++) {
1885  working_space[fSize - i + 1] = working_space[i - 1];
1886  working_space[fSize - i + 1 + fSize] =
1887  -working_space[i - 1 + fSize];
1888  }
1889  working_space[fSize / 2] = 0;
1890  working_space[fSize / 2 + fSize] = 0;
1891  Fourier(working_space, fSize, 0, kTransformInverse, 1);
1892  for (i = 0; i < fSize / 2; i++) {
1893  destVector[i] = working_space[i];
1894  }
1895  break;
1896  case kTransformSin:
1897  fSize = 2 * fSize;
1898  working_space[fSize / 2] =
1899  working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
1900  for (i = fSize / 2 - 1; i > 0; i--) {
1901  a = pi * (Double_t) i / (Double_t) fSize;
1902  working_space[i + fSize] =
1903  -(Double_t) working_space[i - 1] * TMath::Cos(a);
1904  working_space[i] = (Double_t) working_space[i - 1] * TMath::Sin(a);
1905  } for (i = 2; i <= (fSize / 2); i++) {
1906  working_space[fSize - i + 1] = working_space[i - 1];
1907  working_space[fSize - i + 1 + fSize] =
1908  -working_space[i - 1 + fSize];
1909  }
1910  working_space[0] = 0;
1911  working_space[fSize] = 0;
1912  working_space[fSize / 2 + fSize] = 0;
1913  Fourier(working_space, fSize, 0, kTransformInverse, 0);
1914  for (i = 0; i < fSize / 2; i++) {
1915  destVector[i] = working_space[i];
1916  }
1917  break;
1918  case kTransformFourier:
1919  Fourier(working_space, fSize, 0, kTransformInverse, 0);
1920  for (i = 0; i < fSize; i++) {
1921  destVector[i] = working_space[i];
1922  }
1923  break;
1924  case kTransformHartley:
1925  Fourier(working_space, fSize, 1, kTransformInverse, 0);
1926  for (i = 0; i < fSize; i++) {
1927  destVector[i] = working_space[i];
1928  }
1929  break;
1930  case kTransformFourierWalsh:
1931  case kTransformFourierHaar:
1932  case kTransformWalshHaar:
1933  case kTransformCosWalsh:
1934  case kTransformCosHaar:
1935  case kTransformSinWalsh:
1936  case kTransformSinHaar:
1937  if (fTransformType > kTransformWalshHaar)
1938  k = (Int_t) TMath::Power(2, fDegree - 1);
1939 
1940  else
1941  k = (Int_t) TMath::Power(2, fDegree);
1942  j = fSize / k;
1943  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
1944  working_space[fSize + l + i / j] = working_space[i];
1945  working_space[fSize + l + i / j + 2 * fSize] =
1946  working_space[i + 2 * fSize];
1947  }
1948  for (i = 0; i < fSize; i++) {
1949  working_space[i] = working_space[fSize + i];
1950  working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
1951  }
1952  if (fTransformType == kTransformFourierWalsh
1953  || fTransformType == kTransformFourierHaar
1954  || fTransformType == kTransformWalshHaar) {
1955  GeneralInv(working_space, fSize, fDegree, fTransformType);
1956  for (i = 0; i < j; i++)
1957  BitReverseHaar(working_space, fSize, k, i * k);
1958  }
1959 
1960  else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
1961  j = (Int_t) TMath::Power(2, fDegree) / 2;
1962  m = (Int_t) TMath::Power(2, fDegree);
1963  l = 2 * fSize / m;
1964  for (i = 0; i < fSize; i++) {
1965  k = i / j;
1966  k = 2 * k * j;
1967  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
1968  if (i % j == 0) {
1969  working_space[2 * fSize + k + i % j] =
1970  working_space[i] * TMath::Sqrt(2.0);
1971  working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
1972  }
1973 
1974  else {
1975  b = TMath::Sin(a);
1976  a = TMath::Cos(a);
1977  working_space[4 * fSize + 2 * fSize + k + i % j] =
1978  -(Double_t) working_space[i] * b;
1979  working_space[2 * fSize + k + i % j] =
1980  (Double_t) working_space[i] * a;
1981  } } for (i = 0; i < fSize; i++) {
1982  k = i / j;
1983  k = 2 * k * j;
1984  if (i % j == 0) {
1985  working_space[2 * fSize + k + j] = 0;
1986  working_space[4 * fSize + 2 * fSize + k + j] = 0;
1987  }
1988 
1989  else {
1990  working_space[2 * fSize + k + 2 * j - i % j] =
1991  working_space[2 * fSize + k + i % j];
1992  working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
1993  -working_space[4 * fSize + 2 * fSize + k + i % j];
1994  }
1995  }
1996  for (i = 0; i < 2 * fSize; i++) {
1997  working_space[i] = working_space[2 * fSize + i];
1998  working_space[4 * fSize + i] =
1999  working_space[4 * fSize + 2 * fSize + i];
2000  }
2001  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2002  m = (Int_t) TMath::Power(2, fDegree);
2003  l = 2 * fSize / m;
2004  for (i = 0; i < l; i++)
2005  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2006  }
2007 
2008  else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
2009  j = (Int_t) TMath::Power(2, fDegree) / 2;
2010  m = (Int_t) TMath::Power(2, fDegree);
2011  l = 2 * fSize / m;
2012  for (i = 0; i < fSize; i++) {
2013  k = i / j;
2014  k = 2 * k * j;
2015  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2016  if (i % j == 0) {
2017  working_space[2 * fSize + k + j + i % j] =
2018  working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
2019  working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
2020  }
2021 
2022  else {
2023  b = TMath::Sin(a);
2024  a = TMath::Cos(a);
2025  working_space[4 * fSize + 2 * fSize + k + j + i % j] =
2026  -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
2027  working_space[2 * fSize + k + j + i % j] =
2028  (Double_t) working_space[j + k / 2 - i % j - 1] * a;
2029  } } for (i = 0; i < fSize; i++) {
2030  k = i / j;
2031  k = 2 * k * j;
2032  if (i % j == 0) {
2033  working_space[2 * fSize + k] = 0;
2034  working_space[4 * fSize + 2 * fSize + k] = 0;
2035  }
2036 
2037  else {
2038  working_space[2 * fSize + k + i % j] =
2039  working_space[2 * fSize + k + 2 * j - i % j];
2040  working_space[4 * fSize + 2 * fSize + k + i % j] =
2041  -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
2042  }
2043  }
2044  for (i = 0; i < 2 * fSize; i++) {
2045  working_space[i] = working_space[2 * fSize + i];
2046  working_space[4 * fSize + i] =
2047  working_space[4 * fSize + 2 * fSize + i];
2048  }
2049  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2050  for (i = 0; i < l; i++)
2051  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2052  }
2053  for (i = 0; i < fSize; i++) {
2054  if (fTransformType >= kTransformCosWalsh) {
2055  k = i / j;
2056  k = 2 * k * j;
2057  val = working_space[k + i % j];
2058  }
2059 
2060  else
2061  val = working_space[i];
2062  destVector[i] = val;
2063  }
2064  break;
2065  }
2066  delete[]working_space;
2067  return;
2068 }
2069 
2070 //////////ENHANCE FUNCTION - CALCULATES DIFFERENT 1-D ORTHOGONAL TRANSFORMS, MULTIPLIES GIVEN REGION BY ENHANCE COEFFICIENT AND TRANSFORMS IT BACK//////
2071 ////////////////////////////////////////////////////////////////////////////////
2072 /////////////////////////////////////////////////////////////////////////////////
2073 /// ONE-DIMENSIONAL ENHANCE ZONAL FUNCTION
2074 /// This function transforms the source spectrum. The calling program
2075 /// should fill in input parameters. Then it multiplies transformed
2076 /// coefficients in the given region (fXmin, fXmax) by the given
2077 /// fEnhanceCoeff and transforms it back
2078 /// Processed data are written into dest spectrum.
2079 ///
2080 /// Function parameters:
2081 /// source-pointer to the vector of source spectrum, its length should
2082 /// be size except for inverse FOURIER, FOUR-WALSh, FOUR-HAAR
2083 /// transform. These need 2*size length to supply real and
2084 /// imaginary coefficients.
2085 /// destVector-pointer to the vector of dest data, its length should be
2086 /// size except for direct FOURIER, FOUR-WALSh, FOUR-HAAR. These
2087 /// need 2*size length to store real and imaginary coefficients
2088 ///
2089 /////////////////////////////////////////////////////////////////////////////////
2090 
2091 void TSpectrumTransform::Enhance(const Double_t *source, Double_t *destVector)
2092 {
2093 /*
2094 <div class=Section3>
2095 
2096 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of enhancement</span></b></p>
2097 
2098 <p class=MsoNormal><i>&nbsp;</i></p>
2099 
2100 <p class=MsoNormal><i>Function:</i></p>
2101 
2102 <p class=MsoNormal><b>void TSpectrumTransform::Enhance(const <a
2103 href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> *fSource,
2104 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a>
2105 *fDest)</b></p>
2106 
2107 <p class=MsoNormal><b>&nbsp;</b></p>
2108 
2109 <p class=MsoNormal style='text-align:justify'>This function transforms the
2110 source spectrum (for details see Transform function). Before the Enhance
2111 function is called the class must be created by constructor and the type of the
2112 transform as well as some other parameters must be set using a set of setter functions.
2113 The Enhance function multiplies transformed coefficients in the given region
2114 (fXmin, fXmax) by the given fEnhancCoeff and transforms it back. Enhanced data
2115 are written into dest spectrum.</p>
2116 
2117 <p class=MsoNormal>&nbsp;</p>
2118 
2119 <p class=MsoNormal style='text-align:justify'><i>Example  – script Enhance.c:</i></p>
2120 
2121 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
2122 border=0 width=601 height=402 src="gif/spectrumtransform_enhance_image001.jpg"></span></i></p>
2123 
2124 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'>&nbsp;</span></p>
2125 
2126 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum (black
2127 line) and enhanced spectrum (red line) using Cosine transform (channels 0-1024
2128 were multiplied by 2) </b></p>
2129 
2130 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
2131 
2132 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
2133 
2134 <p class=MsoNormal>&nbsp;</p>
2135 
2136 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
2137 Enhance function (class TSpectrumTransform).</span></p>
2138 
2139 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
2140 do</span></p>
2141 
2142 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Enhance.C</span></p>
2143 
2144 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
2145 
2146 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Enhance() {</span></p>
2147 
2148 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
2149 
2150 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
2151 4096;</span></p>
2152 
2153 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
2154 0;</span></p>
2155 
2156 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
2157 (Double_t)nbins;</span></p>
2158 
2159 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2160 style='font-size:10.0pt'>Double_t * source = new Double_t[nbins];</span></p>
2161 
2162 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t * dest = new
2163 Double_t[nbins];   </span></p>
2164 
2165 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new
2166 TH1F(&quot;h&quot;,&quot;Enhancement using Cosine transform&quot;,nbins,xmin,xmax);</span></p>
2167 
2168 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
2169 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);         </span></p>
2170 
2171 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
2172 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
2173 
2174 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
2175 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
2176 
2177 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
2178 i++) source[i]=h-&gt;GetBinContent(i + 1);     </span></p>
2179 
2180 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 = gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
2181 
2182 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
2183 Transform1 = new
2184 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
2185 
2186 <p class=MsoNormal><span style='font-size:10.0pt'>  
2187 h-&gt;SetAxisRange(700,1024);   </span></p>
2188 
2189 <p class=MsoNormal><span style='font-size:10.0pt'>  
2190 h-&gt;Draw(&quot;L&quot;);   </span></p>
2191 
2192 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
2193 TSpectrum();</span></p>
2194 
2195 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
2196 new TSpectrumTransform(4096);</span></p>
2197 
2198 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
2199 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
2200 
2201 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   t-&gt;SetRegion(0,
2202 1024);</span></p>
2203 
2204 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2205 t-&gt;SetEnhanceCoeff(2);</span></p>
2206 
2207 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2208 t-&gt;Enhance(source,dest);        </span></p>
2209 
2210 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2211 style='font-size:10.0pt'>for (i = 0; i &lt; nbins; i++) d-&gt;SetBinContent(i +
2212 1,dest[i]);</span></p>
2213 
2214 <p class=MsoNormal><span style='font-size:10.0pt'>  
2215 d-&gt;SetLineColor(kRed);   </span></p>
2216 
2217 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
2218 L&quot;);</span></p>
2219 
2220 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
2221 
2222 </div>
2223 
2224  */
2225  int i, j=0, k = 1, m, l;
2226  Double_t val;
2227  Double_t *working_space = 0;
2228  Double_t a, b, pi = 3.14159265358979323846, old_area, new_area;
2229  if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
2230  if (fTransformType >= kTransformCosWalsh)
2231  fDegree += 1;
2232  k = (Int_t) TMath::Power(2, fDegree);
2233  j = fSize / k;
2234  }
2235  switch (fTransformType) {
2236  case kTransformHaar:
2237  case kTransformWalsh:
2238  working_space = new Double_t[2 * fSize];
2239  break;
2240  case kTransformCos:
2241  case kTransformSin:
2242  case kTransformFourier:
2243  case kTransformHartley:
2244  case kTransformFourierWalsh:
2245  case kTransformFourierHaar:
2246  case kTransformWalshHaar:
2247  working_space = new Double_t[4 * fSize];
2248  break;
2249  case kTransformCosWalsh:
2250  case kTransformCosHaar:
2251  case kTransformSinWalsh:
2252  case kTransformSinHaar:
2253  working_space = new Double_t[8 * fSize];
2254  break;
2255  }
2256  switch (fTransformType) {
2257  case kTransformHaar:
2258  for (i = 0; i < fSize; i++) {
2259  working_space[i] = source[i];
2260  }
2261  Haar(working_space, fSize, kTransformForward);
2262  break;
2263  case kTransformWalsh:
2264  for (i = 0; i < fSize; i++) {
2265  working_space[i] = source[i];
2266  }
2267  Walsh(working_space, fSize);
2268  BitReverse(working_space, fSize);
2269  break;
2270  case kTransformCos:
2271  fSize = 2 * fSize;
2272  for (i = 1; i <= (fSize / 2); i++) {
2273  val = source[i - 1];
2274  working_space[i - 1] = val;
2275  working_space[fSize - i] = val;
2276  }
2277  Fourier(working_space, fSize, 0, kTransformForward, 0);
2278  for (i = 0; i < fSize / 2; i++) {
2279  a = pi * (Double_t) i / (Double_t) fSize;
2280  a = TMath::Cos(a);
2281  b = working_space[i];
2282  a = b / a;
2283  working_space[i] = a;
2284  working_space[i + fSize] = 0;
2285  } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
2286  fSize = fSize / 2;
2287  break;
2288  case kTransformSin:
2289  fSize = 2 * fSize;
2290  for (i = 1; i <= (fSize / 2); i++) {
2291  val = source[i - 1];
2292  working_space[i - 1] = val;
2293  working_space[fSize - i] = -val;
2294  }
2295  Fourier(working_space, fSize, 0, kTransformForward, 0);
2296  for (i = 0; i < fSize / 2; i++) {
2297  a = pi * (Double_t) i / (Double_t) fSize;
2298  a = TMath::Sin(a);
2299  b = working_space[i];
2300  if (a != 0)
2301  a = b / a;
2302  working_space[i - 1] = a;
2303  working_space[i + fSize] = 0;
2304  }
2305  working_space[fSize / 2 - 1] =
2306  working_space[fSize / 2] / TMath::Sqrt(2.0);
2307  fSize = fSize / 2;
2308  break;
2309  case kTransformFourier:
2310  for (i = 0; i < fSize; i++) {
2311  working_space[i] = source[i];
2312  }
2313  Fourier(working_space, fSize, 0, kTransformForward, 0);
2314  break;
2315  case kTransformHartley:
2316  for (i = 0; i < fSize; i++) {
2317  working_space[i] = source[i];
2318  }
2319  Fourier(working_space, fSize, 1, kTransformForward, 0);
2320  break;
2321  case kTransformFourierWalsh:
2322  case kTransformFourierHaar:
2323  case kTransformWalshHaar:
2324  case kTransformCosWalsh:
2325  case kTransformCosHaar:
2326  case kTransformSinWalsh:
2327  case kTransformSinHaar:
2328  for (i = 0; i < fSize; i++) {
2329  val = source[i];
2330  if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
2331  j = (Int_t) TMath::Power(2, fDegree) / 2;
2332  k = i / j;
2333  k = 2 * k * j;
2334  working_space[k + i % j] = val;
2335  working_space[k + 2 * j - 1 - i % j] = val;
2336  }
2337 
2338  else if (fTransformType == kTransformSinWalsh
2339  || fTransformType == kTransformSinHaar) {
2340  j = (Int_t) TMath::Power(2, fDegree) / 2;
2341  k = i / j;
2342  k = 2 * k * j;
2343  working_space[k + i % j] = val;
2344  working_space[k + 2 * j - 1 - i % j] = -val;
2345  }
2346 
2347  else
2348  working_space[i] = val;
2349  }
2350  if (fTransformType == kTransformFourierWalsh
2351  || fTransformType == kTransformFourierHaar
2352  || fTransformType == kTransformWalshHaar) {
2353  for (i = 0; i < j; i++)
2354  BitReverseHaar(working_space, fSize, k, i * k);
2355  GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
2356  }
2357 
2358  else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
2359  m = (Int_t) TMath::Power(2, fDegree);
2360  l = 2 * fSize / m;
2361  for (i = 0; i < l; i++)
2362  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2363  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
2364  for (i = 0; i < fSize; i++) {
2365  k = i / j;
2366  k = 2 * k * j;
2367  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2368  a = TMath::Cos(a);
2369  b = working_space[k + i % j];
2370  if (i % j == 0)
2371  a = b / TMath::Sqrt(2.0);
2372 
2373  else
2374  a = b / a;
2375  working_space[i] = a;
2376  working_space[i + 2 * fSize] = 0;
2377  }
2378  }
2379 
2380  else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
2381  m = (Int_t) TMath::Power(2, fDegree);
2382  l = 2 * fSize / m;
2383  for (i = 0; i < l; i++)
2384  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2385  GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
2386  for (i = 0; i < fSize; i++) {
2387  k = i / j;
2388  k = 2 * k * j;
2389  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2390  a = TMath::Cos(a);
2391  b = working_space[j + k + i % j];
2392  if (i % j == 0)
2393  a = b / TMath::Sqrt(2.0);
2394 
2395  else
2396  a = b / a;
2397  working_space[j + k / 2 - i % j - 1] = a;
2398  working_space[i + 2 * fSize] = 0;
2399  }
2400  }
2401  if (fTransformType > kTransformWalshHaar)
2402  k = (Int_t) TMath::Power(2, fDegree - 1);
2403 
2404  else
2405  k = (Int_t) TMath::Power(2, fDegree);
2406  j = fSize / k;
2407  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
2408  working_space[fSize + i] = working_space[l + i / j];
2409  working_space[fSize + i + 2 * fSize] =
2410  working_space[l + i / j + 2 * fSize];
2411  }
2412  for (i = 0; i < fSize; i++) {
2413  working_space[i] = working_space[fSize + i];
2414  working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
2415  }
2416  break;
2417  }
2418  if (!working_space) return;
2419  for (i = 0, old_area = 0; i < fSize; i++) {
2420  old_area += working_space[i];
2421  }
2422  for (i = 0, new_area = 0; i < fSize; i++) {
2423  if (i >= fXmin && i <= fXmax)
2424  working_space[i] *= fEnhanceCoeff;
2425  new_area += working_space[i];
2426  }
2427  if (new_area != 0) {
2428  a = old_area / new_area;
2429  for (i = 0; i < fSize; i++) {
2430  working_space[i] *= a;
2431  }
2432  }
2433  if (fTransformType == kTransformFourier) {
2434  for (i = 0, old_area = 0; i < fSize; i++) {
2435  old_area += working_space[i + fSize];
2436  }
2437  for (i = 0, new_area = 0; i < fSize; i++) {
2438  if (i >= fXmin && i <= fXmax)
2439  working_space[i + fSize] *= fEnhanceCoeff;
2440  new_area += working_space[i + fSize];
2441  }
2442  if (new_area != 0) {
2443  a = old_area / new_area;
2444  for (i = 0; i < fSize; i++) {
2445  working_space[i + fSize] *= a;
2446  }
2447  }
2448  }
2449 
2450  else if (fTransformType == kTransformFourierWalsh
2451  || fTransformType == kTransformFourierHaar) {
2452  for (i = 0, old_area = 0; i < fSize; i++) {
2453  old_area += working_space[i + 2 * fSize];
2454  }
2455  for (i = 0, new_area = 0; i < fSize; i++) {
2456  if (i >= fXmin && i <= fXmax)
2457  working_space[i + 2 * fSize] *= fEnhanceCoeff;
2458  new_area += working_space[i + 2 * fSize];
2459  }
2460  if (new_area != 0) {
2461  a = old_area / new_area;
2462  for (i = 0; i < fSize; i++) {
2463  working_space[i + 2 * fSize] *= a;
2464  }
2465  }
2466  }
2467  switch (fTransformType) {
2468  case kTransformHaar:
2469  Haar(working_space, fSize, kTransformInverse);
2470  for (i = 0; i < fSize; i++) {
2471  destVector[i] = working_space[i];
2472  }
2473  break;
2474  case kTransformWalsh:
2475  BitReverse(working_space, fSize);
2476  Walsh(working_space, fSize);
2477  for (i = 0; i < fSize; i++) {
2478  destVector[i] = working_space[i];
2479  }
2480  break;
2481  case kTransformCos:
2482  fSize = 2 * fSize;
2483  working_space[0] = working_space[0] * TMath::Sqrt(2.0);
2484  for (i = 0; i < fSize / 2; i++) {
2485  a = pi * (Double_t) i / (Double_t) fSize;
2486  b = TMath::Sin(a);
2487  a = TMath::Cos(a);
2488  working_space[i + fSize] = (Double_t) working_space[i] * b;
2489  working_space[i] = (Double_t) working_space[i] * a;
2490  } for (i = 2; i <= (fSize / 2); i++) {
2491  working_space[fSize - i + 1] = working_space[i - 1];
2492  working_space[fSize - i + 1 + fSize] =
2493  -working_space[i - 1 + fSize];
2494  }
2495  working_space[fSize / 2] = 0;
2496  working_space[fSize / 2 + fSize] = 0;
2497  Fourier(working_space, fSize, 0, kTransformInverse, 1);
2498  for (i = 0; i < fSize / 2; i++) {
2499  destVector[i] = working_space[i];
2500  }
2501  break;
2502  case kTransformSin:
2503  fSize = 2 * fSize;
2504  working_space[fSize / 2] =
2505  working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
2506  for (i = fSize / 2 - 1; i > 0; i--) {
2507  a = pi * (Double_t) i / (Double_t) fSize;
2508  working_space[i + fSize] =
2509  -(Double_t) working_space[i - 1] * TMath::Cos(a);
2510  working_space[i] = (Double_t) working_space[i - 1] * TMath::Sin(a);
2511  } for (i = 2; i <= (fSize / 2); i++) {
2512  working_space[fSize - i + 1] = working_space[i - 1];
2513  working_space[fSize - i + 1 + fSize] =
2514  -working_space[i - 1 + fSize];
2515  }
2516  working_space[0] = 0;
2517  working_space[fSize] = 0;
2518  working_space[fSize / 2 + fSize] = 0;
2519  Fourier(working_space, fSize, 0, kTransformInverse, 0);
2520  for (i = 0; i < fSize / 2; i++) {
2521  destVector[i] = working_space[i];
2522  }
2523  break;
2524  case kTransformFourier:
2525  Fourier(working_space, fSize, 0, kTransformInverse, 0);
2526  for (i = 0; i < fSize; i++) {
2527  destVector[i] = working_space[i];
2528  }
2529  break;
2530  case kTransformHartley:
2531  Fourier(working_space, fSize, 1, kTransformInverse, 0);
2532  for (i = 0; i < fSize; i++) {
2533  destVector[i] = working_space[i];
2534  }
2535  break;
2536  case kTransformFourierWalsh:
2537  case kTransformFourierHaar:
2538  case kTransformWalshHaar:
2539  case kTransformCosWalsh:
2540  case kTransformCosHaar:
2541  case kTransformSinWalsh:
2542  case kTransformSinHaar:
2543  if (fTransformType > kTransformWalshHaar)
2544  k = (Int_t) TMath::Power(2, fDegree - 1);
2545 
2546  else
2547  k = (Int_t) TMath::Power(2, fDegree);
2548  j = fSize / k;
2549  for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
2550  working_space[fSize + l + i / j] = working_space[i];
2551  working_space[fSize + l + i / j + 2 * fSize] =
2552  working_space[i + 2 * fSize];
2553  }
2554  for (i = 0; i < fSize; i++) {
2555  working_space[i] = working_space[fSize + i];
2556  working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
2557  }
2558  if (fTransformType == kTransformFourierWalsh
2559  || fTransformType == kTransformFourierHaar
2560  || fTransformType == kTransformWalshHaar) {
2561  GeneralInv(working_space, fSize, fDegree, fTransformType);
2562  for (i = 0; i < j; i++)
2563  BitReverseHaar(working_space, fSize, k, i * k);
2564  }
2565 
2566  else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
2567  j = (Int_t) TMath::Power(2, fDegree) / 2;
2568  m = (Int_t) TMath::Power(2, fDegree);
2569  l = 2 * fSize / m;
2570  for (i = 0; i < fSize; i++) {
2571  k = i / j;
2572  k = 2 * k * j;
2573  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2574  if (i % j == 0) {
2575  working_space[2 * fSize + k + i % j] =
2576  working_space[i] * TMath::Sqrt(2.0);
2577  working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
2578  }
2579 
2580  else {
2581  b = TMath::Sin(a);
2582  a = TMath::Cos(a);
2583  working_space[4 * fSize + 2 * fSize + k + i % j] =
2584  -(Double_t) working_space[i] * b;
2585  working_space[2 * fSize + k + i % j] =
2586  (Double_t) working_space[i] * a;
2587  } } for (i = 0; i < fSize; i++) {
2588  k = i / j;
2589  k = 2 * k * j;
2590  if (i % j == 0) {
2591  working_space[2 * fSize + k + j] = 0;
2592  working_space[4 * fSize + 2 * fSize + k + j] = 0;
2593  }
2594 
2595  else {
2596  working_space[2 * fSize + k + 2 * j - i % j] =
2597  working_space[2 * fSize + k + i % j];
2598  working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
2599  -working_space[4 * fSize + 2 * fSize + k + i % j];
2600  }
2601  }
2602  for (i = 0; i < 2 * fSize; i++) {
2603  working_space[i] = working_space[2 * fSize + i];
2604  working_space[4 * fSize + i] =
2605  working_space[4 * fSize + 2 * fSize + i];
2606  }
2607  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2608  m = (Int_t) TMath::Power(2, fDegree);
2609  l = 2 * fSize / m;
2610  for (i = 0; i < l; i++)
2611  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2612  }
2613 
2614  else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
2615  j = (Int_t) TMath::Power(2, fDegree) / 2;
2616  m = (Int_t) TMath::Power(2, fDegree);
2617  l = 2 * fSize / m;
2618  for (i = 0; i < fSize; i++) {
2619  k = i / j;
2620  k = 2 * k * j;
2621  a = pi * (Double_t) (i % j) / (Double_t) (2 * j);
2622  if (i % j == 0) {
2623  working_space[2 * fSize + k + j + i % j] =
2624  working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
2625  working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
2626  }
2627 
2628  else {
2629  b = TMath::Sin(a);
2630  a = TMath::Cos(a);
2631  working_space[4 * fSize + 2 * fSize + k + j + i % j] =
2632  -(Double_t) working_space[j + k / 2 - i % j - 1] * b;
2633  working_space[2 * fSize + k + j + i % j] =
2634  (Double_t) working_space[j + k / 2 - i % j - 1] * a;
2635  } } for (i = 0; i < fSize; i++) {
2636  k = i / j;
2637  k = 2 * k * j;
2638  if (i % j == 0) {
2639  working_space[2 * fSize + k] = 0;
2640  working_space[4 * fSize + 2 * fSize + k] = 0;
2641  }
2642 
2643  else {
2644  working_space[2 * fSize + k + i % j] =
2645  working_space[2 * fSize + k + 2 * j - i % j];
2646  working_space[4 * fSize + 2 * fSize + k + i % j] =
2647  -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
2648  }
2649  }
2650  for (i = 0; i < 2 * fSize; i++) {
2651  working_space[i] = working_space[2 * fSize + i];
2652  working_space[4 * fSize + i] =
2653  working_space[4 * fSize + 2 * fSize + i];
2654  }
2655  GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
2656  for (i = 0; i < l; i++)
2657  BitReverseHaar(working_space, 2 * fSize, m, i * m);
2658  }
2659  for (i = 0; i < fSize; i++) {
2660  if (fTransformType >= kTransformCosWalsh) {
2661  k = i / j;
2662  k = 2 * k * j;
2663  val = working_space[k + i % j];
2664  }
2665 
2666  else
2667  val = working_space[i];
2668  destVector[i] = val;
2669  }
2670  break;
2671  }
2672  delete[]working_space;
2673  return;
2674 }
2675 
2676 ////////////////////////////////////////////////////////////////////////////////
2677 ///////////////////////////////////////////////////////////////////////////////
2678 /// SETTER FUNCION
2679 ///
2680 /// This function sets the following parameters for transform:
2681 /// -transType - type of transform (Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar)
2682 /// -degree - degree of mixed transform, applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar transforms
2683 ///////////////////////////////////////////////////////////////////////////////
2684 
2685 void TSpectrumTransform::SetTransformType(Int_t transType, Int_t degree)
2686 {
2687  Int_t j, n;
2688  j = 0;
2689  n = 1;
2690  for (; n < fSize;) {
2691  j += 1;
2692  n = n * 2;
2693  }
2694  if (transType < kTransformHaar || transType > kTransformSinHaar){
2695  Error ("TSpectrumTransform","Invalid type of transform");
2696  return;
2697  }
2698  if (transType >= kTransformFourierWalsh && transType <= kTransformSinHaar) {
2699  if (degree > j || degree < 1){
2700  Error ("TSpectrumTransform","Invalid degree of mixed transform");
2701  return;
2702  }
2703  }
2704  fTransformType=transType;
2705  fDegree=degree;
2706 }
2707 
2708 
2709 ////////////////////////////////////////////////////////////////////////////////
2710 ///////////////////////////////////////////////////////////////////////////////
2711 /// SETTER FUNCION
2712 ///
2713 /// This function sets the filtering or enhancement region:
2714 /// -xmin, xmax
2715 ///////////////////////////////////////////////////////////////////////////////
2716 
2717 void TSpectrumTransform::SetRegion(Int_t xmin, Int_t xmax)
2718 {
2719  if(xmin<0 || xmax < xmin || xmax >= fSize){
2720  Error("TSpectrumTransform", "Wrong range");
2721  return;
2722  }
2723  fXmin = xmin;
2724  fXmax = xmax;
2725 }
2726 
2727 ////////////////////////////////////////////////////////////////////////////////
2728 ///////////////////////////////////////////////////////////////////////////////
2729 /// SETTER FUNCION
2730 ///
2731 /// This function sets the direction of the transform:
2732 /// -direction (forward or inverse)
2733 ///////////////////////////////////////////////////////////////////////////////
2734 
2735 void TSpectrumTransform::SetDirection(Int_t direction)
2736 {
2737  if(direction != kTransformForward && direction != kTransformInverse){
2738  Error("TSpectrumTransform", "Wrong direction");
2739  return;
2740  }
2741  fDirection = direction;
2742 }
2743 
2744 ////////////////////////////////////////////////////////////////////////////////
2745 ///////////////////////////////////////////////////////////////////////////////
2746 /// SETTER FUNCION
2747 ///
2748 /// This function sets the filter coefficient:
2749 /// -filterCoeff - after the transform the filtered region (xmin, xmax) is replaced by this coefficient. Applies only for filtereng operation.
2750 ///////////////////////////////////////////////////////////////////////////////
2751 
2752 void TSpectrumTransform::SetFilterCoeff(Double_t filterCoeff)
2753 {
2754  fFilterCoeff = filterCoeff;
2755 }
2756 
2757 ////////////////////////////////////////////////////////////////////////////////
2758 ///////////////////////////////////////////////////////////////////////////////
2759 /// SETTER FUNCION
2760 ///
2761 /// This function sets the enhancement coefficient:
2762 /// -enhanceCoeff - after the transform the enhanced region (xmin, xmax) is multiplied by this coefficient. Applies only for enhancement operation.
2763 ///////////////////////////////////////////////////////////////////////////////
2764 
2765 void TSpectrumTransform::SetEnhanceCoeff(Double_t enhanceCoeff)
2766 {
2767  fEnhanceCoeff = enhanceCoeff;
2768 }
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:2685
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
sample_config.c
c
Definition: sample_config.py:3
sample_config.ip
ip
Definition: sample_config.py:8
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
TSpectrumTransform::Transform
void Transform(const Double_t *source, Double_t *destVector)
ONE-DIMENSIONAL TRANSFORM FUNCTION This function transforms the source spectrum.
Definition: TSpectrumTransform.cxx:743
TSpectrumTransform::SetEnhanceCoeff
void SetEnhanceCoeff(Double_t enhanceCoeff)
SETTER FUNCION.
Definition: TSpectrumTransform.cxx:2765
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:1484
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:2717
m
TMarker * m
Definition: textangle.C:8
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
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
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
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:279
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:2752
TSpectrumTransform::TSpectrumTransform
TSpectrumTransform()
default constructor
Definition: TSpectrumTransform.cxx:68
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:2091
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:2735
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