Chris@4: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@4: 
Chris@4: /*
Chris@4:   Vamp feature extraction plugin for the BeatRoot beat tracker.
Chris@4: 
Chris@4:   Centre for Digital Music, Queen Mary, University of London.
Chris@4:   This file copyright 2011 Simon Dixon, Chris Cannam and QMUL.
Chris@4:     
Chris@4:   This program is free software; you can redistribute it and/or
Chris@4:   modify it under the terms of the GNU General Public License as
Chris@4:   published by the Free Software Foundation; either version 2 of the
Chris@4:   License, or (at your option) any later version.  See the file
Chris@4:   COPYING included with this distribution for more information.
Chris@4: */
Chris@4: 
Chris@4: #include "Peaks.h"
Chris@4: 
Chris@9: int Peaks::pre = 3;
Chris@9: int Peaks::post = 1;
Chris@4: 
Chris@15: 
Chris@15: int Peaks::findPeaks(const vector<double> &data, vector<int> peaks, int width) {
Chris@15:     int peakCount = 0;
Chris@15:     int maxp = 0;
Chris@15:     int mid = 0;
Chris@15:     int end = data.size();
Chris@15:     while (mid < end) {
Chris@15:         int i = mid - width;
Chris@15:         if (i < 0)
Chris@15:             i = 0;
Chris@15:         int stop = mid + width + 1;
Chris@22:         if (stop > (int)data.size())
Chris@15:             stop = data.size();
Chris@15:         maxp = i;
Chris@15:         for (i++; i < stop; i++)
Chris@15:             if (data[i] > data[maxp])
Chris@15:                 maxp = i;
Chris@15:         if (maxp == mid) {
Chris@15:             int j;
Chris@15:             for (j = peakCount; j > 0; j--) {
Chris@15:                 if (data[maxp] <= data[peaks[j-1]])
Chris@15:                     break;
Chris@22:                 else if (j < (int)peaks.size())
Chris@15:                     peaks[j] = peaks[j-1];
Chris@15:             }
Chris@22:             if (j != (int)peaks.size())
Chris@15:                 peaks[j] = maxp;
Chris@22:             if (peakCount != (int)peaks.size())
Chris@15:                 peakCount++;
Chris@15:         }
Chris@15:         mid++;
Chris@15:     }
Chris@15:     return peakCount;
Chris@15: } // findPeaks()
Chris@15: 
Chris@15: vector<int> Peaks::findPeaks(const vector<double> &data, int width,
Chris@15:                              double threshold, double decayRate, bool isRelative) {
Chris@15:     vector<int> peaks;
Chris@15:     int maxp = 0;
Chris@15:     int mid = 0;
Chris@18:     if (data.empty()) return peaks;
Chris@15:     int end = data.size();
Chris@15:     double av = data[0];
Chris@15:     while (mid < end) {
Chris@15:         av = decayRate * av + (1 - decayRate) * data[mid];
Chris@15:         if (av < data[mid])
Chris@15:             av = data[mid];
Chris@15:         int i = mid - width;
Chris@15:         if (i < 0)
Chris@15:             i = 0;
Chris@15:         int stop = mid + width + 1;
Chris@22:         if (stop > (int)data.size())
Chris@15:             stop = data.size();
Chris@15:         maxp = i;
Chris@15:         for (i++; i < stop; i++)
Chris@15:             if (data[i] > data[maxp])
Chris@15:                 maxp = i;
Chris@15:         if (maxp == mid) {
Chris@15:             if (overThreshold(data, maxp, width, threshold, isRelative,av)){
Chris@15:                 peaks.push_back(maxp);
Chris@15:             }
Chris@15:         }
Chris@15:         mid++;
Chris@15:     }
Chris@15:     return peaks;
Chris@15: } // findPeaks()
Chris@15: 
Chris@15: double Peaks::expDecayWithHold(double av, double decayRate,
Chris@15:                                const vector<double> &data, int start, int stop) {
Chris@15:     while (start < stop) {
Chris@15:         av = decayRate * av + (1 - decayRate) * data[start];
Chris@15:         if (av < data[start])
Chris@15:             av = data[start];
Chris@15:         start++;
Chris@15:     }
Chris@15:     return av;
Chris@15: } // expDecayWithHold()
Chris@15: 
Chris@15: bool Peaks::overThreshold(const vector<double> &data, int index, int width,
Chris@15:                           double threshold, bool isRelative,
Chris@15:                           double av) {
Chris@15:     if (data[index] < av)
Chris@15:         return false;
Chris@15:     if (isRelative) {
Chris@15:         int iStart = index - pre * width;
Chris@15:         if (iStart < 0)
Chris@15:             iStart = 0;
Chris@15:         int iStop = index + post * width;
Chris@22:         if (iStop > (int)data.size())
Chris@15:             iStop = data.size();
Chris@15:         double sum = 0;
Chris@15:         int count = iStop - iStart;
Chris@15:         while (iStart < iStop)
Chris@15:             sum += data[iStart++];
Chris@15:         return (data[index] > sum / count + threshold);
Chris@15:     } else
Chris@15:         return (data[index] > threshold);
Chris@15: } // overThreshold()
Chris@15: 
Chris@15: void Peaks::normalise(vector<double> &data) {
Chris@15:     double sx = 0;
Chris@15:     double sxx = 0;
Chris@22:     for (int i = 0; i < (int)data.size(); i++) {
Chris@15:         sx += data[i];
Chris@15:         sxx += data[i] * data[i];
Chris@15:     }
Chris@15:     double mean = sx / data.size();
Chris@15:     double sd = sqrt((sxx - sx * mean) / data.size());
Chris@15:     if (sd == 0)
Chris@15:         sd = 1;		// all data[i] == mean  -> 0; avoids div by 0
Chris@22:     for (int i = 0; i < (int)data.size(); i++) {
Chris@15:         data[i] = (data[i] - mean) / sd;
Chris@15:     }
Chris@15: } // normalise()
Chris@15: 
Chris@15: /** Uses an n-point linear regression to estimate the slope of data.
Chris@15:  *  @param data input data
Chris@15:  *  @param hop spacing of data points
Chris@15:  *  @param n length of linear regression
Chris@15:  *  @param slope output data
Chris@15:  */
Chris@15: void Peaks::getSlope(const vector<double> &data, double hop, int n,
Chris@15:                      vector<double> &slope) {
Chris@15:     int i = 0, j = 0;
Chris@15:     double t;
Chris@15:     double sx = 0, sxx = 0, sy = 0, sxy = 0;
Chris@15:     for ( ; i < n; i++) {
Chris@15:         t = i * hop;
Chris@15:         sx += t;
Chris@15:         sxx += t * t;
Chris@15:         sy += data[i];
Chris@15:         sxy += t * data[i];
Chris@15:     }
Chris@15:     double delta = n * sxx - sx * sx;
Chris@15:     for ( ; j < n / 2; j++)
Chris@15:         slope[j] = (n * sxy - sx * sy) / delta;
Chris@22:     for ( ; j < (int)data.size() - (n + 1) / 2; j++, i++) {
Chris@15:         slope[j] = (n * sxy - sx * sy) / delta;
Chris@15:         sy += data[i] - data[i - n];
Chris@15:         sxy += hop * (n * data[i] - sy);
Chris@15:     }
Chris@22:     for ( ; j < (int)data.size(); j++)
Chris@15:         slope[j] = (n * sxy - sx * sy) / delta;
Chris@15: } // getSlope()
Chris@15: 
Chris@15: int Peaks::imin(const vector<double> &arr) {
Chris@15:     int i = 0;
Chris@22:     for (int j = 1; j < (int)arr.size(); j++)
Chris@15:         if (arr[j] < arr[i])
Chris@15:             i = j;
Chris@15:     return i;
Chris@15: } // imin()
Chris@15: 
Chris@15: int Peaks::imax(const vector<double> &arr) {
Chris@15:     int i = 0;
Chris@22:     for (int j = 1; j < (int)arr.size(); j++)
Chris@15:         if (arr[j] > arr[i])
Chris@15:             i = j;
Chris@15:     return i;
Chris@15: } // imax()
Chris@15: