lib/fpa/Generic/PeakDetector.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <fpa/Generic/PeakDetector.h>
   7 #include <cmath>
   8
   9 // -------------------------------------------------------------------------
  10 fpa::Generic::PeakDetector::
  11 PeakDetector( )
  12   : m_K( 3 ),
  13     m_T( 3.5 ),
  14     m_I( 0.5 )
  15 {
  16 }
  17
  18 // -------------------------------------------------------------------------
  19 fpa::Generic::PeakDetector::
  20 ~PeakDetector( )
  21 {
  22 }
  23
  24 // -------------------------------------------------------------------------
  25 unsigned long fpa::Generic::PeakDetector::
  26 GetKernelSize( ) const
  27 {
  28   return( this->m_K );
  29 }
  30
  31 // -------------------------------------------------------------------------
  32 double fpa::Generic::PeakDetector::
  33 GetThreshold( ) const
  34 {
  35   return( this->m_T );
  36 }
  37
  38 // -------------------------------------------------------------------------
  39 double fpa::Generic::PeakDetector::
  40 GetInfluence( ) const
  41 {
  42   return( this->m_I );
  43 }
  44
  45 // -------------------------------------------------------------------------
  46 void fpa::Generic::PeakDetector::
  47 SetKernelSize( unsigned long k )
  48 {
  49   this->m_K = k;
  50   this->Clear( );
  51 }
  52
  53 // -------------------------------------------------------------------------
  54 void fpa::Generic::PeakDetector::
  55 SetThreshold( double t )
  56 {
  57   this->m_T = t;
  58   this->Clear( );
  59 }
  60
  61 // -------------------------------------------------------------------------
  62 void fpa::Generic::PeakDetector::
  63 SetInfluence( double i )
  64 {
  65   this->m_I = i;
  66   this->Clear( );
  67 }
  68
  69 // -------------------------------------------------------------------------
  70 void fpa::Generic::PeakDetector::
  71 Clear( )
  72 {
  73   this->m_X.clear( );
  74   this->m_Y.clear( );
  75   this->m_YF.clear( );
  76   this->m_Avg.clear( );
  77   this->m_STD.clear( );
  78   this->m_Peaks.clear( );
  79   this->m_Mean = 0;
  80   this->m_Vari = 0;
  81 }
  82
  83 // -------------------------------------------------------------------------
  84 unsigned long fpa::Generic::PeakDetector::
  85 GetNumberOfSamples( ) const
  86 {
  87   return( this->m_X.size( ) );
  88 }
  89
  90 // -------------------------------------------------------------------------
  91 fpa::Generic::PeakDetector::
  92 TPeak fpa::Generic::PeakDetector::
  93 AddValue( double x, double y )
  94 {
  95   this->m_X.push_back( x );
  96   this->m_Y.push_back( y );
  97
  98   if( this->m_YF.size( ) < this->m_K )
  99   {
 100     double n = double( this->m_YF.size( ) + 1 );
 101     this->m_YF.push_back( y );
 102     this->m_Avg.push_back( double( 0 ) );
 103     this->m_STD.push_back( double( 0 ) );
 104     this->m_Peaks.push_back( Self::NoPeak );
 105     if( n > double( 1 ) )
 106       this->m_Vari =
 107         ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * this->m_Vari ) +
 108         ( ( ( y - this->m_Mean ) * ( y - this->m_Mean ) ) / n );
 109     this->m_Mean += ( y - this->m_Mean ) / n;
 110     if( this->m_YF.size( ) == this->m_K )
 111     {
 112       this->m_Avg.push_back( this->m_Mean );
 113       this->m_STD.push_back( std::sqrt( this->m_Vari ) );
 114
 115     } // fi
 116   }
 117   else
 118   {
 119     unsigned long i = this->m_X.size( ) - 1;
 120     if(
 121       ( std::fabs( y - this->m_Avg[ i - 1 ] ) ) >
 122       ( this->m_T * this->m_STD[ i - 1 ] )
 123       )
 124     {
 125       this->m_Peaks.push_back(
 126         ( y > this->m_Avg[ i - 1 ] )? Self::PosPeak: Self::NegPeak
 127         );
 128       this->m_YF.push_back(
 129         ( this->m_I * y ) + ( ( 1.0 - this->m_I ) * this->m_YF[ i - 1 ] )
 130         );
 131     }
 132     else
 133     {
 134       this->m_Peaks.push_back( Self::NoPeak );
 135       this->m_YF.push_back( y );
 136
 137     } // fi
 138
 139     double avg = double( 0 );
 140     double var = double( 0 );
 141     unsigned long k = 0;
 142     for( unsigned long j = i - this->m_K; j <= i; ++j, ++k )
 143     {
 144       double v = this->m_YF[ j ];
 145       double n = double( k + 1 );
 146       if( k > 1 )
 147         var =
 148           ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * var ) +
 149           ( ( ( v - avg ) * ( v - avg ) ) / n );
 150       avg += ( v - avg ) / n;
 151
 152     } // rof
 153     this->m_Avg.push_back( avg );
 154     this->m_STD.push_back( std::sqrt( var ) );
 155
 156   } // fi
 157   return( this->m_Peaks.back( ) );
 158 }
 159
 160 // eof - $RCSfile$