hist » spectrum » src

TSpectrumTransform.cxx

Go to the documentation of this file.
00001 // @(#)root/spectrum:$Id: TSpectrumTransform.cxx 36003 2010-10-01 12:54:39Z brun $
00002 // Author: Miroslav Morhac   25/09/06
00003 
00004 //__________________________________________________________________________
00005 //   THIS CLASS CONTAINS ORTHOGONAL TRANSFORM  FUNCTIONS.                  //
00006 //                                                                         //
00007 //   These functions were written by:                                      //
00008 //   Miroslav Morhac                                                       //
00009 //   Institute of Physics                                                  //
00010 //   Slovak Academy of Sciences                                            //
00011 //   Dubravska cesta 9, 842 28 BRATISLAVA                                  //
00012 //   SLOVAKIA                                                              //
00013 //                                                                         //
00014 //   email:fyzimiro@savba.sk,    fax:+421 7 54772479                       //
00015 //                                                                         //
00016 //  The original code in C has been repackaged as a C++ class by R.Brun    //
00017 //                                                                         //
00018 //  The algorithms in this class have been published in the following      //
00019 //  references:                                                            //
00020 //                                                                         //
00021 //  [1] C.V. Hampton, B. Lian, Wm. C. McHarris: Fast-Fourier-transform     //
00022 //      spectral enhancement techniques for gamma-ray spectroscopy.NIM A353//
00023 //      (1994) 280-284.                                                    //
00024 //  [2] Morhac M., Matousek V., New adaptive Cosine-Walsh  transform and   //
00025 //      its application to nuclear data compression, IEEE Transactions on  //
00026 //      Signal Processing 48 (2000) 2693.                                  //  
00027 //  [3] Morhac M., Matousek V., Data compression using new fast adaptive   //
00028 //      Cosine-Haar transforms, Digital Signal Processing 8 (1998) 63.     //
00029 //  [4] Morhac M., Matousek V.: Multidimensional nuclear data compression  //
00030 //      using fast adaptive Walsh-Haar transform. Acta Physica Slovaca 51  //
00031 //     (2001) 307.                                                         //
00032 //____________________________________________________________________________
00033 
00034 #include "TSpectrumTransform.h"
00035 #include "TMath.h"
00036 
00037 ClassImp(TSpectrumTransform) 
00038  
00039 //____________________________________________________________________________    
00040 TSpectrumTransform::TSpectrumTransform() 
00041 {
00042    //default constructor
00043    fSize=0;
00044    fTransformType=kTransformCos;
00045    fDegree=0;
00046    fDirection=kTransformForward;
00047    fXmin=0;
00048    fXmax=0;
00049    fFilterCoeff=0;
00050    fEnhanceCoeff=0.5;
00051 }
00052 
00053 //____________________________________________________________________________    
00054 TSpectrumTransform::TSpectrumTransform(Int_t size):TNamed("SpectrumTransform", "Miroslav Morhac transformer") 
00055 {
00056 //the constructor creates TSpectrumTransform object. Its size must be > than zero and must be power of 2.
00057 //It sets default transform type to be Cosine transform. Transform parameters can be changed using setter functions.
00058    Int_t j,n;
00059    if (size <= 0){
00060       Error ("TSpectrumTransform","Invalid length, must be > than 0");
00061       return;
00062    }    
00063    j = 0;
00064    n = 1;
00065    for (; n < size;) {
00066       j += 1;
00067       n = n * 2;
00068    }
00069    if (n != size){
00070       Error ("TSpectrumTransform","Invalid length, must be power of 2");
00071       return;   
00072    }
00073    fSize=size;
00074    fTransformType=kTransformCos;
00075    fDegree=0;
00076    fDirection=kTransformForward;
00077    fXmin=size/4;
00078    fXmax=size-1;
00079    fFilterCoeff=0;
00080    fEnhanceCoeff=0.5;
00081 }
00082 
00083 
00084 //______________________________________________________________________________
00085 TSpectrumTransform::~TSpectrumTransform() 
00086 {
00087    //destructor
00088 }
00089 
00090 //_____________________________________________________________________________
00091 void TSpectrumTransform::Haar(float *working_space, int num, int direction) 
00092 {   
00093 //////////////////////////////////////////////////////////////////////////////////
00094 //   AUXILIARY FUNCION                                                          //
00095 //                                                                              //
00096 //   This function calculates Haar transform of a part of data                   //
00097 //      Function parameters:                                                    //
00098 //              -working_space-pointer to vector of transformed data            //
00099 //              -num-length of processed data                                   //
00100 //              -direction-forward or inverse transform                         //
00101 //                                                                              //
00102 //////////////////////////////////////////////////////////////////////////////////
00103    int i, ii, li, l2, l3, j, jj, jj1, lj, iter, m, jmin, jmax;
00104    double a, b, c, wlk;
00105    float val;
00106    for (i = 0; i < num; i++)
00107       working_space[i + num] = 0;
00108    i = num;
00109    iter = 0;
00110    for (; i > 1;) {
00111       iter += 1;
00112       i = i / 2;
00113    }
00114    if (direction == kTransformForward) {
00115       for (m = 1; m <= iter; m++) {
00116          li = iter + 1 - m;
00117          l2 = (int) TMath::Power(2, li - 1);
00118          for (i = 0; i < (2 * l2); i++) {
00119             working_space[num + i] = working_space[i];
00120          }
00121          for (j = 0; j < l2; j++) {
00122             l3 = l2 + j;
00123             jj = 2 * j;
00124             val = working_space[jj + num] + working_space[jj + 1 + num];
00125             working_space[j] = val;
00126             val = working_space[jj + num] - working_space[jj + 1 + num];
00127             working_space[l3] = val;
00128          }
00129       }
00130    }
00131    val = working_space[0];
00132    val = val / TMath::Sqrt(TMath::Power(2, iter));
00133    working_space[0] = val;
00134    val = working_space[1];
00135    val = val / TMath::Sqrt(TMath::Power(2, iter));
00136    working_space[1] = val;
00137    for (ii = 2; ii <= iter; ii++) {
00138       i = ii - 1;
00139       wlk = 1 / TMath::Sqrt(TMath::Power(2, iter - i));
00140       jmin = (int) TMath::Power(2, i);
00141       jmax = (int) TMath::Power(2, ii) - 1;
00142       for (j = jmin; j <= jmax; j++) {
00143          val = working_space[j];
00144          a = val;
00145          a = a * wlk;
00146          val = a;
00147          working_space[j] = val;
00148       }
00149    }
00150    if (direction == kTransformInverse) {
00151       for (m = 1; m <= iter; m++) {
00152          a = 2;
00153          b = m - 1;
00154          c = TMath::Power(a, b);
00155          li = (int) c;
00156          for (i = 0; i < (2 * li); i++) {
00157             working_space[i + num] = working_space[i];
00158          }
00159          for (j = 0; j < li; j++) {
00160             lj = li + j;
00161             jj = 2 * (j + 1) - 1;
00162             jj1 = jj - 1;
00163             val = working_space[j + num] - working_space[lj + num];
00164             working_space[jj] = val;
00165             val = working_space[j + num] + working_space[lj + num];
00166             working_space[jj1] = val;
00167          }
00168       }
00169    }
00170    return;
00171 }
00172 
00173 //____________________________________________________________________________    
00174 void TSpectrumTransform::Walsh(float *working_space, int num) 
00175 {
00176 //////////////////////////////////////////////////////////////////////////////////
00177 //   AUXILIARY FUNCION                                                          //
00178 //                                                                              //
00179 //   This function calculates Walsh transform of a part of data                  //
00180 //      Function parameters:                                                    //
00181 //              -working_space-pointer to vector of transformed data            //
00182 //              -num-length of processed data                                   //
00183 //                                                                              //
00184 //////////////////////////////////////////////////////////////////////////////////
00185    int i, m, nump = 1, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter;
00186    double a;
00187    float val1, val2;
00188    for (i = 0; i < num; i++)
00189       working_space[i + num] = 0;
00190    i = num;
00191    iter = 0;
00192    for (; i > 1;) {
00193       iter += 1;
00194       i = i / 2;
00195    }
00196    for (m = 1; m <= iter; m++) {
00197       if (m == 1)
00198          nump = 1;
00199       
00200       else
00201          nump = nump * 2;
00202       mnum = num / nump;
00203       mnum2 = mnum / 2;
00204       for (mp = 0; mp < nump; mp++) {
00205          ib = mp * mnum;
00206          for (mp2 = 0; mp2 < mnum2; mp2++) {
00207             mnum21 = mnum2 + mp2 + ib;
00208             iba = ib + mp2;
00209             val1 = working_space[iba];
00210             val2 = working_space[mnum21];
00211             working_space[iba + num] = val1 + val2;
00212             working_space[mnum21 + num] = val1 - val2;
00213          }
00214       }
00215       for (i = 0; i < num; i++) {
00216          working_space[i] = working_space[i + num];
00217       }
00218    }
00219    a = num;
00220    a = TMath::Sqrt(a);
00221    val2 = a;
00222    for (i = 0; i < num; i++) {
00223       val1 = working_space[i];
00224       val1 = val1 / val2;
00225       working_space[i] = val1;
00226    }
00227    return;
00228 }
00229 
00230 //____________________________________________________________________________    
00231 void TSpectrumTransform::BitReverse(float *working_space, int num) 
00232 {   
00233 //////////////////////////////////////////////////////////////////////////////////
00234 //   AUXILIARY FUNCION                                                          //
00235 //                                                                              //
00236 //   This function carries out bir-reverse reordering of data                    //
00237 //      Function parameters:                                                    //
00238 //              -working_space-pointer to vector of processed data              //
00239 //              -num-length of processed data                                   //
00240 //                                                                              //
00241 //////////////////////////////////////////////////////////////////////////////////
00242    int ipower[26];
00243    int i, ib, il, ibd, ip, ifac, i1;
00244    for (i = 0; i < num; i++) {
00245       working_space[i + num] = working_space[i];
00246    }
00247    for (i = 1; i <= num; i++) {
00248       ib = i - 1;
00249       il = 1;
00250    lab9:ibd = ib / 2;
00251       ipower[il - 1] = 1;
00252       if (ib == (ibd * 2))
00253          ipower[il - 1] = 0;
00254       if (ibd == 0)
00255          goto lab10;
00256       ib = ibd;
00257       il = il + 1;
00258       goto lab9;
00259    lab10:ip = 1;
00260       ifac = num;
00261       for (i1 = 1; i1 <= il; i1++) {
00262          ifac = ifac / 2;
00263          ip = ip + ifac * ipower[i1 - 1];
00264       }
00265       working_space[ip - 1] = working_space[i - 1 + num];
00266    }
00267    return;
00268 }
00269 
00270 //____________________________________________________________________________    
00271 void TSpectrumTransform::Fourier(float *working_space, int num, int hartley,
00272                           int direction, int zt_clear) 
00273 {   
00274 //////////////////////////////////////////////////////////////////////////////////
00275 //   AUXILIARY FUNCION                                                          //
00276 //                                                                              //
00277 //   This function calculates Fourier based transform of a part of data          //
00278 //      Function parameters:                                                    //
00279 //              -working_space-pointer to vector of transformed data            //
00280 //              -num-length of processed data                                   //
00281 //              -hartley-1 if it is Hartley transform, 0 othewise               //
00282 //              -direction-forward or inverse transform                         //
00283 //                                                                              //
00284 //////////////////////////////////////////////////////////////////////////////////
00285    int nxp2, nxp, i, j, k, m, iter, mxp, j1, j2, n1, n2, it;
00286    double a, b, c, d, sign, wpwr, arg, wr, wi, tr, ti, pi =
00287        3.14159265358979323846;
00288    float val1, val2, val3, val4;
00289    if (direction == kTransformForward && zt_clear == 0) {
00290       for (i = 0; i < num; i++)
00291          working_space[i + num] = 0;
00292    }
00293    i = num;
00294    iter = 0;
00295    for (; i > 1;) {
00296       iter += 1;
00297       i = i / 2;
00298    }
00299    sign = -1;
00300    if (direction == kTransformInverse)
00301       sign = 1;
00302    nxp2 = num;
00303    for (it = 1; it <= iter; it++) {
00304       nxp = nxp2;
00305       nxp2 = nxp / 2;
00306       a = nxp2;
00307       wpwr = pi / a;
00308       for (m = 1; m <= nxp2; m++) {
00309          a = m - 1;
00310          arg = a * wpwr;
00311          wr = TMath::Cos(arg);
00312          wi = sign * TMath::Sin(arg);
00313          for (mxp = nxp; mxp <= num; mxp += nxp) {
00314             j1 = mxp - nxp + m;
00315             j2 = j1 + nxp2;
00316             val1 = working_space[j1 - 1];
00317             val2 = working_space[j2 - 1];
00318             val3 = working_space[j1 - 1 + num];
00319             val4 = working_space[j2 - 1 + num];
00320             a = val1;
00321             b = val2;
00322             c = val3;
00323             d = val4;
00324             tr = a - b;
00325             ti = c - d;
00326             a = a + b;
00327             val1 = a;
00328             working_space[j1 - 1] = val1;
00329             c = c + d;
00330             val1 = c;
00331             working_space[j1 - 1 + num] = val1;
00332             a = tr * wr - ti * wi;
00333             val1 = a;
00334             working_space[j2 - 1] = val1;
00335             a = ti * wr + tr * wi;
00336             val1 = a;
00337             working_space[j2 - 1 + num] = val1;
00338          }
00339       }
00340    }
00341    n2 = num / 2;
00342    n1 = num - 1;
00343    j = 1;
00344    for (i = 1; i <= n1; i++) {
00345       if (i >= j)
00346          goto lab55;
00347       val1 = working_space[j - 1];
00348       val2 = working_space[j - 1 + num];
00349       val3 = working_space[i - 1];
00350       working_space[j - 1] = val3;
00351       working_space[j - 1 + num] = working_space[i - 1 + num];
00352       working_space[i - 1] = val1;
00353       working_space[i - 1 + num] = val2;
00354       lab55: k = n2;
00355       lab60: if (k >= j) goto lab65;
00356       j = j - k;
00357       k = k / 2;
00358       goto lab60;
00359       lab65: j = j + k;
00360    }
00361    a = num;
00362    a = TMath::Sqrt(a);
00363    for (i = 0; i < num; i++) {
00364       if (hartley == 0) {
00365          val1 = working_space[i];
00366          b = val1;
00367          b = b / a;
00368          val1 = b;
00369          working_space[i] = val1;
00370          b = working_space[i + num];
00371          b = b / a;
00372          working_space[i + num] = b;
00373       }
00374       
00375       else {
00376          b = working_space[i];
00377          c = working_space[i + num];
00378          b = (b + c) / a;
00379          working_space[i] = b;
00380          working_space[i + num] = 0;
00381       }
00382    }
00383    if (hartley == 1 && direction == kTransformInverse) {
00384       for (i = 1; i < num; i++)
00385          working_space[num - i + num] = working_space[i];
00386       working_space[0 + num] = working_space[0];
00387       for (i = 0; i < num; i++) {
00388          working_space[i] = working_space[i + num];
00389          working_space[i + num] = 0;
00390       }
00391    }
00392    return;
00393 }
00394 
00395 //____________________________________________________________________________    
00396 void TSpectrumTransform::BitReverseHaar(float *working_space, int shift, int num,
00397                                  int start) 
00398 {
00399 //////////////////////////////////////////////////////////////////////////////////
00400 //   AUXILIARY FUNCION                                                          //
00401 //                                                                              //
00402 //   This function carries out bir-reverse reordering for Haar transform         //
00403 //      Function parameters:                                                    //
00404 //              -working_space-pointer to vector of processed data              //
00405 //              -shift-shift of position of processing                          //
00406 //              -start-initial position of processed data                       //
00407 //              -num-length of processed data                                   //
00408 //                                                                              //
00409 //////////////////////////////////////////////////////////////////////////////////
00410    int ipower[26];
00411    int i, ib, il, ibd, ip, ifac, i1;
00412    for (i = 0; i < num; i++) {
00413       working_space[i + shift + start] = working_space[i + start];
00414       working_space[i + shift + start + 2 * shift] =
00415           working_space[i + start + 2 * shift];
00416    }
00417    for (i = 1; i <= num; i++) {
00418       ib = i - 1;
00419       il = 1;
00420       lab9: ibd = ib / 2;
00421       ipower[il - 1] = 1;
00422       if (ib == (ibd * 2))
00423          ipower[il - 1] = 0;
00424       if (ibd == 0)
00425          goto lab10;
00426       ib = ibd;
00427       il = il + 1;
00428       goto lab9;
00429       lab10: ip = 1;
00430       ifac = num;
00431       for (i1 = 1; i1 <= il; i1++) {
00432          ifac = ifac / 2;
00433          ip = ip + ifac * ipower[i1 - 1];
00434       }
00435       working_space[ip - 1 + start] =
00436       working_space[i - 1 + shift + start];
00437       working_space[ip - 1 + start + 2 * shift] =
00438       working_space[i - 1 + shift + start + 2 * shift];
00439    }
00440    return;
00441 }
00442 
00443 //____________________________________________________________________________    
00444 int TSpectrumTransform::GeneralExe(float *working_space, int zt_clear, int num,
00445                             int degree, int type) 
00446 {
00447 //////////////////////////////////////////////////////////////////////////////////
00448 //   AUXILIARY FUNCION                                                          //
00449 //                                                                              //
00450 //   This function calculates generalized (mixed) transforms of different degrees//
00451 //      Function parameters:                                                    //
00452 //              -working_space-pointer to vector of transformed data            //
00453 //              -zt_clear-flag to clear imaginary data before staring           //
00454 //              -num-length of processed data                                   //
00455 //              -degree-degree of transform (see manual)                        //
00456 //              -type-type of mixed transform (see manual)                      //
00457 //                                                                              //
00458 //////////////////////////////////////////////////////////////////////////////////
00459    int i, j, k, m, nump, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter,
00460        mp2step, mppom, ring;
00461    double a, b, c, d, wpwr, arg, wr, wi, tr, ti, pi =
00462        3.14159265358979323846;
00463    float val1, val2, val3, val4, a0oldr = 0, b0oldr = 0, a0r, b0r;
00464    if (zt_clear == 0) {
00465       for (i = 0; i < num; i++)
00466          working_space[i + 2 * num] = 0;
00467    }
00468    i = num;
00469    iter = 0;
00470    for (; i > 1;) {
00471       iter += 1;
00472       i = i / 2;
00473    }
00474    a = num;
00475    wpwr = 2.0 * pi / a;
00476    nump = num;
00477    mp2step = 1;
00478    ring = num;
00479    for (i = 0; i < iter - degree; i++)
00480       ring = ring / 2;
00481    for (m = 1; m <= iter; m++) {
00482       nump = nump / 2;
00483       mnum = num / nump;
00484       mnum2 = mnum / 2;
00485       if (m > degree
00486            && (type == kTransformFourierHaar
00487                || type == kTransformWalshHaar
00488                || type == kTransformCosHaar
00489                || type == kTransformSinHaar))
00490          mp2step *= 2;
00491       if (ring > 1)
00492          ring = ring / 2;
00493       for (mp = 0; mp < nump; mp++) {
00494          if (type != kTransformWalshHaar) {
00495             mppom = mp;
00496             mppom = mppom % ring;
00497             a = 0;
00498             j = 1;
00499             k = num / 4;
00500             for (i = 0; i < (iter - 1); i++) {
00501                if ((mppom & j) != 0)
00502                   a = a + k;
00503                j = j * 2;
00504                k = k / 2;
00505             }
00506             arg = a * wpwr;
00507             wr = TMath::Cos(arg);
00508             wi = TMath::Sin(arg);
00509          }
00510          
00511          else {
00512             wr = 1;
00513             wi = 0;
00514          }
00515          ib = mp * mnum;
00516          for (mp2 = 0; mp2 < mnum2; mp2++) {
00517             mnum21 = mnum2 + mp2 + ib;
00518             iba = ib + mp2;
00519             if (mp2 % mp2step == 0) {
00520                a0r = a0oldr;
00521                b0r = b0oldr;
00522                a0r = 1 / TMath::Sqrt(2.0);
00523                b0r = 1 / TMath::Sqrt(2.0);
00524             }
00525             
00526             else {
00527                a0r = 1;
00528                b0r = 0;
00529             }
00530             val1 = working_space[iba];
00531             val2 = working_space[mnum21];
00532             val3 = working_space[iba + 2 * num];
00533             val4 = working_space[mnum21 + 2 * num];
00534             a = val1;
00535             b = val2;
00536             c = val3;
00537             d = val4;
00538             tr = a * a0r + b * b0r;
00539             val1 = tr;
00540             working_space[num + iba] = val1;
00541             ti = c * a0r + d * b0r;
00542             val1 = ti;
00543             working_space[num + iba + 2 * num] = val1;
00544             tr =
00545                 a * b0r * wr - c * b0r * wi - b * a0r * wr + d * a0r * wi;
00546             val1 = tr;
00547             working_space[num + mnum21] = val1;
00548             ti =
00549                 c * b0r * wr + a * b0r * wi - d * a0r * wr - b * a0r * wi;
00550             val1 = ti;
00551             working_space[num + mnum21 + 2 * num] = val1;
00552          }
00553       }
00554       for (i = 0; i < num; i++) {
00555          val1 = working_space[num + i];
00556          working_space[i] = val1;
00557          val1 = working_space[num + i + 2 * num];
00558          working_space[i + 2 * num] = val1;
00559       }
00560    }
00561    return (0);
00562 }
00563 
00564 //____________________________________________________________________________    
00565 int TSpectrumTransform::GeneralInv(float *working_space, int num, int degree,
00566                             int type) 
00567 {
00568 //////////////////////////////////////////////////////////////////////////////////
00569 //   AUXILIARY FUNCION                                                          //
00570 //                                                                              //
00571 //   This function calculates inverse generalized (mixed) transforms             //
00572 //      Function parameters:                                                    //
00573 //              -working_space-pointer to vector of transformed data            //
00574 //              -num-length of processed data                                   //
00575 //              -degree-degree of transform (see manual)                        //
00576 //              -type-type of mixed transform (see manual)                      //
00577 //                                                                              //
00578 //////////////////////////////////////////////////////////////////////////////////
00579    int i, j, k, m, nump =
00580        1, mnum, mnum2, mp, ib, mp2, mnum21, iba, iter, mp2step, mppom,
00581        ring;
00582    double a, b, c, d, wpwr, arg, wr, wi, tr, ti, pi =
00583        3.14159265358979323846;
00584    float val1, val2, val3, val4, a0oldr = 0, b0oldr = 0, a0r, b0r;
00585    i = num;
00586    iter = 0;
00587    for (; i > 1;) {
00588       iter += 1;
00589       i = i / 2;
00590    }
00591    a = num;
00592    wpwr = 2.0 * pi / a;
00593    mp2step = 1;
00594    if (type == kTransformFourierHaar || type == kTransformWalshHaar
00595         || type == kTransformCosHaar || type == kTransformSinHaar) {
00596       for (i = 0; i < iter - degree; i++)
00597          mp2step *= 2;
00598    }
00599    ring = 1;
00600    for (m = 1; m <= iter; m++) {
00601       if (m == 1)
00602          nump = 1;
00603       
00604       else
00605          nump = nump * 2;
00606       mnum = num / nump;
00607       mnum2 = mnum / 2;
00608       if (m > iter - degree + 1)
00609          ring *= 2;
00610       for (mp = nump - 1; mp >= 0; mp--) {
00611          if (type != kTransformWalshHaar) {
00612             mppom = mp;
00613             mppom = mppom % ring;
00614             a = 0;
00615             j = 1;
00616             k = num / 4;
00617             for (i = 0; i < (iter - 1); i++) {
00618                if ((mppom & j) != 0)
00619                   a = a + k;
00620                j = j * 2;
00621                k = k / 2;
00622             }
00623             arg = a * wpwr;
00624             wr = TMath::Cos(arg);
00625             wi = TMath::Sin(arg);
00626          }
00627          
00628          else {
00629             wr = 1;
00630             wi = 0;
00631          }
00632          ib = mp * mnum;
00633          for (mp2 = 0; mp2 < mnum2; mp2++) {
00634             mnum21 = mnum2 + mp2 + ib;
00635             iba = ib + mp2;
00636             if (mp2 % mp2step == 0) {
00637                a0r = a0oldr;
00638                b0r = b0oldr;
00639                a0r = 1 / TMath::Sqrt(2.0);
00640                b0r = 1 / TMath::Sqrt(2.0);
00641             }
00642             
00643             else {
00644                a0r = 1;
00645                b0r = 0;
00646             }
00647             val1 = working_space[iba];
00648             val2 = working_space[mnum21];
00649             val3 = working_space[iba + 2 * num];
00650             val4 = working_space[mnum21 + 2 * num];
00651             a = val1;
00652             b = val2;
00653             c = val3;
00654             d = val4;
00655             tr = a * a0r + b * wr * b0r + d * wi * b0r;
00656             val1 = tr;
00657             working_space[num + iba] = val1;
00658             ti = c * a0r + d * wr * b0r - b * wi * b0r;
00659             val1 = ti;
00660             working_space[num + iba + 2 * num] = val1;
00661             tr = a * b0r - b * wr * a0r - d * wi * a0r;
00662             val1 = tr;
00663             working_space[num + mnum21] = val1;
00664             ti = c * b0r - d * wr * a0r + b * wi * a0r;
00665             val1 = ti;
00666             working_space[num + mnum21 + 2 * num] = val1;
00667          }
00668       }
00669       if (m <= iter - degree
00670            && (type == kTransformFourierHaar
00671                || type == kTransformWalshHaar
00672                || type == kTransformCosHaar
00673                || type == kTransformSinHaar))
00674          mp2step /= 2;
00675       for (i = 0; i < num; i++) {
00676          val1 = working_space[num + i];
00677          working_space[i] = val1;
00678          val1 = working_space[num + i + 2 * num];
00679          working_space[i + 2 * num] = val1;
00680       }
00681    }
00682    return (0);
00683 }
00684 
00685 
00686 //////////END OF AUXILIARY FUNCTIONS FOR TRANSFORM! FUNCTION////////////////////////
00687 //////////TRANSFORM FUNCTION - CALCULATES DIFFERENT 1-D DIRECT AND INVERSE ORTHOGONAL TRANSFORMS//////
00688 
00689 //____________________________________________________________________________    
00690 void TSpectrumTransform::Transform(const float *source, float *destVector)
00691 {   
00692 ///////////////////////////////////////////////////////////////////////////////
00693 //        ONE-DIMENSIONAL TRANSFORM FUNCTION                                 
00694 //        This function transforms the source spectrum. The calling program 
00695 //        should fill in input parameters.                                    
00696 //        Transformed data are written into dest spectrum.                  
00697 //                                                                           
00698 //        Function parameters:                                               
00699 //        source-pointer to the vector of source spectrum, its length should 
00700 //             be size except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR     
00701 //             transform. These need 2*size length to supply real and        
00702 //             imaginary coefficients.                                       
00703 //        destVector-pointer to the vector of dest data, its length should be
00704 //             size except for direct FOURIER, FOUR-WALSH, FOUR-HAAR. These    
00705 //             need 2*size length to store real and imaginary coefficients    
00706 //                                                                         
00707 ///////////////////////////////////////////////////////////////////////////////
00708 //Begin_Html <!--
00709 /* -->
00710 <div class=Section1>
00711 
00712 <p class=MsoNormal><b><span style='font-size:14.0pt'>Transform methods</span></b></p>
00713 
00714 <p class=MsoNormal style='text-align:justify'><i>&nbsp;</i></p>
00715 
00716 <p class=MsoNormal style='text-align:justify'><i>Goal: to analyze experimental
00717 data using orthogonal transforms</i></p>
00718 
00719 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00720 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00721 </span>orthogonal transforms can be successfully used for the processing of
00722 nuclear spectra (not only) </p>
00723 
00724 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00725 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00726 </span>they can be used to remove high frequency noise, to increase
00727 signal-to-background ratio as well as to enhance low intensity components [1],
00728 to carry out e.g. Fourier analysis etc. </p>
00729 
00730 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00731 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00732 </span>we have implemented the function for the calculation of the commonly
00733 used orthogonal transforms as well as functions for the filtration and
00734 enhancement of experimental data</p>
00735 
00736 <p class=MsoNormal><i>&nbsp;</i></p>
00737 
00738 <p class=MsoNormal><i>Function:</i></p>
00739 
00740 <p class=MsoNormal><b>void TSpectrumTransform::Transform(const <a
00741 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> *fSource,
00742 <a href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a>
00743 *fDest)</b></p>
00744 
00745 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00746 
00747 <p class=MsoNormal style='text-align:justify'>This function transforms the
00748 source spectrum according to the given input parameters. Transformed data are
00749 written into dest spectrum. Before the Transform function is called the class
00750 must be created by constructor and the type of the transform as well as some
00751 other parameters must be set using a set of setter functions.</p>
00752 
00753 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00754 
00755 <p class=MsoNormal><i><span style='color:red'>Member variables of
00756 TSpectrumTransform class:</span></i></p>
00757 
00758 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'> <b>fSource</b>-pointer
00759 to the vector of source spectrum. Its length should be equal to the “fSize”
00760 parameter except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR transforms. These
00761 need “2*fSize” length to supply real and imaginary coefficients.                   </p>
00762 
00763 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify'><b>fDest</b>-pointer
00764 to the vector of destination spectrum. Its length should be equal to the
00765 “fSize” parameter except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR transforms.
00766 These need “2*fSize” length to store real and imaginary coefficients. </p>
00767 
00768 <p class=MsoNormal style='text-align:justify'>        <b>fSize</b>-basic length
00769 of the source and dest spectrum. <span style='color:fuchsia'>It should be power
00770 of 2.</span></p>
00771 
00772 <p class=MsoNormal style='margin-left:25.65pt;text-align:justify;text-indent:
00773 -2.85pt'><b>fType</b>-type of transform</p>
00774 
00775 <p class=MsoNormal style='text-align:justify'>            Classic transforms:</p>
00776 
00777 <p class=MsoNormal style='text-align:justify'>                        kTransformHaar
00778 </p>
00779 
00780 <p class=MsoNormal style='text-align:justify'>                        kTransformWalsh
00781 </p>
00782 
00783 <p class=MsoNormal style='text-align:justify'>                        kTransformCos
00784 </p>
00785 
00786 <p class=MsoNormal style='text-align:justify'>                        kTransformSin
00787 </p>
00788 
00789 <p class=MsoNormal style='text-align:justify'>                        kTransformFourier
00790 </p>
00791 
00792 <p class=MsoNormal style='text-align:justify'>                        kTransformHartey
00793 </p>
00794 
00795 <p class=MsoNormal style='text-align:justify'>            Mixed transforms:</p>
00796 
00797 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierWalsh
00798 </p>
00799 
00800 <p class=MsoNormal style='text-align:justify'>                        kTransformFourierHaar
00801 </p>
00802 
00803 <p class=MsoNormal style='text-align:justify'>                        kTransformWalshHaar
00804 </p>
00805 
00806 <p class=MsoNormal style='text-align:justify'>                        kTransformCosWalsh
00807 </p>
00808 
00809 <p class=MsoNormal style='text-align:justify'>                        kTransformCosHaar
00810 </p>
00811 
00812 <p class=MsoNormal style='text-align:justify'>                        kTransformSinWalsh
00813 </p>
00814 
00815 <p class=MsoNormal style='text-align:justify'>                        kTransformSinHaar
00816 </p>
00817 
00818 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDirection</b>-direction-transform
00819 direction (forward, inverse)</p>
00820 
00821 <p class=MsoNormal style='text-align:justify'>                        kTransformForward
00822 </p>
00823 
00824 <p class=MsoNormal style='text-align:justify'>                        kTransformInverse
00825 </p>
00826 
00827 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'><b>fDegree</b>-applies
00828 only for mixed transforms [2], [3], [4]. </p>
00829 
00830 <p class=MsoNormal style='text-align:justify;text-indent:22.8pt'>                
00831 <span style='color:fuchsia'> Allowed range  <sub><img border=0 width=100
00832 height=27 src="gif/spectrumtransform_transform_image001.gif"></sub>. </span></p>
00833 
00834 <p class=MsoNormal style='text-align:justify'><b><i>References:</i></b></p>
00835 
00836 <p class=MsoNormal style='text-align:justify'>[1] C.V. Hampton, B. Lian, Wm. C.
00837 McHarris: Fast-Fourier-transform spectral enhancement techniques for gamma-ray
00838 spectroscopy. NIM A353 (1994) 280-284. </p>
00839 
00840 <p class=MsoNormal style='text-align:justify'>[2] Morhá&#269; M., Matoušek V.,
00841 New adaptive Cosine-Walsh  transform and its application to nuclear data
00842 compression, IEEE Transactions on Signal Processing 48 (2000) 2693.  </p>
00843 
00844 <p class=MsoNormal style='text-align:justify'>[3] Morhá&#269; M., Matoušek V.,
00845 Data compression using new fast adaptive Cosine-Haar transforms, Digital Signal
00846 Processing 8 (1998) 63. </p>
00847 
00848 <p class=MsoNormal style='text-align:justify'>[4] Morhá&#269; M., Matoušek V.:
00849 Multidimensional nuclear data compression using fast adaptive Walsh-Haar
00850 transform. Acta Physica Slovaca 51 (2001) 307. </p>
00851 
00852 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00853 
00854 <p class=MsoNormal style='text-align:justify'><i>Example  – script Transform.c:</i></p>
00855 
00856 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'><img
00857 width=600 height=324 src="gif/spectrumtransform_transform_image002.jpg"></span></p>
00858 
00859 <p class=MsoNormal><b>Fig. 1 Original gamma-ray spectrum</b></p>
00860 
00861 <p class=MsoNormal><b><span style='font-size:14.0pt'><img border=0 width=601
00862 height=402 src="gif/spectrumtransform_transform_image003.jpg"></span></b></p>
00863 
00864 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'>&nbsp;</span></p>
00865 
00866 <p class=MsoNormal><b>Fig. 2 Transformed spectrum from Fig. 1 using Cosine
00867 transform</b></p>
00868 
00869 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>&nbsp;</span></b></p>
00870 
00871 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
00872 
00873 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
00874 Transform function (class TSpectrumTransform).</span></p>
00875 
00876 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
00877 do</span></p>
00878 
00879 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Transform.C</span></p>
00880 
00881 <p class=MsoNormal><span style='font-size:10.0pt'>   </span></p>
00882 
00883 <p class=MsoNormal><span style='font-size:10.0pt'>#include &lt;TSpectrum&gt;</span></p>
00884 
00885 <p class=MsoNormal><span style='font-size:10.0pt'>#include
00886 &lt;TSpectrumTransform&gt;</span></p>
00887 
00888 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
00889 
00890 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Transform() {</span></p>
00891 
00892 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
00893 
00894 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
00895 4096;</span></p>
00896 
00897 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
00898 0;</span></p>
00899 
00900 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
00901 (Double_t)nbins;</span></p>
00902 
00903 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
00904 style='font-size:10.0pt'>Float_t * source = new float[nbins];</span></p>
00905 
00906 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t * dest = new
00907 float[nbins];   </span></p>
00908 
00909 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new TH1F(&quot;h&quot;,&quot;Transformed
00910 spectrum using Cosine transform&quot;,nbins,xmin,xmax);</span></p>
00911 
00912 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
00913 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
00914 
00915 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
00916 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
00917 
00918 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
00919 i++) source[i]=h-&gt;GetBinContent(i + 1);         </span></p>
00920 
00921 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 =
00922 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
00923 
00924 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
00925 Transform1 = new
00926 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
00927 
00928 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
00929 TSpectrum();</span></p>
00930 
00931 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
00932 new TSpectrumTransform(4096);</span></p>
00933 
00934 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
00935 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
00936 
00937 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
00938 t-&gt;SetDirection(t-&gt;kTransformForward);</span></p>
00939 
00940 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
00941 style='font-size:10.0pt'>t-&gt;Transform(source,dest);</span></p>
00942 
00943 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
00944 i++) h-&gt;SetBinContent(i + 1,dest[i]);   </span></p>
00945 
00946 <p class=MsoNormal><span style='font-size:10.0pt'>  
00947 h-&gt;SetLineColor(kRed);      </span></p>
00948 
00949 <p class=MsoNormal><span style='font-size:10.0pt'>   h-&gt;Draw(&quot;L&quot;);</span></p>
00950 
00951 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
00952 
00953 </div>
00954 
00955 <!-- */
00956 // --> End_Html
00957    int i, j=0, k = 1, m, l;
00958    float val;
00959    double a, b, pi = 3.14159265358979323846;
00960    float *working_space = 0;
00961    if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
00962       if (fTransformType >= kTransformCosWalsh)
00963          fDegree += 1;
00964       k = (int) TMath::Power(2, fDegree);
00965       j = fSize / k;
00966    }
00967    switch (fTransformType) {
00968    case kTransformHaar:
00969    case kTransformWalsh:
00970       working_space = new float[2 * fSize];
00971       break;
00972    case kTransformCos:
00973    case kTransformSin:
00974    case kTransformFourier:
00975    case kTransformHartley:
00976    case kTransformFourierWalsh:
00977    case kTransformFourierHaar:
00978    case kTransformWalshHaar:
00979       working_space = new float[4 * fSize];
00980       break;
00981    case kTransformCosWalsh:
00982    case kTransformCosHaar:
00983    case kTransformSinWalsh:
00984    case kTransformSinHaar:
00985       working_space = new float[8 * fSize];
00986       break;
00987    }
00988    if (fDirection == kTransformForward) {
00989       switch (fTransformType) {
00990       case kTransformHaar:
00991          for (i = 0; i < fSize; i++) {
00992             working_space[i] = source[i];
00993          }
00994          Haar(working_space, fSize, fDirection);
00995          for (i = 0; i < fSize; i++) {
00996             destVector[i] = working_space[i];
00997          }
00998          break;
00999       case kTransformWalsh:
01000          for (i = 0; i < fSize; i++) {
01001             working_space[i] = source[i];
01002          }
01003          Walsh(working_space, fSize);
01004          BitReverse(working_space, fSize);
01005          for (i = 0; i < fSize; i++) {
01006             destVector[i] = working_space[i];
01007          }
01008          break;
01009       case kTransformCos:
01010          fSize = 2 * fSize;
01011          for (i = 1; i <= (fSize / 2); i++) {
01012             val = source[i - 1];
01013             working_space[i - 1] = val;
01014             working_space[fSize - i] = val;
01015          }
01016          Fourier(working_space, fSize, 0, kTransformForward, 0);
01017          for (i = 0; i < fSize / 2; i++) {
01018             a = pi * (double) i / (double) fSize;
01019             a = TMath::Cos(a);
01020             b = working_space[i];
01021             a = b / a;
01022             working_space[i] = a;
01023             working_space[i + fSize] = 0;
01024          } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
01025          for (i = 0; i < fSize / 2; i++) {
01026             destVector[i] = working_space[i];
01027          }
01028          break;
01029       case kTransformSin:
01030          fSize = 2 * fSize;
01031          for (i = 1; i <= (fSize / 2); i++) {
01032             val = source[i - 1];
01033             working_space[i - 1] = val;
01034             working_space[fSize - i] = -val;
01035          }
01036          Fourier(working_space, fSize, 0, kTransformForward, 0);
01037          for (i = 0; i < fSize / 2; i++) {
01038             a = pi * (double) i / (double) fSize;
01039             a = TMath::Sin(a);
01040             b = working_space[i];
01041             if (a != 0)
01042                a = b / a;
01043             working_space[i - 1] = a;
01044             working_space[i + fSize] = 0;
01045          }
01046          working_space[fSize / 2 - 1] =
01047              working_space[fSize / 2] / TMath::Sqrt(2.0);
01048          for (i = 0; i < fSize / 2; i++) {
01049             destVector[i] = working_space[i];
01050          }
01051          break;
01052       case kTransformFourier:
01053          for (i = 0; i < fSize; i++) {
01054             working_space[i] = source[i];
01055          }
01056          Fourier(working_space, fSize, 0, kTransformForward, 0);
01057          for (i = 0; i < 2 * fSize; i++) {
01058             destVector[i] = working_space[i];
01059          }
01060          break;
01061       case kTransformHartley:
01062          for (i = 0; i < fSize; i++) {
01063             working_space[i] = source[i];
01064          }
01065          Fourier(working_space, fSize, 1, kTransformForward, 0);
01066          for (i = 0; i < fSize; i++) {
01067             destVector[i] = working_space[i];
01068          }
01069          break;
01070       case kTransformFourierWalsh:
01071       case kTransformFourierHaar:
01072       case kTransformWalshHaar:
01073       case kTransformCosWalsh:
01074       case kTransformCosHaar:
01075       case kTransformSinWalsh:
01076       case kTransformSinHaar:
01077          for (i = 0; i < fSize; i++) {
01078             val = source[i];
01079             if (fTransformType == kTransformCosWalsh
01080                  || fTransformType == kTransformCosHaar) {
01081                j = (int) TMath::Power(2, fDegree) / 2;
01082                k = i / j;
01083                k = 2 * k * j;
01084                working_space[k + i % j] = val;
01085                working_space[k + 2 * j - 1 - i % j] = val;
01086             }
01087             
01088             else if (fTransformType == kTransformSinWalsh
01089                      || fTransformType == kTransformSinHaar) {
01090                j = (int) TMath::Power(2, fDegree) / 2;
01091                k = i / j;
01092                k = 2 * k * j;
01093                working_space[k + i % j] = val;
01094                working_space[k + 2 * j - 1 - i % j] = -val;
01095             }
01096             
01097             else
01098                working_space[i] = val;
01099          }
01100          if (fTransformType == kTransformFourierWalsh
01101               || fTransformType == kTransformFourierHaar
01102               || fTransformType == kTransformWalshHaar) {
01103             for (i = 0; i < j; i++)
01104                BitReverseHaar(working_space, fSize, k, i * k);
01105             GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
01106          }
01107          
01108          else if (fTransformType == kTransformCosWalsh
01109                   || fTransformType == kTransformCosHaar) {
01110             m = (int) TMath::Power(2, fDegree);
01111             l = 2 * fSize / m;
01112             for (i = 0; i < l; i++)
01113                BitReverseHaar(working_space, 2 * fSize, m, i * m);
01114             GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
01115             for (i = 0; i < fSize; i++) {
01116                k = i / j;
01117                k = 2 * k * j;
01118                a = pi * (double) (i % j) / (double) (2 * j);
01119                a = TMath::Cos(a);
01120                b = working_space[k + i % j];
01121                if (i % j == 0)
01122                   a = b / TMath::Sqrt(2.0);
01123                
01124                else
01125                   a = b / a;
01126                working_space[i] = a;
01127                working_space[i + 2 * fSize] = 0;
01128             }
01129          }
01130          
01131          else if (fTransformType == kTransformSinWalsh
01132                   || fTransformType == kTransformSinHaar) {
01133             m = (int) TMath::Power(2, fDegree);
01134             l = 2 * fSize / m;
01135             for (i = 0; i < l; i++)
01136                BitReverseHaar(working_space, 2 * fSize, m, i * m);
01137             GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
01138             for (i = 0; i < fSize; i++) {
01139                k = i / j;
01140                k = 2 * k * j;
01141                a = pi * (double) (i % j) / (double) (2 * j);
01142                a = TMath::Cos(a);
01143                b = working_space[j + k + i % j];
01144                if (i % j == 0)
01145                   a = b / TMath::Sqrt(2.0);
01146                
01147                else
01148                   a = b / a;
01149                working_space[j + k / 2 - i % j - 1] = a;
01150                working_space[i + 2 * fSize] = 0;
01151             }
01152          }
01153          if (fTransformType > kTransformWalshHaar)
01154             k = (int) TMath::Power(2, fDegree - 1);
01155          
01156          else
01157             k = (int) TMath::Power(2, fDegree);
01158          j = fSize / k;
01159          for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
01160             working_space[fSize + i] = working_space[l + i / j];
01161             working_space[fSize + i + 2 * fSize] =
01162                 working_space[l + i / j + 2 * fSize];
01163          }
01164          for (i = 0; i < fSize; i++) {
01165             working_space[i] = working_space[fSize + i];
01166             working_space[i + 2 * fSize] =
01167                 working_space[fSize + i + 2 * fSize];
01168          }
01169          for (i = 0; i < fSize; i++) {
01170             destVector[i] = working_space[i];
01171          }
01172          if (fTransformType == kTransformFourierWalsh
01173               || fTransformType == kTransformFourierHaar) {
01174             for (i = 0; i < fSize; i++) {
01175                destVector[fSize + i] = working_space[i + 2 * fSize];
01176             }
01177          }
01178          break;
01179       }
01180    }
01181    
01182    else if (fDirection == kTransformInverse) {
01183       switch (fTransformType) {
01184       case kTransformHaar:
01185          for (i = 0; i < fSize; i++) {
01186             working_space[i] = source[i];
01187          }
01188          Haar(working_space, fSize, fDirection);
01189          for (i = 0; i < fSize; i++) {
01190             destVector[i] = working_space[i];
01191          }
01192          break;
01193       case kTransformWalsh:
01194          for (i = 0; i < fSize; i++) {
01195             working_space[i] = source[i];
01196          }
01197          BitReverse(working_space, fSize);
01198          Walsh(working_space, fSize);
01199          for (i = 0; i < fSize; i++) {
01200             destVector[i] = working_space[i];
01201          }
01202          break;
01203       case kTransformCos:
01204          for (i = 0; i < fSize; i++) {
01205             working_space[i] = source[i];
01206          }
01207          fSize = 2 * fSize;
01208          working_space[0] = working_space[0] * TMath::Sqrt(2.0);
01209          for (i = 0; i < fSize / 2; i++) {
01210             a = pi * (double) i / (double) fSize;
01211             b = TMath::Sin(a);
01212             a = TMath::Cos(a);
01213             working_space[i + fSize] = (double) working_space[i] * b;
01214             working_space[i] = (double) working_space[i] * a;
01215          } for (i = 2; i <= (fSize / 2); i++) {
01216             working_space[fSize - i + 1] = working_space[i - 1];
01217             working_space[fSize - i + 1 + fSize] =
01218                 -working_space[i - 1 + fSize];
01219          }
01220          working_space[fSize / 2] = 0;
01221          working_space[fSize / 2 + fSize] = 0;
01222          Fourier(working_space, fSize, 0, kTransformInverse, 1);
01223          for (i = 0; i < fSize / 2; i++) {
01224             destVector[i] = working_space[i];
01225          }
01226          break;
01227       case kTransformSin:
01228          for (i = 0; i < fSize; i++) {
01229             working_space[i] = source[i];
01230          }
01231          fSize = 2 * fSize;
01232          working_space[fSize / 2] =
01233              working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
01234          for (i = fSize / 2 - 1; i > 0; i--) {
01235             a = pi * (double) i / (double) fSize;
01236             working_space[i + fSize] =
01237                 -(double) working_space[i - 1] * TMath::Cos(a);
01238             working_space[i] =
01239                 (double) working_space[i - 1] * TMath::Sin(a);
01240          } for (i = 2; i <= (fSize / 2); i++) {
01241             working_space[fSize - i + 1] = working_space[i - 1];
01242             working_space[fSize - i + 1 + fSize] =
01243                 -working_space[i - 1 + fSize];
01244          }
01245          working_space[0] = 0;
01246          working_space[fSize] = 0;
01247          working_space[fSize / 2 + fSize] = 0;
01248          Fourier(working_space, fSize, 0, kTransformInverse, 0);
01249          for (i = 0; i < fSize / 2; i++) {
01250             destVector[i] = working_space[i];
01251          }
01252          break;
01253       case kTransformFourier:
01254          for (i = 0; i < 2 * fSize; i++) {
01255             working_space[i] = source[i];
01256          }
01257          Fourier(working_space, fSize, 0, kTransformInverse, 0);
01258          for (i = 0; i < fSize; i++) {
01259             destVector[i] = working_space[i];
01260          }
01261          break;
01262       case kTransformHartley:
01263          for (i = 0; i < fSize; i++) {
01264             working_space[i] = source[i];
01265          }
01266          Fourier(working_space, fSize, 1, kTransformInverse, 0);
01267          for (i = 0; i < fSize; i++) {
01268             destVector[i] = working_space[i];
01269          }
01270          break;
01271       case kTransformFourierWalsh:
01272       case kTransformFourierHaar:
01273       case kTransformWalshHaar:
01274       case kTransformCosWalsh:
01275       case kTransformCosHaar:
01276       case kTransformSinWalsh:
01277       case kTransformSinHaar:
01278          for (i = 0; i < fSize; i++) {
01279             working_space[i] = source[i];
01280          }
01281          if (fTransformType == kTransformFourierWalsh
01282               || fTransformType == kTransformFourierHaar) {
01283             for (i = 0; i < fSize; i++) {
01284                working_space[i + 2 * fSize] = source[fSize + i];
01285             }
01286          }
01287          if (fTransformType > kTransformWalshHaar)
01288             k = (int) TMath::Power(2, fDegree - 1);
01289          
01290          else
01291             k = (int) TMath::Power(2, fDegree);
01292          j = fSize / k;
01293          for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
01294             working_space[fSize + l + i / j] = working_space[i];
01295             working_space[fSize + l + i / j + 2 * fSize] =
01296                 working_space[i + 2 * fSize];
01297          }
01298          for (i = 0; i < fSize; i++) {
01299             working_space[i] = working_space[fSize + i];
01300             working_space[i + 2 * fSize] =
01301                 working_space[fSize + i + 2 * fSize];
01302          }
01303          if (fTransformType == kTransformFourierWalsh
01304               || fTransformType == kTransformFourierHaar
01305               || fTransformType == kTransformWalshHaar) {
01306             GeneralInv(working_space, fSize, fDegree, fTransformType);
01307             for (i = 0; i < j; i++)
01308                BitReverseHaar(working_space, fSize, k, i * k);
01309          }
01310          
01311          else if (fTransformType == kTransformCosWalsh
01312                   || fTransformType == kTransformCosHaar) {
01313             j = (int) TMath::Power(2, fDegree) / 2;
01314             m = (int) TMath::Power(2, fDegree);
01315             l = 2 * fSize / m;
01316             for (i = 0; i < fSize; i++) {
01317                k = i / j;
01318                k = 2 * k * j;
01319                a = pi * (double) (i % j) / (double) (2 * j);
01320                if (i % j == 0) {
01321                   working_space[2 * fSize + k + i % j] =
01322                       working_space[i] * TMath::Sqrt(2.0);
01323                   working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
01324                }
01325                
01326                else {
01327                   b = TMath::Sin(a);
01328                   a = TMath::Cos(a);
01329                   working_space[4 * fSize + 2 * fSize + k + i % j] =
01330                       -(double) working_space[i] * b;
01331                   working_space[2 * fSize + k + i % j] =
01332                       (double) working_space[i] * a;
01333             } } for (i = 0; i < fSize; i++) {
01334                k = i / j;
01335                k = 2 * k * j;
01336                if (i % j == 0) {
01337                   working_space[2 * fSize + k + j] = 0;
01338                   working_space[4 * fSize + 2 * fSize + k + j] = 0;
01339                }
01340                
01341                else {
01342                   working_space[2 * fSize + k + 2 * j - i % j] =
01343                       working_space[2 * fSize + k + i % j];
01344                   working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
01345                       -working_space[4 * fSize + 2 * fSize + k + i % j];
01346                }
01347             }
01348             for (i = 0; i < 2 * fSize; i++) {
01349                working_space[i] = working_space[2 * fSize + i];
01350                working_space[4 * fSize + i] =
01351                    working_space[4 * fSize + 2 * fSize + i];
01352             }
01353             GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
01354             m = (int) TMath::Power(2, fDegree);
01355             l = 2 * fSize / m;
01356             for (i = 0; i < l; i++)
01357                BitReverseHaar(working_space, 2 * fSize, m, i * m);
01358          }
01359          
01360          else if (fTransformType == kTransformSinWalsh
01361                   || fTransformType == kTransformSinHaar) {
01362             j = (int) TMath::Power(2, fDegree) / 2;
01363             m = (int) TMath::Power(2, fDegree);
01364             l = 2 * fSize / m;
01365             for (i = 0; i < fSize; i++) {
01366                k = i / j;
01367                k = 2 * k * j;
01368                a = pi * (double) (i % j) / (double) (2 * j);
01369                if (i % j == 0) {
01370                   working_space[2 * fSize + k + j + i % j] =
01371                       working_space[j + k / 2 - i % j -
01372                                     1] * TMath::Sqrt(2.0);
01373                   working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
01374                }
01375                
01376                else {
01377                   b = TMath::Sin(a);
01378                   a = TMath::Cos(a);
01379                   working_space[4 * fSize + 2 * fSize + k + j + i % j] =
01380                       -(double) working_space[j + k / 2 - i % j - 1] * b;
01381                   working_space[2 * fSize + k + j + i % j] =
01382                       (double) working_space[j + k / 2 - i % j - 1] * a;
01383             } } for (i = 0; i < fSize; i++) {
01384                k = i / j;
01385                k = 2 * k * j;
01386                if (i % j == 0) {
01387                   working_space[2 * fSize + k] = 0;
01388                   working_space[4 * fSize + 2 * fSize + k] = 0;
01389                }
01390                
01391                else {
01392                   working_space[2 * fSize + k + i % j] =
01393                       working_space[2 * fSize + k + 2 * j - i % j];
01394                   working_space[4 * fSize + 2 * fSize + k + i % j] =
01395                       -working_space[4 * fSize + 2 * fSize + k + 2 * j -
01396                                      i % j];
01397                }
01398             }
01399             for (i = 0; i < 2 * fSize; i++) {
01400                working_space[i] = working_space[2 * fSize + i];
01401                working_space[4 * fSize + i] =
01402                    working_space[4 * fSize + 2 * fSize + i];
01403             }
01404             GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
01405             for (i = 0; i < l; i++)
01406                BitReverseHaar(working_space, 2 * fSize, m, i * m);
01407          }
01408          for (i = 0; i < fSize; i++) {
01409             if (fTransformType >= kTransformCosWalsh) {
01410                k = i / j;
01411                k = 2 * k * j;
01412                val = working_space[k + i % j];
01413             }
01414             
01415             else
01416                val = working_space[i];
01417             destVector[i] = val;
01418          }
01419          break;
01420       }
01421    }
01422    delete[]working_space;
01423    return;
01424 }
01425 
01426 //////////FilterZonal FUNCTION - CALCULATES DIFFERENT 1-D ORTHOGONAL TRANSFORMS, SETS GIVEN REGION TO FILTER COEFFICIENT AND TRANSFORMS IT BACK//////
01427 
01428 //______________________________________________________________________________
01429 void TSpectrumTransform::FilterZonal(const float *source, float *destVector)
01430 {   
01431 ////////////////////////////////////////////////////////////////////////////////
01432 //        ONE-DIMENSIONAL FILTER ZONAL FUNCTION                               
01433 //        This function transforms the source spectrum. The calling program  
01434 //        should fill in input parameters. Then it sets transformed          
01435 //        coefficients in the given region (fXmin, fXmax) to the given         
01436 //        fFilterCoeff and transforms it back.
01437 //        Filtered data are written into dest spectrum.                     
01438 //                                                                           
01439 //        Function parameters:                                               
01440 //        source-pointer to the vector of source spectrum, its length should 
01441 //             be size except for inverse FOURIER, FOUR-WALSH, FOUR-HAAR    
01442 //             transform. These need 2*size length to supply real and       
01443 //             imaginary coefficients.                                      
01444 //        destVector-pointer to the vector of dest data, its length should be
01445 //           size except for direct FOURIER, FOUR-WALSH, FOUR-HAAR. These  
01446 //           need 2*size length to store real and imaginary coefficients   
01447 //                                                                          
01448 ////////////////////////////////////////////////////////////////////////////////
01449 //       
01450 //Begin_Html <!--
01451 /* -->
01452 <div class=Section2>
01453 
01454 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of zonal filtering</span></b></p>
01455 
01456 <p class=MsoNormal><i>&nbsp;</i></p>
01457 
01458 <p class=MsoNormal><i>Function:</i></p>
01459 
01460 <p class=MsoNormal><b>void TSpectrumTransform::FilterZonal(const <a
01461 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> *fSource,
01462 <a href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> *fDest)</b></p>
01463 
01464 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
01465 
01466 <p class=MsoNormal style='text-align:justify'>This function transforms the
01467 source spectrum (for details see Transform function). Before the FilterZonal
01468 function is called the class must be created by constructor and the type of the
01469 transform as well as some other parameters must be set using a set of setter
01470 functions. The FilterZonal function sets transformed coefficients in the given
01471 region (fXmin, fXmax) to the given fFilterCoeff and transforms it back.
01472 Filtered data are written into dest spectrum. </p>
01473 
01474 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'>&nbsp;</span></i></p>
01475 
01476 <p class=MsoNormal style='text-align:justify'><i>Example – script Filter.c:</i></p>
01477 
01478 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
01479 border=0 width=601 height=402 src="gif/spectrumtransform_filter_image001.jpg"></span></i></p>
01480 
01481 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum
01482 (black line) and filtered spectrum (red line) using Cosine transform and zonal
01483 filtration (channels 2048-4095 were set to 0) </b></p>
01484 
01485 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
01486 
01487 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
01488 
01489 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
01490 FilterZonal function (class TSpectrumTransform).</span></p>
01491 
01492 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
01493 do</span></p>
01494 
01495 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Filter.C</span></p>
01496 
01497 <p class=MsoNormal><span style='font-size:10.0pt'>   </span></p>
01498 
01499 <p class=MsoNormal><span style='font-size:10.0pt'>#include &lt;TSpectrum&gt;</span></p>
01500 
01501 <p class=MsoNormal><span style='font-size:10.0pt'>#include
01502 &lt;TSpectrumTransform&gt;</span></p>
01503 
01504 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
01505 
01506 <p class=MsoNormal><span style='font-size:10.0pt'>void Filter() {</span></p>
01507 
01508 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
01509 
01510 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
01511 4096;</span></p>
01512 
01513 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
01514 0;</span></p>
01515 
01516 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
01517 (Double_t)nbins;</span></p>
01518 
01519 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
01520 style='font-size:10.0pt'>Float_t * source = new float[nbins];</span></p>
01521 
01522 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t * dest = new
01523 float[nbins];   </span></p>
01524 
01525 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new
01526 TH1F(&quot;h&quot;,&quot;Zonal filtering using Cosine
01527 transform&quot;,nbins,xmin,xmax);</span></p>
01528 
01529 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
01530 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);         </span></p>
01531 
01532 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
01533 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
01534 
01535 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
01536 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
01537 
01538 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
01539 i++) source[i]=h-&gt;GetBinContent(i + 1);     </span></p>
01540 
01541 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 =
01542 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
01543 
01544 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
01545 Transform1 = new
01546 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
01547 
01548 <p class=MsoNormal><span style='font-size:10.0pt'>  
01549 h-&gt;SetAxisRange(700,1024);   </span></p>
01550 
01551 <p class=MsoNormal><span style='font-size:10.0pt'>  
01552 h-&gt;Draw(&quot;L&quot;);   </span></p>
01553 
01554 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
01555 TSpectrum();</span></p>
01556 
01557 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
01558 new TSpectrumTransform(4096);</span></p>
01559 
01560 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
01561 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
01562 
01563 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01564 t-&gt;SetRegion(2048, 4095);</span></p>
01565 
01566 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01567 t-&gt;FilterZonal(source,dest);     </span></p>
01568 
01569 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
01570 style='font-size:10.0pt'>for (i = 0; i &lt; nbins; i++) d-&gt;SetBinContent(i +
01571 1,dest[i]);</span></p>
01572 
01573 <p class=MsoNormal><span style='font-size:10.0pt'>  
01574 d-&gt;SetLineColor(kRed);   </span></p>
01575 
01576 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
01577 L&quot;);</span></p>
01578 
01579 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
01580 
01581 </div>
01582 
01583 <!-- */
01584 // --> End_Html
01585    int i, j=0, k = 1, m, l;
01586    float val;
01587    float *working_space = 0;
01588    double a, b, pi = 3.14159265358979323846, old_area, new_area;
01589    if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
01590       if (fTransformType >= kTransformCosWalsh)
01591          fDegree += 1;
01592       k = (int) TMath::Power(2, fDegree);
01593       j = fSize / k;
01594    }
01595    switch (fTransformType) {
01596    case kTransformHaar:
01597    case kTransformWalsh:
01598       working_space = new float[2 * fSize];
01599       break;
01600    case kTransformCos:
01601    case kTransformSin:
01602    case kTransformFourier:
01603    case kTransformHartley:
01604    case kTransformFourierWalsh:
01605    case kTransformFourierHaar:
01606    case kTransformWalshHaar:
01607       working_space = new float[4 * fSize];
01608       break;
01609    case kTransformCosWalsh:
01610    case kTransformCosHaar:
01611    case kTransformSinWalsh:
01612    case kTransformSinHaar:
01613       working_space = new float[8 * fSize];
01614       break;
01615    }
01616    switch (fTransformType) {
01617    case kTransformHaar:
01618       for (i = 0; i < fSize; i++) {
01619          working_space[i] = source[i];
01620       }
01621       Haar(working_space, fSize, kTransformForward);
01622       break;
01623    case kTransformWalsh:
01624       for (i = 0; i < fSize; i++) {
01625          working_space[i] = source[i];
01626       }
01627       Walsh(working_space, fSize);
01628       BitReverse(working_space, fSize);
01629       break;
01630    case kTransformCos:
01631       fSize = 2 * fSize;
01632       for (i = 1; i <= (fSize / 2); i++) {
01633          val = source[i - 1];
01634          working_space[i - 1] = val;
01635          working_space[fSize - i] = val;
01636       }
01637       Fourier(working_space, fSize, 0, kTransformForward, 0);
01638       for (i = 0; i < fSize / 2; i++) {
01639          a = pi * (double) i / (double) fSize;
01640          a = TMath::Cos(a);
01641          b = working_space[i];
01642          a = b / a;
01643          working_space[i] = a;
01644          working_space[i + fSize] = 0;
01645       } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
01646       fSize = fSize / 2;
01647       break;
01648    case kTransformSin:
01649       fSize = 2 * fSize;
01650       for (i = 1; i <= (fSize / 2); i++) {
01651          val = source[i - 1];
01652          working_space[i - 1] = val;
01653          working_space[fSize - i] = -val;
01654       }
01655       Fourier(working_space, fSize, 0, kTransformForward, 0);
01656       for (i = 0; i < fSize / 2; i++) {
01657          a = pi * (double) i / (double) fSize;
01658          a = TMath::Sin(a);
01659          b = working_space[i];
01660          if (a != 0)
01661             a = b / a;
01662          working_space[i - 1] = a;
01663          working_space[i + fSize] = 0;
01664       }
01665       working_space[fSize / 2 - 1] =
01666           working_space[fSize / 2] / TMath::Sqrt(2.0);
01667       fSize = fSize / 2;
01668       break;
01669    case kTransformFourier:
01670       for (i = 0; i < fSize; i++) {
01671          working_space[i] = source[i];
01672       }
01673       Fourier(working_space, fSize, 0, kTransformForward, 0);
01674       break;
01675    case kTransformHartley:
01676       for (i = 0; i < fSize; i++) {
01677          working_space[i] = source[i];
01678       }
01679       Fourier(working_space, fSize, 1, kTransformForward, 0);
01680       break;
01681    case kTransformFourierWalsh:
01682    case kTransformFourierHaar:
01683    case kTransformWalshHaar:
01684    case kTransformCosWalsh:
01685    case kTransformCosHaar:
01686    case kTransformSinWalsh:
01687    case kTransformSinHaar:
01688       for (i = 0; i < fSize; i++) {
01689          val = source[i];
01690          if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
01691             j = (int) TMath::Power(2, fDegree) / 2;
01692             k = i / j;
01693             k = 2 * k * j;
01694             working_space[k + i % j] = val;
01695             working_space[k + 2 * j - 1 - i % j] = val;
01696          }
01697          
01698          else if (fTransformType == kTransformSinWalsh
01699                   || fTransformType == kTransformSinHaar) {
01700             j = (int) TMath::Power(2, fDegree) / 2;
01701             k = i / j;
01702             k = 2 * k * j;
01703             working_space[k + i % j] = val;
01704             working_space[k + 2 * j - 1 - i % j] = -val;
01705          }
01706          
01707          else
01708             working_space[i] = val;
01709       }
01710       if (fTransformType == kTransformFourierWalsh
01711            || fTransformType == kTransformFourierHaar
01712            || fTransformType == kTransformWalshHaar) {
01713          for (i = 0; i < j; i++)
01714             BitReverseHaar(working_space, fSize, k, i * k);
01715          GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
01716       }
01717       
01718       else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
01719          m = (int) TMath::Power(2, fDegree);
01720          l = 2 * fSize / m;
01721          for (i = 0; i < l; i++)
01722             BitReverseHaar(working_space, 2 * fSize, m, i * m);
01723          GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
01724          for (i = 0; i < fSize; i++) {
01725             k = i / j;
01726             k = 2 * k * j;
01727             a = pi * (double) (i % j) / (double) (2 * j);
01728             a = TMath::Cos(a);
01729             b = working_space[k + i % j];
01730             if (i % j == 0)
01731                a = b / TMath::Sqrt(2.0);
01732             
01733             else
01734                a = b / a;
01735             working_space[i] = a;
01736             working_space[i + 2 * fSize] = 0;
01737          }
01738       }
01739       
01740       else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
01741          m = (int) TMath::Power(2, fDegree);
01742          l = 2 * fSize / m;
01743          for (i = 0; i < l; i++)
01744             BitReverseHaar(working_space, 2 * fSize, m, i * m);
01745          GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
01746          for (i = 0; i < fSize; i++) {
01747             k = i / j;
01748             k = 2 * k * j;
01749             a = pi * (double) (i % j) / (double) (2 * j);
01750             a = TMath::Cos(a);
01751             b = working_space[j + k + i % j];
01752             if (i % j == 0)
01753                a = b / TMath::Sqrt(2.0);
01754             
01755             else
01756                a = b / a;
01757             working_space[j + k / 2 - i % j - 1] = a;
01758             working_space[i + 2 * fSize] = 0;
01759          }
01760       }
01761       if (fTransformType > kTransformWalshHaar)
01762          k = (int) TMath::Power(2, fDegree - 1);
01763       
01764       else
01765          k = (int) TMath::Power(2, fDegree);
01766       j = fSize / k;
01767       for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
01768          working_space[fSize + i] = working_space[l + i / j];
01769          working_space[fSize + i + 2 * fSize] =
01770              working_space[l + i / j + 2 * fSize];
01771       }
01772       for (i = 0; i < fSize; i++) {
01773          working_space[i] = working_space[fSize + i];
01774          working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
01775       }
01776       break;
01777    }
01778    if (!working_space) return;
01779    for (i = 0, old_area = 0; i < fSize; i++) {
01780       old_area += working_space[i];
01781    }
01782    for (i = 0, new_area = 0; i < fSize; i++) {
01783       if (i >= fXmin && i <= fXmax)
01784          working_space[i] = fFilterCoeff;
01785       new_area += working_space[i];
01786    }
01787    if (new_area != 0) {
01788       a = old_area / new_area;
01789       for (i = 0; i < fSize; i++) {
01790          working_space[i] *= a;
01791       }
01792    }
01793    if (fTransformType == kTransformFourier) {
01794       for (i = 0, old_area = 0; i < fSize; i++) {
01795          old_area += working_space[i + fSize];
01796       }
01797       for (i = 0, new_area = 0; i < fSize; i++) {
01798          if (i >= fXmin && i <= fXmax)
01799             working_space[i + fSize] = fFilterCoeff;
01800          new_area += working_space[i + fSize];
01801       }
01802       if (new_area != 0) {
01803          a = old_area / new_area;
01804          for (i = 0; i < fSize; i++) {
01805             working_space[i + fSize] *= a;
01806          }
01807       }
01808    }
01809    
01810    else if (fTransformType == kTransformFourierWalsh
01811             || fTransformType == kTransformFourierHaar) {
01812       for (i = 0, old_area = 0; i < fSize; i++) {
01813          old_area += working_space[i + 2 * fSize];
01814       }
01815       for (i = 0, new_area = 0; i < fSize; i++) {
01816          if (i >= fXmin && i <= fXmax)
01817             working_space[i + 2 * fSize] = fFilterCoeff;
01818          new_area += working_space[i + 2 * fSize];
01819       }
01820       if (new_area != 0) {
01821          a = old_area / new_area;
01822          for (i = 0; i < fSize; i++) {
01823             working_space[i + 2 * fSize] *= a;
01824          }
01825       }
01826    }
01827    switch (fTransformType) {
01828    case kTransformHaar:
01829       Haar(working_space, fSize, kTransformInverse);
01830       for (i = 0; i < fSize; i++) {
01831          destVector[i] = working_space[i];
01832       }
01833       break;
01834    case kTransformWalsh:
01835       BitReverse(working_space, fSize);
01836       Walsh(working_space, fSize);
01837       for (i = 0; i < fSize; i++) {
01838          destVector[i] = working_space[i];
01839       }
01840       break;
01841    case kTransformCos:
01842       fSize = 2 * fSize;
01843       working_space[0] = working_space[0] * TMath::Sqrt(2.0);
01844       for (i = 0; i < fSize / 2; i++) {
01845          a = pi * (double) i / (double) fSize;
01846          b = TMath::Sin(a);
01847          a = TMath::Cos(a);
01848          working_space[i + fSize] = (double) working_space[i] * b;
01849          working_space[i] = (double) working_space[i] * a;
01850       } for (i = 2; i <= (fSize / 2); i++) {
01851          working_space[fSize - i + 1] = working_space[i - 1];
01852          working_space[fSize - i + 1 + fSize] =
01853              -working_space[i - 1 + fSize];
01854       }
01855       working_space[fSize / 2] = 0;
01856       working_space[fSize / 2 + fSize] = 0;
01857       Fourier(working_space, fSize, 0, kTransformInverse, 1);
01858       for (i = 0; i < fSize / 2; i++) {
01859          destVector[i] = working_space[i];
01860       }
01861       break;
01862    case kTransformSin:
01863       fSize = 2 * fSize;
01864       working_space[fSize / 2] =
01865           working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
01866       for (i = fSize / 2 - 1; i > 0; i--) {
01867          a = pi * (double) i / (double) fSize;
01868          working_space[i + fSize] =
01869              -(double) working_space[i - 1] * TMath::Cos(a);
01870          working_space[i] = (double) working_space[i - 1] * TMath::Sin(a);
01871       } for (i = 2; i <= (fSize / 2); i++) {
01872          working_space[fSize - i + 1] = working_space[i - 1];
01873          working_space[fSize - i + 1 + fSize] =
01874              -working_space[i - 1 + fSize];
01875       }
01876       working_space[0] = 0;
01877       working_space[fSize] = 0;
01878       working_space[fSize / 2 + fSize] = 0;
01879       Fourier(working_space, fSize, 0, kTransformInverse, 0);
01880       for (i = 0; i < fSize / 2; i++) {
01881          destVector[i] = working_space[i];
01882       }
01883       break;
01884    case kTransformFourier:
01885       Fourier(working_space, fSize, 0, kTransformInverse, 0);
01886       for (i = 0; i < fSize; i++) {
01887          destVector[i] = working_space[i];
01888       }
01889       break;
01890    case kTransformHartley:
01891       Fourier(working_space, fSize, 1, kTransformInverse, 0);
01892       for (i = 0; i < fSize; i++) {
01893          destVector[i] = working_space[i];
01894       }
01895       break;
01896    case kTransformFourierWalsh:
01897    case kTransformFourierHaar:
01898    case kTransformWalshHaar:
01899    case kTransformCosWalsh:
01900    case kTransformCosHaar:
01901    case kTransformSinWalsh:
01902    case kTransformSinHaar:
01903       if (fTransformType > kTransformWalshHaar)
01904          k = (int) TMath::Power(2, fDegree - 1);
01905       
01906       else
01907          k = (int) TMath::Power(2, fDegree);
01908       j = fSize / k;
01909       for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
01910          working_space[fSize + l + i / j] = working_space[i];
01911          working_space[fSize + l + i / j + 2 * fSize] =
01912              working_space[i + 2 * fSize];
01913       }
01914       for (i = 0; i < fSize; i++) {
01915          working_space[i] = working_space[fSize + i];
01916          working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
01917       }
01918       if (fTransformType == kTransformFourierWalsh
01919            || fTransformType == kTransformFourierHaar
01920            || fTransformType == kTransformWalshHaar) {
01921          GeneralInv(working_space, fSize, fDegree, fTransformType);
01922          for (i = 0; i < j; i++)
01923             BitReverseHaar(working_space, fSize, k, i * k);
01924       }
01925       
01926       else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
01927          j = (int) TMath::Power(2, fDegree) / 2;
01928          m = (int) TMath::Power(2, fDegree);
01929          l = 2 * fSize / m;
01930          for (i = 0; i < fSize; i++) {
01931             k = i / j;
01932             k = 2 * k * j;
01933             a = pi * (double) (i % j) / (double) (2 * j);
01934             if (i % j == 0) {
01935                working_space[2 * fSize + k + i % j] =
01936                    working_space[i] * TMath::Sqrt(2.0);
01937                working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
01938             }
01939             
01940             else {
01941                b = TMath::Sin(a);
01942                a = TMath::Cos(a);
01943                working_space[4 * fSize + 2 * fSize + k + i % j] =
01944                    -(double) working_space[i] * b;
01945                working_space[2 * fSize + k + i % j] =
01946                    (double) working_space[i] * a;
01947          } } for (i = 0; i < fSize; i++) {
01948             k = i / j;
01949             k = 2 * k * j;
01950             if (i % j == 0) {
01951                working_space[2 * fSize + k + j] = 0;
01952                working_space[4 * fSize + 2 * fSize + k + j] = 0;
01953             }
01954             
01955             else {
01956                working_space[2 * fSize + k + 2 * j - i % j] =
01957                    working_space[2 * fSize + k + i % j];
01958                working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
01959                    -working_space[4 * fSize + 2 * fSize + k + i % j];
01960             }
01961          }
01962          for (i = 0; i < 2 * fSize; i++) {
01963             working_space[i] = working_space[2 * fSize + i];
01964             working_space[4 * fSize + i] =
01965                 working_space[4 * fSize + 2 * fSize + i];
01966          }
01967          GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
01968          m = (int) TMath::Power(2, fDegree);
01969          l = 2 * fSize / m;
01970          for (i = 0; i < l; i++)
01971             BitReverseHaar(working_space, 2 * fSize, m, i * m);
01972       }
01973       
01974       else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
01975          j = (int) TMath::Power(2, fDegree) / 2;
01976          m = (int) TMath::Power(2, fDegree);
01977          l = 2 * fSize / m;
01978          for (i = 0; i < fSize; i++) {
01979             k = i / j;
01980             k = 2 * k * j;
01981             a = pi * (double) (i % j) / (double) (2 * j);
01982             if (i % j == 0) {
01983                working_space[2 * fSize + k + j + i % j] =
01984                    working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
01985                working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
01986             }
01987             
01988             else {
01989                b = TMath::Sin(a);
01990                a = TMath::Cos(a);
01991                working_space[4 * fSize + 2 * fSize + k + j + i % j] =
01992                    -(double) working_space[j + k / 2 - i % j - 1] * b;
01993                working_space[2 * fSize + k + j + i % j] =
01994                    (double) working_space[j + k / 2 - i % j - 1] * a;
01995          } } for (i = 0; i < fSize; i++) {
01996             k = i / j;
01997             k = 2 * k * j;
01998             if (i % j == 0) {
01999                working_space[2 * fSize + k] = 0;
02000                working_space[4 * fSize + 2 * fSize + k] = 0;
02001             }
02002             
02003             else {
02004                working_space[2 * fSize + k + i % j] =
02005                    working_space[2 * fSize + k + 2 * j - i % j];
02006                working_space[4 * fSize + 2 * fSize + k + i % j] =
02007                    -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
02008             }
02009          }
02010          for (i = 0; i < 2 * fSize; i++) {
02011             working_space[i] = working_space[2 * fSize + i];
02012             working_space[4 * fSize + i] =
02013                 working_space[4 * fSize + 2 * fSize + i];
02014          }
02015          GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
02016          for (i = 0; i < l; i++)
02017             BitReverseHaar(working_space, 2 * fSize, m, i * m);
02018       }
02019       for (i = 0; i < fSize; i++) {
02020          if (fTransformType >= kTransformCosWalsh) {
02021             k = i / j;
02022             k = 2 * k * j;
02023             val = working_space[k + i % j];
02024          }
02025          
02026          else
02027             val = working_space[i];
02028          destVector[i] = val;
02029       }
02030       break;
02031    }
02032    delete[]working_space;
02033    return;
02034 }
02035 
02036 //////////ENHANCE FUNCTION - CALCULATES DIFFERENT 1-D ORTHOGONAL TRANSFORMS, MULTIPLIES GIVEN REGION BY ENHANCE COEFFICIENT AND TRANSFORMS IT BACK//////
02037 //___________________________________________________________________________
02038 void TSpectrumTransform::Enhance(const float *source, float *destVector) 
02039 {   
02040 ////////////////////////////////////////////////////////////////////////////////
02041 //        ONE-DIMENSIONAL ENHANCE ZONAL FUNCTION                             
02042 //        This function transforms the source spectrum. The calling program  
02043 //      should fill in input parameters. Then it multiplies transformed      
02044 //      coefficients in the given region (fXmin, fXmax) by the given          
02045 //      fEnhanceCoeff and transforms it back                                   
02046 //        Processed data are written into dest spectrum.                      
02047 //                                                                            
02048 //        Function parameters:                                                
02049 //        source-pointer to the vector of source spectrum, its length should  
02050 //             be size except for inverse FOURIER, FOUR-WALSh, FOUR-HAAR      
02051 //             transform. These need 2*size length to supply real and         
02052 //             imaginary coefficients.                                        
02053 //        destVector-pointer to the vector of dest data, its length should be 
02054 //           size except for direct FOURIER, FOUR-WALSh, FOUR-HAAR. These     
02055 //           need 2*size length to store real and imaginary coefficients      
02056 //                                                                           
02057 ////////////////////////////////////////////////////////////////////////////////
02058 //Begin_Html <!--
02059 /* -->
02060 <div class=Section3>
02061 
02062 <p class=MsoNormal><b><span style='font-size:14.0pt'>Example of enhancement</span></b></p>
02063 
02064 <p class=MsoNormal><i>&nbsp;</i></p>
02065 
02066 <p class=MsoNormal><i>Function:</i></p>
02067 
02068 <p class=MsoNormal><b>void TSpectrumTransform::Enhance(const <a
02069 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> *fSource,
02070 <a href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a>
02071 *fDest)</b></p>
02072 
02073 <p class=MsoNormal><b>&nbsp;</b></p>
02074 
02075 <p class=MsoNormal style='text-align:justify'>This function transforms the
02076 source spectrum (for details see Transform function). Before the Enhance
02077 function is called the class must be created by constructor and the type of the
02078 transform as well as some other parameters must be set using a set of setter functions.
02079 The Enhance function multiplies transformed coefficients in the given region
02080 (fXmin, fXmax) by the given fEnhancCoeff and transforms it back. Enhanced data
02081 are written into dest spectrum.</p>
02082 
02083 <p class=MsoNormal>&nbsp;</p>
02084 
02085 <p class=MsoNormal style='text-align:justify'><i>Example  – script Enhance.c:</i></p>
02086 
02087 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
02088 border=0 width=601 height=402 src="gif/spectrumtransform_enhance_image001.jpg"></span></i></p>
02089 
02090 <p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'>&nbsp;</span></p>
02091 
02092 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum (black
02093 line) and enhanced spectrum (red line) using Cosine transform (channels 0-1024
02094 were multiplied by 2) </b></p>
02095 
02096 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
02097 
02098 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
02099 
02100 <p class=MsoNormal>&nbsp;</p>
02101 
02102 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
02103 Enhance function (class TSpectrumTransform).</span></p>
02104 
02105 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
02106 do</span></p>
02107 
02108 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x Enhance.C</span></p>
02109 
02110 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
02111 
02112 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void Enhance() {</span></p>
02113 
02114 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t i;</span></p>
02115 
02116 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t nbins =
02117 4096;</span></p>
02118 
02119 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmin  =
02120 0;</span></p>
02121 
02122 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t xmax  =
02123 (Double_t)nbins;</span></p>
02124 
02125 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
02126 style='font-size:10.0pt'>Float_t * source = new float[nbins];</span></p>
02127 
02128 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t * dest = new
02129 float[nbins];   </span></p>
02130 
02131 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new
02132 TH1F(&quot;h&quot;,&quot;Enhancement using Cosine transform&quot;,nbins,xmin,xmax);</span></p>
02133 
02134 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
02135 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);         </span></p>
02136 
02137 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
02138 TFile(&quot;spectra\\TSpectrum.root&quot;);</span></p>
02139 
02140 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
02141 f-&gt;Get(&quot;transform1;1&quot;);   </span></p>
02142 
02143 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
02144 i++) source[i]=h-&gt;GetBinContent(i + 1);     </span></p>
02145 
02146 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Transform1 = gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Transform1&quot;);</span></p>
02147 
02148 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Transform1)
02149 Transform1 = new
02150 TCanvas(&quot;Transform&quot;,&quot;Transform1&quot;,10,10,1000,700);</span></p>
02151 
02152 <p class=MsoNormal><span style='font-size:10.0pt'>  
02153 h-&gt;SetAxisRange(700,1024);   </span></p>
02154 
02155 <p class=MsoNormal><span style='font-size:10.0pt'>  
02156 h-&gt;Draw(&quot;L&quot;);   </span></p>
02157 
02158 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
02159 TSpectrum();</span></p>
02160 
02161 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumTransform *t =
02162 new TSpectrumTransform(4096);</span></p>
02163 
02164 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
02165 style='font-size:10.0pt'>t-&gt;SetTransformType(t-&gt;kTransformCos,0);</span></p>
02166 
02167 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   t-&gt;SetRegion(0,
02168 1024);</span></p>
02169 
02170 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
02171 t-&gt;SetEnhanceCoeff(2);</span></p>
02172 
02173 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
02174 t-&gt;Enhance(source,dest);        </span></p>
02175 
02176 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
02177 style='font-size:10.0pt'>for (i = 0; i &lt; nbins; i++) d-&gt;SetBinContent(i +
02178 1,dest[i]);</span></p>
02179 
02180 <p class=MsoNormal><span style='font-size:10.0pt'>  
02181 d-&gt;SetLineColor(kRed);   </span></p>
02182 
02183 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
02184 L&quot;);</span></p>
02185 
02186 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
02187 
02188 </div>
02189 
02190 <!-- */
02191 // --> End_Html
02192    int i, j=0, k = 1, m, l;
02193    float val;
02194    float *working_space = 0;
02195    double a, b, pi = 3.14159265358979323846, old_area, new_area;
02196    if (fTransformType >= kTransformFourierWalsh && fTransformType <= kTransformSinHaar) {
02197       if (fTransformType >= kTransformCosWalsh)
02198          fDegree += 1;
02199       k = (int) TMath::Power(2, fDegree);
02200       j = fSize / k;
02201    }
02202    switch (fTransformType) {
02203    case kTransformHaar:
02204    case kTransformWalsh:
02205       working_space = new float[2 * fSize];
02206       break;
02207    case kTransformCos:
02208    case kTransformSin:
02209    case kTransformFourier:
02210    case kTransformHartley:
02211    case kTransformFourierWalsh:
02212    case kTransformFourierHaar:
02213    case kTransformWalshHaar:
02214       working_space = new float[4 * fSize];
02215       break;
02216    case kTransformCosWalsh:
02217    case kTransformCosHaar:
02218    case kTransformSinWalsh:
02219    case kTransformSinHaar:
02220       working_space = new float[8 * fSize];
02221       break;
02222    }
02223    switch (fTransformType) {
02224    case kTransformHaar:
02225       for (i = 0; i < fSize; i++) {
02226          working_space[i] = source[i];
02227       }
02228       Haar(working_space, fSize, kTransformForward);
02229       break;
02230    case kTransformWalsh:
02231       for (i = 0; i < fSize; i++) {
02232          working_space[i] = source[i];
02233       }
02234       Walsh(working_space, fSize);
02235       BitReverse(working_space, fSize);
02236       break;
02237    case kTransformCos:
02238       fSize = 2 * fSize;
02239       for (i = 1; i <= (fSize / 2); i++) {
02240          val = source[i - 1];
02241          working_space[i - 1] = val;
02242          working_space[fSize - i] = val;
02243       }
02244       Fourier(working_space, fSize, 0, kTransformForward, 0);
02245       for (i = 0; i < fSize / 2; i++) {
02246          a = pi * (double) i / (double) fSize;
02247          a = TMath::Cos(a);
02248          b = working_space[i];
02249          a = b / a;
02250          working_space[i] = a;
02251          working_space[i + fSize] = 0;
02252       } working_space[0] = working_space[0] / TMath::Sqrt(2.0);
02253       fSize = fSize / 2;
02254       break;
02255    case kTransformSin:
02256       fSize = 2 * fSize;
02257       for (i = 1; i <= (fSize / 2); i++) {
02258          val = source[i - 1];
02259          working_space[i - 1] = val;
02260          working_space[fSize - i] = -val;
02261       }
02262       Fourier(working_space, fSize, 0, kTransformForward, 0);
02263       for (i = 0; i < fSize / 2; i++) {
02264          a = pi * (double) i / (double) fSize;
02265          a = TMath::Sin(a);
02266          b = working_space[i];
02267          if (a != 0)
02268             a = b / a;
02269          working_space[i - 1] = a;
02270          working_space[i + fSize] = 0;
02271       }
02272       working_space[fSize / 2 - 1] =
02273           working_space[fSize / 2] / TMath::Sqrt(2.0);
02274       fSize = fSize / 2;
02275       break;
02276    case kTransformFourier:
02277       for (i = 0; i < fSize; i++) {
02278          working_space[i] = source[i];
02279       }
02280       Fourier(working_space, fSize, 0, kTransformForward, 0);
02281       break;
02282    case kTransformHartley:
02283       for (i = 0; i < fSize; i++) {
02284          working_space[i] = source[i];
02285       }
02286       Fourier(working_space, fSize, 1, kTransformForward, 0);
02287       break;
02288    case kTransformFourierWalsh:
02289    case kTransformFourierHaar:
02290    case kTransformWalshHaar:
02291    case kTransformCosWalsh:
02292    case kTransformCosHaar:
02293    case kTransformSinWalsh:
02294    case kTransformSinHaar:
02295       for (i = 0; i < fSize; i++) {
02296          val = source[i];
02297          if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
02298             j = (int) TMath::Power(2, fDegree) / 2;
02299             k = i / j;
02300             k = 2 * k * j;
02301             working_space[k + i % j] = val;
02302             working_space[k + 2 * j - 1 - i % j] = val;
02303          }
02304          
02305          else if (fTransformType == kTransformSinWalsh
02306                   || fTransformType == kTransformSinHaar) {
02307             j = (int) TMath::Power(2, fDegree) / 2;
02308             k = i / j;
02309             k = 2 * k * j;
02310             working_space[k + i % j] = val;
02311             working_space[k + 2 * j - 1 - i % j] = -val;
02312          }
02313          
02314          else
02315             working_space[i] = val;
02316       }
02317       if (fTransformType == kTransformFourierWalsh
02318            || fTransformType == kTransformFourierHaar
02319            || fTransformType == kTransformWalshHaar) {
02320          for (i = 0; i < j; i++)
02321             BitReverseHaar(working_space, fSize, k, i * k);
02322          GeneralExe(working_space, 0, fSize, fDegree, fTransformType);
02323       }
02324       
02325       else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
02326          m = (int) TMath::Power(2, fDegree);
02327          l = 2 * fSize / m;
02328          for (i = 0; i < l; i++)
02329             BitReverseHaar(working_space, 2 * fSize, m, i * m);
02330          GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
02331          for (i = 0; i < fSize; i++) {
02332             k = i / j;
02333             k = 2 * k * j;
02334             a = pi * (double) (i % j) / (double) (2 * j);
02335             a = TMath::Cos(a);
02336             b = working_space[k + i % j];
02337             if (i % j == 0)
02338                a = b / TMath::Sqrt(2.0);
02339             
02340             else
02341                a = b / a;
02342             working_space[i] = a;
02343             working_space[i + 2 * fSize] = 0;
02344          }
02345       }
02346       
02347       else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
02348          m = (int) TMath::Power(2, fDegree);
02349          l = 2 * fSize / m;
02350          for (i = 0; i < l; i++)
02351             BitReverseHaar(working_space, 2 * fSize, m, i * m);
02352          GeneralExe(working_space, 0, 2 * fSize, fDegree, fTransformType);
02353          for (i = 0; i < fSize; i++) {
02354             k = i / j;
02355             k = 2 * k * j;
02356             a = pi * (double) (i % j) / (double) (2 * j);
02357             a = TMath::Cos(a);
02358             b = working_space[j + k + i % j];
02359             if (i % j == 0)
02360                a = b / TMath::Sqrt(2.0);
02361             
02362             else
02363                a = b / a;
02364             working_space[j + k / 2 - i % j - 1] = a;
02365             working_space[i + 2 * fSize] = 0;
02366          }
02367       }
02368       if (fTransformType > kTransformWalshHaar)
02369          k = (int) TMath::Power(2, fDegree - 1);
02370       
02371       else
02372          k = (int) TMath::Power(2, fDegree);
02373       j = fSize / k;
02374       for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
02375          working_space[fSize + i] = working_space[l + i / j];
02376          working_space[fSize + i + 2 * fSize] =
02377              working_space[l + i / j + 2 * fSize];
02378       }
02379       for (i = 0; i < fSize; i++) {
02380          working_space[i] = working_space[fSize + i];
02381          working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
02382       }
02383       break;
02384    }
02385    if (!working_space) return;
02386    for (i = 0, old_area = 0; i < fSize; i++) {
02387       old_area += working_space[i];
02388    }
02389    for (i = 0, new_area = 0; i < fSize; i++) {
02390       if (i >= fXmin && i <= fXmax)
02391          working_space[i] *= fEnhanceCoeff;
02392       new_area += working_space[i];
02393    }
02394    if (new_area != 0) {
02395       a = old_area / new_area;
02396       for (i = 0; i < fSize; i++) {
02397          working_space[i] *= a;
02398       }
02399    }
02400    if (fTransformType == kTransformFourier) {
02401       for (i = 0, old_area = 0; i < fSize; i++) {
02402          old_area += working_space[i + fSize];
02403       }
02404       for (i = 0, new_area = 0; i < fSize; i++) {
02405          if (i >= fXmin && i <= fXmax)
02406             working_space[i + fSize] *= fEnhanceCoeff;
02407          new_area += working_space[i + fSize];
02408       }
02409       if (new_area != 0) {
02410          a = old_area / new_area;
02411          for (i = 0; i < fSize; i++) {
02412             working_space[i + fSize] *= a;
02413          }
02414       }
02415    }
02416    
02417    else if (fTransformType == kTransformFourierWalsh
02418             || fTransformType == kTransformFourierHaar) {
02419       for (i = 0, old_area = 0; i < fSize; i++) {
02420          old_area += working_space[i + 2 * fSize];
02421       }
02422       for (i = 0, new_area = 0; i < fSize; i++) {
02423          if (i >= fXmin && i <= fXmax)
02424             working_space[i + 2 * fSize] *= fEnhanceCoeff;
02425          new_area += working_space[i + 2 * fSize];
02426       }
02427       if (new_area != 0) {
02428          a = old_area / new_area;
02429          for (i = 0; i < fSize; i++) {
02430             working_space[i + 2 * fSize] *= a;
02431          }
02432       }
02433    }
02434    switch (fTransformType) {
02435    case kTransformHaar:
02436       Haar(working_space, fSize, kTransformInverse);
02437       for (i = 0; i < fSize; i++) {
02438          destVector[i] = working_space[i];
02439       }
02440       break;
02441    case kTransformWalsh:
02442       BitReverse(working_space, fSize);
02443       Walsh(working_space, fSize);
02444       for (i = 0; i < fSize; i++) {
02445          destVector[i] = working_space[i];
02446       }
02447       break;
02448    case kTransformCos:
02449       fSize = 2 * fSize;
02450       working_space[0] = working_space[0] * TMath::Sqrt(2.0);
02451       for (i = 0; i < fSize / 2; i++) {
02452          a = pi * (double) i / (double) fSize;
02453          b = TMath::Sin(a);
02454          a = TMath::Cos(a);
02455          working_space[i + fSize] = (double) working_space[i] * b;
02456          working_space[i] = (double) working_space[i] * a;
02457       } for (i = 2; i <= (fSize / 2); i++) {
02458          working_space[fSize - i + 1] = working_space[i - 1];
02459          working_space[fSize - i + 1 + fSize] =
02460              -working_space[i - 1 + fSize];
02461       }
02462       working_space[fSize / 2] = 0;
02463       working_space[fSize / 2 + fSize] = 0;
02464       Fourier(working_space, fSize, 0, kTransformInverse, 1);
02465       for (i = 0; i < fSize / 2; i++) {
02466          destVector[i] = working_space[i];
02467       }
02468       break;
02469    case kTransformSin:
02470       fSize = 2 * fSize;
02471       working_space[fSize / 2] =
02472           working_space[fSize / 2 - 1] * TMath::Sqrt(2.0);
02473       for (i = fSize / 2 - 1; i > 0; i--) {
02474          a = pi * (double) i / (double) fSize;
02475          working_space[i + fSize] =
02476              -(double) working_space[i - 1] * TMath::Cos(a);
02477          working_space[i] = (double) working_space[i - 1] * TMath::Sin(a);
02478       } for (i = 2; i <= (fSize / 2); i++) {
02479          working_space[fSize - i + 1] = working_space[i - 1];
02480          working_space[fSize - i + 1 + fSize] =
02481              -working_space[i - 1 + fSize];
02482       }
02483       working_space[0] = 0;
02484       working_space[fSize] = 0;
02485       working_space[fSize / 2 + fSize] = 0;
02486       Fourier(working_space, fSize, 0, kTransformInverse, 0);
02487       for (i = 0; i < fSize / 2; i++) {
02488          destVector[i] = working_space[i];
02489       }
02490       break;
02491    case kTransformFourier:
02492       Fourier(working_space, fSize, 0, kTransformInverse, 0);
02493       for (i = 0; i < fSize; i++) {
02494          destVector[i] = working_space[i];
02495       }
02496       break;
02497    case kTransformHartley:
02498       Fourier(working_space, fSize, 1, kTransformInverse, 0);
02499       for (i = 0; i < fSize; i++) {
02500          destVector[i] = working_space[i];
02501       }
02502       break;
02503    case kTransformFourierWalsh:
02504    case kTransformFourierHaar:
02505    case kTransformWalshHaar:
02506    case kTransformCosWalsh:
02507    case kTransformCosHaar:
02508    case kTransformSinWalsh:
02509    case kTransformSinHaar:
02510       if (fTransformType > kTransformWalshHaar)
02511          k = (int) TMath::Power(2, fDegree - 1);
02512       
02513       else
02514          k = (int) TMath::Power(2, fDegree);
02515       j = fSize / k;
02516       for (i = 0, l = 0; i < fSize; i++, l = (l + k) % fSize) {
02517          working_space[fSize + l + i / j] = working_space[i];
02518          working_space[fSize + l + i / j + 2 * fSize] =
02519              working_space[i + 2 * fSize];
02520       }
02521       for (i = 0; i < fSize; i++) {
02522          working_space[i] = working_space[fSize + i];
02523          working_space[i + 2 * fSize] = working_space[fSize + i + 2 * fSize];
02524       }
02525       if (fTransformType == kTransformFourierWalsh
02526            || fTransformType == kTransformFourierHaar
02527            || fTransformType == kTransformWalshHaar) {
02528          GeneralInv(working_space, fSize, fDegree, fTransformType);
02529          for (i = 0; i < j; i++)
02530             BitReverseHaar(working_space, fSize, k, i * k);
02531       }
02532       
02533       else if (fTransformType == kTransformCosWalsh || fTransformType == kTransformCosHaar) {
02534          j = (int) TMath::Power(2, fDegree) / 2;
02535          m = (int) TMath::Power(2, fDegree);
02536          l = 2 * fSize / m;
02537          for (i = 0; i < fSize; i++) {
02538             k = i / j;
02539             k = 2 * k * j;
02540             a = pi * (double) (i % j) / (double) (2 * j);
02541             if (i % j == 0) {
02542                working_space[2 * fSize + k + i % j] =
02543                    working_space[i] * TMath::Sqrt(2.0);
02544                working_space[4 * fSize + 2 * fSize + k + i % j] = 0;
02545             }
02546             
02547             else {
02548                b = TMath::Sin(a);
02549                a = TMath::Cos(a);
02550                working_space[4 * fSize + 2 * fSize + k + i % j] =
02551                    -(double) working_space[i] * b;
02552                working_space[2 * fSize + k + i % j] =
02553                    (double) working_space[i] * a;
02554          } } for (i = 0; i < fSize; i++) {
02555             k = i / j;
02556             k = 2 * k * j;
02557             if (i % j == 0) {
02558                working_space[2 * fSize + k + j] = 0;
02559                working_space[4 * fSize + 2 * fSize + k + j] = 0;
02560             }
02561             
02562             else {
02563                working_space[2 * fSize + k + 2 * j - i % j] =
02564                    working_space[2 * fSize + k + i % j];
02565                working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j] =
02566                    -working_space[4 * fSize + 2 * fSize + k + i % j];
02567             }
02568          }
02569          for (i = 0; i < 2 * fSize; i++) {
02570             working_space[i] = working_space[2 * fSize + i];
02571             working_space[4 * fSize + i] =
02572                 working_space[4 * fSize + 2 * fSize + i];
02573          }
02574          GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
02575          m = (int) TMath::Power(2, fDegree);
02576          l = 2 * fSize / m;
02577          for (i = 0; i < l; i++)
02578             BitReverseHaar(working_space, 2 * fSize, m, i * m);
02579       }
02580       
02581       else if (fTransformType == kTransformSinWalsh || fTransformType == kTransformSinHaar) {
02582          j = (int) TMath::Power(2, fDegree) / 2;
02583          m = (int) TMath::Power(2, fDegree);
02584          l = 2 * fSize / m;
02585          for (i = 0; i < fSize; i++) {
02586             k = i / j;
02587             k = 2 * k * j;
02588             a = pi * (double) (i % j) / (double) (2 * j);
02589             if (i % j == 0) {
02590                working_space[2 * fSize + k + j + i % j] =
02591                    working_space[j + k / 2 - i % j - 1] * TMath::Sqrt(2.0);
02592                working_space[4 * fSize + 2 * fSize + k + j + i % j] = 0;
02593             }
02594             
02595             else {
02596                b = TMath::Sin(a);
02597                a = TMath::Cos(a);
02598                working_space[4 * fSize + 2 * fSize + k + j + i % j] =
02599                    -(double) working_space[j + k / 2 - i % j - 1] * b;
02600                working_space[2 * fSize + k + j + i % j] =
02601                    (double) working_space[j + k / 2 - i % j - 1] * a;
02602          } } for (i = 0; i < fSize; i++) {
02603             k = i / j;
02604             k = 2 * k * j;
02605             if (i % j == 0) {
02606                working_space[2 * fSize + k] = 0;
02607                working_space[4 * fSize + 2 * fSize + k] = 0;
02608             }
02609             
02610             else {
02611                working_space[2 * fSize + k + i % j] =
02612                    working_space[2 * fSize + k + 2 * j - i % j];
02613                working_space[4 * fSize + 2 * fSize + k + i % j] =
02614                    -working_space[4 * fSize + 2 * fSize + k + 2 * j - i % j];
02615             }
02616          }
02617          for (i = 0; i < 2 * fSize; i++) {
02618             working_space[i] = working_space[2 * fSize + i];
02619             working_space[4 * fSize + i] =
02620                 working_space[4 * fSize + 2 * fSize + i];
02621          }
02622          GeneralInv(working_space, 2 * fSize, fDegree, fTransformType);
02623          for (i = 0; i < l; i++)
02624             BitReverseHaar(working_space, 2 * fSize, m, i * m);
02625       }
02626       for (i = 0; i < fSize; i++) {
02627          if (fTransformType >= kTransformCosWalsh) {
02628             k = i / j;
02629             k = 2 * k * j;
02630             val = working_space[k + i % j];
02631          }
02632          
02633          else
02634             val = working_space[i];
02635          destVector[i] = val;
02636       }
02637       break;
02638    }
02639    delete[]working_space;
02640    return;
02641 }
02642 
02643 //___________________________________________________________________________
02644 void TSpectrumTransform::SetTransformType(Int_t transType, Int_t degree)
02645 {
02646 //////////////////////////////////////////////////////////////////////////////
02647 //   SETTER FUNCION                                                      
02648 //                                                     
02649 //   This function sets the following parameters for transform:
02650 //         -transType - type of transform (Haar, Walsh, Cosine, Sine, Fourier, Hartley, Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar)
02651 //         -degree - degree of mixed transform, applies only for Fourier-Walsh, Fourier-Haar, Walsh-Haar, Cosine-Walsh, Cosine-Haar, Sine-Walsh, Sine-Haar transforms
02652 //////////////////////////////////////////////////////////////////////////////      
02653    Int_t j, n;
02654    j = 0;
02655    n = 1;
02656    for (; n < fSize;) {
02657       j += 1;
02658       n = n * 2;
02659    } 
02660    if (transType < kTransformHaar || transType > kTransformSinHaar){
02661       Error ("TSpectrumTransform","Invalid type of transform");
02662       return;       
02663    }
02664    if (transType >= kTransformFourierWalsh && transType <= kTransformSinHaar) {
02665       if (degree > j || degree < 1){
02666          Error ("TSpectrumTransform","Invalid degree of mixed transform");
02667          return;          
02668       }
02669    }
02670    fTransformType=transType;
02671    fDegree=degree;
02672 }
02673 
02674 
02675 //___________________________________________________________________________
02676 void TSpectrumTransform::SetRegion(Int_t xmin, Int_t xmax)
02677 {
02678 //////////////////////////////////////////////////////////////////////////////
02679 //   SETTER FUNCION                                                      
02680 //                                                     
02681 //   This function sets the filtering or enhancement region:
02682 //         -xmin, xmax
02683 //////////////////////////////////////////////////////////////////////////////         
02684    if(xmin<0 || xmax < xmin || xmax >= fSize){ 
02685       Error("TSpectrumTransform", "Wrong range");      
02686       return;
02687    }         
02688    fXmin = xmin;
02689    fXmax = xmax;
02690 }
02691 
02692 //___________________________________________________________________________
02693 void TSpectrumTransform::SetDirection(Int_t direction)
02694 {
02695 //////////////////////////////////////////////////////////////////////////////
02696 //   SETTER FUNCION                                                      
02697 //                                                     
02698 //   This function sets the direction of the transform:
02699 //         -direction (forward or inverse)
02700 //////////////////////////////////////////////////////////////////////////////   
02701    if(direction != kTransformForward && direction != kTransformInverse){ 
02702       Error("TSpectrumTransform", "Wrong direction");      
02703       return;
02704    }         
02705    fDirection = direction;
02706 }
02707 
02708 //___________________________________________________________________________
02709 void TSpectrumTransform::SetFilterCoeff(Float_t filterCoeff)
02710 {
02711 //////////////////////////////////////////////////////////////////////////////
02712 //   SETTER FUNCION                                                      
02713 //                                                     
02714 //   This function sets the filter coefficient:
02715 //         -filterCoeff - after the transform the filtered region (xmin, xmax) is replaced by this coefficient. Applies only for filtereng operation.
02716 //////////////////////////////////////////////////////////////////////////////   
02717    fFilterCoeff = filterCoeff;
02718 }
02719 
02720 //___________________________________________________________________________
02721 void TSpectrumTransform::SetEnhanceCoeff(Float_t enhanceCoeff)
02722 {
02723 //////////////////////////////////////////////////////////////////////////////
02724 //   SETTER FUNCION                                                      
02725 //                                                     
02726 //   This function sets the enhancement coefficient:
02727 //         -enhanceCoeff - after the transform the enhanced region (xmin, xmax) is multiplied by this coefficient. Applies only for enhancement operation.
02728 //////////////////////////////////////////////////////////////////////////////   
02729    fEnhanceCoeff = enhanceCoeff;
02730 }
Generated on Tue Jul 5 14:24:39 2011 for ROOT_528-00b_version by  doxygen 1.5.1