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

Merge branch 'release/1.0.5'

--- a/src/BTrack.cpp	Sat Jun 18 10:50:06 2016 +0100
+++ b/src/BTrack.cpp	Fri Aug 18 20:07:33 2023 +0200
@@ -21,6 +21,7 @@
 
 #include <cmath>
 #include <algorithm>
+#include <numeric>
 #include "BTrack.h"
 #include "samplerate.h"
 #include <iostream>
@@ -29,21 +30,21 @@
 BTrack::BTrack()
  :  odf (512, 1024, ComplexSpectralDifferenceHWR, HanningWindow)
 {
-    initialise (512, 1024);
+    initialise (512);
 }
 
 //=======================================================================
-BTrack::BTrack (int hopSize_)
- :  odf(hopSize_, 2*hopSize_, ComplexSpectralDifferenceHWR, HanningWindow)
-{	
-    initialise (hopSize_, 2*hopSize_);
+BTrack::BTrack (int hop)
+ :  odf (hop, 2 * hop, ComplexSpectralDifferenceHWR, HanningWindow)
+{
+    initialise (hop);
 }
 
 //=======================================================================
-BTrack::BTrack (int hopSize_, int frameSize_)
- : odf (hopSize_, frameSize_, ComplexSpectralDifferenceHWR, HanningWindow)
+BTrack::BTrack (int hop, int frame)
+ : odf (hop, frame, ComplexSpectralDifferenceHWR, HanningWindow)
 {
-    initialise (hopSize_, frameSize_);
+    initialise (hop);
 }
 
 //=======================================================================
@@ -66,82 +67,65 @@
 }
 
 //=======================================================================
-double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int fs)
+double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int samplingFrequency)
 {
-    double hop = (double) hopSize;
-    double samplingFrequency = (double) fs;
-    double frameNum = (double) frameNumber;
-    
-    return ((hop / samplingFrequency) * frameNum);
+    return ((static_cast<double> (hopSize) / static_cast<double> (samplingFrequency)) * static_cast<double> (frameNumber));
 }
 
 //=======================================================================
-double BTrack::getBeatTimeInSeconds (int frameNumber, int hopSize, int fs)
+void BTrack::initialise (int hop)
 {
-    long frameNum = (long) frameNumber;
+    // 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);
     
-    return getBeatTimeInSeconds (frameNum, hopSize, fs);
-}
-
-
-
-//=======================================================================
-void BTrack::initialise (int hopSize_, int frameSize_)
-{
-    double rayparam = 43;
-	double pi = 3.14159265;
-	
+    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)));
-	}
+        weightingVector[n] = ((double) n / pow (rayleighParameter, 2)) * exp((-1 * pow((double) - n, 2)) / (2 * pow (rayleighParameter, 2)));
 	
-	// initialise prev_delta
-	for (int i = 0; i < 41; i++)
-	{
-		prevDelta[i] = 1;
-	}
-	
+    // 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 i = 0; i < 41; i++)
 	{
-		for (int j = 0;j < 41;j++)
+		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)) );
+			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;
@@ -163,12 +147,11 @@
 }
 
 //=======================================================================
-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);
@@ -190,13 +173,13 @@
 }
 
 //=======================================================================
-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);
 }
 
 //=======================================================================
@@ -220,7 +203,7 @@
 //=======================================================================
 double BTrack::getLatestCumulativeScoreValue()
 {
-    return latestCumulativeScoreValue;
+    return cumulativeScore[cumulativeScore.size() - 1];
 }
 
 //=======================================================================
@@ -244,8 +227,8 @@
     // 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
@@ -254,14 +237,12 @@
 	// 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
 		
@@ -273,44 +254,35 @@
 
 //=======================================================================
 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);
+	int tempoIndex = (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;
-	
+    // 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 newBeatPeriod = (int) round (60 / ((((double) hopSize) / 44100) * tempo));
 	
-	int bcounter = 1;
-	// initialise df_buffer to zeros
-	for (int i = (onsetDFBufferSize-1);i >= 0;i--)
+	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 (bcounter == 1)
+		if (k == 1)
 		{
 			cumulativeScore[i] = 150;
 			onsetDF[i] = 150;
@@ -321,21 +293,21 @@
 			onsetDF[i] = 10;
 		}
 		
-		bcounter++;
+		k++;
 		
-		if (bcounter > new_bperiod)
+		if (k > newBeatPeriod)
 		{
-			bcounter = 1;
+			k = 1;
 		}
 	}
 	
 	/////////// INDICATE THAT THIS IS A BEAT //////////////////
 	
 	// beat is now
