view Peaks.h @ 4:c06cf6f7cb04

Bring in Peaks code from BeatRoot

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

/*
  Vamp feature extraction plugin for the BeatRoot beat tracker.

  Centre for Digital Music, Queen Mary, University of London.
  This file copyright 2011 Simon Dixon, Chris Cannam and QMUL.
    
  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.  See the file
  COPYING included with this distribution for more information.
*/

#ifndef _BEATROOT_PEAKS_H_
#define _BEATROOT_PEAKS_H_

#include <vector>
#include <cmath>

using std::vector;

class Peaks
{
protected:
    static int pre;
    static int post;
	
public:
    /** General peak picking method for finding n local maxima in an array
     *  @param data input data
     *  @param peaks list of peak indexes
     *  @param width minimum distance between peaks
     */
    static int findPeaks(const vector<double> &data, vector<int> peaks, int width) {
	int peakCount = 0;
	int maxp = 0;
	int mid = 0;
	int end = data.size();
	while (mid < end) {
	    int i = mid - width;
	    if (i < 0)
		i = 0;
	    int stop = mid + width + 1;
	    if (stop > data.size())
		stop = data.size();
	    maxp = i;
	    for (i++; i < stop; i++)
		if (data[i] > data[maxp])
		    maxp = i;
	    if (maxp == mid) {
		int j;
		for (j = peakCount; j > 0; j--) {
		    if (data[maxp] <= data[peaks[j-1]])
			break;
		    else if (j < peaks.size())
			peaks[j] = peaks[j-1];
		}
		if (j != peaks.size())
		    peaks[j] = maxp;
		if (peakCount != peaks.size())
		    peakCount++;
	    }
	    mid++;
	}
	return peakCount;
    } // findPeaks()

    /** General peak picking method for finding local maxima in an array
     *  @param data input data
     *  @param width minimum distance between peaks
     *  @param threshold minimum value of peaks
     *  @return list of peak indexes
     */ 
    static vector<int> findPeaks(const vector<double> &data, int width,
				 double threshold) {
	return findPeaks(data, width, threshold, 0, false);
    } // findPeaks()
	
    /** General peak picking method for finding local maxima in an array
     *  @param data input data
     *  @param width minimum distance between peaks
     *  @param threshold minimum value of peaks
     *  @param decayRate how quickly previous peaks are forgotten
     *  @param isRelative minimum value of peaks is relative to local average
     *  @return list of peak indexes
     */
    static vector<int> findPeaks(const vector<double> &data, int width,
				 double threshold, double decayRate, bool isRelative) {
	vector<int> peaks;
	int maxp = 0;
	int mid = 0;
	int end = data.size();
	double av = data[0];
	while (mid < end) {
	    av = decayRate * av + (1 - decayRate) * data[mid];
	    if (av < data[mid])
		av = data[mid];
	    int i = mid - width;
	    if (i < 0)
		i = 0;
	    int stop = mid + width + 1;
	    if (stop > data.size())
		stop = data.size();
	    maxp = i;
	    for (i++; i < stop; i++)
		if (data[i] > data[maxp])
		    maxp = i;
	    if (maxp == mid) {
		if (overThreshold(data, maxp, width, threshold, isRelative,av)){
		    peaks.push_back(maxp);
		}
	    }
	    mid++;
	}
	return peaks;
    } // findPeaks()

    static double expDecayWithHold(double av, double decayRate,
				   const vector<double> &data, int start, int stop) {
	while (start < stop) {
	    av = decayRate * av + (1 - decayRate) * data[start];
	    if (av < data[start])
		av = data[start];
	    start++;
	}
	return av;
    } // expDecayWithHold()

    static bool overThreshold(const vector<double> &data, int index, int width,
				 double threshold, bool isRelative,
				 double av) {
	if (data[index] < av)
	    return false;
	if (isRelative) {
	    int iStart = index - pre * width;
	    if (iStart < 0)
		iStart = 0;
	    int iStop = index + post * width;
	    if (iStop > data.size())
		iStop = data.size();
	    double sum = 0;
	    int count = iStop - iStart;
	    while (iStart < iStop)
		sum += data[iStart++];
	    return (data[index] > sum / count + threshold);
	} else
	    return (data[index] > threshold);
    } // overThreshold()

    static void normalise(vector<double> &data) {
	double sx = 0;
	double sxx = 0;
	for (int i = 0; i < data.size(); i++) {
	    sx += data[i];
	    sxx += data[i] * data[i];
	}
	double mean = sx / data.size();
	double sd = sqrt((sxx - sx * mean) / data.size());
	if (sd == 0)
	    sd = 1;		// all data[i] == mean  -> 0; avoids div by 0
	for (int i = 0; i < data.size(); i++) {
	    data[i] = (data[i] - mean) / sd;
	}
    } // normalise()

    /** Uses an n-point linear regression to estimate the slope of data.
     *  @param data input data
     *  @param hop spacing of data points
     *  @param n length of linear regression
     *  @param slope output data
     */
    static void getSlope(const vector<double> &data, double hop, int n,
			 vector<double> &slope) {
	int i = 0, j = 0;
	double t;
	double sx = 0, sxx = 0, sy = 0, sxy = 0;
	for ( ; i < n; i++) {
	    t = i * hop;
	    sx += t;
	    sxx += t * t;
	    sy += data[i];
	    sxy += t * data[i];
	}
	double delta = n * sxx - sx * sx;
	for ( ; j < n / 2; j++)
	    slope[j] = (n * sxy - sx * sy) / delta;
	for ( ; j < data.size() - (n + 1) / 2; j++, i++) {
	    slope[j] = (n * sxy - sx * sy) / delta;
	    sy += data[i] - data[i - n];
	    sxy += hop * (n * data[i] - sy);
	}
	for ( ; j < data.size(); j++)
	    slope[j] = (n * sxy - sx * sy) / delta;
    } // getSlope()

    static double min(const vector<double> &arr) { return arr[imin(arr)]; }

    static double max(const vector<double> &arr) { return arr[imax(arr)]; }

    static int imin(const vector<double> &arr) {
	int i = 0;
	for (int j = 1; j < arr.size(); j++)
	    if (arr[j] < arr[i])
		i = j;
	return i;
    } // imin()

    static int imax(const vector<double> &arr) {
	int i = 0;
	for (int j = 1; j < arr.size(); j++)
	    if (arr[j] > arr[i])
		i = j;
	return i;
    } // imax()

}; // class Peaks

#endif