import com.cycling74.max.*;
import com.cycling74.msp.*;

public class pitch_YIN extends MSPPerformer
{
	private float _gain = 0.f;
	private float pitch;

	// YIN variables
	private static float[] df = new float[512];
	private static float[] cmnd = new float[512];
	private static float[] cand_array = new float[500];
	private static float[] f = new float[3];

	private static final String[] INLET_ASSIST = new String[]{
		"input (sig)"
	};
	private static final String[] OUTLET_ASSIST = new String[]{
		"output (sig)"
	};
	

	public pitch_YIN(float gain)
	{
		declareInlets(new int[]{SIGNAL,DataTypes.ALL});
		declareOutlets(new int[]{SIGNAL,DataTypes.ALL});

		setInletAssist(INLET_ASSIST);
		setOutletAssist(OUTLET_ASSIST);

		_gain = gain;
	}
    
	public void inlet(float f)
	{
		if (getInlet() == 1)
		{	
			_gain = f;
		}
	}

	public void dspsetup(MSPSignal[] ins, MSPSignal[] outs)
	{
		//If you forget the fields of MSPSignal you can select the classname above
		//and choose Open Class Reference For Selected Class.. from the Java menu
	}

	public void perform(MSPSignal[] ins, MSPSignal[] outs)
	{
		
		int i;
		float[] in = ins[0].vec;
		float[] out = outs[0].vec;
		for(i = 0; i < in.length;i++)
		{
			/*do something*/
			out[i] = in[i] * _gain;  	

		}

		pitch = pitch_YIN(in,in.length,44100.f);


		outlet(1,pitch);		
	}

	// YIN pitch detection algorithm - x is the signal, n the frame size and fs the sampling freq
	public static float pitch_YIN(float[] x,int n, float fs)
	{
		float sum = 0; // to hold sum value
		int j;	// array index
		float threshold = 0.1f;		// threshold for finding minima
		int numcands = 0;	// number of candidate lags
		float min = 100000;		// to find minimum value - set to 100000 so that everything will likely be smaller than it
		float mindex = -2.f; // to hold index of minimum point of cmnd function. set at -2 to allow us to check if it hasn't been set
		float a,b,c,interpindex,f0;	// for interpolation. f0 is final fundamental freq

		for (int t = 1;t < n;t++)
		{
			j = t-1;	// t is lag value but we need to store from array value 0
			df[j] = 0;
			
			for (int i=0;i < (n-t);i++)
			{
				df[j] = df[j] + ((float) Math.pow((x[i]-x[i+t]),2));  // difference function calculation
			}
		
			df[j] = df[j] / (n-t);  // weight difference function by number of ops
    
			sum = sum + df[j];   // add differnce function value to sum

			cmnd[j] = df[j] / ((1/((float) t))*sum);  // divide difference function by average so far
		}
		
		
		
		// Absolute Threshold
		// set threshold
		// find all minima below threshold
		// pick minima with lowest lag
		for (int i = 1; i < (n-2);i++)
		{
			if (cmnd[i] < threshold) // if below threshold
			{     
				if ((cmnd[i-1] > cmnd[i]) && (cmnd[i] < cmnd[i+1])) // if at a minimum
				{
					cand_array[numcands] = i;	// add index to candidate array
					numcands++; // increment candidate array counter
				}       
			}
		}
		
		
		
		// now work out index of minimum point
		if (numcands == 0) // if f0 candidate array is empty
		{
			for (int i = 0;i < cmnd.length;i++)
			{
				if (cmnd[i] < min)
				{
					mindex = (float) i; // just take index of minimum point of cmnd
					min = cmnd[i];
				}
			}
		}
		else
		{
			mindex = cand_array[0]; // else take candidate with smallest lag
			
			/* alernate method that searches for smallest - but unnecessary as candidates are in order	
			for (int i = 0;i < numcands;i++)
			{
				if (cand_array[i] < min)
				{
					mindex = cand_array[i]; // else take candidate with smallest lag
					min = cand_array[i];
				}
			}
			*/		
		}
		
	
		// parabolic interpolation between three points of minimum
		// this refines the answer to close to the true minimum
		// removes limitation of sampling frequency
	
		if (mindex > 0 && mindex < cmnd.length) // check index isn't first or last of function
		{
			// a, b and c are the indicies
			a = mindex-1;
			b = mindex;
			c = mindex+1;
    
			// f(0), f(1) and f(2) are the function values at a,b and c
			f[0] = cmnd[(int) a];
			f[1] = cmnd[(int) b];
			f[2] = cmnd[(int) c];
    
			// find index between a and c that has the 
			// lowest val according to  parabolic interpolation
			interpindex = b - (0.5f * ((((float) Math.pow((b-a),2)) * (f[1] - f[2])) - (((float) Math.pow((b-c),2)) * (f[1] - f[0]))) / (((b-a) * (f[1] - f[2])) - ((b-c) * (f[1] - f[0]))));
    
			interpindex = interpindex + 1;	// add 1 to index to account for arrays going from 0..n rather than 1..n
			f0 = fs / interpindex;
		}
		else // if index is first or last of frame, just use it to find freq
		{
			mindex = mindex + 1;	// add 1 to index to account for arrays going from 0..n rather than 1..n
			f0 = fs / mindex;
		}
	
		return f0;
	}
}