UnivTools.cc

Go to the documentation of this file.

#include <src/UnivTools.h>
#include <src/Functions.h>


namespace un2 {

  // https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data
  std::vector<std::vector<double>>
  FindPeaks(const std::vector<double>& input, const double totalThreshold)
  {   
    // lag = n time bins, over which the local average will be calculated to detect peaks
    const unsigned int lag = 5;
    // threshold in standard deviations from local average to detect peaks
    const double threshold = 2.5;
    // between 0 and 1, where 1 is normal influence, 0.5 is half
    const double influence = 0.25;

    if (input.size() <= lag + 2)
      return std::vector<std::vector<double>>();

    // Initialise variables
    std::vector<std::vector<double>> peaks;
    std::vector<double> filteredY(input.size());
    std::vector<double> avgFilter(input.size());
    std::vector<double> stdFilter(input.size());
    std::vector<double> subVecStart(input.begin(), input.begin() + lag);
    const double mean = Mean(subVecStart);
    avgFilter[lag] = mean;
    stdFilter[lag] = StandardDeviation(subVecStart, mean);

    for (size_t i = lag + 1, n = input.size(); i < n; ++i) {
      // if current value is 3.5 standard deviations away from the mean of the last 5 time bins
      // and above the total signal threshold, check further
      const auto& inp = input[i];
      if (std::abs(inp - avgFilter[i - 1]) > threshold * stdFilter[i - 1] && inp > totalThreshold) {
        // positive peak is detected if current value is above absolutely larger than the average
        // of the last 5 time bins and if there is no further upwards trend
        if (inp > avgFilter[i - 1] && inp > input[i - 1] && inp > input[i + 1]) {
          peaks.push_back({ double(i), input[i] }); // Positive signal
        } else {
          // negative peaks could be stored here
          // peaks.push_back({i*1., input[i]}); // Negative signal
        }
        // Make influence lower
        // larger influence (closer to one) reduces the likelihood of peak detection based
        // when there were many peaks in the current area
        filteredY[i] = influence * inp + (1 - influence) * filteredY[i - 1];
      } else
        filteredY[i] = inp;

      //Adjust the filters
      const std::vector<double> subVec(filteredY.begin() + (i - lag), filteredY.begin() + i);
      const double mean = Mean(subVec);
      avgFilter[i] = mean;
      stdFilter[i] = StandardDeviation(subVec, mean);
    }

    return peaks;
  }


  double
  MIPtoMIPpersqm(const double mIPval, const double theta)
  {
    return mIPval / (4 * cos(theta));
  }


  std::vector<double>
  MIPtoMIPpersqm(const std::vector<double>& mIPvals, const double theta)
  {
    std::vector<double> returnVec;
    for (const auto m : mIPvals)
      returnVec.push_back(MIPtoMIPpersqm(m, theta));
    return returnVec;
  }


  double
  VEMtoVEMpersqm(const double VEMval, const double theta)
  {
    return VEMval / (10 * cos(theta) + (1.2 + 2*1.8) * sin(theta));
  }


  std::vector<double>
  VEMtoVEMpersqm(const std::vector<double>& mIPvals, const double theta)
  {
    std::vector<double> returnVec;
    for (const auto m : mIPvals)
      returnVec.push_back(VEMtoVEMpersqm(m, theta));
    return returnVec;
  }

}