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@46: #include "BTrack.h"
adamstark@46: #include "samplerate.h"
adamstark@89: #include <iostream>
adamstark@46: 
adamstark@55: //=======================================================================
adamstark@57: BTrack::BTrack() : odf(512,1024,ComplexSpectralDifferenceHWR,HanningWindow)
adamstark@55: {
adamstark@55:     initialise(512, 1024);
adamstark@55: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: BTrack::BTrack(int hopSize_) : odf(hopSize_,2*hopSize_,ComplexSpectralDifferenceHWR,HanningWindow)
adamstark@46: {	
adamstark@57:     initialise(hopSize_, 2*hopSize_);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@57: BTrack::BTrack(int hopSize_,int frameSize_) : odf(hopSize_,frameSize_,ComplexSpectralDifferenceHWR,HanningWindow)
adamstark@55: {
adamstark@57:     initialise(hopSize_, frameSize_);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@88: BTrack::~BTrack()
adamstark@88: {
adamstark@88:     // destroy fft plan
adamstark@88:     fftw_destroy_plan(acfForwardFFT);
adamstark@88:     fftw_destroy_plan(acfBackwardFFT);
adamstark@88:     fftw_free(complexIn);
adamstark@88:     fftw_free(complexOut);
adamstark@88: }
adamstark@88: 
adamstark@88: //=======================================================================
adamstark@55: double BTrack::getBeatTimeInSeconds(long frameNumber,int hopSize,int fs)
adamstark@55: {
adamstark@55:     double hop = (double) hopSize;
adamstark@55:     double samplingFrequency = (double) fs;
adamstark@55:     double frameNum = (double) frameNumber;
adamstark@55:     
adamstark@55:     return ((hop / samplingFrequency) * frameNum);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@55: double BTrack::getBeatTimeInSeconds(int frameNumber,int hopSize,int fs)
adamstark@55: {
adamstark@55:     long frameNum = (long) frameNumber;
adamstark@55:     
adamstark@55:     return getBeatTimeInSeconds(frameNum, hopSize, fs);
adamstark@55: }
adamstark@55: 
adamstark@55: 
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@57: void BTrack::initialise(int hopSize_, int frameSize_)
adamstark@55: {
adamstark@55:     double rayparam = 43;
adamstark@54: 	double pi = 3.14159265;
adamstark@46: 	
adamstark@46: 	
adamstark@46: 	// initialise parameters
adamstark@46: 	tightness = 5;
adamstark@46: 	alpha = 0.9;
adamstark@46: 	tempo = 120;
adamstark@58: 	estimatedTempo = 120.0;
adamstark@59: 	tempoToLagFactor = 60.*44100./512.;
adamstark@46: 	
adamstark@46: 	m0 = 10;
adamstark@58: 	beatCounter = -1;
adamstark@46: 	
adamstark@57: 	beatDueInFrame = false;
adamstark@46: 	
adamstark@58: 
adamstark@46: 	// create rayleigh weighting vector
adamstark@46: 	for (int n = 0;n < 128;n++)
adamstark@46: 	{
adamstark@58: 		weightingVector[n] = ((double) n / pow(rayparam,2)) * exp((-1*pow((double)-n,2)) / (2*pow(rayparam,2)));
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// initialise prev_delta
adamstark@46: 	for (int i = 0;i < 41;i++)
adamstark@46: 	{
adamstark@58: 		prevDelta[i] = 1;
adamstark@46: 	}
adamstark@46: 	
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@46: 	for (int i = 0;i < 41;i++)
adamstark@46: 	{
adamstark@46: 		for (int j = 0;j < 41;j++)
adamstark@46: 		{
adamstark@46: 			x = j+1;
adamstark@46: 			t_mu = i+1;
adamstark@58: 			tempoTransitionMatrix[i][j] = (1 / (m_sig * sqrt(2*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@58:     // initialise latest cumulative score value
adamstark@58:     // in case it is requested before any processing takes place
adamstark@58:     latestCumulativeScoreValue = 0;
adamstark@55:     
adamstark@55:     // initialise algorithm given the hopsize
adamstark@57:     setHopSize(hopSize_);
adamstark@88:     
adamstark@88:     
adamstark@88:     // Set up FFT for calculating the auto-correlation function
adamstark@88:     FFTLengthForACFCalculation = 1024;
adamstark@88:     
adamstark@88:     complexIn = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFTLengthForACFCalculation);		// complex array to hold fft data
adamstark@88:     complexOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFTLengthForACFCalculation);	// complex array to hold fft data
adamstark@88:     
adamstark@88:     acfForwardFFT = fftw_plan_dft_1d(FFTLengthForACFCalculation, complexIn, complexOut, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
adamstark@88:     acfBackwardFFT = fftw_plan_dft_1d(FFTLengthForACFCalculation, complexOut, complexIn, FFTW_BACKWARD, FFTW_ESTIMATE);	// FFT plan initialisation
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::setHopSize(int hopSize_)
adamstark@46: {	
adamstark@57: 	hopSize = hopSize_;
adamstark@58: 	onsetDFBufferSize = (512*512)/hopSize;		// calculate df buffer size
adamstark@46: 	
adamstark@57: 	beatPeriod = round(60/((((double) hopSize)/44100)*tempo));
adamstark@63: 
adamstark@63:     // set size of onset detection function buffer
adamstark@63:     onsetDF.resize(onsetDFBufferSize);
adamstark@63:     
adamstark@63:     // set size of cumulative score buffer
adamstark@63:     cumulativeScore.resize(onsetDFBufferSize);
adamstark@46: 	
adamstark@46: 	// initialise df_buffer to zeros
adamstark@58: 	for (int i = 0;i < onsetDFBufferSize;i++)
adamstark@46: 	{
adamstark@58: 		onsetDF[i] = 0;
adamstark@58: 		cumulativeScore[i] = 0;
adamstark@46: 		
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@65: void BTrack::updateHopAndFrameSize(int hopSize_,int frameSize_)
adamstark@65: {
adamstark@65:     // update the onset detection function object
adamstark@66:     odf.initialise(hopSize_, frameSize_);
adamstark@65:     
adamstark@65:     // update the hop size being used by the beat tracker
adamstark@65:     setHopSize(hopSize_);
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@58:     return latestCumulativeScoreValue;
adamstark@58: }
adamstark@58: 
adamstark@58: //=======================================================================
adamstark@55: void BTrack::processAudioFrame(double *frame)
adamstark@55: {
adamstark@55:     // calculate the onset detection function sample for the frame
adamstark@59:     double sample = odf.calculateOnsetDetectionFunctionSample(frame);
adamstark@55:     
adamstark@56:     
adamstark@55:     
adamstark@55:     // process the new onset detection function sample in the beat tracking algorithm
adamstark@55:     processOnsetDetectionFunctionSample(sample);
adamstark@55: }
adamstark@55: 
adamstark@55: //=======================================================================
adamstark@55: void BTrack::processOnsetDetectionFunctionSample(double newSample)
adamstark@56: {
adamstark@56:     // we need to ensure that the onset
adamstark@56:     // detection function sample is positive
adamstark@56:     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@46: 	m0--;
adamstark@58: 	beatCounter--;
adamstark@57: 	beatDueInFrame = false;
adamstark@90: 		
adamstark@46: 	// add new sample at the end
adamstark@89:     onsetDF.addSampleToEnd(newSample);
adamstark@46: 	
adamstark@46: 	// update cumulative score
adamstark@57: 	updateCumulativeScore(newSample);
adamstark@46: 	
adamstark@46: 	// if we are halfway between beats
adamstark@46: 	if (m0 == 0)
adamstark@46: 	{
adamstark@57: 		predictBeat();
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// if we are at a beat
adamstark@58: 	if (beatCounter == 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@57: void BTrack::setTempo(double tempo)
adamstark@46: {	 
adamstark@46: 	
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@46: 	{
adamstark@46: 		tempo = tempo/2;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	while (tempo < 80)
adamstark@46: 	{
adamstark@46: 		tempo = tempo * 2;
adamstark@46: 	}
adamstark@46: 		
adamstark@46: 	// convert tempo from bpm value to integer index of tempo probability 
adamstark@46: 	int tempo_index = (int) round((tempo - 80)/2);
adamstark@46: 	
adamstark@46: 	// now set previous tempo observations to zero
adamstark@46: 	for (int i=0;i < 41;i++)
adamstark@46: 	{
adamstark@58: 		prevDelta[i] = 0;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// set desired tempo index to 1
adamstark@58: 	prevDelta[tempo_index] = 1;
adamstark@46: 	
adamstark@46: 	
adamstark@46: 	/////////// CUMULATIVE SCORE ARTIFICAL TEMPO UPDATE //////////////////
adamstark@46: 	
adamstark@46: 	// calculate new beat period
adamstark@57: 	int new_bperiod = (int) round(60/((((double) hopSize)/44100)*tempo));
adamstark@46: 	
adamstark@46: 	int bcounter = 1;
adamstark@46: 	// initialise df_buffer to zeros
adamstark@58: 	for (int i = (onsetDFBufferSize-1);i >= 0;i--)
adamstark@46: 	{
adamstark@46: 		if (bcounter == 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@46: 		bcounter++;
adamstark@46: 		
adamstark@46: 		if (bcounter > new_bperiod)
adamstark@46: 		{
adamstark@46: 			bcounter = 1;
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	/////////// INDICATE THAT THIS IS A BEAT //////////////////
adamstark@46: 	
adamstark@46: 	// beat is now
adamstark@58: 	beatCounter = 0;
adamstark@46: 	
adamstark@46: 	// offbeat is half of new beat period away
adamstark@54: 	m0 = (int) round(((double) new_bperiod)/2);
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: 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@46: 	{
adamstark@46: 		tempo = tempo/2;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	while (tempo < 80)
adamstark@46: 	{
adamstark@46: 		tempo = tempo * 2;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// convert tempo from bpm value to integer index of tempo probability 
adamstark@46: 	int tempo_index = (int) round((tempo - 80)/2);
adamstark@46: 	
adamstark@46: 	// now set previous fixed previous tempo observation values to zero
adamstark@46: 	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@58: 	prevDeltaFixed[tempo_index] = 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@89:     
adamstark@89: 
adamstark@58:     float input[onsetDFBufferSize];
adamstark@54:     
adamstark@58:     for (int i = 0;i < onsetDFBufferSize;i++)
adamstark@54:     {
adamstark@58:         input[i] = (float) onsetDF[i];
adamstark@54:     }
adamstark@89:         
adamstark@89:     double src_ratio = 512.0/((double) onsetDFBufferSize);
adamstark@89:     int BUFFER_LEN = onsetDFBufferSize;
adamstark@89:     int output_len;
adamstark@89:     SRC_DATA	src_data ;
adamstark@89:     
adamstark@89:     //output_len = (int) floor (((double) BUFFER_LEN) * src_ratio) ;
adamstark@89:     output_len = 512;
adamstark@89:     
adamstark@89:     src_data.data_in = input;
adamstark@89:     src_data.input_frames = BUFFER_LEN;
adamstark@89:     
adamstark@89:     src_data.src_ratio = src_ratio;
adamstark@89:     
adamstark@89:     src_data.data_out = output;
adamstark@89:     src_data.output_frames = output_len;
adamstark@89:     
adamstark@89:     src_simple (&src_data, SRC_SINC_BEST_QUALITY, 1);
adamstark@89:             
adamstark@89:     for (int i = 0;i < output_len;i++)
adamstark@89:     {
adamstark@89:         resampledOnsetDF[i] = (double) src_data.data_out[i];
adamstark@89:     }
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::calculateTempo()
adamstark@46: {
adamstark@46: 	// adaptive threshold on input
adamstark@58: 	adaptiveThreshold(resampledOnsetDF,512);
adamstark@46: 		
adamstark@46: 	// calculate auto-correlation function of detection function
adamstark@58: 	calculateBalancedACF(resampledOnsetDF);
adamstark@46: 	
adamstark@46: 	// calculate output of comb filterbank
adamstark@57: 	calculateOutputOfCombFilterBank();
adamstark@46: 	
adamstark@46: 	
adamstark@46: 	// adaptive threshold on rcf
adamstark@58: 	adaptiveThreshold(combFilterBankOutput,128);
adamstark@46: 
adamstark@46: 	
adamstark@46: 	int t_index;
adamstark@46: 	int t_index2;
adamstark@59: 	// calculate tempo observation vector from beat period observation vector
adamstark@46: 	for (int i = 0;i < 41;i++)
adamstark@46: 	{
adamstark@59: 		t_index = (int) round(tempoToLagFactor / ((double) ((2*i)+80)));
adamstark@59: 		t_index2 = (int) round(tempoToLagFactor / ((double) ((4*i)+160)));
adamstark@46: 
adamstark@46: 		
adamstark@58: 		tempoObservationVector[i] = combFilterBankOutput[t_index-1] + combFilterBankOutput[t_index2-1];
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	
adamstark@54: 	double maxval;
adamstark@54: 	double maxind;
adamstark@54: 	double curval;
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@46: 		for (int k = 0;k < 41;k++)
adamstark@46: 		{
adamstark@58: 			prevDelta[k] = prevDeltaFixed[k];
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: 		
adamstark@46: 	for (int j=0;j < 41;j++)
adamstark@46: 	{
adamstark@46: 		maxval = -1;
adamstark@46: 		for (int i = 0;i < 41;i++)
adamstark@46: 		{
adamstark@58: 			curval = prevDelta[i]*tempoTransitionMatrix[i][j];
adamstark@46: 			
adamstark@46: 			if (curval > maxval)
adamstark@46: 			{
adamstark@46: 				maxval = curval;
adamstark@46: 			}
adamstark@46: 		}
adamstark@46: 		
adamstark@58: 		delta[j] = maxval*tempoObservationVector[j];
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 
adamstark@57: 	normaliseArray(delta,41);
adamstark@46: 	
adamstark@46: 	maxind = -1;
adamstark@46: 	maxval = -1;
adamstark@46: 	
adamstark@46: 	for (int j=0;j < 41;j++)
adamstark@46: 	{
adamstark@46: 		if (delta[j] > maxval)
adamstark@46: 		{
adamstark@46: 			maxval = delta[j];
adamstark@46: 			maxind = j;
adamstark@46: 		}
adamstark@46: 		
adamstark@58: 		prevDelta[j] = delta[j];
adamstark@46: 	}
adamstark@46: 	
adamstark@57: 	beatPeriod = round((60.0*44100.0)/(((2*maxind)+80)*((double) hopSize)));
adamstark@46: 	
adamstark@57: 	if (beatPeriod > 0)
adamstark@46: 	{
adamstark@58: 		estimatedTempo = 60.0/((((double) hopSize) / 44100.0)*beatPeriod);
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::adaptiveThreshold(double *x,int N)
adamstark@46: {
adamstark@46: 	int i = 0;
adamstark@46: 	int k,t = 0;
adamstark@54: 	double x_thresh[N];
adamstark@46: 	
adamstark@46: 	int p_post = 7;
adamstark@46: 	int p_pre = 8;
adamstark@46: 	
adamstark@52: 	t = std::min(N,p_post);	// what is smaller, p_post of 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@46: 	for (i = 0;i <= t;i++)
adamstark@46: 	{	
adamstark@52: 		k = std::min((i+p_pre),N);
adamstark@57: 		x_thresh[i] = calculateMeanOfArray(x,1,k);
adamstark@46: 	}
adamstark@46: 	// find threshold for bulk of samples across a moving average from [i-p_pre,i+p_post]
adamstark@46: 	for (i = t+1;i < N-p_post;i++)
adamstark@46: 	{
adamstark@57: 		x_thresh[i] = calculateMeanOfArray(x,i-p_pre,i+p_post);
adamstark@46: 	}
adamstark@46: 	// for last few samples calculate threshold, again, not enough samples to do as above
adamstark@46: 	for (i = N-p_post;i < N;i++)
adamstark@46: 	{
adamstark@52: 		k = std::max((i-p_post),1);
adamstark@57: 		x_thresh[i] = calculateMeanOfArray(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@46: 	for (i = 0;i < N;i++)
adamstark@46: 	{
adamstark@46: 		x[i] = x[i] - x_thresh[i];
adamstark@46: 		if (x[i] < 0)
adamstark@46: 		{
adamstark@46: 			x[i] = 0;
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::calculateOutputOfCombFilterBank()
adamstark@46: {
adamstark@46: 	int numelem;
adamstark@46: 	
adamstark@46: 	for (int i = 0;i < 128;i++)
adamstark@46: 	{
adamstark@58: 		combFilterBankOutput[i] = 0;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	numelem = 4;
adamstark@46: 	
adamstark@46: 	for (int i = 2;i <= 127;i++) // max beat period
adamstark@46: 	{
adamstark@46: 		for (int a = 1;a <= numelem;a++) // number of comb elements
adamstark@46: 		{
adamstark@46: 			for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements
adamstark@46: 			{
adamstark@58: 				combFilterBankOutput[i-1] = 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@60: void BTrack::calculateBalancedACF(double *onsetDetectionFunction)
adamstark@46: {
adamstark@88:     int onsetDetectionFunctionLength = 512;
adamstark@88:     
adamstark@88:     // copy into complex array and zero pad
adamstark@88:     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@88:     fftw_execute(acfForwardFFT);
adamstark@88:     
adamstark@88:     // multiply by complex conjugate
adamstark@88:     for (int i = 0;i < FFTLengthForACFCalculation;i++)
adamstark@88:     {
adamstark@88:         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@88:     fftw_execute(acfBackwardFFT);
adamstark@88:     
adamstark@88:     
adamstark@88:     double lag = 512;
adamstark@88:     
adamstark@88:     for (int i = 0;i < 512;i++)
adamstark@88:     {
adamstark@88:         // calculate absolute value of result
adamstark@88:         double absValue = sqrt(complexIn[i][0]*complexIn[i][0] + complexIn[i][1]*complexIn[i][1]);
adamstark@88:         
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@59: double BTrack::calculateMeanOfArray(double *array,int startIndex,int endIndex)
adamstark@46: {
adamstark@46: 	int i;
adamstark@47: 	double sum = 0;
adamstark@47: 
adamstark@59:     int length = endIndex - startIndex;
adamstark@46: 	
adamstark@46: 	// find sum
adamstark@59: 	for (i = startIndex;i < endIndex;i++)
adamstark@46: 	{
adamstark@46: 		sum = sum + array[i];
adamstark@46: 	}
adamstark@46: 	
adamstark@47:     if (length > 0)
adamstark@47:     {
adamstark@47:         return sum / length;	// average and return
adamstark@47:     }
adamstark@47:     else
adamstark@47:     {
adamstark@47:         return 0;
adamstark@47:     }
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::normaliseArray(double *array,int N)
adamstark@46: {
adamstark@46: 	double sum = 0;
adamstark@46: 	
adamstark@46: 	for (int i = 0;i < N;i++)
adamstark@46: 	{
adamstark@46: 		if (array[i] > 0)
adamstark@46: 		{
adamstark@46: 			sum = sum + array[i];
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	if (sum > 0)
adamstark@46: 	{
adamstark@46: 		for (int i = 0;i < N;i++)
adamstark@46: 		{
adamstark@46: 			array[i] = array[i] / sum;
adamstark@46: 		}
adamstark@46: 	}
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@59: void BTrack::updateCumulativeScore(double odfSample)
adamstark@46: {	 
adamstark@46: 	int start, end, winsize;
adamstark@54: 	double max;
adamstark@46: 	
adamstark@58: 	start = onsetDFBufferSize - round(2*beatPeriod);
adamstark@58: 	end = onsetDFBufferSize - round(beatPeriod/2);
adamstark@46: 	winsize = end-start+1;
adamstark@46: 	
adamstark@54: 	double w1[winsize];
adamstark@57: 	double v = -2*beatPeriod;
adamstark@54: 	double wcumscore;
adamstark@46: 	
adamstark@46: 	
adamstark@46: 	// create window
adamstark@46: 	for (int i = 0;i < winsize;i++)
adamstark@46: 	{
adamstark@57: 		w1[i] = exp((-1*pow(tightness*log(-v/beatPeriod),2))/2);
adamstark@46: 		v = v+1;
adamstark@46: 	}	
adamstark@46: 	
adamstark@46: 	// calculate new cumulative score value
adamstark@46: 	max = 0;
adamstark@46: 	int n = 0;
adamstark@46: 	for (int i=start;i <= end;i++)
adamstark@46: 	{
adamstark@58: 			wcumscore = cumulativeScore[i]*w1[n];
adamstark@46: 		
adamstark@46: 			if (wcumscore > max)
adamstark@46: 			{
adamstark@46: 				max = wcumscore;
adamstark@46: 			}
adamstark@46: 		n++;
adamstark@46: 	}
adamstark@46: 	
adamstark@90: 		
adamstark@89:     latestCumulativeScoreValue = ((1-alpha)*odfSample) + (alpha*max);
adamstark@89:     
adamstark@89:     cumulativeScore.addSampleToEnd(latestCumulativeScoreValue);
adamstark@46: }
adamstark@46: 
adamstark@51: //=======================================================================
adamstark@57: void BTrack::predictBeat()
adamstark@46: {	 
adamstark@58: 	int windowSize = (int) beatPeriod;
adamstark@58: 	double futureCumulativeScore[onsetDFBufferSize + windowSize];
adamstark@58: 	double w2[windowSize];
adamstark@46: 	// copy cumscore to first part of fcumscore
adamstark@58: 	for (int i = 0;i < onsetDFBufferSize;i++)
adamstark@46: 	{
adamstark@58: 		futureCumulativeScore[i] = cumulativeScore[i];
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// create future window
adamstark@54: 	double v = 1;
adamstark@58: 	for (int i = 0;i < windowSize;i++)
adamstark@46: 	{
adamstark@57: 		w2[i] = exp((-1*pow((v - (beatPeriod/2)),2))   /  (2*pow((beatPeriod/2) ,2)));
adamstark@46: 		v++;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	// create past window
adamstark@57: 	v = -2*beatPeriod;
adamstark@58: 	int start = onsetDFBufferSize - round(2*beatPeriod);
adamstark@58: 	int end = onsetDFBufferSize - round(beatPeriod/2);
adamstark@46: 	int pastwinsize = end-start+1;
adamstark@54: 	double w1[pastwinsize];
adamstark@46: 
adamstark@46: 	for (int i = 0;i < pastwinsize;i++)
adamstark@46: 	{
adamstark@57: 		w1[i] = exp((-1*pow(tightness*log(-v/beatPeriod),2))/2);
adamstark@46: 		v = v+1;
adamstark@46: 	}
adamstark@46: 
adamstark@46: 	
adamstark@46: 
adamstark@46: 	// calculate future cumulative score
adamstark@54: 	double max;
adamstark@46: 	int n;
adamstark@54: 	double wcumscore;
adamstark@58: 	for (int i = onsetDFBufferSize;i < (onsetDFBufferSize+windowSize);i++)
adamstark@46: 	{
adamstark@57: 		start = i - round(2*beatPeriod);
adamstark@57: 		end = i - round(beatPeriod/2);
adamstark@46: 		
adamstark@46: 		max = 0;
adamstark@46: 		n = 0;
adamstark@46: 		for (int k=start;k <= end;k++)
adamstark@46: 		{
adamstark@58: 			wcumscore = futureCumulativeScore[k]*w1[n];
adamstark@46: 			
adamstark@46: 			if (wcumscore > max)
adamstark@46: 			{
adamstark@46: 				max = wcumscore;
adamstark@46: 			}
adamstark@46: 			n++;
adamstark@46: 		}
adamstark@46: 		
adamstark@58: 		futureCumulativeScore[i] = max;
adamstark@46: 	}
adamstark@46: 	
adamstark@46: 	
adamstark@46: 	// predict beat
adamstark@46: 	max = 0;
adamstark@46: 	n = 0;
adamstark@46: 	
adamstark@58: 	for (int i = onsetDFBufferSize;i < (onsetDFBufferSize+windowSize);i++)
adamstark@46: 	{
adamstark@58: 		wcumscore = futureCumulativeScore[i]*w2[n];
adamstark@46: 		
adamstark@46: 		if (wcumscore > max)
adamstark@46: 		{
adamstark@46: 			max = wcumscore;
adamstark@58: 			beatCounter = n;
adamstark@46: 		}	
adamstark@46: 		
adamstark@46: 		n++;
adamstark@46: 	}
adamstark@46: 		
adamstark@46: 	// set next prediction time
adamstark@58: 	m0 = beatCounter+round(beatPeriod/2);
adamstark@46: 	
adamstark@46: 
adamstark@46: }