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: #ifndef _BEATROOT_PEAKS_H_
Chris@4: #define _BEATROOT_PEAKS_H_
Chris@4: 
Chris@4: #include <vector>
Chris@4: #include <cmath>
Chris@4: 
Chris@4: using std::vector;
Chris@4: 
Chris@4: class Peaks
Chris@4: {
Chris@4: protected:
Chris@4:     static int pre;
Chris@4:     static int post;
Chris@4: 	
Chris@4: public:
Chris@4:     /** General peak picking method for finding n local maxima in an array
Chris@4:      *  @param data input data
Chris@4:      *  @param peaks list of peak indexes
Chris@4:      *  @param width minimum distance between peaks
Chris@4:      */
Chris@4:     static int findPeaks(const vector<double> &data, vector<int> peaks, int width) {
Chris@4: 	int peakCount = 0;
Chris@4: 	int maxp = 0;
Chris@4: 	int mid = 0;
Chris@4: 	int end = data.size();
Chris@4: 	while (mid < end) {
Chris@4: 	    int i = mid - width;
Chris@4: 	    if (i < 0)
Chris@4: 		i = 0;
Chris@4: 	    int stop = mid + width + 1;
Chris@4: 	    if (stop > data.size())
Chris@4: 		stop = data.size();
Chris@4: 	    maxp = i;
Chris@4: 	    for (i++; i < stop; i++)
Chris@4: 		if (data[i] > data[maxp])
Chris@4: 		    maxp = i;
Chris@4: 	    if (maxp == mid) {
Chris@4: 		int j;
Chris@4: 		for (j = peakCount; j > 0; j--) {
Chris@4: 		    if (data[maxp] <= data[peaks[j-1]])
Chris@4: 			break;
Chris@4: 		    else if (j < peaks.size())
Chris@4: 			peaks[j] = peaks[j-1];
Chris@4: 		}
Chris@4: 		if (j != peaks.size())
Chris@4: 		    peaks[j] = maxp;
Chris@4: 		if (peakCount != peaks.size())
Chris@4: 		    peakCount++;
Chris@4: 	    }
Chris@4: 	    mid++;
Chris@4: 	}
Chris@4: 	return peakCount;
Chris@4:     } // findPeaks()
Chris@4: 
Chris@4:     /** General peak picking method for finding local maxima in an array
Chris@4:      *  @param data input data
Chris@4:      *  @param width minimum distance between peaks
Chris@4:      *  @param threshold minimum value of peaks
Chris@4:      *  @return list of peak indexes
Chris@4:      */ 
Chris@4:     static vector<int> findPeaks(const vector<double> &data, int width,
Chris@4: 				 double threshold) {
Chris@4: 	return findPeaks(data, width, threshold, 0, false);
Chris@4:     } // findPeaks()
Chris@4: 	
Chris@4:     /** General peak picking method for finding local maxima in an array
Chris@4:      *  @param data input data
Chris@4:      *  @param width minimum distance between peaks
Chris@4:      *  @param threshold minimum value of peaks
Chris@4:      *  @param decayRate how quickly previous peaks are forgotten
Chris@4:      *  @param isRelative minimum value of peaks is relative to local average
Chris@4:      *  @return list of peak indexes
Chris@4:      */
Chris@4:     static vector<int> findPeaks(const vector<double> &data, int width,
Chris@4: 				 double threshold, double decayRate, bool isRelative) {
Chris@4: 	vector<int> peaks;
Chris@4: 	int maxp = 0;
Chris@4: 	int mid = 0;
Chris@4: 	int end = data.size();
Chris@4: 	double av = data[0];
Chris@4: 	while (mid < end) {
Chris@4: 	    av = decayRate * av + (1 - decayRate) * data[mid];
Chris@4: 	    if (av < data[mid])
Chris@4: 		av = data[mid];
Chris@4: 	    int i = mid - width;
Chris@4: 	    if (i < 0)
Chris@4: 		i = 0;
Chris@4: 	    int stop = mid + width + 1;
Chris@4: 	    if (stop > data.size())
Chris@4: 		stop = data.size();
Chris@4: 	    maxp = i;
Chris@4: 	    for (i++; i < stop; i++)
Chris@4: 		if (data[i] > data[maxp])
Chris@4: 		    maxp = i;
Chris@4: 	    if (maxp == mid) {
Chris@4: 		if (overThreshold(data, maxp, width, threshold, isRelative,av)){
Chris@4: 		    peaks.push_back(maxp);
Chris@4: 		}
Chris@4: 	    }
Chris@4: 	    mid++;
Chris@4: 	}
Chris@4: 	return peaks;
Chris@4:     } // findPeaks()
Chris@4: 
Chris@4:     static double expDecayWithHold(double av, double decayRate,
Chris@4: 				   const vector<double> &data, int start, int stop) {
Chris@4: 	while (start < stop) {
Chris@4: 	    av = decayRate * av + (1 - decayRate) * data[start];
Chris@4: 	    if (av < data[start])
Chris@4: 		av = data[start];
Chris@4: 	    start++;
Chris@4: 	}
Chris@4: 	return av;
Chris@4:     } // expDecayWithHold()
Chris@4: 
Chris@4:     static bool overThreshold(const vector<double> &data, int index, int width,
Chris@4: 				 double threshold, bool isRelative,
Chris@4: 				 double av) {
Chris@4: 	if (data[index] < av)
Chris@4: 	    return false;
Chris@4: 	if (isRelative) {
Chris@4: 	    int iStart = index - pre * width;
Chris@4: 	    if (iStart < 0)
Chris@4: 		iStart = 0;
Chris@4: 	    int iStop = index + post * width;
Chris@4: 	    if (iStop > data.size())
Chris@4: 		iStop = data.size();
Chris@4: 	    double sum = 0;
Chris@4: 	    int count = iStop - iStart;
Chris@4: 	    while (iStart < iStop)
Chris@4: 		sum += data[iStart++];
Chris@4: 	    return (data[index] > sum / count + threshold);
Chris@4: 	} else
Chris@4: 	    return (data[index] > threshold);
Chris@4:     } // overThreshold()
Chris@4: 
Chris@4:     static void normalise(vector<double> &data) {
Chris@4: 	double sx = 0;
Chris@4: 	double sxx = 0;
Chris@4: 	for (int i = 0; i < data.size(); i++) {
Chris@4: 	    sx += data[i];
Chris@4: 	    sxx += data[i] * data[i];
Chris@4: 	}
Chris@4: 	double mean = sx / data.size();
Chris@4: 	double sd = sqrt((sxx - sx * mean) / data.size());
Chris@4: 	if (sd == 0)
Chris@4: 	    sd = 1;		// all data[i] == mean  -> 0; avoids div by 0
Chris@4: 	for (int i = 0; i < data.size(); i++) {
Chris@4: 	    data[i] = (data[i] - mean) / sd;
Chris@4: 	}
Chris@4:     } // normalise()
Chris@4: 
Chris@4:     /** Uses an n-point linear regression to estimate the slope of data.
Chris@4:      *  @param data input data
Chris@4:      *  @param hop spacing of data points
Chris@4:      *  @param n length of linear regression
Chris@4:      *  @param slope output data
Chris@4:      */
Chris@4:     static void getSlope(const vector<double> &data, double hop, int n,
Chris@4: 			 vector<double> &slope) {
Chris@4: 	int i = 0, j = 0;
Chris@4: 	double t;
Chris@4: 	double sx = 0, sxx = 0, sy = 0, sxy = 0;
Chris@4: 	for ( ; i < n; i++) {
Chris@4: 	    t = i * hop;
Chris@4: 	    sx += t;
Chris@4: 	    sxx += t * t;
Chris@4: 	    sy += data[i];
Chris@4: 	    sxy += t * data[i];
Chris@4: 	}
Chris@4: 	double delta = n * sxx - sx * sx;
Chris@4: 	for ( ; j < n / 2; j++)
Chris@4: 	    slope[j] = (n * sxy - sx * sy) / delta;
Chris@4: 	for ( ; j < data.size() - (n + 1) / 2; j++, i++) {
Chris@4: 	    slope[j] = (n * sxy - sx * sy) / delta;
Chris@4: 	    sy += data[i] - data[i - n];
Chris@4: 	    sxy += hop * (n * data[i] - sy);
Chris@4: 	}
Chris@4: 	for ( ; j < data.size(); j++)
Chris@4: 	    slope[j] = (n * sxy - sx * sy) / delta;
Chris@4:     } // getSlope()
Chris@4: 
Chris@4:     static double min(const vector<double> &arr) { return arr[imin(arr)]; }
Chris@4: 
Chris@4:     static double max(const vector<double> &arr) { return arr[imax(arr)]; }
Chris@4: 
Chris@4:     static int imin(const vector<double> &arr) {
Chris@4: 	int i = 0;
Chris@4: 	for (int j = 1; j < arr.size(); j++)
Chris@4: 	    if (arr[j] < arr[i])
Chris@4: 		i = j;
Chris@4: 	return i;
Chris@4:     } // imin()
Chris@4: 
Chris@4:     static int imax(const vector<double> &arr) {
Chris@4: 	int i = 0;
Chris@4: 	for (int j = 1; j < arr.size(); j++)
Chris@4: 	    if (arr[j] > arr[i])
Chris@4: 		i = j;
Chris@4: 	return i;
Chris@4:     } // imax()
Chris@4: 
Chris@4: }; // class Peaks
Chris@4: 
Chris@4: #endif