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