TPC/AliTPCMonitorFFT.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.2  2007/10/12 13:36:27  cvetan
  19 Coding convention fixes from Stefan
  20
  21 Revision 1.1  2007/09/17 10:23:31  cvetan
  22 New TPC monitoring package from Stefan Kniege. The monitoring package can be started by running TPCMonitor.C macro located in macros folder.
  23
  24 */
  25
  26 ////////////////////////////////////////////////////////////////////////
  27 ////
  28 //// AliTPCMonitorFFT class
  29 ////
  30 //// Wrapper class to perform Fast Fourier Transformations.
  31 //// The code is based on the Gnu Scientific Library.
  32 //// See documentation of gsl for further details.
  33 ////
  34 //// Author: Stefan Kniege, IKF, Frankfurt
  35 ////
  36 ////
  37 /////////////////////////////////////////////////////////////////////////
  38
  39
  40 #include "AliTPCMonitorFFT.h"
  41
  42
  43 ClassImp(AliTPCMonitorFFT)
  44
  45 //__________________________________________________________________
  46   AliTPCMonitorFFT::AliTPCMonitorFFT()
  47 {
  48   // Constuctor
  49 }
  50
  51 //__________________________________________________________________
  52 AliTPCMonitorFFT::~AliTPCMonitorFFT()
  53 {
  54   // Destructor
  55 }
  56
  57
  58 //__________________________________________________________________
  59 Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
  60 {
  61
  62   // Wrapper function for forward fourier transformation
  63   Int_t direction = 1;
  64   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  65   return ret ;
  66 }
  67
  68 //__________________________________________________________________
  69 Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
  70 {
  71   // Wrapper function for backward  fourier transformation
  72   Int_t direction = -1;
  73   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  74   return ret ;
  75 }
  76
  77 //__________________________________________________________________
  78 Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
  79 {
  80   // Wrapper function for inverse fourier transformation
  81   Int_t direction = -1;
  82   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  83
  84   if (ret)
  85     {
  86       return ret;
  87     }
  88
  89   const Double_t knorm = 1.0 / n;
  90   size_t i;
  91   for (i = 0; i < n; i++)
  92     {
  93       REAL(data,stride,i) *= knorm;
  94       IMAG(data,stride,i) *= knorm;
  95     }
  96
  97   return ret;
  98 }
  99
 100 //__________________________________________________________________
 101 Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data,   Int_t stride, size_t n, Int_t sign)
 102 {
 103   // Wrapper function for fourier transformation depending on sign  ( forward 1;backward -1 )
 104   Int_t result ;
 105   size_t dual;
 106   size_t bit;
 107   size_t logn = 0;
 108   Int_t status;
 109
 110   if (n == 1) /* identity operation */
 111     {
 112       return 0 ;
 113     }
 114   /* make sure that n is a power of 2 */
 115
 116   result = FFTBinaryLogn(n) ;
 117
 118   if (result == -1)
 119     {
 120       cout << " not a power of 2 " << endl;
 121     }
 122   else
 123     {
 124       logn = result ;
 125     }
 126
 127   /* apply fft recursion */
 128
 129   dual = n / 2;
 130
 131   for (bit = 0; bit < logn; bit++)
 132     {
 133       Double_t wreal = 1.0;
 134       Double_t wimag = 0.0;
 135
 136       const Double_t ktheta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
 137
 138       const Double_t ks = sin (ktheta);
 139       const Double_t kt = sin (ktheta / 2.0);
 140       const Double_t ks2 = 2.0 * kt * kt;
 141
 142       size_t a, b;
 143
 144       for (b = 0; b < dual; b++)
 145         {
 146           for (a = 0; a < n; a+= 2 * dual)
 147             {
 148               const size_t ki = b + a;
 149               const size_t kj = b + a + dual;
 150
 151               const Double_t kt1real = REAL(data,stride,ki) + REAL(data,stride,kj);
 152               const Double_t kt1imag = IMAG(data,stride,ki) + IMAG(data,stride,kj);
 153               const Double_t kt2real = REAL(data,stride,ki) - REAL(data,stride,kj);
 154               const Double_t kt2imag = IMAG(data,stride,ki) - IMAG(data,stride,kj);
 155
 156               REAL(data,stride,ki) = kt1real;
 157               IMAG(data,stride,ki) = kt1imag;
 158               REAL(data,stride,kj) = wreal*kt2real - wimag * kt2imag;
 159               IMAG(data,stride,kj) = wreal*kt2imag + wimag * kt2real;
 160             }
 161
 162           /* trignometric recurrence for w-> exp(i ktheta) w */
 163
 164           {
 165             const Double_t ktmpreal = wreal - ks * wimag - ks2 * wreal;
 166             const Double_t ktmpimag = wimag + ks * wreal - ks2 * wimag;
 167             wreal = ktmpreal;
 168             wimag = ktmpimag;
 169           }
 170         }
 171       dual /= 2;
 172     }
 173
 174   /* bit reverse the ordering of output data for decimation in
 175      frequency algorithm */
 176
 177   //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
 178   status  = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
 179   return 0;
 180 }
 181
 182 //__________________________________________________________________
 183 Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn) const
 184 {
 185   // Wrapper function from gnu scientific library
 186   /* This is the Goldrader bit-reversal algorithm */
 187
 188   size_t i;
 189   size_t j = 0;
 190
 191   logn = 0 ; /* not needed for this algorithm */
 192
 193   for (i = 0; i < n - 1; i++)
 194     {
 195       size_t k = n / 2 ;
 196
 197       if (i < j)
 198         {
 199           const Double_t ktmpreal = REAL(data,stride,i);
 200           const Double_t ktmpimag = IMAG(data,stride,i);
 201           REAL(data,stride,i) = REAL(data,stride,j);
 202           IMAG(data,stride,i) = IMAG(data,stride,j);
 203           REAL(data,stride,j) = ktmpreal;
 204           IMAG(data,stride,j) = ktmpimag;
 205         }
 206
 207       while (k <= j)
 208         {
 209           j = j - k ;
 210           k = k / 2 ;
 211         }
 212
 213       j += k ;
 214     }
 215
 216   return 0;
 217
 218 }
 219
 220 //__________________________________________________________________
 221 Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n) const
 222 {
 223
 224   // Return log on base 2
 225   size_t ntest ;
 226   size_t binarylogn = 0 ;
 227   size_t k = 1;
 228
 229   while (k < n)
 230     {
 231       k *= 2;
 232       binarylogn++;
 233     }
 234
 235   ntest = (1 << binarylogn) ;
 236
 237   if (n != ntest )
 238     {
 239       return -1 ; /* n is not a power of 2 */
 240     }
 241
 242   return binarylogn;
 243 }
 244