plectrum

McLeodPitchMethod.java (12005B)

---

 /*
 *      _______                       _____   _____ _____  
 *     |__   __|                     |  __ \ / ____|  __ \ 
 *        | | __ _ _ __ ___  ___  ___| |  | | (___ | |__) |
 *        | |/ _` | '__/ __|/ _ \/ __| |  | |\___ \|  ___/ 
 *        | | (_| | |  \__ \ (_) \__ \ |__| |____) | |     
 *        |_|\__,_|_|  |___/\___/|___/_____/|_____/|_|     
 *                                                         
 * -------------------------------------------------------------
 *
 * TarsosDSP is developed by Joren Six at IPEM, University Ghent
 *  
 * -------------------------------------------------------------
 *
 *  Info: http://0110.be/tag/TarsosDSP
 *  Github: https://github.com/JorenSix/TarsosDSP
 *  Releases: http://0110.be/releases/TarsosDSP/
 *  
 *  TarsosDSP includes modified source code by various authors,
 *  for credits and info, see README.
 * 
 */
 
 
 /**
  */
 package be.tarsos.dsp.pitch;
 
 import java.util.ArrayList;
 import java.util.List;
 
 /**
  * <p>
  * Implementation of The McLeod Pitch Method (MPM). It is described in the
  * article <a href=
  * "http://miracle.otago.ac.nz/tartini/papers/A_Smarter_Way_to_Find_Pitch.pdf"
  * >A Smarter Way to Find Pitch</a>. According to the article:
  * </p>
  * <blockquote> <bufferCount>
  * <p>
  * A fast, accurate and robust method for finding the continuous pitch in
  * monophonic musical sounds. [It uses] a special normalized version of the
  * Squared Difference Function (SDF) coupled with a peak picking algorithm.
  * </p>
  * <p>
  * MPM runs in real time with a standard 44.1 kHz sampling rate. It operates
  * without using low-pass filtering so it can work on sound with high harmonic
  * frequencies such as a violin and it can display pitch changes of one cent
  * reliably. MPM works well without any post processing to correct the pitch.
  * </p>
  * </bufferCount> </blockquote>
  * <p>
  * For the moment this implementation uses the inefficient way of calculating
  * the pitch. It uses <code>O(Ww)</code> with W the window size in samples and w
  * the desired number of ACF coefficients. The implementation can be optimized
  * to <code>O((W+w)log(W+w))</code> by using an <abbr
  * title="Fast Fourier Transform">FFT</abbr> to calculate the <abbr
  * title="Auto-Correlation Function">ACF</abbr>. But I am still afraid of the
  * dark magic of the FFT and clinging to the familiar, friendly, laggard time
  * domain.
  * </p>
  * 
  * @author Phillip McLeod
  * @author Joren Six
  */
 public final class McLeodPitchMethod implements PitchDetector {
 
 	/**
 	 * The expected size of an audio buffer (in samples).
 	 */
 	public static final int DEFAULT_BUFFER_SIZE = 1024;
 
 	/**
 	 * Overlap defines how much two audio buffers following each other should
 	 * overlap (in samples). 75% overlap is advised in the MPM article.
 	 */
 	public static final int DEFAULT_OVERLAP = 768;
 
 	/**
 	 * Defines the relative size the chosen peak (pitch) has. 0.93 means: choose
 	 * the first peak that is higher than 93% of the highest peak detected. 93%
 	 * is the default value used in the Tartini user interface.
 	 */
 	private static final double DEFAULT_CUTOFF = 0.97;
 	/**
 	 * For performance reasons, peaks below this cutoff are not even considered.
 	 */
 	private static final double SMALL_CUTOFF = 0.5;
 
 	/**
 	 * Pitch annotations below this threshold are considered invalid, they are
 	 * ignored.
 	 */
 	private static final double LOWER_PITCH_CUTOFF = 80.0; // Hz
 
 	/**
 	 * Defines the relative size the chosen peak (pitch) has.
 	 */
 	private final double cutoff;
 
 	/**
 	 * The audio sample rate. Most audio has a sample rate of 44.1kHz.
 	 */
 	private final float sampleRate;
 
 	/**
 	 * Contains a normalized square difference function value for each delay
 	 * (tau).
 	 */
 	private final float[] nsdf;
 
 	/**
 	 * The x and y coordinate of the top of the curve (nsdf).
 	 */
 	private float turningPointX, turningPointY;
 
 	/**
 	 * A list with minimum and maximum values of the nsdf curve.
 	 */
 	private final List<Integer> maxPositions = new ArrayList<Integer>();
 
 	/**
 	 * A list of estimates of the period of the signal (in samples).
 	 */
 	private final List<Float> periodEstimates = new ArrayList<Float>();
 
 	/**
 	 * A list of estimates of the amplitudes corresponding with the period
 	 * estimates.
 	 */
 	private final List<Float> ampEstimates = new ArrayList<Float>();
 
 	/**
 	 * The result of the pitch detection iteration.
 	 */
 	private final PitchDetectionResult result;
 
 	/**
 	 * Initializes the normalized square difference value array and stores the
 	 * sample rate.
 	 * 
 	 * @param audioSampleRate
 	 *            The sample rate of the audio to check.
 	 */
 	public McLeodPitchMethod(final float audioSampleRate) {
 		this(audioSampleRate, DEFAULT_BUFFER_SIZE, DEFAULT_CUTOFF);
 	}
 
 	/**
 	 * Create a new pitch detector.
 	 * 
 	 * @param audioSampleRate
 	 *            The sample rate of the audio.
 	 * @param audioBufferSize
 	 *            The size of one audio buffer 1024 samples is common.
 	 */
 	public McLeodPitchMethod(final float audioSampleRate, final int audioBufferSize) {
 		this(audioSampleRate, audioBufferSize, DEFAULT_CUTOFF);
 	}
 
 	/**
 	 * Create a new pitch detector.
 	 * 
 	 * @param audioSampleRate
 	 *            The sample rate of the audio.
 	 * @param audioBufferSize
 	 *            The size of one audio buffer 1024 samples is common.
 	 * @param cutoffMPM
 	 *            The cutoff (similar to the YIN threshold). In the Tartini
 	 *            paper 0.93 is used.
 	 */
 	public McLeodPitchMethod(final float audioSampleRate, final int audioBufferSize, final double cutoffMPM) {
 		this.sampleRate = audioSampleRate;
 		nsdf = new float[audioBufferSize];
 		this.cutoff = cutoffMPM;
 		result = new PitchDetectionResult();
 	}
 
 	/**
 	 * Implements the normalized square difference function. See section 4 (and
 	 * the explanation before) in the MPM article. This calculation can be
 	 * optimized by using an FFT. The results should remain the same.
 	 * 
 	 * @param audioBuffer
 	 *            The buffer with audio information.
 	 */
 	private void normalizedSquareDifference(final float[] audioBuffer) {
 		for (int tau = 0; tau < audioBuffer.length; tau++) {
 			float acf = 0;
 			float divisorM = 0;
 			for (int i = 0; i < audioBuffer.length - tau; i++) {
 				acf += audioBuffer[i] * audioBuffer[i + tau];
 				divisorM += audioBuffer[i] * audioBuffer[i] + audioBuffer[i + tau] * audioBuffer[i + tau];
 			}
 			nsdf[tau] = 2 * acf / divisorM;
 		}
 	}
 
 	/*
 	 * (non-Javadoc)
 	 * 
 	 * @see be.tarsos.pitch.pure.PurePitchDetector#getPitch(float[])
 	 */
 	public PitchDetectionResult getPitch(final float[] audioBuffer) {
 		final float pitch;
 
 		// 0. Clear previous results (Is this faster than initializing a list
 		// again and again?)
 		maxPositions.clear();
 		periodEstimates.clear();
 		ampEstimates.clear();
 
 		// 1. Calculate the normalized square difference for each Tau value.
 		normalizedSquareDifference(audioBuffer);
 		// 2. Peak picking time: time to pick some peaks.
 		peakPicking();
 
 		double highestAmplitude = Double.NEGATIVE_INFINITY;
 
 		for (final Integer tau : maxPositions) {
 			// make sure every annotation has a probability attached
 			highestAmplitude = Math.max(highestAmplitude, nsdf[tau]);
 
 			if (nsdf[tau] > SMALL_CUTOFF) {
 				// calculates turningPointX and Y
 				parabolicInterpolation(tau);
 				// store the turning points
 				ampEstimates.add(turningPointY);
 				periodEstimates.add(turningPointX);
 				// remember the highest amplitude
 				highestAmplitude = Math.max(highestAmplitude, turningPointY);
 			}
 		}
 
 		if (periodEstimates.isEmpty()) {
 			pitch = -1;
 		} else {
 			// use the overall maximum to calculate a cutoff.
 			// The cutoff value is based on the highest value and a relative
 			// threshold.
 			final double actualCutoff = cutoff * highestAmplitude;
 
 			// find first period above or equal to cutoff
 			int periodIndex = 0;
 			for (int i = 0; i < ampEstimates.size(); i++) {
 				if (ampEstimates.get(i) >= actualCutoff) {
 					periodIndex = i;
 					break;
 				}
 			}
 
 			final double period = periodEstimates.get(periodIndex);
 			final float pitchEstimate = (float) (sampleRate / period);
 			if (pitchEstimate > LOWER_PITCH_CUTOFF) {
 				pitch = pitchEstimate;
 			} else {
 				pitch = -1;
 			}
 
 		}
 		result.setProbability((float) highestAmplitude);
 		result.setPitch(pitch);
 		result.setPitched(pitch != -1);
 		
 		return result;
 	}
 
 	/**
 	 * <p>
 	 * Finds the x value corresponding with the peak of a parabola.
 	 * </p>
 	 * <p>
 	 * a,b,c are three samples that follow each other. E.g. a is at 511, b at
 	 * 512 and c at 513; f(a), f(b) and f(c) are the normalized square
 	 * difference values for those samples; x is the peak of the parabola and is
 	 * what we are looking for. Because the samples follow each other
 	 * <code>b - a = 1</code> the formula for <a
 	 * href="http://fizyka.umk.pl/nrbook/c10-2.pdf">parabolic interpolation</a>
 	 * can be simplified a lot.
 	 * </p>
 	 * <p>
 	 * The following ASCII ART shows it a bit more clear, imagine this to be a
 	 * bit more curvaceous.
 	 * </p>
 	 * 
 	 * <pre>
 	 *     nsdf(x)
 	 *       ^
 	 *       |
 	 * f(x)  |------ ^
 	 * f(b)  |     / |\
 	 * f(a)  |    /  | \
 	 *       |   /   |  \
 	 *       |  /    |   \
 	 * f(c)  | /     |    \
 	 *       |_____________________> x
 	 *            a  x b  c
 	 * </pre>
 	 * 
 	 * @param tau
 	 *            The delay tau, b value in the drawing is the tau value.
 	 */
 	private void parabolicInterpolation(final int tau) {
 		final float nsdfa = nsdf[tau - 1];
 		final float nsdfb = nsdf[tau];
 		final float nsdfc = nsdf[tau + 1];
 		final float bValue = tau;
 		final float bottom = nsdfc + nsdfa - 2 * nsdfb;
 		if (bottom == 0.0) {
 			turningPointX = bValue;
 			turningPointY = nsdfb;
 		} else {
 			final float delta = nsdfa - nsdfc;
 			turningPointX = bValue + delta / (2 * bottom);
 			turningPointY = nsdfb - delta * delta / (8 * bottom);
 		}
 	}
 
 	/**
 	 * <p>
 	 * Implementation based on the GPL'ED code of <a
 	 * href="http://tartini.net">Tartini</a> This code can be found in the file
 	 * <code>general/mytransforms.cpp</code>.
 	 * </p>
 	 * <p>
 	 * Finds the highest value between each pair of positive zero crossings.
 	 * Including the highest value between the last positive zero crossing and
 	 * the end (if any). Ignoring the first maximum (which is at zero). In this
 	 * diagram the desired values are marked with a +
 	 * </p>
 	 * 
 	 * <pre>
 	 *  f(x)
 	 *   ^
 	 *   |
 	 *  1|               +
 	 *   | \      +     /\      +     /\
 	 *  0| _\____/\____/__\/\__/\____/_______> x
 	 *   |   \  /  \  /      \/  \  /
 	 * -1|    \/    \/            \/
 	 *   |
 	 * </pre>
 	 * 
 	 * @param nsdf
 	 *            The array to look for maximum values in. It should contain
 	 *            values between -1 and 1
 	 * @author Phillip McLeod
 	 */
 	private void peakPicking() {
 
 		int pos = 0;
 		int curMaxPos = 0;
 
 		// find the first negative zero crossing
 		while (pos < (nsdf.length - 1) / 3 && nsdf[pos] > 0) {
 			pos++;
 		}
 
 		// loop over all the values below zero
 		while (pos < nsdf.length - 1 && nsdf[pos] <= 0.0) {
 			pos++;
 		}
 
 		// can happen if output[0] is NAN
 		if (pos == 0) {
 			pos = 1;
 		}
 
 		while (pos < nsdf.length - 1) {
 			assert nsdf[pos] >= 0;
 			if (nsdf[pos] > nsdf[pos - 1] && nsdf[pos] >= nsdf[pos + 1]) {
 				if (curMaxPos == 0) {
 					// the first max (between zero crossings)
 					curMaxPos = pos;
 				} else if (nsdf[pos] > nsdf[curMaxPos]) {
 					// a higher max (between the zero crossings)
 					curMaxPos = pos;
 				}
 			}
 			pos++;
 			// a negative zero crossing
 			if (pos < nsdf.length - 1 && nsdf[pos] <= 0) {
 				// if there was a maximum add it to the list of maxima
 				if (curMaxPos > 0) {
 					maxPositions.add(curMaxPos);
 					curMaxPos = 0; // clear the maximum position, so we start
 					// looking for a new ones
 				}
 				while (pos < nsdf.length - 1 && nsdf[pos] <= 0.0f) {
 					pos++; // loop over all the values below zero
 				}
 			}
 		}
 		if (curMaxPos > 0) { // if there was a maximum in the last part
 			maxPositions.add(curMaxPos); // add it to the vector of maxima
 		}
 	}
 }