-	beatCounter = 0;
+	timeToNextBeat = 0;
 	
-	// offbeat is half of new beat period away
-	m0 = (int) round(((double) new_bperiod)/2);
+	// next prediction is on the offbeat, so half of new beat period away
+	timeToNextPrediction = (int) round (((double) newBeatPeriod) / 2);
 }
 
 //=======================================================================
@@ -343,26 +315,22 @@
 {	
 	// 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;
@@ -379,43 +347,35 @@
 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 ;
+    double ratio = 512.0 / ((double) onsetDFBufferSize);
+    int bufferLength = onsetDFBufferSize;
+    int outputLength = 512;
     
-    //output_len = (int) floor (((double) BUFFER_LEN) * src_ratio) ;
-    output_len = 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);
@@ -424,149 +384,127 @@
 	calculateOutputOfCombFilterBank();
 	
 	// adaptive threshold on rcf
-	adaptiveThreshold (combFilterBankOutput,128);
+	adaptiveThreshold (combFilterBankOutput);
 
-	
-	int t_index;
-	int t_index2;
 	// calculate tempo observation vector from beat period observation vector
-	for (int i = 0;i < 41;i++)
+	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];
+		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];
 	}
 	
-	
-	double maxval;
-	double maxind;
-	double curval;
-	
 	// 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++)
+	for (int j = 0; j < 41; j++)
 	{
-		maxval = -1;
-		for (int i = 0;i < 41;i++)
+		double maxValue = -1;
+        
+		for (int i = 0; i < 41; i++)
 		{
-			curval = prevDelta[i] * tempoTransitionMatrix[i][j];
+			double currentValue = prevDelta[i] * tempoTransitionMatrix[i][j];
 			
-			if (curval > maxval)
-			{
-				maxval = curval;
-			}
+			if (currentValue > maxValue)
+                maxValue = currentValue;
 		}
 		
-		delta[j] = maxval * tempoObservationVector[j];
+		delta[j] = maxValue * tempoObservationVector[j];
 	}
 	
-
-	normaliseArray(delta,41);
+	normaliseVector (delta);
 	
-	maxind = -1;
-	maxval = -1;
+	double maxIndex = -1;
+	double maxValue = -1;
 	
-	for (int j=0;j < 41;j++)
+	for (int j = 0; j < 41; j++)
 	{
-		if (delta[j] > maxval)
+		if (delta[j] > maxValue)
 		{
-			maxval = delta[j];
-			maxind = j;
+			maxValue = delta[j];
+			maxIndex = j;
 		}
 		
 		prevDelta[j] = delta[j];
 	}
 	
-	beatPeriod = round ((60.0*44100.0)/(((2*maxind)+80)*((double) hopSize)));
+	beatPeriod = round ((60.0 * 44100.0) / (((2 * maxIndex) + 80) * ((double) hopSize)));
 	
 	if (beatPeriod > 0)
-	{
-		estimatedTempo = 60.0/((((double) hopSize) / 44100.0) * beatPeriod);
-	}
+        estimatedTempo = 60.0 / ((((double) hopSize) / 44100.0) * beatPeriod);
 }
 
 //=======================================================================
