ROOT  6.06/08
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 /*
1752 <div class=Section1>
1753 
1754 <p class=MsoNormal><b><span style='font-size:14.0pt'>Transform methods</span></b></p>
1755 
1756 <p class=MsoNormal style='text-align:justify'><i>&nbsp;</i></p>
1757 
1758 <p class=MsoNormal style='text-align:justify'><i>Goal: to analyze experimental
1759 data using orthogonal transforms</i></p>
1760 
1761 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
1762 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
1763 </span>orthogonal transforms can be successfully used for the processing of
1764 nuclear spectra (not only) </p>
1765 
1766 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
1767 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
1768 </span>they can be used to remove high frequency noise, to increase
1769 signal-to-background ratio as well as to enhance low intensity components [1],
1770 to carry out e.g. Fourier analysis etc. </p>
1771 
1772 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
1773 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
1774 </span>we have implemented the function for the calculation of the commonly
1775 used orthogonal transforms as well as functions for the filtration and
1776 enhancement of experimental data</p>
1777 
1778 <p class=MsoNormal><i>&nbsp;</i></p>
1779 
1780 <p class=MsoNormal><i>Function:</i></p>
1781 
1782 <p class=MsoNormal>void <a
1783 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumTransform2::Transform</b></a><b>(const
1784 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> **fSource,
1785 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> **fDest)</b></p>
1786 
1787 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
1788 
1789 <p class=MsoNormal style='text-align:justify'>This function transforms the
1790 source spectrum according to the given input parameters. Transformed data are
1791 written into dest spectrum. Before the Transform function is called the class
1792 must be created by constructor and the type of the transform as well as some
1793 other parameters must be set using a set of setter functions:</p>
1794 
1795 <p class=MsoNormal><span lang=FR>&nbsp;</span></p>
1796 
1797 <p class=MsoNormal><i><span style='color:red'>Member variables of
1798 TSpectrumTransform2 class:</span></i></p>
1799 
1800 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'><b>fSource</b>-pointer
1801 to the matrix of source spectrum. Its lengths should be equal to the “fSizex,
1802 fSizey” parameters except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR
1803 transforms. These need “2*fSizex*fSizey” length to supply real and imaginary
1804 coefficients.                   </p>
1805 
1806 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'><b>fDest</b>-pointer
1807 to the matrix of destination spectrum. Its lengths should be equal to the
1808 “fSizex, fSizey” parameters except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR
1809 transforms. These need “2*fSizex*fSizey” length to store real and imaginary
1810 coefficients. </p>
1811 
1812 <p class=MsoNormal style='text-align:justify'>        <b>fSizeX,fSizeY</b>-basic
1813 lengths of the source and dest spectra. They<span style='color:fuchsia'> should
1814 be power  </span></p>
1815 
1816 <p class=MsoNormal style='text-align:justify'><span style='color:fuchsia'>     
1817 of 2.</span></p>
1818 
1819 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify;text-indent:
1820 -2.85pt'><b>fType</b>-type of transform</p>
1821 
1822 <p class=MsoNormal style='text-align:justify'>            Classic transforms:</p>
1823 
1824 <p class=MsoNormal style='text-align:justify'>                        kTransformHaar
1825 </p>
1826 
1827 <p class=MsoNormal style='text-align:justify'>                        kTransformWalsh
1828 </p>
1829 
1830 <p class=MsoNormal style='text-align:justify'>                        kTransformCos
1831 </p>
1832 
1833 <p class=MsoNormal style='text-align:justify'>                        kTransformSin
1834 </p>
1835 
1836 <p class=MsoNormal style='text-align:justify'>                        kTransformFourier
1837 </p>
1838 
1839 <p class=MsoNormal style='text-align:justify'>                        kTransformHartley
1840 </p>
1841 
1842 <p class=MsoNormal style='text-align:justify'>            Mixed transforms:</p>
1843 
1844 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierWalsh
1845 </p>
1846 
1847 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierHaar
1848 </p>
1849 
1850 <p class=MsoNormal style='text-align:justify'>                        kTransformWalshHaar
1851 </p>
1852 
1853 <p class=MsoNormal style='text-align:justify'>                        kTransformCosWalsh
1854 </p>
1855 
1856 <p class=MsoNormal style='text-align:justify'>                        kTransformCosHaar
1857 </p>
1858 
1859 <p class=MsoNormal style='text-align:justify'>                        kTransformSinWalsh
1860 </p>
1861 
1862 <p class=MsoNormal style='text-align:justify'>                        kTransformSinHaar
1863 </p>
1864 
1865 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDirection</b>-direction-transform
1866 direction (forward, inverse)</p>
1867 
1868 <p class=MsoNormal style='text-align:justify'>                        kTransformForward
1869 </p>
1870 
1871 <p class=MsoNormal style='text-align:justify'>                        kTransformInverse
1872 </p>
1873 
1874 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDegree</b>-applies
1875 only for mixed transforms [2], [3], [4]. </p>
1876 
1877 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'>                
1878 <span style='color:fuchsia'> Allowed range  <sub><img border=0 width=100
1879 height=27 src="gif/spectrum2transform_transform_image001.gif"></sub>. </span></p>
1880 
1881 <p class=MsoNormal style='text-align:justify'><b><i>References:</i></b></p>
1882 
1883 <p class=MsoNormal style='text-align:justify'>[1] C.V. Hampton, B. Lian, Wm. C.
1884 McHarris: Fast-Fourier-transform spectral enhancement techniques for gamma-ray
1885 spectroscopy. NIM A353 (1994) 280-284. </p>
1886 
1887 <p class=MsoNormal style='text-align:justify'>[2] Morhá&#269; M., Matoušek V.,
1888 New adaptive Cosine-Walsh  transform and its application to nuclear data
1889 compression, IEEE Transactions on Signal Processing 48 (2000) 2693.  </p>
1890 
1891 <p class=MsoNormal style='text-align:justify'>[3] Morhá&#269; M., Matoušek V.,
1892 Data compression using new fast adaptive Cosine-Haar transforms, Digital Signal
1893 Processing 8 (1998) 63. </p>
1894 
1895 <p class=MsoNormal style='text-align:justify'>[4] Morhá&#269; M., Matoušek V.:
1896 Multidimensional nuclear data compression using fast adaptive Walsh-Haar
1897 transform. Acta Physica Slovaca 51 (2001) 307. </p>
1898 
1899 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
1900 
1901 <p class=MsoNormal style='text-align:justify'><i>Example 1 – script Transform2.c:</i></p>
1902 
1903 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><img
1904 border=0 width=602 height=455 src="gif/spectrum2transform_transform_image002.jpg"></span></p>
1905 
1906 <p class=MsoNormal><b>Fig. 1 Original two-dimensional noisy spectrum</b></p>
1907 
1908 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><img
1909 border=0 width=602 height=455 src="gif/spectrum2transform_transform_image003.jpg"></span></p>
1910 
1911 <p class=MsoNormal style='text-align:justify'><b>Fig. 2 Transformed spectrum
1912 from Fig. 1 using Cosine transform. Energy of the trasnsformed data is
1913 concentrated around the beginning of the coordinate system</b></p>
1914 
1915 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>&nbsp;</span></b></p>
1916 
1917 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
1918 
1919 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
1920 Transform function (class TSpectrumTransform2).</span></p>
1921 
1922 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
1923 do</span></p>
1924 
1925 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Transform2.C</span></p>
1926 
1927 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Transform2() {</span></p>
1928 
1929 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i, j;</span></p>
1930 
1931 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbinsx =
1932 256;</span></p>
1933 
1934 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbinsy =
1935 256;   </span></p>
1936 
1937 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmin  = 0;</span></p>
1938 
1939 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
1940 nbinsx;</span></p>
1941 
1942 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymin  = 0;</span></p>
1943 
1944 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymax  = nbinsy;</span></p>
1945 
1946 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
1947 style='font-size:10.0pt'>Double_t ** source = new Double_t *[nbinsx];   </span></p>
1948 
1949 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t ** dest = new
1950 Double_t *[nbinsx];      </span></p>
1951 
1952 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i=0;i&lt;nbinsx;i++)</span></p>
1953 
1954 <p class=MsoNormal><span style='font-size:10.0pt'>                                                source[i]=new
1955 Double_t[nbinsy];</span></p>
1956 
1957 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i=0;i&lt;nbinsx;i++)</span></p>
1958 
1959 <p class=MsoNormal><span style='font-size:10.0pt'>                                                dest[i]=new
1960 Double_t[nbinsy];   </span></p>
1961 
1962 <p class=MsoNormal><span style='font-size:10.0pt'>   TH2F *trans = new
1963 TH2F(&quot;trans&quot;,&quot;Background
1964 estimation&quot;,nbinsx,xmin,xmax,nbinsy,ymin,ymax);</span></p>
1965 
1966 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
1967 TFile(&quot;TSpectrum2.root&quot;);</span></p>
1968 
1969 <p class=MsoNormal><span style='font-size:10.0pt'>   trans=(TH2F*)
1970 f-&gt;Get(&quot;back3;1&quot;);</span></p>
1971 
1972 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Tr = new
1973 TCanvas(&quot;Transform&quot;,&quot;Illustation of transform
1974 function&quot;,10,10,1000,700);</span></p>
1975 
1976 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbinsx;
1977 i++){</span></p>
1978 
1979 <p class=MsoNormal><span style='font-size:10.0pt'>     for (j = 0; j &lt;
1980 nbinsy; j++){</span></p>
1981 
1982 <p class=MsoNormal><span style='font-size:10.0pt'>                    source[i][j]
1983 = trans-&gt;GetBinContent(i + 1,j + 1); </span></p>
1984 
1985 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                 }</span></p>
1986 
1987 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }           </span></p>
1988 
1989 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1990 TSpectrumTransform2 *t = new TSpectrumTransform2(256,256);   </span></p>
1991 
1992 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1993 t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
1994 
1995 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   t-&gt;SetDirection(t-&gt;kTransformForward);</span></p>
1996 
1997 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1998 t-&gt;Transform(source,dest);</span></p>
1999 
2000 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2001 style='font-size:10.0pt'>for (i = 0; i &lt; nbinsx; i++){</span></p>
2002 
2003 <p class=MsoNormal><span style='font-size:10.0pt'>     for (j = 0; j &lt;
2004 nbinsy; j++){</span></p>
2005 
2006 <p class=MsoNormal><span style='font-size:10.0pt'>                  trans-&gt;SetBinContent(i
2007 + 1, j + 1,dest[i][j]);</span></p>
2008 
2009 <p class=MsoNormal><span style='font-size:10.0pt'>                 }</span></p>
2010 
2011 <p class=MsoNormal><span style='font-size:10.0pt'>   }   </span></p>
2012 
2013 <p class=MsoNormal><span style='font-size:10.0pt'>  
2014 trans-&gt;Draw(&quot;SURF&quot;);      </span></p>
2015 
2016 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
2017 
2018 </div>
2019 
2020  */
2021  Int_t i, j;
2022  Int_t size;
2023  Double_t *working_vector = 0, **working_matrix = 0;
2024  size = (Int_t) TMath::Max(fSizeX, fSizeY);
2025  switch (fTransformType) {
2026  case kTransformHaar:
2027  case kTransformWalsh:
2028  working_vector = new Double_t[2 * size];
2029  working_matrix = new Double_t *[fSizeX];
2030  for (i = 0; i < fSizeX; i++)
2031  working_matrix[i] = new Double_t[fSizeY];
2032  break;
2033  case kTransformCos:
2034  case kTransformSin:
2035  case kTransformFourier:
2036  case kTransformHartley:
2037  case kTransformFourierWalsh:
2038  case kTransformFourierHaar:
2039  case kTransformWalshHaar:
2040  working_vector = new Double_t[4 * size];
2041  working_matrix = new Double_t *[fSizeX];
2042  for (i = 0; i < fSizeX; i++)
2043  working_matrix[i] = new Double_t[2 * fSizeY];
2044  break;
2045  case kTransformCosWalsh:
2046  case kTransformCosHaar:
2047  case kTransformSinWalsh:
2048  case kTransformSinHaar:
2049  working_vector = new Double_t[8 * size];
2050  working_matrix = new Double_t *[fSizeX];
2051  for (i = 0; i < fSizeX; i++)
2052  working_matrix[i] = new Double_t[2 * fSizeY];
2053  break;
2054  }
2055  if (fDirection == kTransformForward) {
2056  switch (fTransformType) {
2057  case kTransformHaar:
2058  for (i = 0; i < fSizeX; i++) {
2059  for (j = 0; j < fSizeY; j++) {
2060  working_matrix[i][j] = fSource[i][j];
2061  }
2062  }
2063  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2064  fDirection, kTransformHaar);
2065  for (i = 0; i < fSizeX; i++) {
2066  for (j = 0; j < fSizeY; j++) {
2067  fDest[i][j] = working_matrix[i][j];
2068  }
2069  }
2070  break;
2071  case kTransformWalsh:
2072  for (i = 0; i < fSizeX; i++) {
2073  for (j = 0; j < fSizeY; j++) {
2074  working_matrix[i][j] = fSource[i][j];
2075  }
2076  }
2077  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2078  fDirection, kTransformWalsh);
2079  for (i = 0; i < fSizeX; i++) {
2080  for (j = 0; j < fSizeY; j++) {
2081  fDest[i][j] = working_matrix[i][j];
2082  }
2083  }
2084  break;
2085  case kTransformCos:
2086  for (i = 0; i < fSizeX; i++) {
2087  for (j = 0; j < fSizeY; j++) {
2088  working_matrix[i][j] = fSource[i][j];
2089  }
2090  }
2091  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2092  kTransformCos);
2093  for (i = 0; i < fSizeX; i++) {
2094  for (j = 0; j < fSizeY; j++) {
2095  fDest[i][j] = working_matrix[i][j];
2096  }
2097  }
2098  break;
2099  case kTransformSin:
2100  for (i = 0; i < fSizeX; i++) {
2101  for (j = 0; j < fSizeY; j++) {
2102  working_matrix[i][j] = fSource[i][j];
2103  }
2104  }
2105  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2106  kTransformSin);
2107  for (i = 0; i < fSizeX; i++) {
2108  for (j = 0; j < fSizeY; j++) {
2109  fDest[i][j] = working_matrix[i][j];
2110  }
2111  }
2112  break;
2113  case kTransformFourier:
2114  for (i = 0; i < fSizeX; i++) {
2115  for (j = 0; j < fSizeY; j++) {
2116  working_matrix[i][j] = fSource[i][j];
2117  }
2118  }
2119  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2120  kTransformFourier);
2121  for (i = 0; i < fSizeX; i++) {
2122  for (j = 0; j < fSizeY; j++) {
2123  fDest[i][j] = working_matrix[i][j];
2124  }
2125  }
2126  for (i = 0; i < fSizeX; i++) {
2127  for (j = 0; j < fSizeY; j++) {
2128  fDest[i][j + fSizeY] = working_matrix[i][j + fSizeY];
2129  }
2130  }
2131  break;
2132  case kTransformHartley:
2133  for (i = 0; i < fSizeX; i++) {
2134  for (j = 0; j < fSizeY; j++) {
2135  working_matrix[i][j] = fSource[i][j];
2136  }
2137  }
2138  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2139  kTransformHartley);
2140  for (i = 0; i < fSizeX; i++) {
2141  for (j = 0; j < fSizeY; j++) {
2142  fDest[i][j] = working_matrix[i][j];
2143  }
2144  }
2145  break;
2146  case kTransformFourierWalsh:
2147  case kTransformFourierHaar:
2148  case kTransformWalshHaar:
2149  case kTransformCosWalsh:
2150  case kTransformCosHaar:
2151  case kTransformSinWalsh:
2152  case kTransformSinHaar:
2153  for (i = 0; i < fSizeX; i++) {
2154  for (j = 0; j < fSizeY; j++) {
2155  working_matrix[i][j] = fSource[i][j];
2156  }
2157  }
2158  General2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2159  fTransformType, fDegree);
2160  for (i = 0; i < fSizeX; i++) {
2161  for (j = 0; j < fSizeY; j++) {
2162  fDest[i][j] = working_matrix[i][j];
2163  }
2164  }
2165  if (fTransformType == kTransformFourierWalsh
2166  || fTransformType == kTransformFourierHaar) {
2167  for (i = 0; i < fSizeX; i++) {
2168  for (j = 0; j < fSizeY; j++) {
2169  fDest[i][j + fSizeY] = working_matrix[i][j + fSizeY];
2170  }
2171  }
2172  }
2173  break;
2174  }
2175  }
2176 
2177  else if (fDirection == kTransformInverse) {
2178  switch (fTransformType) {
2179  case kTransformHaar:
2180  for (i = 0; i < fSizeX; i++) {
2181  for (j = 0; j < fSizeY; j++) {
2182  working_matrix[i][j] = fSource[i][j];
2183  }
2184  }
2185  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2186  fDirection, kTransformHaar);
2187  for (i = 0; i < fSizeX; i++) {
2188  for (j = 0; j < fSizeY; j++) {
2189  fDest[i][j] = working_matrix[i][j];
2190  }
2191  }
2192  break;
2193  case kTransformWalsh:
2194  for (i = 0; i < fSizeX; i++) {
2195  for (j = 0; j < fSizeY; j++) {
2196  working_matrix[i][j] = fSource[i][j];
2197  }
2198  }
2199  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2200  fDirection, kTransformWalsh);
2201  for (i = 0; i < fSizeX; i++) {
2202  for (j = 0; j < fSizeY; j++) {
2203  fDest[i][j] = working_matrix[i][j];
2204  }
2205  }
2206  break;
2207  case kTransformCos:
2208  for (i = 0; i < fSizeX; i++) {
2209  for (j = 0; j < fSizeY; j++) {
2210  working_matrix[i][j] = fSource[i][j];
2211  }
2212  }
2213  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2214  kTransformCos);
2215  for (i = 0; i < fSizeX; i++) {
2216  for (j = 0; j < fSizeY; j++) {
2217  fDest[i][j] = working_matrix[i][j];
2218  }
2219  }
2220  break;
2221  case kTransformSin:
2222  for (i = 0; i < fSizeX; i++) {
2223  for (j = 0; j < fSizeY; j++) {
2224  working_matrix[i][j] = fSource[i][j];
2225  }
2226  }
2227  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2228  kTransformSin);
2229  for (i = 0; i < fSizeX; i++) {
2230  for (j = 0; j < fSizeY; j++) {
2231  fDest[i][j] = working_matrix[i][j];
2232  }
2233  }
2234  break;
2235  case kTransformFourier:
2236  for (i = 0; i < fSizeX; i++) {
2237  for (j = 0; j < fSizeY; j++) {
2238  working_matrix[i][j] = fSource[i][j];
2239  }
2240  }
2241  for (i = 0; i < fSizeX; i++) {
2242  for (j = 0; j < fSizeY; j++) {
2243  working_matrix[i][j + fSizeY] = fSource[i][j + fSizeY];
2244  }
2245  }
2246  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2247  kTransformFourier);
2248  for (i = 0; i < fSizeX; i++) {
2249  for (j = 0; j < fSizeY; j++) {
2250  fDest[i][j] = working_matrix[i][j];
2251  }
2252  }
2253  break;
2254  case kTransformHartley:
2255  for (i = 0; i < fSizeX; i++) {
2256  for (j = 0; j < fSizeY; j++) {
2257  working_matrix[i][j] = fSource[i][j];
2258  }
2259  }
2260  FourCos2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2261  kTransformHartley);
2262  for (i = 0; i < fSizeX; i++) {
2263  for (j = 0; j < fSizeY; j++) {
2264  fDest[i][j] = working_matrix[i][j];
2265  }
2266  }
2267  break;
2268  case kTransformFourierWalsh:
2269  case kTransformFourierHaar:
2270  case kTransformWalshHaar:
2271  case kTransformCosWalsh:
2272  case kTransformCosHaar:
2273  case kTransformSinWalsh:
2274  case kTransformSinHaar:
2275  for (i = 0; i < fSizeX; i++) {
2276  for (j = 0; j < fSizeY; j++) {
2277  working_matrix[i][j] = fSource[i][j];
2278  }
2279  }
2280  if (fTransformType == kTransformFourierWalsh
2281  || fTransformType == kTransformFourierHaar) {
2282  for (i = 0; i < fSizeX; i++) {
2283  for (j = 0; j < fSizeY; j++) {
2284  working_matrix[i][j + fSizeY] = fSource[i][j + fSizeY];
2285  }
2286  }
2287  }
2288  General2(working_matrix, working_vector, fSizeX, fSizeY, fDirection,
2289  fTransformType, fDegree);
2290  for (i = 0; i < fSizeX; i++) {
2291  for (j = 0; j < fSizeY; j++) {
2292  fDest[i][j] = working_matrix[i][j];
2293  }
2294  }
2295  break;
2296  }
2297  }
2298  for (i = 0; i < fSizeX; i++) {
2299  if (working_matrix) delete[]working_matrix[i];
2300  }
2301  delete[]working_matrix;
2302  delete[]working_vector;
2303  return;
2304 }
2305 //////////END OF TRANSFORM2 FUNCTION/////////////////////////////////
2306 //_______________________________________________________________________________________
2307 //////////FILTER2_ZONAL FUNCTION - CALCULATES DIFFERENT 2-D ORTHOGONAL TRANSFORMS, SETS GIVEN REGION TO FILTER COEFFICIENT AND TRANSFORMS IT BACK//////
2308 void TSpectrum2Transform::FilterZonal(const Double_t **fSource, Double_t **fDest)
2309 {
2310 //////////////////////////////////////////////////////////////////////////////////////////
2311 /* TWO-DIMENSIONAL FILTER ZONAL FUNCTION */
2312 /* This function transforms the source spectrum. The calling program */
2313 /* should fill in input parameters. Then it sets transformed */
2314 /* coefficients in the given region to the given */
2315 /* filter_coeff and transforms it back */
2316 /* Filtered data are written into dest spectrum. */
2317 /* */
2318 /* Function parameters: */
2319 /* fSource-pointer to the matrix of source spectrum, its size should */
2320 /* be fSizeX*fSizeY */
2321 /* fDest-pointer to the matrix of destination data, its size should be */
2322 /* fSizeX*fSizeY */
2323 /* */
2324 //////////////////////////////////////////////////////////////////////////////////////////
2325 /*
2326 <div class=Section2>
2327 
2328 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of zonal filtering</span></b></p>
2329 
2330 <p class=MsoNormal><i>&nbsp;</i></p>
2331 
2332 <p class=MsoNormal><i>Function:</i></p>
2333 
2334 <p class=MsoNormal>void <a
2335 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumTransform2::FilterZonal</b></a><b>(const
2336 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> **fSource,
2337 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> **fDest)</b></p>
2338 
2339 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
2340 
2341 <p class=MsoNormal style='text-align:justify'>This function transforms the
2342 source spectrum (for details see Transform function).  Before the FilterZonal
2343 function is called the class must be created by constructor and the type of the
2344 transform as well as some other parameters must be set using a set of setter
2345 functions. The FilterZonal function sets transformed coefficients in the given
2346 region (fXmin, fXmax) to the given fFilterCoeff and transforms it back. Filtered
2347 data are written into dest spectrum. </p>
2348 
2349 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
2350 
2351 <p class=MsoNormal style='text-align:justify'><i>Example  – script Fitler2.c:</i></p>
2352 
2353 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><img
2354 border=0 width=602 height=455 src="gif/spectrum2transform_filter_image001.jpg"></span></p>
2355 
2356 <p class=MsoNormal><b>Fig. 1 Original two-dimensional noisy spectrum</b></p>
2357 
2358 <p class=MsoNormal><b><span style='font-size:14.0pt'><img border=0 width=602
2359 height=455 src="gif/spectrum2transform_filter_image002.jpg"></span></b></p>
2360 
2361 <p class=MsoNormal style='text-align:justify'><b>Fig. 2 Filtered spectrum using
2362 Cosine transform and zonal filtration (channels in regions (128-255)x(0-255)
2363 and (0-255)x(128-255) were set to 0).  </b></p>
2364 
2365 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
2366 
2367 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
2368 
2369 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
2370 zonal filtration (class TSpectrumTransform2).</span></p>
2371 
2372 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
2373 do</span></p>
2374 
2375 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Filter2.C</span></p>
2376 
2377 <p class=MsoNormal><span style='font-size:10.0pt'> </span></p>
2378 
2379 <p class=MsoNormal><span style='font-size:10.0pt'>void Filter2() {</span></p>
2380 
2381 <p class=MsoNormal><span style='font-size:10.0pt'>   Int_t i, j;</span></p>
2382 
2383 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
2384 style='font-size:10.0pt'>Int_t nbinsx = 256;</span></p>
2385 
2386 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbinsy =
2387 256;   </span></p>
2388 
2389 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmin  = 0;</span></p>
2390 
2391 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
2392 nbinsx;</span></p>
2393 
2394 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymin  = 0;</span></p>
2395 
2396 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymax  =
2397 nbinsy;</span></p>
2398 
2399 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2400 style='font-size:10.0pt'>Double_t ** source = new Double_t *[nbinsx];   </span></p>
2401 
2402 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t ** dest = new
2403 Double_t *[nbinsx];      </span></p>
2404 
2405 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i=0;i&lt;nbinsx;i++)</span></p>
2406 
2407 <p class=MsoNormal><span style='font-size:10.0pt'>                                                source[i]=new
2408 Double_t[nbinsy];</span></p>
2409 
2410 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i=0;i&lt;nbinsx;i++)</span></p>
2411 
2412 <p class=MsoNormal><span style='font-size:10.0pt'>                                                dest[i]=new
2413 Double_t[nbinsy];   </span></p>
2414 
2415 <p class=MsoNormal><span style='font-size:10.0pt'>   TH2F *trans = new
2416 TH2F(&quot;trans&quot;,&quot;Background
2417 estimation&quot;,nbinsx,xmin,xmax,nbinsy,ymin,ymax);</span></p>
2418 
2419 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
2420 TFile(&quot;TSpectrum2.root&quot;);</span></p>
2421 
2422 <p class=MsoNormal><span style='font-size:10.0pt'>   trans=(TH2F*)
2423 f-&gt;Get(&quot;back3;1&quot;);</span></p>
2424 
2425 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Tr = new
2426 TCanvas(&quot;Transform&quot;,&quot;Illustation of transform
2427 function&quot;,10,10,1000,700);</span></p>
2428 
2429 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbinsx;
2430 i++){</span></p>
2431 
2432 <p class=MsoNormal><span style='font-size:10.0pt'>     for (j = 0; j &lt;
2433 nbinsy; j++){</span></p>
2434 
2435 <p class=MsoNormal><span style='font-size:10.0pt'>                    </span><span
2436 lang=FR style='font-size:10.0pt'>source[i][j] = trans-&gt;GetBinContent(i + 1,j
2437 + 1); </span></p>
2438 
2439 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                 }</span></p>
2440 
2441 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }           </span></p>
2442 
2443 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2444 TSpectrumTransform2 *t = new TSpectrumTransform2(256,256);   </span></p>
2445 
2446 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   t-&gt;SetTransformType(t-&gt;kTransformCos,0);  
2447 </span></p>
2448 
2449 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2450 t-&gt;SetRegion(0,255,128,255);</span></p>
2451 
2452 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2453 t-&gt;FilterZonal(source,dest);     </span></p>
2454 
2455 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2456 style='font-size:10.0pt'>for (i = 0; i &lt; nbinsx; i++){</span></p>
2457 
2458 <p class=MsoNormal><span style='font-size:10.0pt'>     for (j = 0; j &lt;
2459 nbinsy; j++){</span></p>
2460 
2461 <p class=MsoNormal><span style='font-size:10.0pt'>                    </span><span
2462 lang=FR style='font-size:10.0pt'>source[i][j] = dest[i][j]; </span></p>
2463 
2464 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                 }</span></p>
2465 
2466 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }   </span></p>
2467 
2468 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2469 t-&gt;SetRegion(128,255,0,255);</span></p>
2470 
2471 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   t-&gt;FilterZonal(source,dest);       
2472 </span></p>
2473 
2474 <p class=MsoNormal><span style='font-size:10.0pt'>  
2475 trans-&gt;Draw(&quot;SURF&quot;);     </span></p>
2476 
2477 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
2478 
2479 </div>
2480 
2481  */
2482 
2483  Int_t i, j;
2484  Double_t a, old_area = 0, new_area = 0;
2485  Int_t size;
2486  Double_t *working_vector = 0, **working_matrix = 0;
2487  size = (Int_t) TMath::Max(fSizeX, fSizeY);
2488  switch (fTransformType) {
2489  case kTransformHaar:
2490  case kTransformWalsh:
2491  working_vector = new Double_t[2 * size];
2492  working_matrix = new Double_t *[fSizeX];
2493  for (i = 0; i < fSizeX; i++)
2494  working_matrix[i] = new Double_t[fSizeY];
2495  break;
2496  case kTransformCos:
2497  case kTransformSin:
2498  case kTransformFourier:
2499  case kTransformHartley:
2500  case kTransformFourierWalsh:
2501  case kTransformFourierHaar:
2502  case kTransformWalshHaar:
2503  working_vector = new Double_t[4 * size];
2504  working_matrix = new Double_t *[fSizeX];
2505  for (i = 0; i < fSizeX; i++)
2506  working_matrix[i] = new Double_t[2 * fSizeY];
2507  break;
2508  case kTransformCosWalsh:
2509  case kTransformCosHaar:
2510  case kTransformSinWalsh:
2511  case kTransformSinHaar:
2512  working_vector = new Double_t[8 * size];
2513  working_matrix = new Double_t *[fSizeX];
2514  for (i = 0; i < fSizeX; i++)
2515  working_matrix[i] = new Double_t[2 * fSizeY];
2516  break;
2517  }
2518  switch (fTransformType) {
2519  case kTransformHaar:
2520  for (i = 0; i < fSizeX; i++) {
2521  for (j = 0; j < fSizeY; j++) {
2522  working_matrix[i][j] = fSource[i][j];
2523  old_area = old_area + fSource[i][j];
2524  }
2525  }
2526  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2527  kTransformForward, kTransformHaar);
2528  break;
2529  case kTransformWalsh:
2530  for (i = 0; i < fSizeX; i++) {
2531  for (j = 0; j < fSizeY; j++) {
2532  working_matrix[i][j] = fSource[i][j];
2533  old_area = old_area + fSource[i][j];
2534  }
2535  }
2536  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2537  kTransformForward, kTransformWalsh);
2538  break;
2539  case kTransformCos:
2540  for (i = 0; i < fSizeX; i++) {
2541  for (j = 0; j < fSizeY; j++) {
2542  working_matrix[i][j] = fSource[i][j];
2543  old_area = old_area + fSource[i][j];
2544  }
2545  }
2546  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2547  kTransformForward, kTransformCos);
2548  break;
2549  case kTransformSin:
2550  for (i = 0; i < fSizeX; i++) {
2551  for (j = 0; j < fSizeY; j++) {
2552  working_matrix[i][j] = fSource[i][j];
2553  old_area = old_area + fSource[i][j];
2554  }
2555  }
2556  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2557  kTransformForward, kTransformSin);
2558  break;
2559  case kTransformFourier:
2560  for (i = 0; i < fSizeX; i++) {
2561  for (j = 0; j < fSizeY; j++) {
2562  working_matrix[i][j] = fSource[i][j];
2563  old_area = old_area + fSource[i][j];
2564  }
2565  }
2566  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2567  kTransformForward, kTransformFourier);
2568  break;
2569  case kTransformHartley:
2570  for (i = 0; i < fSizeX; i++) {
2571  for (j = 0; j < fSizeY; j++) {
2572  working_matrix[i][j] = fSource[i][j];
2573  old_area = old_area + fSource[i][j];
2574  }
2575  }
2576  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2577  kTransformForward, kTransformHartley);
2578  break;
2579  case kTransformFourierWalsh:
2580  case kTransformFourierHaar:
2581  case kTransformWalshHaar:
2582  case kTransformCosWalsh:
2583  case kTransformCosHaar:
2584  case kTransformSinWalsh:
2585  case kTransformSinHaar:
2586  for (i = 0; i < fSizeX; i++) {
2587  for (j = 0; j < fSizeY; j++) {
2588  working_matrix[i][j] = fSource[i][j];
2589  old_area = old_area + fSource[i][j];
2590  }
2591  }
2592  General2(working_matrix, working_vector, fSizeX, fSizeY,
2593  kTransformForward, fTransformType, fDegree);
2594  break;
2595  }
2596  for (i = 0; i < fSizeX; i++) {
2597  for (j = 0; j < fSizeY; j++) {
2598  if (i >= fXmin && i <= fXmax && j >= fYmin && j <= fYmax)
2599  if (working_matrix) working_matrix[i][j] = fFilterCoeff;
2600  }
2601  }
2602  if (fTransformType == kTransformFourier || fTransformType == kTransformFourierWalsh
2603  || fTransformType == kTransformFourierHaar) {
2604  for (i = 0; i < fSizeX; i++) {
2605  for (j = 0; j < fSizeY; j++) {
2606  if (i >= fXmin && i <= fXmax && j >= fYmin && j <= fYmax)
2607  if (working_matrix) working_matrix[i][j + fSizeY] = fFilterCoeff;
2608  }
2609  }
2610  }
2611  switch (fTransformType) {
2612  case kTransformHaar:
2613  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2614  kTransformInverse, kTransformHaar);
2615  for (i = 0; i < fSizeX; i++) {
2616  for (j = 0; j < fSizeY; j++) {
2617  new_area = new_area + working_matrix[i][j];
2618  }
2619  }
2620  if (new_area != 0) {
2621  a = old_area / new_area;
2622  for (i = 0; i < fSizeX; i++) {
2623  for (j = 0; j < fSizeY; j++) {
2624  fDest[i][j] = working_matrix[i][j] * a;
2625  }
2626  }
2627  }
2628  break;
2629  case kTransformWalsh:
2630  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2631  kTransformInverse, kTransformWalsh);
2632  for (i = 0; i < fSizeX; i++) {
2633  for (j = 0; j < fSizeY; j++) {
2634  new_area = new_area + working_matrix[i][j];
2635  }
2636  }
2637  if (new_area != 0) {
2638  a = old_area / new_area;
2639  for (i = 0; i < fSizeX; i++) {
2640  for (j = 0; j < fSizeY; j++) {
2641  fDest[i][j] = working_matrix[i][j] * a;
2642  }
2643  }
2644  }
2645  break;
2646  case kTransformCos:
2647  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2648  kTransformInverse, kTransformCos);
2649  for (i = 0; i < fSizeX; i++) {
2650  for (j = 0; j < fSizeY; j++) {
2651  new_area = new_area + working_matrix[i][j];
2652  }
2653  }
2654  if (new_area != 0) {
2655  a = old_area / new_area;
2656  for (i = 0; i < fSizeX; i++) {
2657  for (j = 0; j < fSizeY; j++) {
2658  fDest[i][j] = working_matrix[i][j] * a;
2659  }
2660  }
2661  }
2662  break;
2663  case kTransformSin:
2664  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2665  kTransformInverse, kTransformSin);
2666  for (i = 0; i < fSizeX; i++) {
2667  for (j = 0; j < fSizeY; j++) {
2668  new_area = new_area + working_matrix[i][j];
2669  }
2670  }
2671  if (new_area != 0) {
2672  a = old_area / new_area;
2673  for (i = 0; i < fSizeX; i++) {
2674  for (j = 0; j < fSizeY; j++) {
2675  fDest[i][j] = working_matrix[i][j] * a;
2676  }
2677  }
2678  }
2679  break;
2680  case kTransformFourier:
2681  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2682  kTransformInverse, kTransformFourier);
2683  for (i = 0; i < fSizeX; i++) {
2684  for (j = 0; j < fSizeY; j++) {
2685  new_area = new_area + working_matrix[i][j];
2686  }
2687  }
2688  if (new_area != 0) {
2689  a = old_area / new_area;
2690  for (i = 0; i < fSizeX; i++) {
2691  for (j = 0; j < fSizeY; j++) {
2692  fDest[i][j] = working_matrix[i][j] * a;
2693  }
2694  }
2695  }
2696  break;
2697  case kTransformHartley:
2698  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2699  kTransformInverse, kTransformHartley);
2700  for (i = 0; i < fSizeX; i++) {
2701  for (j = 0; j < fSizeY; j++) {
2702  new_area = new_area + working_matrix[i][j];
2703  }
2704  }
2705  if (new_area != 0) {
2706  a = old_area / new_area;
2707  for (i = 0; i < fSizeX; i++) {
2708  for (j = 0; j < fSizeY; j++) {
2709  fDest[i][j] = working_matrix[i][j] * a;
2710  }
2711  }
2712  }
2713  break;
2714  case kTransformFourierWalsh:
2715  case kTransformFourierHaar:
2716  case kTransformWalshHaar:
2717  case kTransformCosWalsh:
2718  case kTransformCosHaar:
2719  case kTransformSinWalsh:
2720  case kTransformSinHaar:
2721  General2(working_matrix, working_vector, fSizeX, fSizeY,
2722  kTransformInverse, fTransformType, fDegree);
2723  for (i = 0; i < fSizeX; i++) {
2724  for (j = 0; j < fSizeY; j++) {
2725  new_area = new_area + working_matrix[i][j];
2726  }
2727  }
2728  if (new_area != 0) {
2729  a = old_area / new_area;
2730  for (i = 0; i < fSizeX; i++) {
2731  for (j = 0; j < fSizeY; j++) {
2732  fDest[i][j] = working_matrix[i][j] * a;
2733  }
2734  }
2735  }
2736  break;
2737  }
2738  for (i = 0; i < fSizeX; i++) {
2739  if (working_matrix) delete[]working_matrix[i];
2740  }
2741  delete[]working_matrix;
2742  delete[]working_vector;
2743  return;
2744 }
2745 
2746 
2747 ////////// END OF FILTER2_ZONAL FUNCTION/////////////////////////////////
2748 //////////ENHANCE2 FUNCTION - CALCULATES DIFFERENT 2-D ORTHOGONAL TRANSFORMS, MULTIPLIES GIVEN REGION BY ENHANCE COEFFICIENT AND TRANSFORMS IT BACK//////
2749 ////////////////////////////////////////////////////////////////////////////////
2750 ///////////////////////////////////////////////////////////////////////////////////////////
2751 
2752 void TSpectrum2Transform::Enhance(const Double_t **fSource, Double_t **fDest)
2753 {
2754 /* TWO-DIMENSIONAL ENHANCE ZONAL FUNCTION */
2755 /* This function transforms the source spectrum. The calling program */
2756 /* should fill in input parameters. Then it multiplies transformed */
2757 /* coefficients in the given region by the given */
2758 /* enhance_coeff and transforms it back */
2759 /* */
2760 /* Function parameters: */
2761 /* fSource-pointer to the matrix of source spectrum, its size should */
2762 /* be fSizeX*fSizeY */
2763 /* fDest-pointer to the matrix of destination data, its size should be */
2764 /* fSizeX*fSizeY */
2765 /* */
2766 //////////////////////////////////////////////////////////////////////////////////////////
2767 /*
2768 <div class=Section3>
2769 
2770 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of enhancement</span></b></p>
2771 
2772 <p class=MsoNormal><i>&nbsp;</i></p>
2773 
2774 <p class=MsoNormal><i>Function:</i></p>
2775 
2776 <p class=MsoNormal>void <a
2777 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumTransform2::Enhance</b></a><b>(const
2778 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a>
2779 **fSource, <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a>
2780 **fDest)</b></p>
2781 
2782 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
2783 
2784 <p class=MsoNormal style='text-align:justify'>This function transforms the
2785 source spectrum (for details see Transform function).  Before the Enhance
2786 function is called the class must be created by constructor and the type of the
2787 transform as well as some other parameters must be set using a set of setter
2788 functions. The Enhance function multiplies transformed coefficients in the given
2789 region (fXmin, fXmax, fYmin, fYmax) by the given fEnhancCoeff and transforms it
2790 back. Enhanced data are written into dest spectrum.</p>
2791 
2792 <p class=MsoNormal style='text-align:justify'><i>Example – script Enhance2.c:</i></p>
2793 
2794 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><img
2795 border=0 width=602 height=455 src="gif/spectrum2transform_enhance_image001.jpg"></span></p>
2796 
2797 <p class=MsoNormal><b>Fig. 1 Original two-dimensional noisy spectrum</b></p>
2798 
2799 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
2800 border=0 width=602 height=455 src="gif/spectrum2transform_enhance_image002.jpg"></span></i></p>
2801 
2802 <p class=MsoNormal style='text-align:justify'><b>Fig. 2 Enhanced spectrum of
2803 the data from Fig. 1 using Cosine transform (channels in region (0-63)x(0-63)
2804 were multiplied by 5) </b></p>
2805 
2806 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>&nbsp;</span></b></p>
2807 
2808 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
2809 
2810 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
2811 enhancement (class TSpectrumTransform2).</span></p>
2812 
2813 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
2814 do</span></p>
2815 
2816 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Enhance2.C</span></p>
2817 
2818 <p class=MsoNormal><span style='font-size:10.0pt'> </span></p>
2819 
2820 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Enhance2() {</span></p>
2821 
2822 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i, j;</span></p>
2823 
2824 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbinsx =
2825 256;</span></p>
2826 
2827 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbinsy =
2828 256;   </span></p>
2829 
2830 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmin  = 0;</span></p>
2831 
2832 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
2833 nbinsx;</span></p>
2834 
2835 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymin  = 0;</span></p>
2836 
2837 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t ymax  =
2838 nbinsy;</span></p>
2839 
2840 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2841 style='font-size:10.0pt'>Double_t ** source = new Double_t *[nbinsx];   </span></p>
2842 
2843 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t ** dest = new
2844 Double_t *[nbinsx];      </span></p>
2845 
2846 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i=0;i&lt;nbinsx;i++)</span></p>
2847 
2848 <p class=MsoNormal><span style='font-size:10.0pt'>                                                source[i]=new
2849 Double_t[nbinsy];</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'>                                                dest[i]=new
2854 Double_t[nbinsy];   </span></p>
2855 
2856 <p class=MsoNormal><span style='font-size:10.0pt'>   TH2F *trans = new
2857 TH2F(&quot;trans&quot;,&quot;Background
2858 estimation&quot;,nbinsx,xmin,xmax,nbinsy,ymin,ymax);</span></p>
2859 
2860 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
2861 TFile(&quot;TSpectrum2.root&quot;);</span></p>
2862 
2863 <p class=MsoNormal><span style='font-size:10.0pt'>   trans=(TH2F*)
2864 f-&gt;Get(&quot;back3;1&quot;);</span></p>
2865 
2866 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Tr = new
2867 TCanvas(&quot;Transform&quot;,&quot;Illustation of transform
2868 function&quot;,10,10,1000,700);</span></p>
2869 
2870 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbinsx;
2871 i++){</span></p>
2872 
2873 <p class=MsoNormal><span style='font-size:10.0pt'>     for (j = 0; j &lt;
2874 nbinsy; j++){</span></p>
2875 
2876 <p class=MsoNormal><span style='font-size:10.0pt'>                    </span><span
2877 lang=FR style='font-size:10.0pt'>source[i][j] = trans-&gt;GetBinContent(i + 1,j
2878 + 1); </span></p>
2879 
2880 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                 }</span></p>
2881 
2882 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }           </span></p>
2883 
2884 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2885 TSpectrumTransform2 *t = new TSpectrumTransform2(256,256);   </span></p>
2886 
2887 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2888 t-&gt;SetTransformType(t-&gt;kTransformCos,0);   </span></p>
2889 
2890 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2891 t-&gt;SetRegion(0,63,0,63);   </span></p>
2892 
2893 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2894 t-&gt;SetEnhanceCoeff(5);</span></p>
2895 
2896 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
2897 t-&gt;Enhance(source,dest);   </span></p>
2898 
2899 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'> </span><span
2900 style='font-size:10.0pt'>  trans-&gt;Draw(&quot;SURF&quot;);     </span></p>
2901 
2902 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
2903 
2904 </div>
2905 
2906  */
2907 
2908  Int_t i, j;
2909  Double_t a, old_area = 0, new_area = 0;
2910  Int_t size;
2911  Double_t *working_vector = 0, **working_matrix = 0;
2912  size = (Int_t) TMath::Max(fSizeX, fSizeY);
2913  switch (fTransformType) {
2914  case kTransformHaar:
2915  case kTransformWalsh:
2916  working_vector = new Double_t[2 * size];
2917  working_matrix = new Double_t *[fSizeX];
2918  for (i = 0; i < fSizeX; i++)
2919  working_matrix[i] = new Double_t[fSizeY];
2920  break;
2921  case kTransformCos:
2922  case kTransformSin:
2923  case kTransformFourier:
2924  case kTransformHartley:
2925  case kTransformFourierWalsh:
2926  case kTransformFourierHaar:
2927  case kTransformWalshHaar:
2928  working_vector = new Double_t[4 * size];
2929  working_matrix = new Double_t *[fSizeX];
2930  for (i = 0; i < fSizeX; i++)
2931  working_matrix[i] = new Double_t[2 * fSizeY];
2932  break;
2933  case kTransformCosWalsh:
2934  case kTransformCosHaar:
2935  case kTransformSinWalsh:
2936  case kTransformSinHaar:
2937  working_vector = new Double_t[8 * size];
2938  working_matrix = new Double_t *[fSizeX];
2939  for (i = 0; i < fSizeX; i++)
2940  working_matrix[i] = new Double_t[2 * fSizeY];
2941  break;
2942  }
2943  switch (fTransformType) {
2944  case kTransformHaar:
2945  for (i = 0; i < fSizeX; i++) {
2946  for (j = 0; j < fSizeY; j++) {
2947  working_matrix[i][j] = fSource[i][j];
2948  old_area = old_area + fSource[i][j];
2949  }
2950  }
2951  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2952  kTransformForward, kTransformHaar);
2953  break;
2954  case kTransformWalsh:
2955  for (i = 0; i < fSizeX; i++) {
2956  for (j = 0; j < fSizeY; j++) {
2957  working_matrix[i][j] = fSource[i][j];
2958  old_area = old_area + fSource[i][j];
2959  }
2960  }
2961  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
2962  kTransformForward, kTransformWalsh);
2963  break;
2964  case kTransformCos:
2965  for (i = 0; i < fSizeX; i++) {
2966  for (j = 0; j < fSizeY; j++) {
2967  working_matrix[i][j] = fSource[i][j];
2968  old_area = old_area + fSource[i][j];
2969  }
2970  }
2971  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2972  kTransformForward, kTransformCos);
2973  break;
2974  case kTransformSin:
2975  for (i = 0; i < fSizeX; i++) {
2976  for (j = 0; j < fSizeY; j++) {
2977  working_matrix[i][j] = fSource[i][j];
2978  old_area = old_area + fSource[i][j];
2979  }
2980  }
2981  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2982  kTransformForward, kTransformSin);
2983  break;
2984  case kTransformFourier:
2985  for (i = 0; i < fSizeX; i++) {
2986  for (j = 0; j < fSizeY; j++) {
2987  working_matrix[i][j] = fSource[i][j];
2988  old_area = old_area + fSource[i][j];
2989  }
2990  }
2991  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
2992  kTransformForward, kTransformFourier);
2993  break;
2994  case kTransformHartley:
2995  for (i = 0; i < fSizeX; i++) {
2996  for (j = 0; j < fSizeY; j++) {
2997  working_matrix[i][j] = fSource[i][j];
2998  old_area = old_area + fSource[i][j];
2999  }
3000  }
3001  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
3002  kTransformForward, kTransformHartley);
3003  break;
3004  case kTransformFourierWalsh:
3005  case kTransformFourierHaar:
3006  case kTransformWalshHaar:
3007  case kTransformCosWalsh:
3008  case kTransformCosHaar:
3009  case kTransformSinWalsh:
3010  case kTransformSinHaar:
3011  for (i = 0; i < fSizeX; i++) {
3012  for (j = 0; j < fSizeY; j++) {
3013  working_matrix[i][j] = fSource[i][j];
3014  old_area = old_area + fSource[i][j];
3015  }
3016  }
3017  General2(working_matrix, working_vector, fSizeX, fSizeY,
3018  kTransformForward, fTransformType, fDegree);
3019  break;
3020  }
3021  for (i = 0; i < fSizeX; i++) {
3022  for (j = 0; j < fSizeY; j++) {
3023  if (i >= fXmin && i <= fXmax && j >= fYmin && j <= fYmax)
3024  if (working_matrix) working_matrix[i][j] *= fEnhanceCoeff;
3025  }
3026  }
3027  if (fTransformType == kTransformFourier || fTransformType == kTransformFourierWalsh
3028  || fTransformType == kTransformFourierHaar) {
3029  for (i = 0; i < fSizeX; i++) {
3030  for (j = 0; j < fSizeY; j++) {
3031  if (i >= fXmin && i <= fXmax && j >= fYmin && j <= fYmax)
3032  working_matrix[i][j + fSizeY] *= fEnhanceCoeff;
3033  }
3034  }
3035  }
3036  switch (fTransformType) {
3037  case kTransformHaar:
3038  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
3039  kTransformInverse, kTransformHaar);
3040  for (i = 0; i < fSizeX; i++) {
3041  for (j = 0; j < fSizeY; j++) {
3042  new_area = new_area + working_matrix[i][j];
3043  }
3044  }
3045  if (new_area != 0) {
3046  a = old_area / new_area;
3047  for (i = 0; i < fSizeX; i++) {
3048  for (j = 0; j < fSizeY; j++) {
3049  fDest[i][j] = working_matrix[i][j] * a;
3050  }
3051  }
3052  }
3053  break;
3054  case kTransformWalsh:
3055  HaarWalsh2(working_matrix, working_vector, fSizeX, fSizeY,
3056  kTransformInverse, kTransformWalsh);
3057  for (i = 0; i < fSizeX; i++) {
3058  for (j = 0; j < fSizeY; j++) {
3059  new_area = new_area + working_matrix[i][j];
3060  }
3061  }
3062  if (new_area != 0) {
3063  a = old_area / new_area;
3064  for (i = 0; i < fSizeX; i++) {
3065  for (j = 0; j < fSizeY; j++) {
3066  fDest[i][j] = working_matrix[i][j] * a;
3067  }
3068  }
3069  }
3070  break;
3071  case kTransformCos:
3072  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
3073  kTransformInverse, kTransformCos);
3074  for (i = 0; i < fSizeX; i++) {
3075  for (j = 0; j < fSizeY; j++) {
3076  new_area = new_area + working_matrix[i][j];
3077  }
3078  }
3079  if (new_area != 0) {
3080  a = old_area / new_area;
3081  for (i = 0; i < fSizeX; i++) {
3082  for (j = 0; j < fSizeY; j++) {
3083  fDest[i][j] = working_matrix[i][j] * a;
3084  }
3085  }
3086  }
3087  break;
3088  case kTransformSin:
3089  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
3090  kTransformInverse, kTransformSin);
3091  for (i = 0; i < fSizeX; i++) {
3092  for (j = 0; j < fSizeY; j++) {
3093  new_area = new_area + working_matrix[i][j];
3094  }
3095  }
3096  if (new_area != 0) {
3097  a = old_area / new_area;
3098  for (i = 0; i < fSizeX; i++) {
3099  for (j = 0; j < fSizeY; j++) {
3100  fDest[i][j] = working_matrix[i][j] * a;
3101  }
3102  }
3103  }
3104  break;
3105  case kTransformFourier:
3106  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
3107  kTransformInverse, kTransformFourier);
3108  for (i = 0; i < fSizeX; i++) {
3109  for (j = 0; j < fSizeY; j++) {
3110  new_area = new_area + working_matrix[i][j];
3111  }
3112  }
3113  if (new_area != 0) {
3114  a = old_area / new_area;
3115  for (i = 0; i < fSizeX; i++) {
3116  for (j = 0; j < fSizeY; j++) {
3117  fDest[i][j] = working_matrix[i][j] * a;
3118  }
3119  }
3120  }
3121  break;
3122  case kTransformHartley:
3123  FourCos2(working_matrix, working_vector, fSizeX, fSizeY,
3124  kTransformInverse, kTransformHartley);
3125  for (i = 0; i < fSizeX; i++) {
3126  for (j = 0; j < fSizeY; j++) {
3127  new_area = new_area + working_matrix[i][j];
3128  }
3129  }
3130  if (new_area != 0) {
3131  a = old_area / new_area;
3132  for (i = 0; i < fSizeX; i++) {
3133  for (j = 0; j < fSizeY; j++) {
3134  fDest[i][j] = working_matrix[i][j] * a;
3135  }
3136  }
3137  }
3138  break;
3139  case kTransformFourierWalsh:
3140  case kTransformFourierHaar:
3141  case kTransformWalshHaar:
3142  case kTransformCosWalsh:
3143  case kTransformCosHaar:
3144  case kTransformSinWalsh:
3145  case kTransformSinHaar:
3146  General2(working_matrix, working_vector, fSizeX, fSizeY,
3147  kTransformInverse, fTransformType, fDegree);
3148  for (i = 0; i < fSizeX; i++) {
3149  for (j = 0; j < fSizeY; j++) {
3150  new_area = new_area + working_matrix[i][j];
3151  }
3152  }
3153  if (new_area != 0) {
3154  a = old_area / new_area;
3155  for (i = 0; i < fSizeX; i++) {
3156  for (j = 0; j < fSizeY; j++) {
3157  fDest[i][j] = working_matrix[i][j] * a;
3158  }
3159  }
3160  }
3161  break;
3162  }
3163  for (i = 0; i < fSizeX; i++) {
3164  if (working_matrix) delete[]working_matrix[i];
3165  }
3166  delete[]working_matrix;
3167  delete[]working_vector;
3168  return;
3169 }
3170 
3171 
3172 ////////// END OF ENHANCE2 FUNCTION/////////////////////////////////
3173 
3174 ////////////////////////////////////////////////////////////////////////////////
3175 ///////////////////////////////////////////////////////////////////////////////
3176 /// SETTER FUNCION
3177 ///
3178 /// This function sets the following parameters for transform:
3179 /// -transType - type of transform (Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar)
3180 /// -degree - degree of mixed transform, applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar transforms
3181 ///////////////////////////////////////////////////////////////////////////////
3182 
3183 void TSpectrum2Transform::SetTransformType(Int_t transType, Int_t degree)
3184 {
3185  Int_t j1, j2, n;
3186  j1 = 0;
3187  n = 1;
3188  for (; n < fSizeX;) {
3189  j1 += 1;
3190  n = n * 2;
3191  }
3192  j2 = 0;
3193  n = 1;
3194  for (; n < fSizeY;) {
3195  j2 += 1;
3196  n = n * 2;
3197  }
3198  if (transType < kTransformHaar || transType > kTransformSinHaar){
3199  Error ("TSpectrumTransform","Invalid type of transform");
3200  return;
3201  }
3202  if (transType >= kTransformFourierWalsh && transType <= kTransformSinHaar) {
3203  if (degree > j1 || degree > j2 || degree < 1){
3204  Error ("TSpectrumTransform","Invalid degree of mixed transform");
3205  return;
3206  }
3207  }
3208  fTransformType = transType;
3209  fDegree = degree;
3210 }
3211 
3212 ////////////////////////////////////////////////////////////////////////////////
3213 ///////////////////////////////////////////////////////////////////////////////
3214 /// SETTER FUNCION
3215 ///
3216 /// This function sets the filtering or enhancement region:
3217 /// -xmin, xmax, ymin, ymax
3218 ///////////////////////////////////////////////////////////////////////////////
3219 
3220 void TSpectrum2Transform::SetRegion(Int_t xmin, Int_t xmax, Int_t ymin, Int_t ymax)
3221 {
3222  if(xmin<0 || xmax < xmin || xmax >= fSizeX){
3223  Error("TSpectrumTransform", "Wrong range");
3224  return;
3225  }
3226  if(ymin<0 || ymax < ymin || ymax >= fSizeY){
3227  Error("TSpectrumTransform", "Wrong range");
3228  return;
3229  }
3230  fXmin = xmin;
3231  fXmax = xmax;
3232  fYmin = ymin;
3233  fYmax = ymax;
3234 }
3235 
3236 ////////////////////////////////////////////////////////////////////////////////
3237 ///////////////////////////////////////////////////////////////////////////////
3238 /// SETTER FUNCION
3239 ///
3240 /// This function sets the direction of the transform:
3241 /// -direction (forward or inverse)
3242 ///////////////////////////////////////////////////////////////////////////////
3243 
3244 void TSpectrum2Transform::SetDirection(Int_t direction)
3245 {
3246  if(direction != kTransformForward && direction != kTransformInverse){
3247  Error("TSpectrumTransform", "Wrong direction");
3248  return;
3249  }
3250  fDirection = direction;
3251 }
3252 
3253 ////////////////////////////////////////////////////////////////////////////////
3254 ///////////////////////////////////////////////////////////////////////////////
3255 /// SETTER FUNCION
3256 ///
3257 /// This function sets the filter coefficient:
3258 /// -filterCoeff - after the transform the filtered region (xmin, xmax, ymin, ymax) is replaced by this coefficient. Applies only for filtereng operation.
3259 ///////////////////////////////////////////////////////////////////////////////
3260 
3261 void TSpectrum2Transform::SetFilterCoeff(Double_t filterCoeff)
3262 {
3263  fFilterCoeff = filterCoeff;
3264 }
3265 
3266 ////////////////////////////////////////////////////////////////////////////////
3267 ///////////////////////////////////////////////////////////////////////////////
3268 /// SETTER FUNCION
3269 ///
3270 /// This function sets the enhancement coefficient:
3271 /// -enhanceCoeff - after the transform the enhanced region (xmin, xmax, ymin, ymax) is multiplied by this coefficient. Applies only for enhancement operation.
3272 ///////////////////////////////////////////////////////////////////////////////
3273 
3274 void TSpectrum2Transform::SetEnhanceCoeff(Double_t enhanceCoeff)
3275 {
3276  fEnhanceCoeff = enhanceCoeff;
3277 }
3278 
TSpectrum2Transform::SetRegion
void SetRegion(Int_t xmin, Int_t xmax, Int_t ymin, Int_t ymax)
SETTER FUNCION.
Definition: TSpectrum2Transform.cxx:3220
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
TSpectrum2Transform::TSpectrum2Transform
TSpectrum2Transform()
default constructor
Definition: TSpectrum2Transform.cxx:66
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:3244
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
sample_config.c
c
Definition: sample_config.py:3
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
sample_config.ip
ip
Definition: sample_config.py:8
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
TSpectrum2Transform::FilterZonal
void FilterZonal(const Double_t **fSource, Double_t **fDest)
Definition: TSpectrum2Transform.cxx:2308
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:3261
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
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
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:2752
TMath::Cos
Double_t Cos(Double_t)
Definition: TMath.h:424
TSpectrum2Transform::SetEnhanceCoeff
void SetEnhanceCoeff(Double_t enhanceCoeff)
SETTER FUNCION.
Definition: TSpectrum2Transform.cxx:3274
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:279
TSpectrum2Transform::SetTransformType
void SetTransformType(Int_t transType, Int_t degree)
SETTER FUNCION.
Definition: TSpectrum2Transform.cxx:3183
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
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
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