view src/DynamicBayesianArray.cpp @ 1:1a32ce016bb9

Changed bestEstimate timing to work via time sent from Max not the elapsed time. This had caused some problems, but this version now working surprisingly well on MIDI files with variable timing.

/*
 *  DynamicDynamicBayesianArray.cpp
 *  midiCannamReader
 *
 *  Created by Andrew on 17/07/2011.
 *  Copyright 2011 QMUL. All rights reserved.
 *
 */

#include "DynamicBayesianArray.h"
#include "math.h"
#include "ofMain.h"

DynamicBayesianArray::DynamicBayesianArray(){
	
//	prior.createVector(240);
//	likelihood.createVector(240);
//	posterior.createVector(240);	
	testVector.createVector(240);
	testVector.addGaussianShape(100,10, 0.1);
	
	
	likelihoodNoise = 0.5;
	likelihoodMean = ARRAY_SIZE/2;
	likelihoodStdDev = ARRAY_SIZE / 12;
	initialiseArray();
}

void DynamicBayesianArray::initialiseArray(){
	
	//maximumIndex = 12;//change this
	setGaussianPrior(ARRAY_SIZE/2, ARRAY_SIZE/1);	
	setGaussianLikelihood(ARRAY_SIZE/2, ARRAY_SIZE/1);//likelihoodMean, likelihoodStdDev);
	
	calculatePosterior();
	renormalisePosterior();
	posteriorDecayRate = 0.06;
	
	eighthNoteProportion = 0.35;//must be less than 0.5 to discriminate - was 0.4
	earlySixteenthNoteProportion = 0;
	lateSixteenthNoteProportion = 0;
	decayNoiseAmount = 0.1;
	decayNoiseStdDev = ARRAY_SIZE/24;
	standardDeviation = likelihoodStdDev;
	setDecayNoiseGaussian(ARRAY_SIZE/2, decayNoiseStdDev);
	
	setGaussianLikelihood(likelihoodMean, likelihoodStdDev);
}


void DynamicBayesianArray::setGaussianPrior(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		prior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
		//posterior[i] = prior[i];
	}
}

void DynamicBayesianArray::setGaussianPosterior(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}


void DynamicBayesianArray::setGaussianLikelihood(float mean, float StdDev){
	if (mean >= 0 && mean <= ARRAY_SIZE){
		int i;	float eighthDifference;
		int eighthPosition = ((int)mean + ARRAY_SIZE/2)%ARRAY_SIZE;
		float mainDifference;
		float gaussianProportion = 1 - likelihoodNoise;
		
		for (i=0;i < ARRAY_SIZE;i++){
			
			mainDifference = min( double(fabs(i-mean)) , (double)(i + ARRAY_SIZE - mean));
			//without * (1 - eighthNoteProportion)
			likelihood[i] = gaussianProportion  * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
			
			likelihood[i] += (likelihoodNoise / ARRAY_SIZE);
			//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
			//(double) (likelihoodNoise / ARRAY_SIZE) );
		}
		//	renormaliseArray(&likelihood[0], ARRAY_SIZE);
	}//end if mean within limits
}

void DynamicBayesianArray::calculatePosterior(){
	int i;
	for (i=0;i < ARRAY_SIZE;i++){
		posterior[i] = likelihood[i] * prior[i];
	}
	//renormalisePosterior();
}


float DynamicBayesianArray::getMaximum(float *ptr, int length){
	int i;
	float max = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max)
			max = *(ptr+i);
	}
	maximumValue = max;
	return max;
}

float* DynamicBayesianArray::getMaximumEstimate(float *ptr, int length){
	float returnArray[2];
	int i;
	float max = 0;
	maximumIndex = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max){
			max = *(ptr+i);
			maximumIndex = i;
		}
	}
	returnArray[0] = max;
	returnArray[1] = maximumIndex;
	maximumValue = max;
	return &returnArray[0];
}



double DynamicBayesianArray::getIntegratedEstimateIndex(){
	int i;
	float integratedQuantity = 0;
	float integratedTotal = 0;
	double integratedIndex = 0;
	for (i=0;i < ARRAY_SIZE;i++){	
		integratedQuantity += posterior[i];//the values of the probability distribution
		integratedTotal += i*posterior[i];
	}
	if (integratedQuantity > 0){
		integratedIndex = integratedTotal / integratedQuantity;
	}
	integratedEstimate = (float) integratedIndex;
	return integratedIndex;
}


double DynamicBayesianArray::calculateStandardDeviation(){
	
	double total = 0;
	double pdfSum;
	double variance = 0;
	for (int i=0;i < ARRAY_SIZE;i++){	
		//*posterior[i] * 		
		total += posterior[i] * (i - integratedEstimate) * (i - integratedEstimate);//the values of the probability distribution
		pdfSum += posterior[i];
	}
	
	if (pdfSum > 0)
		variance = total / pdfSum;
	else 
		variance = ARRAY_SIZE;
	
	standardDeviation = sqrt(variance);
	return standardDeviation;
}



void DynamicBayesianArray::renormaliseArray(float *ptr, int length){
	int i;
	float totalArea = 0;
	for (i=0;i < length;i++){
		totalArea += *(ptr+i);
	}
	
	for (i=0;i < length;i++){
		*(ptr+i) /= totalArea;
	}
	
}

void DynamicBayesianArray::resetPrior(){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		prior[i] = posterior[i];
	}
}

