diff src/BTrack.cpp @ 96:c58f01834337

Merge branch 'release/1.0.4'

--- a/src/BTrack.cpp	Sun Jan 10 11:36:52 2016 +0000
+++ b/src/BTrack.cpp	Sat Jun 18 10:50:06 2016 +0100
@@ -23,27 +23,50 @@
 #include <algorithm>
 #include "BTrack.h"
 #include "samplerate.h"
+#include <iostream>
 
 //=======================================================================
-BTrack::BTrack() : odf(512,1024,ComplexSpectralDifferenceHWR,HanningWindow)
+BTrack::BTrack()
+ :  odf (512, 1024, ComplexSpectralDifferenceHWR, HanningWindow)
 {
-    initialise(512, 1024);
+    initialise (512, 1024);
 }
 
 //=======================================================================
-BTrack::BTrack(int hopSize_) : odf(hopSize_,2*hopSize_,ComplexSpectralDifferenceHWR,HanningWindow)
+BTrack::BTrack (int hopSize_)
+ :  odf(hopSize_, 2*hopSize_, ComplexSpectralDifferenceHWR, HanningWindow)
 {	
-    initialise(hopSize_, 2*hopSize_);
+    initialise (hopSize_, 2*hopSize_);
 }
 
 //=======================================================================
-BTrack::BTrack(int hopSize_,int frameSize_) : odf(hopSize_,frameSize_,ComplexSpectralDifferenceHWR,HanningWindow)
+BTrack::BTrack (int hopSize_, int frameSize_)
+ : odf (hopSize_, frameSize_, ComplexSpectralDifferenceHWR, HanningWindow)
 {
-    initialise(hopSize_, frameSize_);
+    initialise (hopSize_, frameSize_);
 }
 
 //=======================================================================
-double BTrack::getBeatTimeInSeconds(long frameNumber,int hopSize,int fs)
+BTrack::~BTrack()
+{
+#ifdef USE_FFTW
+    // destroy fft plan
+    fftw_destroy_plan (acfForwardFFT);
+    fftw_destroy_plan (acfBackwardFFT);
+    fftw_free (complexIn);
+    fftw_free (complexOut);
+#endif
+    
+#ifdef USE_KISS_FFT
+    free (cfgForwards);
+    free (cfgBackwards);
+    delete [] fftIn;
+    delete [] fftOut;
+#endif
+}
+
+//=======================================================================
+double BTrack::getBeatTimeInSeconds (long frameNumber, int hopSize, int fs)
 {
     double hop = (double) hopSize;
     double samplingFrequency = (double) fs;
@@ -53,17 +76,17 @@
 }
 
 //=======================================================================
-double BTrack::getBeatTimeInSeconds(int frameNumber,int hopSize,int fs)
+double BTrack::getBeatTimeInSeconds (int frameNumber, int hopSize, int fs)
 {
     long frameNum = (long) frameNumber;
     
-    return getBeatTimeInSeconds(frameNum, hopSize, fs);
+    return getBeatTimeInSeconds (frameNum, hopSize, fs);
 }
 
 
 
 //=======================================================================
-void BTrack::initialise(int hopSize_, int frameSize_)
+void BTrack::initialise (int hopSize_, int frameSize_)
 {
     double rayparam = 43;
 	double pi = 3.14159265;
@@ -83,13 +106,13 @@
 	
 
 	// create rayleigh weighting vector
-	for (int n = 0;n < 128;n++)
+	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++)
+	for (int i = 0; i < 41; i++)
 	{
 		prevDelta[i] = 1;
 	}
@@ -118,10 +141,29 @@
     
     // initialise algorithm given the hopsize
     setHopSize(hopSize_);
+    
+    
+    // Set up FFT for calculating the auto-correlation function
+    FFTLengthForACFCalculation = 1024;
+    
+#ifdef USE_FFTW
+    complexIn = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * FFTLengthForACFCalculation);		// complex array to hold fft data
+    complexOut = (fftw_complex*) fftw_malloc (sizeof(fftw_complex) * FFTLengthForACFCalculation);	// complex array to hold fft data
+    
+    acfForwardFFT = fftw_plan_dft_1d (FFTLengthForACFCalculation, complexIn, complexOut, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
+    acfBackwardFFT = fftw_plan_dft_1d (FFTLengthForACFCalculation, complexOut, complexIn, FFTW_BACKWARD, FFTW_ESTIMATE);	// FFT plan initialisation
+#endif
+    
+#ifdef USE_KISS_FFT
+    fftIn = new kiss_fft_cpx[FFTLengthForACFCalculation];
+    fftOut = new kiss_fft_cpx[FFTLengthForACFCalculation];
+    cfgForwards = kiss_fft_alloc (FFTLengthForACFCalculation, 0, 0, 0);
+    cfgBackwards = kiss_fft_alloc (FFTLengthForACFCalculation, 1, 0, 0);
+#endif
 }
 
 //=======================================================================