-void BTrack::adaptiveThreshold (double* x, int N)
+void BTrack::adaptiveThreshold (std::vector<double>& x)
 {
-	int i = 0;
-	int k,t = 0;
-	double x_thresh[N];
+    int N = static_cast<int> (x.size());
+	double threshold[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
+	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 (i = 0;i <= t;i++)
+	for (int i = 0; i <= t; i++)
 	{	
-		k = std::min ((i+p_pre),N);
-		x_thresh[i] = calculateMeanOfArray (x,1,k);
+		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 (i = t+1;i < N-p_post;i++)
+	for (int i = t + 1; i < N - p_post; i++)
 	{
-		x_thresh[i] = calculateMeanOfArray (x,i-p_pre,i+p_post);
+		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 (i = N-p_post;i < N;i++)
+	for (int i = N - p_post; i < N; i++)
 	{
-		k = std::max ((i-p_post),1);
-		x_thresh[i] = calculateMeanOfArray (x,k,N);
+		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 (i = 0; i < N; i++)
+	for (int i = 0; i < N; i++)
 	{
-		x[i] = x[i] - x_thresh[i];
+		x[i] = x[i] - threshold[i];
+        
 		if (x[i] < 0)
-		{
-			x[i] = 0;
-		}
+            x[i] = 0;
 	}
 }
 
 //=======================================================================
 void BTrack::calculateOutputOfCombFilterBank()
 {
-	int numelem;
-	
-	for (int i = 0;i < 128;i++)
-	{
-		combFilterBankOutput[i] = 0;
-	}
-	
-	numelem = 4;
+    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 <= numelem; a++) // number of comb elements
+		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
+			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
+				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)
         {
@@ -584,9 +522,9 @@
     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;
     }
     
@@ -597,7 +535,7 @@
     
 #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)
         {
@@ -615,7 +553,7 @@
     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;
@@ -632,7 +570,7 @@
     {
 #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
@@ -652,167 +590,155 @@
 }
 
 //=======================================================================
-double BTrack::calculateMeanOfArray (double* array, int startIndex, int endIndex)
+double BTrack::calculateMeanOfVector (std::vector<double>& vector, int startIndex, int endIndex)
 {
-	int i;
-	double sum = 0;
+    int length = endIndex - startIndex;
+    double sum = std::accumulate (vector.begin() + startIndex, vector.begin() + endIndex, 0.0);
 
-    int length = endIndex - startIndex;
+    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);
 	
-	// find sum
-	for (i = startIndex; i < endIndex; i++)
-	{
-		sum = sum + array[i];
-	}
-	
-    if (length > 0)
+	if (sum > 0)
     {
-        return sum / length;	// average and return
-    }
-    else
-    {
-        return 0;
+        for (int i = 0; i < vector.size(); i++)
+            vector[i] = vector[i] / sum;
     }
 }
 
 //=======================================================================
-void BTrack::normaliseArray (double* array, int N)
+void BTrack::updateCumulativeScore (double onsetDetectionFunctionSample)
 {
-	double sum = 0;
+	int windowStart = onsetDFBufferSize - round (2. * beatPeriod);
+	int windowEnd = onsetDFBufferSize - round (beatPeriod / 2.);
+	int windowSize = windowEnd - windowStart + 1;
 	
-	for (int i = 0; i < N; i++)
-	{
-		if (array[i] > 0)
-		{
-			sum = sum + array[i];
-		}
-	}
+    // create log gaussian transition window
+	double logGaussianTransitionWeighting[windowSize];
+    createLogGaussianTransitionWeighting (logGaussianTransitionWeighting, windowSize, beatPeriod);
 	
-	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);
+    // calculate the new cumulative score value
+    double cumulativeScoreValue = calculateNewCumulativeScoreValue (cumulativeScore, logGaussianTransitionWeighting, windowStart, windowEnd, onsetDetectionFunctionSample, alpha);
     
-    cumulativeScore.addSampleToEnd (latestCumulativeScoreValue);
+    // 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];
+	int beatExpectationWindowSize = static_cast<int> (beatPeriod);
+	double futureCumulativeScore[onsetDFBufferSize + beatExpectationWindowSize];
+	double beatExpectationWindow[beatExpectationWindowSize];
     
-	// copy cumscore to first part of fcumscore
-	for (int i = 0;i < onsetDFBufferSize;i++)
+	// 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 < beatExpectationWindowSize; i++)
 	{
-		futureCumulativeScore[i] = cumulativeScore[i];
-	}
-	
-	// create future window
-	double v = 1;
-	for (int i = 0; i < windowSize; i++)
-	{
-		w2[i] = exp((-1*pow((v - (beatPeriod/2)),2))   /  (2*pow((beatPeriod/2) ,2)));
+		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];
+	// 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);
 
-	for (int i = 0;i < pastwinsize;i++)
+	// 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++)
 	{
-		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;
+        // 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 beat
-	max = 0;
-	n = 0;
+	// 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 + windowSize); i++)
+	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + beatExpectationWindowSize); i++)
 	{
-		wcumscore = futureCumulativeScore[i]*w2[n];
+		double weightedCumulativeScore = futureCumulativeScore[i] * beatExpectationWindow[n];
 		
-		if (wcumscore > max)
+		if (weightedCumulativeScore > maxValue)
 		{
-			max = wcumscore;
-			beatCounter = n;
+			maxValue = weightedCumulativeScore;
+			timeToNextBeat = n;
 		}	
 		
 		n++;
 	}
 		
-	// set next prediction time
-	m0 = beatCounter + round (beatPeriod / 2);
-}
\ No newline at end of file
+	// 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;
+}