void DynamicBayesianArray::renormalisePosterior(){
	int i;
	float totalArea = 0;
	for (i=0;i < ARRAY_SIZE;i++){
		totalArea += posterior[i];
	}
	
	for (i=0;i < ARRAY_SIZE;i++){
		posterior[i] /= totalArea;
	}
	
}

void DynamicBayesianArray::decayPosterior(){
	float *pointer;
	pointer = getMaximumEstimate(&posterior[0], ARRAY_SIZE);	
	float maximum;
	maximum = *pointer;
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] += (maximum - posterior[i]) * posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}
	maximumIndex = *(pointer+1);
}

void DynamicBayesianArray::setDecayNoiseGaussian(float mean, float StdDev){
	int i;
	for (i=0;i<ARRAY_SIZE;i++){
		decayNoiseArray[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}

void DynamicBayesianArray::decayPosteriorWithGaussianNoise(){
	
	int i;
	float currentMaximum = getMaximum(&posterior[0], ARRAY_SIZE);
	for (i=0;i<ARRAY_SIZE;i++){
		posterior[i] += decayNoiseArray[(i - (int)maximumIndex + ((3*ARRAY_SIZE)/2)) % ARRAY_SIZE] * currentMaximum * decayNoiseAmount;
		//posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}
	
}

void DynamicBayesianArray::resetMaximumPosterior(){
	int i;
	float max = 0;
	for (i=0;i < ARRAY_SIZE;i++){	
		if (posterior[i]>max){
			maximumIndex = i;
			max = posterior[i];
		}
	}
}

void DynamicBayesianArray::translateDistribution(int translationIndex){
	int tmpIndex;
	//copy array
	int i;
	for (i=0;i < ARRAY_SIZE;i++){
		tempPosteriorArray[i] =	posterior[i] ;
	}
	//translate values
	for (i=0;i < ARRAY_SIZE;i++){
		tmpIndex = (i + translationIndex + ARRAY_SIZE)%ARRAY_SIZE;
		posterior[tmpIndex] = tempPosteriorArray[i]; 
	}
	//now delete tmp array
}


void DynamicBayesianArray::drawFloatArray(float* arrayToDraw, const int& minIndex, const int& maxIndex){
	
	if (minIndex >= 0){
	
		double stepSize = ofGetWidth() / (double)(maxIndex - minIndex);
		double screenHeight = ofGetHeight();
		double maxVal = getMaximum(&arrayToDraw[0], maxIndex);

		for (int i = minIndex+1;i < maxIndex;i++){
		
			ofLine (stepSize*(i-1), screenHeight * (1 - arrayToDraw[i-1] / maxVal), stepSize*i, screenHeight * (1 - arrayToDraw[i] / maxVal) );
		}
	
	}
		
		
}

 

/*
void DynamicBayesianArray::drawDoubleArray(double[]& arrayToDraw, const int& minIndex, const int& maxIndex){

	if (minIndex >= 0 && maxIndex <= arrayToDraw.size(0))

}
*/
/*
 void DynamicBayesianArray::setGaussianLikelihoodForBeats(float mean, float StdDev){
 //this has eighth and sixteenth positions included
 
 if (mean >= 0 && mean <= ARRAY_SIZE){
 int i;	float eighthDifference;
 int eighthPosition = ((int)mean + ARRAY_SIZE/2)%ARRAY_SIZE;
 int earlySixteenthPosition = ((int)mean + (3*ARRAY_SIZE/4))%ARRAY_SIZE;;
 int lateSixteenthPosition = ((int)mean + (ARRAY_SIZE/4))%ARRAY_SIZE;;
 
 float mainDifference, sixteenthDifference;
 float gaussianProportion = 1 - likelihoodNoise;
 float mainProportion = (1 - eighthNoteProportion - earlySixteenthNoteProportion - lateSixteenthNoteProportion);
 
 for (i=0;i < ARRAY_SIZE;i++){
 
 mainDifference = min( fabs(i-mean) , (double)(i + ARRAY_SIZE - mean));
 likelihood[i] = gaussianProportion * mainProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
 
 eighthDifference = min( abs(i - eighthPosition) , i + ARRAY_SIZE - eighthPosition);
 eighthDifference = min(eighthDifference , (float)(ARRAY_SIZE + eighthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * eighthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(eighthDifference)*(eighthDifference)/(2*StdDev*StdDev)) ;
 
 sixteenthDifference = min( abs(i - earlySixteenthPosition) , i + ARRAY_SIZE - earlySixteenthPosition);
 sixteenthDifference = min(sixteenthDifference , (float)(ARRAY_SIZE + earlySixteenthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * earlySixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
 
 sixteenthDifference = min( abs(i - lateSixteenthPosition) , i + ARRAY_SIZE - lateSixteenthPosition);
 sixteenthDifference = min(sixteenthDifference , (float)(ARRAY_SIZE + lateSixteenthPosition - i ));
 //for e.g. +0.43, or -0.47 we require the gaussian around the half note too
 likelihood[i] += gaussianProportion * lateSixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
 
 
 
 likelihood[i] += (likelihoodNoise / ARRAY_SIZE);
 //likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
 //(double) (likelihoodNoise / ARRAY_SIZE) );
 }
 //	renormaliseArray(&likelihood[0], ARRAY_SIZE);
 }//end if mean within limits
 }
 */