adamstark@46: //=======================================================================
adamstark@46: /** @file BTrack.cpp
adamstark@47:  *  @brief BTrack - a real-time beat tracker
adamstark@46:  *  @author Adam Stark
adamstark@46:  *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
adamstark@46:  *
adamstark@46:  * This program is free software: you can redistribute it and/or modify
adamstark@46:  * it under the terms of the GNU General Public License as published by
adamstark@46:  * the Free Software Foundation, either version 3 of the License, or
adamstark@46:  * (at your option) any later version.
adamstark@46:  *
adamstark@46:  * This program is distributed in the hope that it will be useful,
adamstark@46:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
adamstark@46:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
adamstark@46:  * GNU General Public License for more details.
adamstark@46:  *
adamstark@46:  * You should have received a copy of the GNU General Public License
adamstark@46:  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
adamstark@46:  */
adamstark@46: //=======================================================================
adamstark@46: 
adamstark@46: #include <cmath>
adamstark@52: #include <algorithm>
adamstark@97: #include <numeric>
adamstark@46: #include "BTrack.h"
adamstark@46: #include "samplerate.h"
adamstark@89: #include <iostream>
adamstark@46: 
adamstark@55: //=======================================================================
adamstark@91: BTrack::BTrack()
adamstark@91:  :  odf (512, 1024, ComplexSpectralDifferenceHWR, HanningWindow)
adamstark@55: {
adamstark@108:     initialise (512);
adamstark@55: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@108: BTrack::BTrack (int hop)
adamstark@108:  :  odf (hop, 2 * hop, ComplexSpectralDifferenceHWR, HanningWindow)
adamstark@111: {
adamstark@108:     initialise (hop);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@108: BTrack::BTrack (int hop, int frame)
adamstark@108:  : odf (hop, frame, ComplexSpectralDifferenceHWR, HanningWindow)
adamstark@55: {
adamstark@108:     initialise (hop);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@88: BTrack::~BTrack()
adamstark@88: {
adamstark@93: #ifdef USE_FFTW
adamstark@88:     // destroy fft plan
adamstark@91:     fftw_destroy_plan (acfForwardFFT);
adamstark@91:     fftw_destroy_plan (acfBackwardFFT);
adamstark@91:     fftw_free (complexIn);
adamstark@91:     fftw_free (complexOut);
adamstark@93: #endif
adamstark@93:     
adamstark@93: #ifdef USE_KISS_FFT
adamstark@93:     free (cfgForwards);
adamstark@93:     free (cfgBackwards);
adamstark@93:     delete [] fftIn;
adamstark@93:     delete [] fftOut;
adamstark@93: #endif
adamstark@88: }
adamstark@88: 
adamstark@88: //=======================================================================
adamstark@108: double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int samplingFrequency)
adamstark@55: {
adamstark@108:     return ((static_cast<double> (hopSize) / static_cast<double> (samplingFrequency)) * static_cast<double> (frameNumber));
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@108: void BTrack::initialise (int hop)
adamstark@55: {
adamstark@97:     // set vector sizes
adamstark@97:     resampledOnsetDF.resize (512);
adamstark@97:     acf.resize (512);
adamstark@97:     weightingVector.resize (128);
adamstark@97:     combFilterBankOutput.resize (128);
adamstark@97:     tempoObservationVector.resize (41);
adamstark@97:     delta.resize (41);
adamstark@97:     prevDelta.resize (41);
adamstark@97:     prevDeltaFixed.resize (41);
adamstark@97:     
adamstark@98:     double rayleighParameter = 43;
adamstark@46: 	
adamstark@46: 	// initialise parameters
adamstark@46: 	tightness = 5;
adamstark@46: 	alpha = 0.9;
adamstark@58: 	estimatedTempo = 120.0;
adamstark@46: 	
adamstark@105: 	timeToNextPrediction = 10;
adamstark@105: 	timeToNextBeat = -1;
adamstark@46: 	
adamstark@57: 	beatDueInFrame = false;
adamstark@46: 	
adamstark@58: 
adamstark@46: 	// create rayleigh weighting vector
adamstark@91: 	for (int n = 0; n < 128; n++)
adamstark@98:         weightingVector[n] = ((double) n / pow (rayleighParameter, 2)) * exp((-1 * pow((double) - n, 2)) / (2 * pow (rayleighParameter, 2)));
adamstark@46: 	
adamstark@100:     // initialise prevDelta
adamstark@97:     std::fill (prevDelta.begin(), prevDelta.end(), 1);
adamstark@97:     
adamstark@54: 	double t_mu = 41/2;
adamstark@54: 	double m_sig;
adamstark@54: 	double x;
adamstark@46: 	// create tempo transition matrix
adamstark@46: 	m_sig = 41/8;
adamstark@111:     
adamstark@111: 	for (int i = 0; i < 41; i++)
adamstark@46: 	{
adamstark@111: 		for (int j = 0; j < 41; j++)
adamstark@46: 		{
adamstark@111: 			x = j + 1;
adamstark@111: 			t_mu = i + 1;
adamstark@111: 			tempoTransitionMatrix[i][j] = (1 / (m_sig * sqrt (2 * M_PI))) * exp((-1 * pow ((x - t_mu), 2)) / (2 * pow (m_sig, 2)) );
adamstark@46: 		}
adamstark@55: 	}
adamstark@46: 	
adamstark@46: 	// tempo is not fixed
adamstark@58: 	tempoFixed = false;
adamstark@58:     
adamstark@55:     // initialise algorithm given the hopsize
adamstark@108:     setHopSize (hop);
adamstark@88:     
adamstark@88:     // Set up FFT for calculating the auto-correlation function
adamstark@88:     FFTLengthForACFCalculation = 1024;
adamstark@88:     
adamstark@93: #ifdef USE_FFTW
adamstark@91:     complexIn = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * FFTLengthForACFCalculation);		// complex array to hold fft data
adamstark@91:     complexOut = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * FFTLengthForACFCalculation);	// complex array to hold fft data
adamstark@88:     
adamstark@91:     acfForwardFFT = fftw_plan_dft_1d (FFTLengthForACFCalculation, complexIn, complexOut, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
adamstark@91:     acfBackwardFFT = fftw_plan_dft_1d (FFTLengthForACFCalculation, complexOut, complexIn, FFTW_BACKWARD, FFTW_ESTIMATE);	// FFT plan initialisation
adamstark@93: #endif
adamstark@93:     
adamstark@93: #ifdef USE_KISS_FFT
adamstark@93:     fftIn = new kiss_fft_cpx[FFTLengthForACFCalculation];
adamstark@93:     fftOut = new kiss_fft_cpx[FFTLengthForACFCalculation];
adamstark@93:     cfgForwards = kiss_fft_alloc (FFTLengthForACFCalculation, 0, 0, 0);
adamstark@93:     cfgBackwards = kiss_fft_alloc (FFTLengthForACFCalculation, 1, 0, 0);
adamstark@93: #endif
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@108: void BTrack::setHopSize (int hop)
adamstark@46: {	
adamstark@108: 	hopSize = hop;
adamstark@97: 	onsetDFBufferSize = (512 * 512) / hopSize;		// calculate df buffer size
adamstark@115: 	beatPeriod = round (60 / ((((double) hopSize) / 44100) * 120.));
adamstark@63: 
adamstark@63:     // set size of onset detection function buffer
adamstark@91:     onsetDF.resize (onsetDFBufferSize);
adamstark@63:     
adamstark@63:     // set size of cumulative score buffer
adamstark@91:     cumulativeScore.resize (onsetDFBufferSize);
adamstark@46: 	
adamstark@46: 	// initialise df_buffer to zeros
adamstark@91: 	for (int i = 0; i < onsetDFBufferSize; i++)
adamstark@46: 	{
adamstark@58: 		onsetDF[i] = 0;
adamstark@58: 		cumulativeScore[i] = 0;
adamstark@46: 		
adamstark@57: 		if ((i %  ((int) round(beatPeriod))) == 0)
adamstark@46: 		{
adamstark@58: 			onsetDF[i] = 1;
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@108: void BTrack::updateHopAndFrameSize (int hop, int frame)
adamstark@65: {
adamstark@65:     // update the onset detection function object
adamstark@108:     odf.initialise (hop, frame);
adamstark@65:     
adamstark@65:     // update the hop size being used by the beat tracker
adamstark@108:     setHopSize (hop);
adamstark@65: }
adamstark@65: 
adamstark@65: //=======================================================================
adamstark@57: bool BTrack::beatDueInCurrentFrame()
adamstark@57: {
adamstark@57:     return beatDueInFrame;
adamstark@57: }
adamstark@57: 
adamstark@57: //=======================================================================
adamstark@78: double BTrack::getCurrentTempoEstimate()
adamstark@78: {
adamstark@78:     return estimatedTempo;
adamstark@78: }
adamstark@78: 
adamstark@78: //=======================================================================
adamstark@57: int BTrack::getHopSize()
adamstark@57: {
adamstark@57:     return hopSize;
adamstark@57: }
adamstark@57: 
adamstark@57: //=======================================================================
adamstark@58: double BTrack::getLatestCumulativeScoreValue()
adamstark@58: {
adamstark@105:     return cumulativeScore[cumulativeScore.size() - 1];
adamstark@58: }
adamstark@58: 
adamstark@58: //=======================================================================
adamstark@91: void BTrack::processAudioFrame (double* frame)
adamstark@55: {
adamstark@55:     // calculate the onset detection function sample for the frame
adamstark@91:     double sample = odf.calculateOnsetDetectionFunctionSample (frame);
adamstark@55:     
adamstark@55:     // process the new onset detection function sample in the beat tracking algorithm
adamstark@91:     processOnsetDetectionFunctionSample (sample);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@91: void BTrack::processOnsetDetectionFunctionSample (double newSample)
adamstark@56: {
adamstark@56:     // we need to ensure that the onset
adamstark@56:     // detection function sample is positive
adamstark@91:     newSample = fabs (newSample);
adamstark@56:     
adamstark@56:     // add a tiny constant to the sample to stop it from ever going
adamstark@56:     // to zero. this is to avoid problems further down the line
adamstark@56:     newSample = newSample + 0.0001;
adamstark@56:     
adamstark@105: 	timeToNextPrediction--;
adamstark@105: 	timeToNextBeat--;
adamstark@57: 	beatDueInFrame = false;
adamstark@90: 		
adamstark@46: 	// add new sample at the end
adamstark@91:     onsetDF.addSampleToEnd (newSample);
adamstark@46: 	
adamstark@46: 	// update cumulative score
adamstark@91: 	updateCumulativeScore (newSample);
adamstark@46: 	
adamstark@97: 	// if we are halfway between beats, predict a beat
adamstark@105: 	if (timeToNextPrediction == 0)
adamstark@97:         predictBeat();
adamstark@46: 	
adamstark@46: 	// if we are at a beat
adamstark@105: 	if (timeToNextBeat == 0)
adamstark@46: 	{
adamstark@57: 		beatDueInFrame = true;	// indicate a beat should be output
adamstark@46: 		
adamstark@46: 		// recalculate the tempo
adamstark@57: 		resampleOnsetDetectionFunction();
adamstark@57: 		calculateTempo();
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@91: void BTrack::setTempo (double tempo)
adamstark@97: {
adamstark@46: 	/////////// TEMPO INDICATION RESET //////////////////
adamstark@46: 	
adamstark@46: 	// firstly make sure tempo is between 80 and 160 bpm..
adamstark@46: 	while (tempo > 160)
adamstark@97:         tempo = tempo / 2;
adamstark@46: 	
adamstark@46: 	while (tempo < 80)
adamstark@97:         tempo = tempo * 2;
adamstark@46: 		
adamstark@46: 	// convert tempo from bpm value to integer index of tempo probability 
adamstark@105: 	int tempoIndex = (int) round ((tempo - 80.) / 2);
adamstark@46: 	
adamstark@97:     // now set previous tempo observations to zero and set desired tempo index to 1
adamstark@97:     std::fill (prevDelta.begin(), prevDelta.end(), 0);
adamstark@105: 	prevDelta[tempoIndex] = 1;
adamstark@46: 	
adamstark@46: 	/////////// CUMULATIVE SCORE ARTIFICAL TEMPO UPDATE //////////////////
adamstark@46: 	
adamstark@46: 	// calculate new beat period
adamstark@97: 	int newBeatPeriod = (int) round (60 / ((((double) hopSize) / 44100) * tempo));
adamstark@46: 	
adamstark@97: 	int k = 1;
adamstark@97:     
adamstark@97:     // initialise onset detection function with delta functions spaced
adamstark@97:     // at the new beat period
adamstark@97: 	for (int i = onsetDFBufferSize - 1; i >= 0; i--)
adamstark@46: 	{
adamstark@97: 		if (k == 1)
adamstark@46: 		{
adamstark@58: 			cumulativeScore[i] = 150;
adamstark@58: 			onsetDF[i] = 150;
adamstark@46: 		}
adamstark@46: 		else
adamstark@46: 		{
adamstark@58: 			cumulativeScore[i] = 10;
adamstark@58: 			onsetDF[i] = 10;
adamstark@46: 		}
adamstark@46: 		
adamstark@97: 		k++;
adamstark@46: 		
adamstark@97: 		if (k > newBeatPeriod)
adamstark@46: 		{
adamstark@97: 			k = 1;
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	/////////// INDICATE THAT THIS IS A BEAT //////////////////
adamstark@46: 	
adamstark@46: 	// beat is now
adamstark@105: 	timeToNextBeat = 0;
adamstark@46: 	
adamstark@105: 	// next prediction is on the offbeat, so half of new beat period away
adamstark@105: 	timeToNextPrediction = (int) round (((double) newBeatPeriod) / 2);
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@91: void BTrack::fixTempo (double tempo)
adamstark@46: {	
adamstark@46: 	// firstly make sure tempo is between 80 and 160 bpm..
adamstark@46: 	while (tempo > 160)
adamstark@100:         tempo = tempo / 2;
adamstark@46: 	
adamstark@46: 	while (tempo < 80)
adamstark@100:         tempo = tempo * 2;
adamstark@46: 	
adamstark@46: 	// convert tempo from bpm value to integer index of tempo probability 
adamstark@100: 	int tempoIndex = (int) round((tempo - 80) / 2);
adamstark@46: 	
adamstark@46: 	// now set previous fixed previous tempo observation values to zero
adamstark@115: 	for (int i = 0; i < 41; i++)
adamstark@46: 	{
adamstark@58: 		prevDeltaFixed[i] = 0;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// set desired tempo index to 1
adamstark@100: 	prevDeltaFixed[tempoIndex] = 1;
adamstark@46: 		
adamstark@46: 	// set the tempo fix flag
adamstark@58: 	tempoFixed = true;
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::doNotFixTempo()
adamstark@46: {	
adamstark@46: 	// set the tempo fix flag
adamstark@58: 	tempoFixed = false;
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::resampleOnsetDetectionFunction()
adamstark@46: {
adamstark@46: 	float output[512];
adamstark@58:     float input[onsetDFBufferSize];
adamstark@54:     
adamstark@115:     for (int i = 0; i < onsetDFBufferSize; i++)
adamstark@58:         input[i] = (float) onsetDF[i];
adamstark@89:         
adamstark@97:     double ratio = 512.0 / ((double) onsetDFBufferSize);
adamstark@97:     int bufferLength = onsetDFBufferSize;
adamstark@97:     int outputLength = 512;
adamstark@89:     
adamstark@97:     SRC_DATA src_data;
adamstark@89:     src_data.data_in = input;
adamstark@97:     src_data.input_frames = bufferLength;
adamstark@97:     src_data.src_ratio = ratio;
adamstark@89:     src_data.data_out = output;
adamstark@97:     src_data.output_frames = outputLength;
adamstark@89:     
adamstark@89:     src_simple (&src_data, SRC_SINC_BEST_QUALITY, 1);
adamstark@89:             
adamstark@97:     for (int i = 0; i < outputLength; i++)
adamstark@89:         resampledOnsetDF[i] = (double) src_data.data_out[i];
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::calculateTempo()
adamstark@46: {
adamstark@105:     double tempoToLagFactor = 60. * 44100. / 512.;
adamstark@105:     
adamstark@46: 	// adaptive threshold on input
adamstark@100: 	adaptiveThreshold (resampledOnsetDF);
adamstark@46: 		
adamstark@46: 	// calculate auto-correlation function of detection function
adamstark@100: 	calculateBalancedACF (resampledOnsetDF);
adamstark@46: 	
adamstark@46: 	// calculate output of comb filterbank
adamstark@57: 	calculateOutputOfCombFilterBank();
adamstark@46: 	
adamstark@46: 	// adaptive threshold on rcf
adamstark@100: 	adaptiveThreshold (combFilterBankOutput);
adamstark@46: 
adamstark@59: 	// calculate tempo observation vector from beat period observation vector
adamstark@100: 	for (int i = 0; i < 41; i++)
adamstark@46: 	{
adamstark@105: 		int tempoIndex1 = (int) round (tempoToLagFactor / ((double) ((2 * i) + 80)));
adamstark@105: 		int tempoIndex2 = (int) round (tempoToLagFactor / ((double) ((4 * i) + 160)));
adamstark@100: 		tempoObservationVector[i] = combFilterBankOutput[tempoIndex1 - 1] + combFilterBankOutput[tempoIndex2 - 1];
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// if tempo is fixed then always use a fixed set of tempi as the previous observation probability function
adamstark@58: 	if (tempoFixed)
adamstark@46: 	{
adamstark@100: 		for (int k = 0; k < 41; k++)
adamstark@100:             prevDelta[k] = prevDeltaFixed[k];
adamstark@46: 	}
adamstark@46: 		
adamstark@100: 	for (int j = 0; j < 41; j++)
adamstark@46: 	{
adamstark@100: 		double maxValue = -1;
adamstark@100:         
adamstark@100: 		for (int i = 0; i < 41; i++)
adamstark@46: 		{
adamstark@100: 			double currentValue = prevDelta[i] * tempoTransitionMatrix[i][j];
adamstark@46: 			
adamstark@100: 			if (currentValue > maxValue)
adamstark@100:                 maxValue = currentValue;
adamstark@46: 		}
adamstark@46: 		
adamstark@100: 		delta[j] = maxValue * tempoObservationVector[j];
adamstark@46: 	}
adamstark@46: 	
adamstark@100: 	normaliseVector (delta);
adamstark@46: 	
adamstark@100: 	double maxIndex = -1;
adamstark@100: 	double maxValue = -1;
adamstark@46: 	
adamstark@100: 	for (int j = 0; j < 41; j++)
adamstark@46: 	{
adamstark@100: 		if (delta[j] > maxValue)
adamstark@46: 		{
adamstark@100: 			maxValue = delta[j];
adamstark@100: 			maxIndex = j;
adamstark@46: 		}
adamstark@46: 		
adamstark@58: 		prevDelta[j] = delta[j];
adamstark@46: 	}
adamstark@46: 	
adamstark@100: 	beatPeriod = round ((60.0 * 44100.0) / (((2 * maxIndex) + 80) * ((double) hopSize)));
adamstark@46: 	
adamstark@57: 	if (beatPeriod > 0)
adamstark@115:         estimatedTempo = 60.0 / ((((double) hopSize) / 44100.0) * beatPeriod);
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@100: void BTrack::adaptiveThreshold (std::vector<double>& x)
adamstark@46: {
adamstark@100:     int N = static_cast<int> (x.size());
adamstark@100: 	double threshold[N];
adamstark@46: 	
adamstark@46: 	int p_post = 7;
adamstark@46: 	int p_pre = 8;
adamstark@46: 	
adamstark@100: 	int t = std::min (N, p_post);	// what is smaller, p_post or df size. This is to avoid accessing outside of arrays
adamstark@46: 	
adamstark@46: 	// find threshold for first 't' samples, where a full average cannot be computed yet 
adamstark@100: 	for (int i = 0; i <= t; i++)
adamstark@46: 	{	
adamstark@100: 		int k = std::min ((i + p_pre), N);
adamstark@100: 		threshold[i] = calculateMeanOfVector (x, 1, k);
adamstark@46: 	}
adamstark@100:     
adamstark@46: 	// find threshold for bulk of samples across a moving average from [i-p_pre,i+p_post]
adamstark@100: 	for (int i = t + 1; i < N - p_post; i++)
adamstark@46: 	{
adamstark@100: 		threshold[i] = calculateMeanOfVector (x, i - p_pre, i + p_post);
adamstark@46: 	}
adamstark@100:     
adamstark@46: 	// for last few samples calculate threshold, again, not enough samples to do as above
adamstark@100: 	for (int i = N - p_post; i < N; i++)
adamstark@46: 	{
adamstark@100: 		int k = std::max ((i - p_post), 1);
adamstark@100: 		threshold[i] = calculateMeanOfVector (x, k, N);
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// subtract the threshold from the detection function and check that it is not less than 0
adamstark@100: 	for (int i = 0; i < N; i++)
adamstark@46: 	{
adamstark@100: 		x[i] = x[i] - threshold[i];
adamstark@100:         
adamstark@46: 		if (x[i] < 0)
adamstark@100:             x[i] = 0;
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::calculateOutputOfCombFilterBank()
adamstark@46: {
adamstark@100:     std::fill (combFilterBankOutput.begin(), combFilterBankOutput.end(), 0.0);
adamstark@100: 	int numCombElements = 4;
adamstark@46: 	
adamstark@91: 	for (int i = 2; i <= 127; i++) // max beat period
adamstark@46: 	{
adamstark@100: 		for (int a = 1; a <= numCombElements; a++) // number of comb elements
adamstark@46: 		{
adamstark@100: 			for (int b = 1 - a; b <= a - 1; b++) // general state using normalisation of comb elements
adamstark@46: 			{
adamstark@115: 				combFilterBankOutput[i - 1] += (acf[(a * i + b) - 1] * weightingVector[i - 1]) / (2 * a - 1);	// calculate value for comb filter row
adamstark@46: 			}
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@100: void BTrack::calculateBalancedACF (std::vector<double>& onsetDetectionFunction)
adamstark@46: {
adamstark@88:     int onsetDetectionFunctionLength = 512;
adamstark@88:     
adamstark@93: #ifdef USE_FFTW
adamstark@88:     // copy into complex array and zero pad
adamstark@111:     for (int i = 0; i < FFTLengthForACFCalculation; i++)
adamstark@88:     {
adamstark@88:         if (i < onsetDetectionFunctionLength)
adamstark@88:         {
adamstark@88:             complexIn[i][0] = onsetDetectionFunction[i];
adamstark@88:             complexIn[i][1] = 0.0;
adamstark@88:         }
adamstark@88:         else
adamstark@88:         {
adamstark@88:             complexIn[i][0] = 0.0;
adamstark@88:             complexIn[i][1] = 0.0;
adamstark@88:         }
adamstark@88:     }
adamstark@88:     
adamstark@88:     // perform the fft
adamstark@91:     fftw_execute (acfForwardFFT);
adamstark@88:     
adamstark@88:     // multiply by complex conjugate
adamstark@115:     for (int i = 0; i < FFTLengthForACFCalculation; i++)
adamstark@88:     {
adamstark@111:         complexOut[i][0] = complexOut[i][0] * complexOut[i][0] + complexOut[i][1] * complexOut[i][1];
adamstark@88:         complexOut[i][1] = 0.0;
adamstark@88:     }
adamstark@88:     
adamstark@88:     // perform the ifft
adamstark@91:     fftw_execute (acfBackwardFFT);
adamstark@88:     
adamstark@93: #endif
adamstark@93:     
adamstark@93: #ifdef USE_KISS_FFT
adamstark@93:     // copy into complex array and zero pad
adamstark@111:     for (int i = 0; i < FFTLengthForACFCalculation; i++)
adamstark@93:     {
adamstark@93:         if (i < onsetDetectionFunctionLength)
adamstark@93:         {
adamstark@93:             fftIn[i].r = onsetDetectionFunction[i];
adamstark@93:             fftIn[i].i = 0.0;
adamstark@93:         }
adamstark@93:         else
adamstark@93:         {
adamstark@93:             fftIn[i].r = 0.0;
adamstark@93:             fftIn[i].i = 0.0;
adamstark@93:         }
adamstark@93:     }
adamstark@93:     
adamstark@93:     // execute kiss fft
adamstark@93:     kiss_fft (cfgForwards, fftIn, fftOut);
adamstark@93:     
adamstark@93:     // multiply by complex conjugate
adamstark@115:     for (int i = 0; i < FFTLengthForACFCalculation; i++)
adamstark@93:     {
adamstark@93:         fftOut[i].r = fftOut[i].r * fftOut[i].r + fftOut[i].i * fftOut[i].i;
adamstark@93:         fftOut[i].i = 0.0;
adamstark@93:     }
adamstark@93:     
adamstark@93:     // perform the ifft
adamstark@93:     kiss_fft (cfgBackwards, fftOut, fftIn);
adamstark@93:     
adamstark@93: #endif
adamstark@88:     
adamstark@88:     double lag = 512;
adamstark@88:     
adamstark@91:     for (int i = 0; i < 512; i++)
adamstark@88:     {
adamstark@93: #ifdef USE_FFTW
adamstark@88:         // calculate absolute value of result
adamstark@111:         double absValue = sqrt (complexIn[i][0] * complexIn[i][0] + complexIn[i][1] * complexIn[i][1]);
adamstark@93: #endif
adamstark@88:         
adamstark@93: #ifdef USE_KISS_FFT
adamstark@93:         // calculate absolute value of result
adamstark@93:         double absValue = sqrt (fftIn[i].r * fftIn[i].r + fftIn[i].i * fftIn[i].i);
adamstark@93: #endif
adamstark@88:         // divide by inverse lad to deal with scale bias towards small lags
adamstark@88:         acf[i] = absValue / lag;
adamstark@88:         
adamstark@88:         // this division by 1024 is technically unnecessary but it ensures the algorithm produces
adamstark@88:         // exactly the same ACF output as the old time domain implementation. The time difference is
adamstark@88:         // minimal so I decided to keep it
adamstark@88:         acf[i] = acf[i] / 1024.;
adamstark@88:         
adamstark@88:         lag = lag - 1.;
adamstark@88:     }
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@100: double BTrack::calculateMeanOfVector (std::vector<double>& vector, int startIndex, int endIndex)
adamstark@46: {
adamstark@97:     int length = endIndex - startIndex;
adamstark@100:     double sum = std::accumulate (vector.begin() + startIndex, vector.begin() + endIndex, 0.0);
adamstark@47: 
adamstark@47:     if (length > 0)
adamstark@97:         return sum / static_cast<double> (length);	// average and return
adamstark@47:     else
adamstark@47:         return 0;
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@100: void BTrack::normaliseVector (std::vector<double>& vector)
adamstark@46: {
adamstark@100:     double sum = std::accumulate (vector.begin(), vector.end(), 0.0);
adamstark@46: 	
adamstark@46: 	if (sum > 0)
adamstark@97:     {
adamstark@100:         for (int i = 0; i < vector.size(); i++)
adamstark@100:             vector[i] = vector[i] / sum;
adamstark@97:     }
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@100: void BTrack::updateCumulativeScore (double onsetDetectionFunctionSample)
adamstark@98: {
adamstark@100: 	int windowStart = onsetDFBufferSize - round (2. * beatPeriod);
adamstark@100: 	int windowEnd = onsetDFBufferSize - round (beatPeriod / 2.);
adamstark@100: 	int windowSize = windowEnd - windowStart + 1;
adamstark@46: 	
adamstark@104:     // create log gaussian transition window
adamstark@103: 	double logGaussianTransitionWeighting[windowSize];
adamstark@103:     createLogGaussianTransitionWeighting (logGaussianTransitionWeighting, windowSize, beatPeriod);
adamstark@46: 	
adamstark@104:     // calculate the new cumulative score value
adamstark@105:     double cumulativeScoreValue = calculateNewCumulativeScoreValue (cumulativeScore, logGaussianTransitionWeighting, windowStart, windowEnd, onsetDetectionFunctionSample, alpha);
adamstark@104:     
adamstark@104:     // add the new cumulative score value to the buffer
adamstark@105:     cumulativeScore.addSampleToEnd (cumulativeScoreValue);
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::predictBeat()
adamstark@46: {	 
adamstark@104: 	int beatExpectationWindowSize = static_cast<int> (beatPeriod);
adamstark@102: 	double futureCumulativeScore[onsetDFBufferSize + beatExpectationWindowSize];
adamstark@102: 	double beatExpectationWindow[beatExpectationWindowSize];
adamstark@93:     
adamstark@102: 	// copy cumulativeScore to first part of futureCumulativeScore
adamstark@115: 	for (int i = 0; i < onsetDFBufferSize; i++)
adamstark@102:         futureCumulativeScore[i] = cumulativeScore[i];
adamstark@102: 	
adamstark@102: 	// Create a beat expectation window for predicting future beats from the "future" of the cumulative score.
adamstark@102:     // We are making this beat prediction at the midpoint between beats, and so we make a Gaussian
adamstark@102:     // weighting centred on the most likely beat position (half a beat period into the future)
adamstark@102:     // This is W2 in Adam Stark's PhD thesis, equation 3.6, page 62
adamstark@102:     
adamstark@102: 	double v = 1;
adamstark@102: 	for (int i = 0; i < beatExpectationWindowSize; i++)
adamstark@46: 	{
adamstark@102: 		beatExpectationWindow[i] = exp((-1 * pow ((v - (beatPeriod / 2)), 2))   /  (2 * pow (beatPeriod / 2, 2)));
adamstark@46: 		v++;
adamstark@46: 	}
adamstark@46: 	
adamstark@102: 	// Create window for "synthesizing" the cumulative score into the future
adamstark@102:     // It is a log-Gaussian transition weighting running from from 2 beat periods
adamstark@102:     // in the past to half a beat period in the past. It favours the time exactly
adamstark@102:     // one beat period in the past
adamstark@102:     
adamstark@102: 	int startIndex = onsetDFBufferSize - round (2 * beatPeriod);
adamstark@102: 	int endIndex = onsetDFBufferSize - round (beatPeriod / 2);
adamstark@102: 	int pastWindowSize = endIndex - startIndex + 1;
adamstark@103:     
adamstark@102: 	double logGaussianTransitionWeighting[pastWindowSize];
adamstark@103:     createLogGaussianTransitionWeighting (logGaussianTransitionWeighting, pastWindowSize, beatPeriod);
adamstark@46: 
adamstark@104: 	// Calculate the future cumulative score, by shifting the log Gaussian transition weighting from its
adamstark@106:     // start position of [-2 beat periods, - 0.5 beat periods] forwards over the size of the beat
adamstark@104:     // expectation window, calculating a new cumulative score where the onset detection function sample
adamstark@104:     // is zero. This uses the "momentum" of the function to generate itself into the future.
adamstark@102: 	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + beatExpectationWindowSize); i++)
adamstark@46: 	{
adamstark@106:         // note here that we pass 0.0 in for the onset detection function sample and 1.0 for the alpha weighting factor
adamstark@104:         // see equation 3.4 and page 60 - 62 of Adam Stark's PhD thesis for details
adamstark@104:         futureCumulativeScore[i] = calculateNewCumulativeScoreValue (futureCumulativeScore, logGaussianTransitionWeighting, startIndex, endIndex, 0.0, 1.0);
adamstark@104:         
adamstark@104:         startIndex++;
adamstark@104:         endIndex++;
adamstark@46: 	}
adamstark@46: 	
adamstark@102: 	// Predict the next beat, finding the maximum point of the future cumulative score
adamstark@102:     // over the next beat, after being weighted by the beat expectation window
adamstark@102:     
adamstark@102: 	double maxValue = 0;
adamstark@102: 	int n = 0;
adamstark@46: 	
adamstark@102: 	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + beatExpectationWindowSize); i++)
adamstark@46: 	{
adamstark@102: 		double weightedCumulativeScore = futureCumulativeScore[i] * beatExpectationWindow[n];
adamstark@46: 		
adamstark@102: 		if (weightedCumulativeScore > maxValue)
adamstark@46: 		{
adamstark@102: 			maxValue = weightedCumulativeScore;
adamstark@105: 			timeToNextBeat = n;
adamstark@46: 		}	
adamstark@46: 		
adamstark@46: 		n++;
adamstark@46: 	}
adamstark@46: 		
adamstark@105: 	// set next prediction time as on the offbeat after the next beat
adamstark@105: 	timeToNextPrediction = timeToNextBeat + round (beatPeriod / 2);
adamstark@97: }
adamstark@103: 
adamstark@103: //=======================================================================
adamstark@103: void BTrack::createLogGaussianTransitionWeighting (double* weightingArray, int numSamples, double beatPeriod)
adamstark@103: {
adamstark@105:     // (This is W1 in Adam Stark's PhD thesis, equation 3.2, page 60)
adamstark@105:     
adamstark@103:     double v = -2. * beatPeriod;
adamstark@103:     
adamstark@103:     for (int i = 0; i < numSamples; i++)
adamstark@103:     {
adamstark@103:         double a = tightness * log (-v / beatPeriod);
adamstark@103:         weightingArray[i] = exp ((-1. * a * a) / 2.);
adamstark@103:         v++;
adamstark@103:     }
adamstark@103: }
adamstark@104: 
adamstark@104: //=======================================================================
adamstark@104: template <typename T>
adamstark@104: double BTrack::calculateNewCumulativeScoreValue (T cumulativeScoreArray, double* logGaussianTransitionWeighting, int startIndex, int endIndex, double onsetDetectionFunctionSample, double alphaWeightingFactor)
adamstark@104: {
adamstark@104:     // calculate new cumulative score value by weighting the cumulative score between
adamstark@104:     // startIndex and endIndex and finding the maximum value
adamstark@104:     double maxValue = 0;
adamstark@104:     int n = 0;
adamstark@104:     for (int i = startIndex; i <= endIndex; i++)
adamstark@104:     {
adamstark@104:         double weightedCumulativeScore = cumulativeScoreArray[i] * logGaussianTransitionWeighting[n];
adamstark@104:         
adamstark@104:         if (weightedCumulativeScore > maxValue)
adamstark@104:             maxValue = weightedCumulativeScore;
adamstark@104:         
adamstark@104:         n++;
adamstark@104:     }
adamstark@104:     
adamstark@104:     // now mix with the incoming onset detection function sample
adamstark@104:     // (equation 3.4 on page 60 of Adam Stark's PhD thesis)
adamstark@104:     double cumulativeScoreValue = ((1. - alphaWeightingFactor) * onsetDetectionFunctionSample) + (alphaWeightingFactor * maxValue);
adamstark@104:     
adamstark@104:     return cumulativeScoreValue;
adamstark@104: }