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