comparison src/BTrack.cpp @ 117:ca2d83d29814 tip master

Merge branch 'release/1.0.5'

  */
 //=======================================================================
 
 #include <cmath>
 #include <algorithm>
+#include <numeric>
 #include "BTrack.h"
 #include "samplerate.h"
 #include <iostream>
 
 //=======================================================================
 BTrack::BTrack()
  :  odf (512, 1024, ComplexSpectralDifferenceHWR, HanningWindow)
 {
-    initialise (512, 1024);
-}
-
-//=======================================================================
-BTrack::BTrack (int hopSize_)
- :  odf(hopSize_, 2*hopSize_, ComplexSpectralDifferenceHWR, HanningWindow)
-{	
-    initialise (hopSize_, 2*hopSize_);
-}
-
-//=======================================================================
-BTrack::BTrack (int hopSize_, int frameSize_)
- : odf (hopSize_, frameSize_, ComplexSpectralDifferenceHWR, HanningWindow)
-{
-    initialise (hopSize_, frameSize_);
+    initialise (512);
+}
+
+//=======================================================================
+BTrack::BTrack (int hop)
+ :  odf (hop, 2 * hop, ComplexSpectralDifferenceHWR, HanningWindow)
+{
+    initialise (hop);
+}
+
+//=======================================================================
+BTrack::BTrack (int hop, int frame)
+ : odf (hop, frame, ComplexSpectralDifferenceHWR, HanningWindow)
+{
+    initialise (hop);
 }
 
 //=======================================================================
 BTrack::~BTrack()
 {

     delete [] fftOut;
 #endif
 }
 
 //=======================================================================
-double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int fs)
-{
-    double hop = (double) hopSize;
-    double samplingFrequency = (double) fs;
-    double frameNum = (double) frameNumber;
-    
-    return ((hop / samplingFrequency) * frameNum);
-}
-
-//=======================================================================
-double BTrack::getBeatTimeInSeconds (int frameNumber, int hopSize, int fs)
-{
-    long frameNum = (long) frameNumber;
-    
-    return getBeatTimeInSeconds (frameNum, hopSize, fs);
-}
-
-
-
-//=======================================================================
-void BTrack::initialise (int hopSize_, int frameSize_)
-{
-    double rayparam = 43;
-	double pi = 3.14159265;
-	
+double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int samplingFrequency)
+{
+    return ((static_cast<double> (hopSize) / static_cast<double> (samplingFrequency)) * static_cast<double> (frameNumber));
+}
+
+//=======================================================================
+void BTrack::initialise (int hop)
+{
+    // set vector sizes
+    resampledOnsetDF.resize (512);
+    acf.resize (512);
+    weightingVector.resize (128);
+    combFilterBankOutput.resize (128);
+    tempoObservationVector.resize (41);
+    delta.resize (41);
+    prevDelta.resize (41);
+    prevDeltaFixed.resize (41);
+    
+    double rayleighParameter = 43;
 	
 	// initialise parameters
 	tightness = 5;
 	alpha = 0.9;
-	tempo = 120;
 	estimatedTempo = 120.0;
-	tempoToLagFactor = 60.*44100./512.;
-	
-	m0 = 10;
-	beatCounter = -1;
+	
+	timeToNextPrediction = 10;
+	timeToNextBeat = -1;
 	
 	beatDueInFrame = false;
 	
 
 	// create rayleigh weighting vector
 	for (int n = 0; n < 128; n++)
-	{
-		weightingVector[n] = ((double) n / pow(rayparam,2)) * exp((-1*pow((double)-n,2)) / (2*pow(rayparam,2)));
-	}
-	
-	// initialise prev_delta
-	for (int i = 0; i < 41; i++)
-	{
-		prevDelta[i] = 1;
-	}
-	
+        weightingVector[n] = ((double) n / pow (rayleighParameter, 2)) * exp((-1 * pow((double) - n, 2)) / (2 * pow (rayleighParameter, 2)));
+	
+    // initialise prevDelta
+    std::fill (prevDelta.begin(), prevDelta.end(), 1);
+    
 	double t_mu = 41/2;
 	double m_sig;
 	double x;
 	// create tempo transition matrix
 	m_sig = 41/8;
-	for (int i = 0;i < 41;i++)
-	{
-		for (int j = 0;j < 41;j++)
-		{
-			x = j+1;
-			t_mu = i+1;
-			tempoTransitionMatrix[i][j] = (1 / (m_sig * sqrt(2*pi))) * exp( (-1*pow((x-t_mu),2)) / (2*pow(m_sig,2)) );
+    
+	for (int i = 0; i < 41; i++)
+	{
+		for (int j = 0; j < 41; j++)
+		{
+			x = j + 1;
+			t_mu = i + 1;
+			tempoTransitionMatrix[i][j] = (1 / (m_sig * sqrt (2 * M_PI))) * exp((-1 * pow ((x - t_mu), 2)) / (2 * pow (m_sig, 2)) );
 		}
 	}
 	
 	// tempo is not fixed
 	tempoFixed = false;
     
-    // initialise latest cumulative score value
-    // in case it is requested before any processing takes place
-    latestCumulativeScoreValue = 0;
-    
     // initialise algorithm given the hopsize
-    setHopSize(hopSize_);
-    
+    setHopSize (hop);
     
     // Set up FFT for calculating the auto-correlation function
     FFTLengthForACFCalculation = 1024;
     
 #ifdef USE_FFTW

     cfgBackwards = kiss_fft_alloc (FFTLengthForACFCalculation, 1, 0, 0);
 #endif
 }
 
 //=======================================================================
-void BTrack::setHopSize (int hopSize_)
+void BTrack::setHopSize (int hop)
 {	
-	hopSize = hopSize_;
-	onsetDFBufferSize = (512*512)/hopSize;		// calculate df buffer size
-	
-	beatPeriod = round(60/((((double) hopSize)/44100)*tempo));
+	hopSize = hop;
+	onsetDFBufferSize = (512 * 512) / hopSize;		// calculate df buffer size
+	beatPeriod = round (60 / ((((double) hopSize) / 44100) * 120.));
 
     // set size of onset detection function buffer
     onsetDF.resize (onsetDFBufferSize);
     
     // set size of cumulative score buffer

 		}
 	}
 }
 
 //=======================================================================
