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 */
  19
  20 #include "AliTPCMonitorFFT.h"
  21
  22
  23 ClassImp(AliTPCMonitorFFT)
  24
  25 //__________________________________________________________________
  26   AliTPCMonitorFFT::AliTPCMonitorFFT()
  27 {
  28   // Constuctor
  29 }
  30
  31 //__________________________________________________________________
  32 AliTPCMonitorFFT::~AliTPCMonitorFFT()
  33 {
  34   // Destructor
  35 }
  36
  37
  38 //__________________________________________________________________
  39 Int_t AliTPCMonitorFFT::ComplexRadix2ForwardWrap(Double_t* data, Int_t stride, size_t n)
  40 {
  41
  42   // Wrapper function for forward fourier transformation
  43   Int_t direction = 1;
  44   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  45   return ret ;
  46 }
  47
  48 //__________________________________________________________________
  49 Int_t AliTPCMonitorFFT::ComplexRadix2BackwardWrap(Double_t* data,Int_t stride, size_t n)
  50 {
  51   // Wrapper function for backward  fourier transformation
  52   Int_t direction = -1;
  53   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  54   return ret ;
  55 }
  56
  57 //__________________________________________________________________
  58 Int_t AliTPCMonitorFFT::ComplexRadix2InverseWrap(Double_t* data,Int_t stride, size_t n)
  59 {
  60   // Wrapper function for inverse fourier transformation
  61   Int_t direction = -1;
  62   Int_t ret = ComplexRadix2TransformWrap(data,stride,n,direction);
  63
  64   if (ret)
  65     {
  66       return ret;
  67     }
  68
  69   const Double_t norm = 1.0 / n;
  70   size_t i;
  71   for (i = 0; i < n; i++)
  72     {
  73       REAL(data,stride,i) *= norm;
  74       IMAG(data,stride,i) *= norm;
  75     }
  76
  77   return ret;
  78 }
  79
  80 //__________________________________________________________________
  81 Int_t AliTPCMonitorFFT::ComplexRadix2TransformWrap(Double_t* data,   Int_t stride, size_t n, Int_t sign)
  82 {
  83   // Wrapper function for fourier transformation depending on sign  ( forward 1;backward -1 )
  84   Int_t result ;
  85   size_t dual;
  86   size_t bit;
  87   size_t logn = 0;
  88   Int_t status;
  89
  90   if (n == 1) /* identity operation */
  91     {
  92       return 0 ;
  93     }
  94   /* make sure that n is a power of 2 */
  95
  96   result = FFTBinaryLogn(n) ;
  97
  98   if (result == -1)
  99     {
 100       cout << " not a power of 2 " << endl;
 101     }
 102   else
 103     {
 104       logn = result ;
 105     }
 106
 107   /* apply fft recursion */
 108
 109   dual = n / 2;
 110
 111   for (bit = 0; bit < logn; bit++)
 112     {
 113       Double_t w_real = 1.0;
 114       Double_t w_imag = 0.0;
 115
 116       const Double_t theta = 2.0 * ((Int_t) sign) * M_PI / ((Double_t) (2 * dual));
 117
 118       const Double_t s = sin (theta);
 119       const Double_t t = sin (theta / 2.0);
 120       const Double_t s2 = 2.0 * t * t;
 121
 122       size_t a, b;
 123
 124       for (b = 0; b < dual; b++)
 125         {
 126           for (a = 0; a < n; a+= 2 * dual)
 127             {
 128               const size_t i = b + a;
 129               const size_t j = b + a + dual;
 130
 131               const Double_t t1_real = REAL(data,stride,i) + REAL(data,stride,j);
 132               const Double_t t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
 133               const Double_t t2_real = REAL(data,stride,i) - REAL(data,stride,j);
 134               const Double_t t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
 135
 136               REAL(data,stride,i) = t1_real;
 137               IMAG(data,stride,i) = t1_imag;
 138               REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
 139               IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
 140             }
 141
 142           /* trignometric recurrence for w-> exp(i theta) w */
 143
 144           {
 145             const Double_t tmp_real = w_real - s * w_imag - s2 * w_real;
 146             const Double_t tmp_imag = w_imag + s * w_real - s2 * w_imag;
 147             w_real = tmp_real;
 148             w_imag = tmp_imag;
 149           }
 150         }
 151       dual /= 2;
 152     }
 153
 154   /* bit reverse the ordering of output data for decimation in
 155      frequency algorithm */
 156
 157   //status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
 158   status  = ComplexBitReverseOrderWrap(data, stride, n, logn) ;
 159   return 0;
 160 }
 161
 162 //__________________________________________________________________
 163 Int_t AliTPCMonitorFFT::ComplexBitReverseOrderWrap(Double_t* data, Int_t stride, size_t n, Int_t logn)
 164 {
 165   // Wrapper function from gnu scientific library
 166   /* This is the Goldrader bit-reversal algorithm */
 167
 168   size_t i;
 169   size_t j = 0;
 170
 171   logn = 0 ; /* not needed for this algorithm */
 172
 173   for (i = 0; i < n - 1; i++)
 174     {
 175       size_t k = n / 2 ;
 176
 177       if (i < j)
 178         {
 179           const Double_t tmp_real = REAL(data,stride,i);
 180           const Double_t tmp_imag = IMAG(data,stride,i);
 181           REAL(data,stride,i) = REAL(data,stride,j);
 182           IMAG(data,stride,i) = IMAG(data,stride,j);
 183           REAL(data,stride,j) = tmp_real;
 184           IMAG(data,stride,j) = tmp_imag;
 185         }
 186
 187       while (k <= j)
 188         {
 189           j = j - k ;
 190           k = k / 2 ;
 191         }
 192
 193       j += k ;
 194     }
 195
 196   return 0;
 197
 198 }
 199
 200 //__________________________________________________________________
 201 Int_t AliTPCMonitorFFT::FFTBinaryLogn(size_t n)
 202 {
 203
 204   // Return log on base 2
 205   size_t ntest ;
 206   size_t binary_logn = 0 ;
 207   size_t k = 1;
 208
 209   while (k < n)
 210     {
 211       k *= 2;
 212       binary_logn++;
 213     }
 214
 215   ntest = (1 << binary_logn) ;
 216
 217   if (n != ntest )
 218     {
 219       return -1 ; /* n is not a power of 2 */
 220     }
 221
 222   return binary_logn;
 223 }
 224