-void BTrack::setHopSize(int hopSize_)
+void BTrack::setHopSize (int hopSize_)
 {	
 	hopSize = hopSize_;
 	onsetDFBufferSize = (512*512)/hopSize;		// calculate df buffer size
@@ -129,18 +171,17 @@
 	beatPeriod = round(60/((((double) hopSize)/44100)*tempo));
 
     // set size of onset detection function buffer
-    onsetDF.resize(onsetDFBufferSize);
+    onsetDF.resize (onsetDFBufferSize);
     
     // set size of cumulative score buffer
-    cumulativeScore.resize(onsetDFBufferSize);
+    cumulativeScore.resize (onsetDFBufferSize);
 	
 	// initialise df_buffer to zeros
-	for (int i = 0;i < onsetDFBufferSize;i++)
+	for (int i = 0; i < onsetDFBufferSize; i++)
 	{
 		onsetDF[i] = 0;
 		cumulativeScore[i] = 0;
 		
-		
 		if ((i %  ((int) round(beatPeriod))) == 0)
 		{
 			onsetDF[i] = 1;
@@ -149,13 +190,13 @@
 }
 
 //=======================================================================
-void BTrack::updateHopAndFrameSize(int hopSize_,int frameSize_)
+void BTrack::updateHopAndFrameSize (int hopSize_, int frameSize_)
 {
     // update the onset detection function object
-    odf.initialise(hopSize_, frameSize_);
+    odf.initialise (hopSize_, frameSize_);
     
     // update the hop size being used by the beat tracker
-    setHopSize(hopSize_);
+    setHopSize (hopSize_);
 }
 
 //=======================================================================
@@ -183,23 +224,21 @@
 }
 
 //=======================================================================
-void BTrack::processAudioFrame(double *frame)
+void BTrack::processAudioFrame (double* frame)
 {
     // calculate the onset detection function sample for the frame
-    double sample = odf.calculateOnsetDetectionFunctionSample(frame);
-    
-    
+    double sample = odf.calculateOnsetDetectionFunctionSample (frame);
     
     // process the new onset detection function sample in the beat tracking algorithm
-    processOnsetDetectionFunctionSample(sample);
+    processOnsetDetectionFunctionSample (sample);
 }
 
 //=======================================================================
-void BTrack::processOnsetDetectionFunctionSample(double newSample)
+void BTrack::processOnsetDetectionFunctionSample (double newSample)
 {
     // we need to ensure that the onset
     // detection function sample is positive
-    newSample = fabs(newSample);
+    newSample = fabs (newSample);
     
     // add a tiny constant to the sample to stop it from ever going
     // to zero. this is to avoid problems further down the line
@@ -208,18 +247,12 @@
 	m0--;
 	beatCounter--;
 	beatDueInFrame = false;
-	
-	// move all samples back one step
-	for (int i=0;i < (onsetDFBufferSize-1);i++)
-	{
-		onsetDF[i] = onsetDF[i+1];
-	}
-	
+		
 	// add new sample at the end
-	onsetDF[onsetDFBufferSize-1] = newSample;
+    onsetDF.addSampleToEnd (newSample);
 	
 	// update cumulative score
-	updateCumulativeScore(newSample);
+	updateCumulativeScore (newSample);
 	
 	// if we are halfway between beats
 	if (m0 == 0)
@@ -239,7 +272,7 @@
 }
 
 //=======================================================================
-void BTrack::setTempo(double tempo)
+void BTrack::setTempo (double tempo)
 {	 
 	
 	/////////// TEMPO INDICATION RESET //////////////////
@@ -306,7 +339,7 @@
 }
 
 //=======================================================================
-void BTrack::fixTempo(double tempo)
+void BTrack::fixTempo (double tempo)
 {	
 	// firstly make sure tempo is between 80 and 160 bpm..
 	while (tempo > 160)
@@ -346,52 +379,52 @@
 void BTrack::resampleOnsetDetectionFunction()
 {
 	float output[512];
+    
     float input[onsetDFBufferSize];
     
     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;
-	
-	src_data.data_in = input;
-	src_data.input_frames = BUFFER_LEN;
-	
-	src_data.src_ratio = src_ratio;
-	
-	src_data.data_out = output;
-	src_data.output_frames = output_len;
-	
-	src_simple (&src_data, SRC_SINC_BEST_QUALITY, 1);
-			
-	for (int i = 0;i < output_len;i++)
-	{
-		resampledOnsetDF[i] = (double) src_data.data_out[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;
+    
+    src_data.data_in = input;
+    src_data.input_frames = BUFFER_LEN;
+    
+    src_data.src_ratio = src_ratio;
+    
+    src_data.data_out = output;
+    src_data.output_frames = output_len;
+    
+    src_simple (&src_data, SRC_SINC_BEST_QUALITY, 1);
+            
+    for (int i = 0;i < output_len;i++)
+    {
+        resampledOnsetDF[i] = (double) src_data.data_out[i];
+    }
 }
 
 //=======================================================================
 void BTrack::calculateTempo()
 {
 	// adaptive threshold on input
-	adaptiveThreshold(resampledOnsetDF,512);
+	adaptiveThreshold (resampledOnsetDF,512);
 		
 	// calculate auto-correlation function of detection function
-	calculateBalancedACF(resampledOnsetDF);
+	calculateBalancedACF (resampledOnsetDF);
 	
 	// calculate output of comb filterbank
 	calculateOutputOfCombFilterBank();
 	
-	
 	// adaptive threshold on rcf
-	adaptiveThreshold(combFilterBankOutput,128);
+	adaptiveThreshold (combFilterBankOutput,128);
 
 	
 	int t_index;
@@ -399,8 +432,8 @@
 	// 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)));
+		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];
@@ -425,7 +458,7 @@
 		maxval = -1;
 		for (int i = 0;i < 41;i++)
 		{
-			curval = prevDelta[i]*tempoTransitionMatrix[i][j];
+			curval = prevDelta[i] * tempoTransitionMatrix[i][j];
 			
 			if (curval > maxval)
 			{
@@ -433,7 +466,7 @@
 			}
 		}
 		
-		delta[j] = maxval*tempoObservationVector[j];
+		delta[j] = maxval * tempoObservationVector[j];
 	}
 	
 
@@ -453,16 +486,16 @@
 		prevDelta[j] = delta[j];
 	}
 	
-	beatPeriod = round((60.0*44100.0)/(((2*maxind)+80)*((double) hopSize)));
+	beatPeriod = round ((60.0*44100.0)/(((2*maxind)+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 (double* x, int N)
 {
 	int i = 0;
 	int k,t = 0;
@@ -476,23 +509,23 @@
 	// 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);
+		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);
+		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);
+		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++)
+	for (i = 0; i < N; i++)
 	{
 		x[i] = x[i] - x_thresh[i];
 		if (x[i] < 0)
@@ -514,11 +547,11 @@
 	
 	numelem = 4;
 	
-	for (int i = 2;i <= 127;i++) // max beat period
+	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 <= numelem; 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
 			}
@@ -527,29 +560,99 @@
 }
 
 //=======================================================================
-void BTrack::calculateBalancedACF(double *onsetDetectionFunction)
+void BTrack::calculateBalancedACF (double* onsetDetectionFunction)
 {
-	int l, n = 0;
-	double sum, tmp;
-	
-	// for l lags from 0-511
-	for (l = 0;l < 512;l++)
-	{
-		sum = 0;	
-		
-		// for n samples from 0 - (512-lag)
-		for (n = 0;n < (512-l);n++)
-		{
-			tmp = onsetDetectionFunction[n] * onsetDetectionFunction[n+l];	// multiply current sample n by sample (n+l)
-			sum = sum + tmp;	// add to sum
-		}
-		
-		acf[l] = sum / (512-l);		// weight by number of mults and add to acf buffer
-	}
+    int onsetDetectionFunctionLength = 512;
+    
+#ifdef USE_FFTW
+    // copy into complex array and zero pad
+    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    {
+        if (i < onsetDetectionFunctionLength)
+        {
+            complexIn[i][0] = onsetDetectionFunction[i];
+            complexIn[i][1] = 0.0;
+        }
+        else
+        {
+            complexIn[i][0] = 0.0;
+            complexIn[i][1] = 0.0;
+        }
+    }
+    
+    // perform the fft
+    fftw_execute (acfForwardFFT);
+    
+    // multiply by complex conjugate
+    for (int i = 0;i < FFTLengthForACFCalculation;i++)
+    {
+        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++)
+    {
+        if (i < onsetDetectionFunctionLength)
+        {
+            fftIn[i].r = onsetDetectionFunction[i];
+            fftIn[i].i = 0.0;
+        }
+        else
+        {
+            fftIn[i].r = 0.0;
+            fftIn[i].i = 0.0;
+        }
+    }
+    
+    // execute kiss fft
+    kiss_fft (cfgForwards, fftIn, fftOut);
+    
+    // multiply by complex conjugate
+    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;
+    }
+    
+    // perform the ifft
+    kiss_fft (cfgBackwards, fftOut, fftIn);
+    
+#endif
+    
+    double lag = 512;
+    
+    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]);
+#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);
+#endif
+        // divide by inverse lad to deal with scale bias towards small lags
+        acf[i] = absValue / lag;
+        
+        // this division by 1024 is technically unnecessary but it ensures the algorithm produces
+        // exactly the same ACF output as the old time domain implementation. The time difference is
+        // minimal so I decided to keep it
+        acf[i] = acf[i] / 1024.;
+        
+        lag = lag - 1.;
+    }
 }
 
 //=======================================================================
