/*
*
      
_______
                       
_____
   
_____ _____
  

*
     
|__
   
__|
                     
|
  
__ \ / ____|__ \
 

*
        
| | __ _ _ __ ___
  
______| || | (___ | |__) |
*
        
| |/ _` | '__/ __|/ _ \/ __| |
  
| |\___ \|___/
 

*
        
| | (_| | |
  
\__ \ (_) \__ \ |__| |____) | |
     

*
        
|_|\__,_|_|
  
|___/\___/|___/_____/|_____/|_|
     

*
                                                         

* -------------------------------------------------------------
*
* 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 bc
	 
* </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
		
}
	
}
}