SignalProcessing.hpp

Go to the documentation of this file.

 
 
template <class T> 
double* SignalProcessing<T>::simpleMovingAverage(T* a_signal, 
                                unsigned int a_length,
                                unsigned int a_windowSize
                                )
{
  //get a variable representing the new length
  int numRtnSamples = a_length;
  //allocate new storage for the filtered signal
  double* rtrnPtr = new double[numRtnSamples];
  //straigh forward implementation would be nested loop recalculating the 
  //average for each out put sample. Instead of the loop subtracting the first
  //of the last element and adding the next element past the window will 
  //reducing the complexity
 
  //initialize the average to the first a_length samples
  double avg =0;
  for(int iS =0;iS < a_windowSize;iS++)
  {
    avg += static_cast<double>(a_signal[iS]);
  }
  //set the first element windowsize elements
  for(int iS =0;iS < a_windowSize;iS++)
  {
    rtrnPtr[iS] = avg/static_cast<double>(a_windowSize);
  }
  for(int iS =a_windowSize; iS < numRtnSamples;iS++)
  {
    //subtract the first element of the last window
    avg -= static_cast<double>(a_signal[iS-a_windowSize]);
    //add the next element past the last window
    avg += static_cast<double>(a_signal[iS]);
    //set the next element of the array
    rtrnPtr[iS] = avg/static_cast<double>(a_windowSize);
  }
  return rtrnPtr;
}
 
template <class T> double* SignalProcessing<T>::
                    iterativeMovingAverage(T* a_signal, 
                              unsigned int a_length,
                              unsigned int a_windowSize,
                              unsigned int a_numIter
                              )
{
  double* smoothedInput = simpleMovingAverage(a_signal,
                                              a_length,
                                              a_windowSize);
  for(int iI = 0; iI < a_numIter;iI++)
  {
    double* tmpPtr = smoothedInput;
    smoothedInput = 
     SignalProcessing<double>::simpleMovingAverage(smoothedInput,
                                                   a_length,
                                                   a_windowSize);
     //delete the last array
     delete[] tmpPtr;
  }
  return smoothedInput;
}
 
template <class T> double* SignalProcessing<T>::
                    singlePoleLowPassFilter(T* a_signal,
                    unsigned int a_length,
                    double a_xDecay
                   )
{
  //get a variable representing the length
  int numRtnSamples = a_length;
  //allocate new storage for the filtered signal
  double* rtrnPtr = new double[numRtnSamples];
  rtrnPtr[0] = static_cast<double>(a_signal[0]);
  for(int iS =1;iS < a_length;iS++)
  {
    rtrnPtr[iS] = (1-a_xDecay)*static_cast<double>(a_signal[iS])
                  +a_xDecay*rtrnPtr[iS-1];
  }
  return rtrnPtr;
} 
 
template <class T> double* SignalProcessing<T>::
                    singlePoleHighPassFilter(T* a_signal,
                    unsigned int a_length,
                    double a_xDecay
                   )
{
  //get a variable representing the length
  int numRtnSamples = a_length;
  //allocate new storage for the filtered signal
  double* rtrnPtr = new double[numRtnSamples];
  rtrnPtr[0] = static_cast<double>(a_signal[0]);
  for(int iS =1;iS < a_length;iS++)
  {
    rtrnPtr[iS] = 0.5*(1+a_xDecay)*static_cast<double>(a_signal[iS])
                  -0.5*(1+a_xDecay)*static_cast<double>(a_signal[iS-1])
                  +a_xDecay*rtrnPtr[iS-1];
  }
  return rtrnPtr;
} 
 
template <class T> double* SignalProcessing<T>::
                    fourStageLowPassFilter(T* a_signal,
                    unsigned int a_length,
                    double a_xDecay
                   )
{
  //get a variable representing the length
  int numRtnSamples = a_length;
  //allocate new storage for the filtered signal
  double* rtrnPtr = new double[numRtnSamples];
  double av4FirstSamples = static_cast<double>((a_signal[0]+a_signal[1]
     +a_signal[2]+a_signal[3])/4.);
  rtrnPtr[0] = av4FirstSamples;//static_cast<double>(a_signal[0]);
  rtrnPtr[1] = av4FirstSamples;//static_cast<double>(a_signal[1]);
  rtrnPtr[2] = av4FirstSamples;//static_cast<double>(a_signal[2]);
  rtrnPtr[3] = av4FirstSamples;//static_cast<double>(a_signal[3]);
  double coeff1 = pow(1-a_xDecay,4);
  double coeff2 = 4*a_xDecay;
  double coeff3 = -6*a_xDecay*a_xDecay;
  double coeff4 = 4*pow(a_xDecay,3);
  double coeff5 = -pow(a_xDecay,4);
  for(int iS =4;iS < a_length;iS++)
  {
    rtrnPtr[iS] = coeff1*static_cast<double>(a_signal[iS])
                  +coeff2*rtrnPtr[iS-1]
                  +coeff3*rtrnPtr[iS-2]
                  +coeff4*rtrnPtr[iS-3]
                  +coeff5*rtrnPtr[iS-4];
  }
  return rtrnPtr;
} 
 
template <class T> double* SignalProcessing<T>::
                      trapezoidalFilter(T* a_signal,
                      unsigned int a_length,
                      int a_k, 
                      int a_m  
                     )
{
  //get a variable representing the length
  int numRtnSamples = a_length;
  //allocate new storage for the filtered signal
  double* rtrnPtr = new double[numRtnSamples];
  for(int iS =0;iS < a_length;iS++)
  {
    if(iS<2*a_k+a_m)
    {
      rtrnPtr[iS] = static_cast<double>(a_signal[iS]);
    }
    else
    {
      rtrnPtr[iS] = static_cast<double>(a_signal[iS])
                    + rtrnPtr[iS-1]
                    - rtrnPtr[iS-a_k]
                    - rtrnPtr[iS-a_k-a_m]
                    + rtrnPtr[iS-2*a_k-a_m];
    }
  }
  return rtrnPtr;
}
 
template <class T> double SignalProcessing<T>::
                      calculateTIPS(T* a_signal,
                                    unsigned int a_length,
                                    int a_kFast,
                                    int a_mFast,
                                    int a_kSlow,
                                    int a_mSlow
                                    )
{
  double* slowTrap = trapezoidalFilter(a_signal,a_length,a_kSlow,a_mSlow);
  double* fastTrap = trapezoidalFilter(a_signal,a_length,a_kFast,a_mFast);
  double maxSlow = 0;
  double maxFast = 0;
  for(int iS=0;iS<a_length;iS++)
  {
    if(slowTrap[iS]>maxSlow)
    {
      maxSlow = slowTrap[iS];
    }
    if(fastTrap[iS]>maxFast)
    {
      maxFast = fastTrap[iS];
    }
  }
  return maxSlow/maxFast;
}
 
#ifdef _CLING_
template class SignalProcessing<short>;
template class SignalProcessing<double>;
#endif