-double BTrack::calculateMeanOfArray(double *array,int startIndex,int endIndex)
+double BTrack::calculateMeanOfArray (double* array, int startIndex, int endIndex)
 {
 	int i;
 	double sum = 0;
@@ -557,7 +660,7 @@
     int length = endIndex - startIndex;
 	
 	// find sum
-	for (i = startIndex;i < endIndex;i++)
+	for (i = startIndex; i < endIndex; i++)
 	{
 		sum = sum + array[i];
 	}
@@ -573,11 +676,11 @@
 }
 
 //=======================================================================
-void BTrack::normaliseArray(double *array,int N)
+void BTrack::normaliseArray (double* array, int N)
 {
 	double sum = 0;
 	
-	for (int i = 0;i < N;i++)
+	for (int i = 0; i < N; i++)
 	{
 		if (array[i] > 0)
 		{
@@ -587,7 +690,7 @@
 	
 	if (sum > 0)
 	{
-		for (int i = 0;i < N;i++)
+		for (int i = 0; i < N; i++)
 		{
 			array[i] = array[i] / sum;
 		}
@@ -595,31 +698,30 @@
 }
 
 //=======================================================================
-void BTrack::updateCumulativeScore(double odfSample)
+void BTrack::updateCumulativeScore (double odfSample)
 {	 
 	int start, end, winsize;
 	double max;
 	
-	start = onsetDFBufferSize - round(2*beatPeriod);
-	end = onsetDFBufferSize - round(beatPeriod/2);
+	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++)
+	for (int i = 0; i < winsize; i++)
 	{
-		w1[i] = exp((-1*pow(tightness*log(-v/beatPeriod),2))/2);
+		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++)
+	for (int i=start; i <= end; i++)
 	{
 			wcumscore = cumulativeScore[i]*w1[n];
 		
@@ -630,18 +732,9 @@
 		n++;
 	}
 	
-	
-	// shift cumulative score back one
-	for (int i = 0;i < (onsetDFBufferSize-1);i++)
-	{
-		cumulativeScore[i] = cumulativeScore[i+1];
-	}
-	
-	// add new value to cumulative score
-	cumulativeScore[onsetDFBufferSize-1] = ((1-alpha)*odfSample) + (alpha*max);
-	
-	latestCumulativeScoreValue = cumulativeScore[onsetDFBufferSize-1];
-			
+    latestCumulativeScoreValue = ((1 - alpha) * odfSample) + (alpha * max);
+    
+    cumulativeScore.addSampleToEnd (latestCumulativeScoreValue);
 }
 
 //=======================================================================
@@ -650,6 +743,7 @@
 	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++)
 	{
@@ -658,7 +752,7 @@
 	
 	// create future window
 	double v = 1;
-	for (int i = 0;i < windowSize;i++)
+	for (int i = 0; i < windowSize; i++)
 	{
 		w2[i] = exp((-1*pow((v - (beatPeriod/2)),2))   /  (2*pow((beatPeriod/2) ,2)));
 		v++;
@@ -677,16 +771,14 @@
 		v = v+1;
 	}
 
-	
-
 	// calculate future cumulative score
 	double max;
 	int n;
 	double wcumscore;
-	for (int i = onsetDFBufferSize;i < (onsetDFBufferSize+windowSize);i++)
+	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + windowSize); i++)
 	{
-		start = i - round(2*beatPeriod);
-		end = i - round(beatPeriod/2);
+		start = i - round (2*beatPeriod);
+		end = i - round (beatPeriod/2);
 		
 		max = 0;
 		n = 0;
@@ -704,12 +796,11 @@
 		futureCumulativeScore[i] = max;
 	}
 	
-	
 	// predict beat
 	max = 0;
 	n = 0;
 	
-	for (int i = onsetDFBufferSize;i < (onsetDFBufferSize+windowSize);i++)
+	for (int i = onsetDFBufferSize; i < (onsetDFBufferSize + windowSize); i++)
 	{
 		wcumscore = futureCumulativeScore[i]*w2[n];
 		
@@ -723,7 +814,5 @@
 	}
 		
 	// set next prediction time
-	m0 = beatCounter+round(beatPeriod/2);
-	
-
+	m0 = beatCounter + round (beatPeriod / 2);
 }
\ No newline at end of file