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