-void BTrack::updateHopAndFrameSize (int hopSize_, int frameSize_)
+void BTrack::updateHopAndFrameSize (int hop, int frame)
 {
     // update the onset detection function object
-    odf.initialise (hopSize_, frameSize_);
+    odf.initialise (hop, frame);
     
     // update the hop size being used by the beat tracker
-    setHopSize (hopSize_);
+    setHopSize (hop);
 }
 
 //=======================================================================
 bool BTrack::beatDueInCurrentFrame()
 {

 }
 
 //=======================================================================
 double BTrack::getLatestCumulativeScoreValue()
 {
-    return latestCumulativeScoreValue;
+    return cumulativeScore[cumulativeScore.size() - 1];
 }
 
 //=======================================================================
 void BTrack::processAudioFrame (double* frame)
 {

     
     // add a tiny constant to the sample to stop it from ever going
     // to zero. this is to avoid problems further down the line
     newSample = newSample + 0.0001;
     
-	m0--;
-	beatCounter--;
+	timeToNextPrediction--;
+	timeToNextBeat--;
 	beatDueInFrame = false;
 		
 	// add new sample at the end
     onsetDF.addSampleToEnd (newSample);
 	
 	// update cumulative score
 	updateCumulativeScore (newSample);
 	
-	// if we are halfway between beats
-	if (m0 == 0)
-	{
-		predictBeat();
-	}
+	// if we are halfway between beats, predict a beat
+	if (timeToNextPrediction == 0)
+        predictBeat();
 	
 	// if we are at a beat
-	if (beatCounter == 0)
+	if (timeToNextBeat == 0)
 	{
 		beatDueInFrame = true;	// indicate a beat should be output
 		
 		// recalculate the tempo
 		resampleOnsetDetectionFunction();

 	}
 }
 
 //=======================================================================
 void BTrack::setTempo (double tempo)
