Chris@2: /*
Chris@2: 	Copyright (C) 2001, 2006 by Simon Dixon
Chris@2: 
Chris@2: 	This program is free software; you can redistribute it and/or modify
Chris@2: 	it under the terms of the GNU General Public License as published by
Chris@2: 	the Free Software Foundation; either version 2 of the License, or
Chris@2: 	(at your option) any later version.
Chris@2: 
Chris@2: 	This program is distributed in the hope that it will be useful,
Chris@2: 	but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@2: 	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@2: 	GNU General Public License for more details.
Chris@2: 
Chris@2: 	You should have received a copy of the GNU General Public License along
Chris@2: 	with this program (the file gpl.txt); if not, download it from
Chris@2: 	http://www.gnu.org/licenses/gpl.txt or write to the
Chris@2: 	Free Software Foundation, Inc.,
Chris@2: 	51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
Chris@2: */
Chris@2: 
Chris@2: package at.ofai.music.util;
Chris@2: 
Chris@2: import java.util.LinkedList;
Chris@2: 
Chris@2: public class Peaks {
Chris@2: 
Chris@2: 	public static boolean debug = false;
Chris@2: 	public static int pre = 3;
Chris@2: 	public static int post = 1;
Chris@2: 	
Chris@2: 	/** General peak picking method for finding n local maxima in an array
Chris@2: 	 *  @param data input data
Chris@2: 	 *  @param peaks list of peak indexes
Chris@2: 	 *  @param width minimum distance between peaks
Chris@2: 	 */
Chris@2: 	public static int findPeaks(double[] data, int[] peaks, int width) {
Chris@2: 		int peakCount = 0;
Chris@2: 		int maxp = 0;
Chris@2: 		int mid = 0;
Chris@2: 		int end = data.length;
Chris@2: 		while (mid < end) {
Chris@2: 			int i = mid - width;
Chris@2: 			if (i < 0)
Chris@2: 				i = 0;
Chris@2: 			int stop = mid + width + 1;
Chris@2: 			if (stop > data.length)
Chris@2: 				stop = data.length;
Chris@2: 			maxp = i;
Chris@2: 			for (i++; i < stop; i++)
Chris@2: 				if (data[i] > data[maxp])
Chris@2: 					maxp = i;
Chris@2: 			if (maxp == mid) {
Chris@2: 				int j;
Chris@2: 				for (j = peakCount; j > 0; j--) {
Chris@2: 					if (data[maxp] <= data[peaks[j-1]])
Chris@2: 						break;
Chris@2: 					else if (j < peaks.length)
Chris@2: 						peaks[j] = peaks[j-1];
Chris@2: 				}
Chris@2: 				if (j != peaks.length)
Chris@2: 					peaks[j] = maxp;
Chris@2: 				if (peakCount != peaks.length)
Chris@2: 					peakCount++;
Chris@2: 			}
Chris@2: 			mid++;
Chris@2: 		}
Chris@2: 		return peakCount;
Chris@2: 	} // findPeaks()
Chris@2: 
Chris@2: 	/** General peak picking method for finding local maxima in an array
Chris@2: 	 *  @param data input data
Chris@2: 	 *  @param width minimum distance between peaks
Chris@2: 	 *  @param threshold minimum value of peaks
Chris@2: 	 *  @return list of peak indexes
Chris@2: 	 */
Chris@2: 	public static LinkedList<Integer> findPeaks(double[] data, int width,
Chris@2: 												double threshold) {
Chris@2: 		return findPeaks(data, width, threshold, 0, false);
Chris@2: 	} // findPeaks()
Chris@2: 	
Chris@2: 	/** General peak picking method for finding local maxima in an array
Chris@2: 	 *  @param data input data
Chris@2: 	 *  @param width minimum distance between peaks
Chris@2: 	 *  @param threshold minimum value of peaks
Chris@2: 	 *  @param decayRate how quickly previous peaks are forgotten
Chris@2: 	 *  @param isRelative minimum value of peaks is relative to local average
Chris@2: 	 *  @return list of peak indexes
Chris@2: 	 */
Chris@2: 	public static LinkedList<Integer> findPeaks(double[] data, int width,
Chris@2: 				double threshold, double decayRate, boolean isRelative) {
Chris@2: 		LinkedList<Integer> peaks = new LinkedList<Integer>();
Chris@2: 		int maxp = 0;
Chris@2: 		int mid = 0;
Chris@2: 		int end = data.length;
Chris@2: 		double av = data[0];
Chris@2: 		while (mid < end) {
Chris@2: 			av = decayRate * av + (1 - decayRate) * data[mid];
Chris@2: 			if (av < data[mid])
Chris@2: 				av = data[mid];
Chris@2: 			int i = mid - width;
Chris@2: 			if (i < 0)
Chris@2: 				i = 0;
Chris@2: 			int stop = mid + width + 1;
Chris@2: 			if (stop > data.length)
Chris@2: 				stop = data.length;
Chris@2: 			maxp = i;
Chris@2: 			for (i++; i < stop; i++)
Chris@2: 				if (data[i] > data[maxp])
Chris@2: 					maxp = i;
Chris@2: 			if (maxp == mid) {
Chris@2: 				if (overThreshold(data, maxp, width, threshold, isRelative,av)){
Chris@2: 					if (debug)
Chris@2: 						System.out.println(" peak");
Chris@2: 					peaks.add(new Integer(maxp));
Chris@2: 				} else if (debug)
Chris@2: 					System.out.println();
Chris@2: 			}
Chris@2: 			mid++;
Chris@2: 		}
Chris@2: 		return peaks;
Chris@2: 	} // findPeaks()
Chris@2: 
Chris@2: 	public static double expDecayWithHold(double av, double decayRate,
Chris@2: 										  double[] data, int start, int stop) {
Chris@2: 		while (start < stop) {
Chris@2: 			av = decayRate * av + (1 - decayRate) * data[start];
Chris@2: 			if (av < data[start])
Chris@2: 				av = data[start];
Chris@2: 			start++;
Chris@2: 		}
Chris@2: 		return av;
Chris@2: 	} // expDecayWithHold()
Chris@2: 
Chris@2: 	public static boolean overThreshold(double[] data, int index, int width,
Chris@2: 										double threshold, boolean isRelative,
Chris@2: 										double av) {
Chris@2: 		if (debug)
Chris@2: 			System.out.printf("%4d : %6.3f     Av1: %6.3f    ",
Chris@2: 								index, data[index], av);
Chris@2: 		if (data[index] < av)
Chris@2: 			return false;
Chris@2: 		if (isRelative) {
Chris@2: 			int iStart = index - pre * width;
Chris@2: 			if (iStart < 0)
Chris@2: 				iStart = 0;
Chris@2: 			int iStop = index + post * width;
Chris@2: 			if (iStop > data.length)
Chris@2: 				iStop = data.length;
Chris@2: 			double sum = 0;
Chris@2: 			int count = iStop - iStart;
Chris@2: 			while (iStart < iStop)
Chris@2: 				sum += data[iStart++];
Chris@2: 			if (debug)
Chris@2: 				System.out.printf("    %6.3f    %6.3f   ", sum / count,
Chris@2: 							data[index] - sum / count - threshold);
Chris@2: 			return (data[index] > sum / count + threshold);
Chris@2: 		} else
Chris@2: 			return (data[index] > threshold);
Chris@2: 	} // overThreshold()
Chris@2: 
Chris@2: 	public static void normalise(double[] data) {
Chris@2: 		double sx = 0;
Chris@2: 		double sxx = 0;
Chris@2: 		for (int i = 0; i < data.length; i++) {
Chris@2: 			sx += data[i];
Chris@2: 			sxx += data[i] * data[i];
Chris@2: 		}
Chris@2: 		double mean = sx / data.length;
Chris@2: 		double sd = Math.sqrt((sxx - sx * mean) / data.length);
Chris@2: 		if (sd == 0)
Chris@2: 			sd = 1;		// all data[i] == mean  -> 0; avoids div by 0
Chris@2: 		for (int i = 0; i < data.length; i++) {
Chris@2: 			data[i] = (data[i] - mean) / sd;
Chris@2: 		}
Chris@2: 	} // normalise()
Chris@2: 
Chris@2: 	/** Uses an n-point linear regression to estimate the slope of data.
Chris@2: 	 *  @param data input data
Chris@2: 	 *  @param hop spacing of data points
Chris@2: 	 *  @param n length of linear regression
Chris@2: 	 *  @param slope output data
Chris@2: 	 */
Chris@2: 	public static void getSlope(double[] data, double hop, int n,
Chris@2: 								double[] slope) {
Chris@2: 		int i = 0, j = 0;
Chris@2: 		double t;
Chris@2: 		double sx = 0, sxx = 0, sy = 0, sxy = 0;
Chris@2: 		for ( ; i < n; i++) {
Chris@2: 			t = i * hop;
Chris@2: 			sx += t;
Chris@2: 			sxx += t * t;
Chris@2: 			sy += data[i];
Chris@2: 			sxy += t * data[i];
Chris@2: 		}
Chris@2: 		double delta = n * sxx - sx * sx;
Chris@2: 		for ( ; j < n / 2; j++)
Chris@2: 			slope[j] = (n * sxy - sx * sy) / delta;
Chris@2: 		for ( ; j < data.length - (n + 1) / 2; j++, i++) {
Chris@2: 			slope[j] = (n * sxy - sx * sy) / delta;
Chris@2: 			sy += data[i] - data[i - n];
Chris@2: 			sxy += hop * (n * data[i] - sy);
Chris@2: 		}
Chris@2: 		for ( ; j < data.length; j++)
Chris@2: 			slope[j] = (n * sxy - sx * sy) / delta;
Chris@2: 	} // getSlope()
Chris@2: 
Chris@2: 	public static double min(double[] arr) { return arr[imin(arr)]; }
Chris@2: 
Chris@2: 	public static double max(double[] arr) { return arr[imax(arr)]; }
Chris@2: 
Chris@2: 	public static int imin(double[] arr) {
Chris@2: 		int i = 0;
Chris@2: 		for (int j = 1; j < arr.length; j++)
Chris@2: 			if (arr[j] < arr[i])
Chris@2: 				i = j;
Chris@2: 		return i;
Chris@2: 	} // imin()
Chris@2: 
Chris@2: 	public static int imax(double[] arr) {
Chris@2: 		int i = 0;
Chris@2: 		for (int j = 1; j < arr.length; j++)
Chris@2: 			if (arr[j] > arr[i])
Chris@2: 				i = j;
Chris@2: 		return i;
Chris@2: 	} // imax()
Chris@2: 
Chris@2: } // class Peaks