-{	 
-	
+{
 	/////////// TEMPO INDICATION RESET //////////////////
 	
 	// firstly make sure tempo is between 80 and 160 bpm..
 	while (tempo > 160)
-	{
-		tempo = tempo/2;
-	}
+        tempo = tempo / 2;
 	
 	while (tempo < 80)
-	{
-		tempo = tempo * 2;
-	}
+        tempo = tempo * 2;
 		
 	// convert tempo from bpm value to integer index of tempo probability 
-	int tempo_index = (int) round((tempo - 80)/2);
-	
-	// now set previous tempo observations to zero
-	for (int i=0;i < 41;i++)
-	{
-		prevDelta[i] = 0;
-	}
-	
-	// set desired tempo index to 1
-	prevDelta[tempo_index] = 1;
-	
+	int tempoIndex = (int) round ((tempo - 80.) / 2);
+	
+    // now set previous tempo observations to zero and set desired tempo index to 1
+    std::fill (prevDelta.begin(), prevDelta.end(), 0);
+	prevDelta[tempoIndex] = 1;
 	
 	/////////// CUMULATIVE SCORE ARTIFICAL TEMPO UPDATE //////////////////
 	
 	// calculate new beat period
-	int new_bperiod = (int) round(60/((((double) hopSize)/44100)*tempo));
-	
-	int bcounter = 1;
-	// initialise df_buffer to zeros
-	for (int i = (onsetDFBufferSize-1);i >= 0;i--)
-	{
-		if (bcounter == 1)
+	int newBeatPeriod = (int) round (60 / ((((double) hopSize) / 44100) * tempo));
+	
+	int k = 1;
+    
+    // initialise onset detection function with delta functions spaced
+    // at the new beat period
+	for (int i = onsetDFBufferSize - 1; i >= 0; i--)
+	{
+		if (k == 1)
 		{
 			cumulativeScore[i] = 150;
 			onsetDF[i] = 150;
 		}
 		else
 		{
 			cumulativeScore[i] = 10;
 			onsetDF[i] = 10;
 		}
 		
-		bcounter++;
-		
-		if (bcounter > new_bperiod)
-		{
-			bcounter = 1;
+		k++;
+		
+		if (k > newBeatPeriod)
+		{
+			k = 1;
 		}
 	}
 	
 	/////////// INDICATE THAT THIS IS A BEAT //////////////////
 	
 	// beat is now
-	beatCounter = 0;
-	
-	// offbeat is half of new beat period away
-	m0 = (int) round(((double) new_bperiod)/2);
+	timeToNextBeat = 0;
+	
+	// next prediction is on the offbeat, so half of new beat period away
+	timeToNextPrediction = (int) round (((double) newBeatPeriod) / 2);
 }
 
 //=======================================================================
 void BTrack::fixTempo (double tempo)
 {	
 	// firstly make sure tempo is between 80 and 160 bpm..
 	while (tempo > 160)
-	{
-		tempo = tempo/2;
-	}
+        tempo = tempo / 2;
 	
 	while (tempo < 80)
-	{
-		tempo = tempo * 2;
-	}
+        tempo = tempo * 2;
 	
 	// convert tempo from bpm value to integer index of tempo probability 
-	int tempo_index = (int) round((tempo - 80)/2);
+	int tempoIndex = (int) round((tempo - 80) / 2);
 	
 	// now set previous fixed previous tempo observation values to zero
-	for (int i=0;i < 41;i++)
+	for (int i = 0; i < 41; i++)
 	{
 		prevDeltaFixed[i] = 0;
 	}
 	
 	// set desired tempo index to 1
-	prevDeltaFixed[tempo_index] = 1;
+	prevDeltaFixed[tempoIndex] = 1;
 		
 	// set the tempo fix flag
 	tempoFixed = true;
 }
 

 
 //=======================================================================
 void BTrack::resampleOnsetDetectionFunction()
 {
 	float output[512];
-    
     float input[onsetDFBufferSize];
     
-    for (int i = 0;i < onsetDFBufferSize;i++)
-    {
+    for (int i = 0; i < onsetDFBufferSize; i++)
         input[i] = (float) onsetDF[i];
-    }
-        
-    double src_ratio = 512.0/((double) onsetDFBufferSize);
-    int BUFFER_LEN = onsetDFBufferSize;
-    int output_len;
-    SRC_DATA	src_data ;
-    
-    //output_len = (int) floor (((double) BUFFER_LEN) * src_ratio) ;
-    output_len = 512;
-    
+        
+    double ratio = 512.0 / ((double) onsetDFBufferSize);
+    int bufferLength = onsetDFBufferSize;
+    int outputLength = 512;
+    
+    SRC_DATA src_data;
     src_data.data_in = input;
-    src_data.input_frames = BUFFER_LEN;
-    
-    src_data.src_ratio = src_ratio;
-    
+    src_data.input_frames = bufferLength;
+    src_data.src_ratio = ratio;
     src_data.data_out = output;
-    src_data.output_frames = output_len;
+    src_data.output_frames = outputLength;
     
     src_simple (&src_data, SRC_SINC_BEST_QUALITY, 1);
             
-    for (int i = 0;i < output_len;i++)
-    {
+    for (int i = 0; i < outputLength; i++)
         resampledOnsetDF[i] = (double) src_data.data_out[i];
-    }
 }
 
 //=======================================================================
 void BTrack::calculateTempo()
 {
+    double tempoToLagFactor = 60. * 44100. / 512.;
+    
 	// adaptive threshold on input
-	adaptiveThreshold (resampledOnsetDF,512);
+	adaptiveThreshold (resampledOnsetDF);
 		
 	// calculate auto-correlation function of detection function
 	calculateBalancedACF (resampledOnsetDF);
 	
 	// calculate output of comb filterbank
 	calculateOutputOfCombFilterBank();
 	
 	// adaptive threshold on rcf
-	adaptiveThreshold (combFilterBankOutput,128);
-
-	
-	int t_index;
-	int t_index2;
+	adaptiveThreshold (combFilterBankOutput);
+
 	// calculate tempo observation vector from beat period observation vector
-	for (int i = 0;i < 41;i++)
-	{
-		t_index = (int) round (tempoToLagFactor / ((double) ((2*i)+80)));
-		t_index2 = (int) round (tempoToLagFactor / ((double) ((4*i)+160)));
-
-		
-		tempoObservationVector[i] = combFilterBankOutput[t_index-1] + combFilterBankOutput[t_index2-1];
-	}
-	
-	
-	double maxval;
-	double maxind;
-	double curval;
+	for (int i = 0; i < 41; i++)
+	{
+		int tempoIndex1 = (int) round (tempoToLagFactor / ((double) ((2 * i) + 80)));
+		int tempoIndex2 = (int) round (tempoToLagFactor / ((double) ((4 * i) + 160)));
+		tempoObservationVector[i] = combFilterBankOutput[tempoIndex1 - 1] + combFilterBankOutput[tempoIndex2 - 1];
+	}
 	
 	// if tempo is fixed then always use a fixed set of tempi as the previous observation probability function
 	if (tempoFixed)
 	{
-		for (int k = 0;k < 41;k++)
-		{
-			prevDelta[k] = prevDeltaFixed[k];
+		for (int k = 0; k < 41; k++)
+            prevDelta[k] = prevDeltaFixed[k];
+	}
+		
+	for (int j = 0; j < 41; j++)
+	{
+		double maxValue = -1;
+        
+		for (int i = 0; i < 41; i++)
+		{
+			double currentValue = prevDelta[i] * tempoTransitionMatrix[i][j];
+			
+			if (currentValue > maxValue)
+                maxValue = currentValue;
 		}
-	}
-		
-	for (int j=0;j < 41;j++)
-	{
-		maxval = -1;
-		for (int i = 0;i < 41;i++)
-		{
-			curval = prevDelta[i] * tempoTransitionMatrix[i][j];
-			
-			if (curval > maxval)
+		
+		delta[j] = maxValue * tempoObservationVector[j];
+	}
+	
+	normaliseVector (delta);
+	
+	double maxIndex = -1;
+	double maxValue = -1;
+	
+	for (int j = 0; j < 41; j++)
+	{
+		if (delta[j] > maxValue)
+		{
+			maxValue = delta[j];
+			maxIndex = j;
+		}
+		
+		prevDelta[j] = delta[j];
+	}
+	
+	beatPeriod = round ((60.0 * 44100.0) / (((2 * maxIndex) + 80) * ((double) hopSize)));
+	
+	if (beatPeriod > 0)
+        estimatedTempo = 60.0 / ((((double) hopSize) / 44100.0) * beatPeriod);
+}
+
+//=======================================================================
+void BTrack::adaptiveThreshold (std::vector<double>& x)
+{
+    int N = static_cast<int> (x.size());
+	double threshold[N];
+	
+	int p_post = 7;
+	int p_pre = 8;
+	
+	int t = std::min (N, p_post);	// what is smaller, p_post or df size. This is to avoid accessing outside of arrays
+	
+	// find threshold for first 't' samples, where a full average cannot be computed yet 
+	for (int i = 0; i <= t; i++)
+	{	
+		int k = std::min ((i + p_pre), N);
+		threshold[i] = calculateMeanOfVector (x, 1, k);
+	}
+    
+	// find threshold for bulk of samples across a moving average from [i-p_pre,i+p_post]
+	for (int i = t + 1; i < N - p_post; i++)
+	{
+		threshold[i] = calculateMeanOfVector (x, i - p_pre, i + p_post);
+	}
+    
+	// for last few samples calculate threshold, again, not enough samples to do as above
+	for (int i = N - p_post; i < N; i++)
+	{
+		int k = std::max ((i - p_post), 1);
+		threshold[i] = calculateMeanOfVector (x, k, N);
+	}
+	
+	// subtract the threshold from the detection function and check that it is not less than 0
+	for (int i = 0; i < N; i++)
+	{
+		x[i] = x[i] - threshold[i];
+        
+		if (x[i] < 0)
+            x[i] = 0;
+	}
+}
+
+//=======================================================================
+void BTrack::calculateOutputOfCombFilterBank()
+{
+    std::fill (combFilterBankOutput.begin(), combFilterBankOutput.end(), 0.0);
+	int numCombElements = 4;
+	
+	for (int i = 2; i <= 127; i++) // max beat period
+	{
+		for (int a = 1; a <= numCombElements; a++) // number of comb elements
+		{
+			for (int b = 1 - a; b <= a - 1; b++) // general state using normalisation of comb elements
 			{
-				maxval = curval;
+				combFilterBankOutput[i - 1] += (acf[(a * i + b) - 1] * weightingVector[i - 1]) / (2 * a - 1);	// calculate value for comb filter row
 			}
 		}
-		
-		delta[j] = maxval * tempoObservationVector[j];
-	}
-	
-
-	normaliseArray(delta,41);
-	
-	maxind = -1;
-	maxval = -1;
-	
-	for (int j=0;j < 41;j++)
-	{
-		if (delta[j] > maxval)
-		{
-			maxval = delta[j];
-			maxind = j;
-		}
-		
-		prevDelta[j] = delta[j];
-	}
-	
-	beatPeriod = round ((60.0*44100.0)/(((2*maxind)+80)*((double) hopSize)));
-	
-	if (beatPeriod > 0)
-	{
-		estimatedTempo = 60.0/((((double) hopSize) / 44100.0) * beatPeriod);
-	}
-}
-
-//=======================================================================
-void BTrack::adaptiveThreshold (double* x, int N)
-{
-	int i = 0;
-	int k,t = 0;
-	double x_thresh[N];
-	
-	int p_post = 7;
-	int p_pre = 8;
-	
-	t = std::min(N,p_post);	// what is smaller, p_post of df size. This is to avoid accessing outside of arrays
-	
-	// find threshold for first 't' samples, where a full average cannot be computed yet 
-	for (i = 0;i <= t;i++)
-	{	
-		k = std::min ((i+p_pre),N);
-		x_thresh[i] = calculateMeanOfArray (x,1,k);
-	}
-	// find threshold for bulk of samples across a moving average from [i-p_pre,i+p_post]
-	for (i = t+1;i < N-p_post;i++)
-	{
-		x_thresh[i] = calculateMeanOfArray (x,i-p_pre,i+p_post);
-	}
-	// for last few samples calculate threshold, again, not enough samples to do as above
-	for (i = N-p_post;i < N;i++)
-	{
-		k = std::max ((i-p_post),1);
-		x_thresh[i] = calculateMeanOfArray (x,k,N);
-	}
-	
-	// subtract the threshold from the detection function and check that it is not less than 0
-	for (i = 0; i < N; i++)
-	{
-		x[i] = x[i] - x_thresh[i];
-		if (x[i] < 0)
-		{
-			x[i] = 0;
-		}
-	}
-}
-
-//=======================================================================
-void BTrack::calculateOutputOfCombFilterBank()
-{
-	int numelem;
-	
-	for (int i = 0;i < 128;i++)
-	{
-		combFilterBankOutput[i] = 0;
-	}
-	
-	numelem = 4;
-	
-	for (int i = 2; i <= 127; i++) // max beat period
-	{
-		for (int a = 1; a <= numelem; a++) // number of comb elements
-		{
-			for (int b = 1-a; b <= a-1; b++) // general state using normalisation of comb elements
-			{
-				combFilterBankOutput[i-1] = combFilterBankOutput[i-1] + (acf[(a*i+b)-1]*weightingVector[i-1])/(2*a-1);	// calculate value for comb filter row
-			}
-		}
-	}
-}
-
-//=======================================================================
-void BTrack::calculateBalancedACF (double* onsetDetectionFunction)
+	}
+}
+
+//=======================================================================
+void BTrack::calculateBalancedACF (std::vector<double>& onsetDetectionFunction)
 {
     int onsetDetectionFunctionLength = 512;
     
 #ifdef USE_FFTW
     // copy into complex array and zero pad
-    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    for (int i = 0; i < FFTLengthForACFCalculation; i++)
     {
         if (i < onsetDetectionFunctionLength)
         {
             complexIn[i][0] = onsetDetectionFunction[i];
             complexIn[i][1] = 0.0;

     
     // perform the fft
     fftw_execute (acfForwardFFT);
     
     // multiply by complex conjugate
-    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    for (int i = 0; i < FFTLengthForACFCalculation; i++)
     {
-        complexOut[i][0] = complexOut[i][0]*complexOut[i][0] + complexOut[i][1]*complexOut[i][1];
+        complexOut[i][0] = complexOut[i][0] * complexOut[i][0] + complexOut[i][1] * complexOut[i][1];
         complexOut[i][1] = 0.0;
     }
     
     // perform the ifft
     fftw_execute (acfBackwardFFT);
     
 #endif
     
 #ifdef USE_KISS_FFT
     // copy into complex array and zero pad
-    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    for (int i = 0; i < FFTLengthForACFCalculation; i++)
     {
         if (i < onsetDetectionFunctionLength)
         {
             fftIn[i].r = onsetDetectionFunction[i];
             fftIn[i].i = 0.0;

     
     // execute kiss fft
     kiss_fft (cfgForwards, fftIn, fftOut);
     
     // multiply by complex conjugate
-    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    for (int i = 0; i < FFTLengthForACFCalculation; i++)
     {
         fftOut[i].r = fftOut[i].r * fftOut[i].r + fftOut[i].i * fftOut[i].i;
         fftOut[i].i = 0.0;
     }
     

     
     for (int i = 0; i < 512; i++)
     {
 #ifdef USE_FFTW
         // calculate absolute value of result
-        double absValue = sqrt (complexIn[i][0]*complexIn[i][0] + complexIn[i][1]*complexIn[i][1]);
+        double absValue = sqrt (complexIn[i][0] * complexIn[i][0] + complexIn[i][1] * complexIn[i][1]);
 #endif
         
 #ifdef USE_KISS_FFT
         // calculate absolute value of result
         double absValue = sqrt (fftIn[i].r * fftIn[i].r + fftIn[i].i * fftIn[i].i);

         lag = lag - 1.;
     }
 }
 
 //=======================================================================
-double BTrack::calculateMeanOfArray (double* array, int startIndex, int endIndex)
-{
-	int i;
-	double sum = 0;
-
+double BTrack::calculateMeanOfVector (std::vector<double>& vector, int startIndex, int endIndex)
+{
     int length = endIndex - startIndex;
-	
-	// find sum
-	for (i = startIndex; i < endIndex; i++)
-	{
-		sum = sum + array[i];
-	}
-	
+    double sum = std::accumulate (vector.begin() + startIndex, vector.begin() + endIndex, 0.0);
+
     if (length > 0)
+        return sum / static_cast<double> (length);	// average and return
+    else
+        return 0;
+}
+
+//=======================================================================
+void BTrack::normaliseVector (std::vector<double>& vector)
+{
+    double sum = std::accumulate (vector.begin(), vector.end(), 0.0);
+	
+	if (sum > 0)
     {
-        return sum / length;	// average and return
+        for (int i = 0; i < vector.size(); i++)
+            vector[i] = vector[i] / sum;
     }
-    else
-    {
-        return 0;
-    }
-}
-
-//=======================================================================
-void BTrack::normaliseArray (double* array, int N)
-{
-	double sum = 0;
-	
-	for (int i = 0; i < N; i++)
-	{
-		if (array[i] > 0)
-		{
-			sum = sum + array[i];
-		}
-	}
-	
-	if (sum > 0)
-	{
-		for (int i = 0; i < N; i++)
-		{
-			array[i] = array[i] / sum;
-		}
-	}
-}
-
-//=======================================================================
-void BTrack::updateCumulativeScore (double odfSample)
-{	 
-	int start, end, winsize;
-	double max;
-	
-	start = onsetDFBufferSize - round (2 * beatPeriod);
-	end = onsetDFBufferSize - round (beatPeriod / 2);
-	winsize = end-start+1;
-	
-	double w1[winsize];
-	double v = -2*beatPeriod;
-	double wcumscore;
-	
-	// create window
-	for (int i = 0; i < winsize; i++)
-	{
-		w1[i] = exp((-1 * pow (tightness * log (-v / beatPeriod), 2)) / 2);
-		v = v+1;
-	}	
-	
-	// calculate new cumulative score value
-	max = 0;
-	int n = 0;
-	for (int i=start; i <= end; i++)
-	{
-			wcumscore = cumulativeScore[i]*w1[n];
-		
-			if (wcumscore > max)
-			{
-				max = wcumscore;
-			}
-		n++;
-	}
-	
-    latestCumulativeScoreValue = ((1 - alpha) * odfSample) + (alpha * max);
-    
-    cumulativeScore.addSampleToEnd (latestCumulativeScoreValue);
+}
+
+//=======================================================================
+void BTrack::updateCumulativeScore (double onsetDetectionFunctionSample)
+{
+	int windowStart = onsetDFBufferSize - round (2. * beatPeriod);
+	int windowEnd = onsetDFBufferSize - round (beatPeriod / 2.);
+	int windowSize = windowEnd - windowStart + 1;
+	
+    // create log gaussian transition window
+	double logGaussianTransitionWeighting[windowSize];
+    createLogGaussianTransitionWeighting (logGaussianTransitionWeighting, windowSize, beatPeriod);
+	
+    // calculate the new cumulative score value
+    double cumulativeScoreValue = calculateNewCumulativeScoreValue (cumulativeScore, logGaussianTransitionWeighting, windowStart, windowEnd, onsetDetectionFunctionSample, alpha);
+    
+    // add the new cumulative score value to the buffer
+    cumulativeScore.addSampleToEnd (cumulativeScoreValue);
 }
 
 //=======================================================================
 void BTrack::predictBeat()
 {	 
-	int windowSize = (int) beatPeriod;
-	double futureCumulativeScore[onsetDFBufferSize + windowSize];
-	double w2[windowSize];
-    
-	// copy cumscore to first part of fcumscore
-	for (int i = 0;i < onsetDFBufferSize;i++)
-	{
-		futureCumulativeScore[i] = cumulativeScore[i];
-	}
-	
-	// create future window
+	int beatExpectationWindowSize = static_cast<int> (beatPeriod);
+	double futureCumulativeScore[onsetDFBufferSize + beatExpectationWindowSize];
+	double beatExpectationWindow[beatExpectationWindowSize];
+    
+	// copy cumulativeScore to first part of futureCumulativeScore
+	for (int i = 0; i < onsetDFBufferSize; i++)
+        futureCumulativeScore[i] = cumulativeScore[i];
+	
+	// Create a beat expectation window for predicting future beats from the "future" of the cumulative score.
+    // We are making this beat prediction at the midpoint between beats, and so we make a Gaussian
+    // weighting centred on the most likely beat position (half a beat period into the future)
+    // This is W2 in Adam Stark's PhD thesis, equation 3.6, page 62
+    
 	double v = 1;
-	for (int i = 0; i < windowSize; i++)
-	{
-		w2[i] = exp((-1*pow((v - (beatPeriod/2)),2))   /  (2*pow((beatPeriod/2) ,2)));
+	for (int i = 0; i < beatExpectationWindowSize; i++)
+	{
+		beatExpectationWindow[i] = exp((-1 * pow ((v - (beatPeriod / 2)), 2))   /  (2 * pow (beatPeriod / 2, 2)));
 		v++;
 	}
 	
-	// create past window
-	v = -2*beatPeriod;
-	int start = onsetDFBufferSize - round(2*beatPeriod);
-	int end = onsetDFBufferSize - round(beatPeriod/2);
-	int pastwinsize = end-start+1;
-	double w1[pastwinsize];
-
-	for (int i = 0;i < pastwinsize;i++)
-	{
-		w1[i] = exp((-1*pow(tightness*log(-v/beatPeriod),2))/2);
-		v = v+1;
-	}
-
-	// calculate future cumulative score
-	double max;
-	int n;
-	double wcumscore;
-	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + windowSize); i++)
-	{
-		start = i - round (2*beatPeriod);
-		end = i - round (beatPeriod/2);
-		
-		max = 0;
-		n = 0;
-		for (int k=start;k <= end;k++)
-		{
-			wcumscore = futureCumulativeScore[k]*w1[n];
-			
-			if (wcumscore > max)
-			{
-				max = wcumscore;
-			}
-			n++;
-		}
-		
-		futureCumulativeScore[i] = max;
-	}
-	
-	// predict beat
-	max = 0;
-	n = 0;
-	
-	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + windowSize); i++)
-	{
-		wcumscore = futureCumulativeScore[i]*w2[n];
-		
-		if (wcumscore > max)
-		{
-			max = wcumscore;
-			beatCounter = n;
+	// Create window for "synthesizing" the cumulative score into the future
+    // It is a log-Gaussian transition weighting running from from 2 beat periods
+    // in the past to half a beat period in the past. It favours the time exactly
+    // one beat period in the past
+    
+	int startIndex = onsetDFBufferSize - round (2 * beatPeriod);
+	int endIndex = onsetDFBufferSize - round (beatPeriod / 2);
+	int pastWindowSize = endIndex - startIndex + 1;
+    
+	double logGaussianTransitionWeighting[pastWindowSize];
+    createLogGaussianTransitionWeighting (logGaussianTransitionWeighting, pastWindowSize, beatPeriod);
+
+	// Calculate the future cumulative score, by shifting the log Gaussian transition weighting from its
+    // start position of [-2 beat periods, - 0.5 beat periods] forwards over the size of the beat
+    // expectation window, calculating a new cumulative score where the onset detection function sample
+    // is zero. This uses the "momentum" of the function to generate itself into the future.
+	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + beatExpectationWindowSize); i++)
+	{
+        // note here that we pass 0.0 in for the onset detection function sample and 1.0 for the alpha weighting factor
+        // see equation 3.4 and page 60 - 62 of Adam Stark's PhD thesis for details
+        futureCumulativeScore[i] = calculateNewCumulativeScoreValue (futureCumulativeScore, logGaussianTransitionWeighting, startIndex, endIndex, 0.0, 1.0);
+        
+        startIndex++;
+        endIndex++;
+	}
+	
+	// Predict the next beat, finding the maximum point of the future cumulative score
+    // over the next beat, after being weighted by the beat expectation window
+    
+	double maxValue = 0;
+	int n = 0;
+	
+	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + beatExpectationWindowSize); i++)
+	{
+		double weightedCumulativeScore = futureCumulativeScore[i] * beatExpectationWindow[n];
+		
+		if (weightedCumulativeScore > maxValue)
+		{
+			maxValue = weightedCumulativeScore;
+			timeToNextBeat = n;
 		}	
 		
 		n++;
 	}
 		
-	// set next prediction time
-	m0 = beatCounter + round (beatPeriod / 2);
-}
+	// set next prediction time as on the offbeat after the next beat
+	timeToNextPrediction = timeToNextBeat + round (beatPeriod / 2);
+}
+
+//=======================================================================
+void BTrack::createLogGaussianTransitionWeighting (double* weightingArray, int numSamples, double beatPeriod)
+{
+    // (This is W1 in Adam Stark's PhD thesis, equation 3.2, page 60)
+    
+    double v = -2. * beatPeriod;
+    
+    for (int i = 0; i < numSamples; i++)
+    {
+        double a = tightness * log (-v / beatPeriod);
+        weightingArray[i] = exp ((-1. * a * a) / 2.);
+        v++;
+    }
+}
+
+//=======================================================================
+template <typename T>
+double BTrack::calculateNewCumulativeScoreValue (T cumulativeScoreArray, double* logGaussianTransitionWeighting, int startIndex, int endIndex, double onsetDetectionFunctionSample, double alphaWeightingFactor)
+{
+    // calculate new cumulative score value by weighting the cumulative score between
+    // startIndex and endIndex and finding the maximum value
+    double maxValue = 0;
+    int n = 0;
+    for (int i = startIndex; i <= endIndex; i++)
+    {
+        double weightedCumulativeScore = cumulativeScoreArray[i] * logGaussianTransitionWeighting[n];
+        
+        if (weightedCumulativeScore > maxValue)
+            maxValue = weightedCumulativeScore;
+        
+        n++;
+    }
+    
+    // now mix with the incoming onset detection function sample
+    // (equation 3.4 on page 60 of Adam Stark's PhD thesis)
+    double cumulativeScoreValue = ((1. - alphaWeightingFactor) * onsetDetectionFunctionSample) + (alphaWeightingFactor * maxValue);
+    
+    return cumulativeScoreValue;
+}