Mercurial > hg > btrack

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/.hgignore	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,1 @@
+python-module/build
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/python-module/btrack_python_module.cpp	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,362 @@
+#include <iostream>
+#include <Python.h>
+#include "../src/OnsetDetectionFunction.h"
+#include "../src/BTrack.h"
+#include <numpy/arrayobject.h>
+
+static PyObject * btrack_onsetdf(PyObject *dummy, PyObject *args)
+{
+    PyObject *arg1=NULL;
+    PyObject *arr1=NULL;
+
+    if (!PyArg_ParseTuple(args, "O", &arg1))
+    {
+        return NULL;
+    }
+
+    arr1 = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (arr1 == NULL)
+    {
+        return NULL;
+    }
+
+
+
+    ////////// GET INPUT DATA ///////////////////
+
+    // get data as array
+    double* data = (double*) PyArray_DATA(arr1);
+
+    // get array size
+    int signal_length = PyArray_Size((PyObject*)arr1);
+    //int k = (int) theSize;
+
+    // get data type
+    char type = PyArray_DESCR(arr1)->type;
+
+    ////////// BEGIN PROCESS ///////////////////
+    int hsize = 512;
+    int fsize = 1024;
+    int df_type = 6;
+    int numframes;
+    double buffer[hsize];	// buffer to hold one hopsize worth of audio samples
+
+
+    // get number of audio frames, given the hop size and signal length
+	numframes = (int) floor(((double) signal_length) / ((double) hsize));
+
+    OnsetDetectionFunction onset(hsize,fsize,df_type,1);
+
+    double df[numframes];
+
+
+
+    ///////////////////////////////////////////
+	//////// Begin Processing Loop ////////////
+
+	for (int i=0;i < numframes;i++)
+	{
+		// add new samples to frame
+		for (int n = 0;n < hsize;n++)
+		{
+			buffer[n] = data[(i*hsize)+n];
+		}
+
+		df[i] = onset.getDFsample(buffer);
+
+	}
+
+	///////// End Processing Loop /////////////
+	///////////////////////////////////////////
+
+
+    ////////// END PROCESS ///////////////////
+
+
+
+    ////////// CREATE ARRAY AND RETURN IT ///////////////////
+    int nd=1;
+    npy_intp m= numframes;
+    //double fArray[5] = {0,1,2,3,4};
+
+    PyObject* c=PyArray_SimpleNew(nd, &m, NPY_DOUBLE);
+
+    void *arr_data = PyArray_DATA((PyArrayObject*)c);
+
+    memcpy(arr_data, df, PyArray_ITEMSIZE((PyArrayObject*) c) * m);
+
+
+    Py_DECREF(arr1);
+    Py_INCREF(Py_None);
+    //return Py_None;
+
+    return (PyObject *)c;
+
+    //return Py_BuildValue("c", type);
+    //return Py_BuildValue("d", sum);
+    //return Py_BuildValue("i", k);
+/*
+fail:
+    Py_XDECREF(arr1);
+    Py_XDECREF(arr2);
+    PyArray_XDECREF_ERR(oarr);
+    return NULL;*/
+}
+
+
+static PyObject * btrack_btrack(PyObject *dummy, PyObject *args)
+{
+    PyObject *arg1=NULL;
+    PyObject *arr1=NULL;
+
+    if (!PyArg_ParseTuple(args, "O", &arg1))
+    {
+        return NULL;
+    }
+
+    arr1 = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (arr1 == NULL)
+    {
+        return NULL;
+    }
+
+
+
+    ////////// GET INPUT DATA ///////////////////
+
+    // get data as array
+    double* data = (double*) PyArray_DATA(arr1);
+
+    // get array size
+    int signal_length = PyArray_Size((PyObject*)arr1);
+    //int k = (int) theSize;
+
+    // get data type
+    char type = PyArray_DESCR(arr1)->type;
+
+    ////////// BEGIN PROCESS ///////////////////
+    int hsize = 512;
+    int fsize = 1024;
+    int df_type = 6;
+    int numframes;
+    double buffer[hsize];	// buffer to hold one hopsize worth of audio samples
+
+
+    // get number of audio frames, given the hop size and signal length
+	numframes = (int) floor(((double) signal_length) / ((double) hsize));
+
+    OnsetDetectionFunction onset(hsize,fsize,df_type,1);
+    BTrack b;
+
+    b.initialise((int) hsize);	// initialise beat tracker
+
+	// set parameters
+    //b.setparams(0.9,5);
+
+    double df[numframes];
+    double beats[5000];
+    int beatnum = 0;
+    float df_val;
+
+    ///////////////////////////////////////////
+	//////// Begin Processing Loop ////////////
+
+	for (int i=0;i < numframes;i++)
+	{
+		// add new samples to frame
+		for (int n = 0;n < hsize;n++)
+		{
+			buffer[n] = data[(i*hsize)+n];
+		}
+
+		df[i] = onset.getDFsample(buffer);
+
+        df_val = (float) (df[i] + 0.0001);
+
+		b.process(df_val);				// process df sample in beat tracker
+
+		if (b.playbeat == 1)
+		{
+			beats[beatnum] = (((double) hsize) / 44100) * ((double) i);
+			beatnum = beatnum + 1;
+		}
+
+	}
+
+	///////// End Processing Loop /////////////
+	///////////////////////////////////////////
+
+
+    ////////// END PROCESS ///////////////////
+
+    double beats_out[beatnum];          // create output array
+
+    // copy beats into output array
+    for (int i = 0;i < beatnum;i++)
+    {
+        beats_out[i] = beats[i];
+    }
+
+
+
+    ////////// CREATE ARRAY AND RETURN IT ///////////////////
+    int nd=1;
+    npy_intp m= beatnum;
+    //double fArray[5] = {0,1,2,3,4};
+
+    PyObject* c=PyArray_SimpleNew(nd, &m, NPY_DOUBLE);
+
+    void *arr_data = PyArray_DATA((PyArrayObject*)c);
+
+    memcpy(arr_data, beats_out, PyArray_ITEMSIZE((PyArrayObject*) c) * m);
+
+
+    Py_DECREF(arr1);
+    Py_INCREF(Py_None);
+    //return Py_None;
+
+    return (PyObject *)c;
+
+    //return Py_BuildValue("c", type);
+    //return Py_BuildValue("d", sum);
+    //return Py_BuildValue("i", k);
+    /*
+     fail:
+     Py_XDECREF(arr1);
+     Py_XDECREF(arr2);
+     PyArray_XDECREF_ERR(oarr);
+     return NULL;*/
+}
+
+static PyObject * btrack_btrack_df(PyObject *dummy, PyObject *args)
+{
+    PyObject *arg1=NULL;
+    PyObject *arr1=NULL;
+
+    if (!PyArg_ParseTuple(args, "O", &arg1))
+    {
+        return NULL;
+    }
+
+    arr1 = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
+    if (arr1 == NULL)
+    {
+        return NULL;
+    }
+
+
+
+    ////////// GET INPUT DATA ///////////////////
+
+    // get data as array
+    double* data = (double*) PyArray_DATA(arr1);
+
+    // get array size
+    int numframes = PyArray_Size((PyObject*)arr1);
+    //int k = (int) theSize;
+
+    // get data type
+    char type = PyArray_DESCR(arr1)->type;
+
+    ////////// BEGIN PROCESS ///////////////////
+    int hsize = 512;
+
+    BTrack b;
+
+    b.initialise((int) hsize);	// initialise beat tracker
+
+	// set parameters
+    //b.setparams(0.9,5);
+
+    double beats[5000];
+    int beatnum = 0;
+    float df_val;
+
+    ///////////////////////////////////////////
+	//////// Begin Processing Loop ////////////
+
+	for (int i=0;i < numframes;i++)
+	{
+        df_val = (float) (data[i] + 0.0001);
+
+		b.process(df_val);				// process df sample in beat tracker
+
+		if (b.playbeat == 1)
+		{
+			beats[beatnum] = (((double) hsize) / 44100) * ((double) i);
+			beatnum = beatnum + 1;
+		}
+
+	}
+
+	///////// End Processing Loop /////////////
+	///////////////////////////////////////////
+
+
+    ////////// END PROCESS ///////////////////
+
+    double beats_out[beatnum];          // create output array
+
+
+    // copy beats into output array
+    for (int i = 0;i < beatnum;i++)
+    {
+        beats_out[i] = beats[i];
+    }
+
+
+    ////////// CREATE ARRAY AND RETURN IT ///////////////////
+    int nd=1;
+    npy_intp m= beatnum;
+    //double fArray[5] = {0,1,2,3,4};
+
+    PyObject* c=PyArray_SimpleNew(nd, &m, NPY_DOUBLE);
+
+    void *arr_data = PyArray_DATA((PyArrayObject*)c);
+
+    memcpy(arr_data, beats_out, PyArray_ITEMSIZE((PyArrayObject*) c) * m);
+
+
+    Py_DECREF(arr1);
+    Py_INCREF(Py_None);
+    //return Py_None;
+
+    return (PyObject *)c;
+
+    //return Py_BuildValue("c", type);
+    //return Py_BuildValue("d", sum);
+    //return Py_BuildValue("i", k);
+    /*
+     fail:
+     Py_XDECREF(arr1);
+     Py_XDECREF(arr2);
+     PyArray_XDECREF_ERR(oarr);
+     return NULL;*/
+}
+
+
+
+static PyMethodDef btrack_methods[] = {
+    { "onsetdf",btrack_onsetdf,METH_VARARGS,"onset detection function"},
+    { "btrack",btrack_btrack,METH_VARARGS,"beat tracker"},
+    { "btrack_df",btrack_btrack_df,METH_VARARGS,"beat tracker with detection function input"},
+    {NULL, NULL, 0, NULL} /* Sentinel */
+};
+
+PyMODINIT_FUNC initbtrack(void)
+{
+    (void)Py_InitModule("btrack", btrack_methods);
+    import_array();
+}
+
+int main(int argc, char *argv[])
+{
+    /* Pass argv[0] to the Python interpreter */
+    Py_SetProgramName(argv[0]);
+
+    /* Initialize the Python interpreter.  Required. */
+    Py_Initialize();
+
+    /* Add a static module */
+    initbtrack();
+}
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/python-module/setup.py	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,16 @@
+# setup.py
+# build command : python setup.py build build_ext --inplace
+from numpy.distutils.core import setup, Extension
+import os, numpy
+
+name = 'btrack'
+sources = ['btrack_python_module.cpp','../src/OnsetDetectionFunction.cpp','../src/BTrack.cpp']
+
+include_dirs = [
+                numpy.get_include(),'/usr/local/include'
+                ]
+
+setup( name = 'BTtrack',
+      include_dirs = include_dirs,
+      ext_modules = [Extension(name, sources,libraries = ['fftw3','samplerate'])]
+      )
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BTrack.cpp	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,665 @@
+//=======================================================================
+/** @file BTrack.cpp
+ *  @brief Implementation file for the BTrack beat tracker
+ *  @author Adam Stark
+ *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+//=======================================================================
+
+#include <iostream>
+#include <cmath>
+#include "BTrack.h"
+#include "samplerate.h"
+using namespace std;
+
+
+
+
+//-------------------------------------------------------------------------------
+// Constructor
+BTrack :: BTrack()
+{
+	float rayparam = 43;
+	float pi = 3.14159265;
+
+
+	// initialise parameters
+	tightness = 5;
+	alpha = 0.9;
+	tempo = 120;
+	est_tempo = 120;
+	p_fact = 60.*44100./512.;
+
+	m0 = 10;
+	beat = -1;
+
+	playbeat = 0;
+
+
+
+
+	// create rayleigh weighting vector
+	for (int n = 0;n < 128;n++)
+	{
+		wv[n] = ((float) n / pow(rayparam,2)) * exp((-1*pow((float)-n,2)) / (2*pow(rayparam,2)));
+	}
+
+	// initialise prev_delta
+	for (int i = 0;i < 41;i++)
+	{
+		prev_delta[i] = 1;
+	}
+
+	float t_mu = 41/2;
+	float m_sig;
+	float x;
+	// create tempo transition matrix
+	m_sig = 41/8;
+	for (int i = 0;i < 41;i++)
+	{
+		for (int j = 0;j < 41;j++)
+		{
+			x = j+1;
+			t_mu = i+1;
+			t_tmat[i][j] = (1 / (m_sig * sqrt(2*pi))) * exp( (-1*pow((x-t_mu),2)) / (2*pow(m_sig,2)) );
+		}
+	}
+
+	// tempo is not fixed
+	tempofix = 0;
+}
+
+//-------------------------------------------------------------------------------
+// Destructor
+BTrack :: ~BTrack()
+{
+
+}
+
+//-------------------------------------------------------------------------------
+// Initialise with frame size and set all frame sizes accordingly
+void BTrack :: initialise(int fsize)
+{
+	framesize = fsize;
+	dfbuffer_size = (512*512)/fsize;		// calculate df buffer size
+
+	bperiod = round(60/((((float) fsize)/44100)*tempo));
+
+	dfbuffer = new float[dfbuffer_size];	// create df_buffer
+	cumscore = new float[dfbuffer_size];	// create cumscore
+
+
+	// initialise df_buffer to zeros
+	for (int i = 0;i < dfbuffer_size;i++)
+	{
+		dfbuffer[i] = 0;
+		cumscore[i] = 0;
+
+
+		if ((i %  ((int) round(bperiod))) == 0)
+		{
+			dfbuffer[i] = 1;
+		}
+	}
+}
+
+//-------------------------------------------------------------------------------
+// Add new sample to buffer and apply beat tracking
+void BTrack :: process(float df_sample)
+{
+	m0--;
+	beat--;
+	playbeat = 0;
+
+	// move all samples back one step
+	for (int i=0;i < (dfbuffer_size-1);i++)
+	{
+		dfbuffer[i] = dfbuffer[i+1];
+	}
+
+	// add new sample at the end
+	dfbuffer[dfbuffer_size-1] = df_sample;
+
+	// update cumulative score
+	updatecumscore(df_sample);
+
+	// if we are halfway between beats
+	if (m0 == 0)
+	{
+		predictbeat();
+	}
+
+	// if we are at a beat
+	if (beat == 0)
+	{
+		playbeat = 1;	// indicate a beat should be output
+
+		// recalculate the tempo
+		dfconvert();
+		calcTempo();
+	}
+}
+
+//-------------------------------------------------------------------------------
+// Set the tempo of the beat tracker
+void BTrack :: settempo(float tempo)
+{
+
+	/////////// TEMPO INDICATION RESET //////////////////
+
+	// firstly make sure tempo is between 80 and 160 bpm..
+	while (tempo > 160)
+	{
+		tempo = tempo/2;
+	}
+
+	while (tempo < 80)
+	{
+		tempo = tempo * 2;
+	}
+
+	// convert tempo from bpm value to integer index of tempo probability
+	int tempo_index = (int) round((tempo - 80)/2);
+
+	// now set previous tempo observations to zero
+	for (int i=0;i < 41;i++)
+	{
+		prev_delta[i] = 0;
+	}
+
+	// set desired tempo index to 1
+	prev_delta[tempo_index] = 1;
+
+
+	/////////// CUMULATIVE SCORE ARTIFICAL TEMPO UPDATE //////////////////
+
+	// calculate new beat period
+	int new_bperiod = (int) round(60/((((float) framesize)/44100)*tempo));
+
+	int bcounter = 1;
+	// initialise df_buffer to zeros
+	for (int i = (dfbuffer_size-1);i >= 0;i--)
+	{
+		if (bcounter == 1)
+		{
+			cumscore[i] = 150;
+			dfbuffer[i] = 150;
+		}
+		else
+		{
+			cumscore[i] = 10;
+			dfbuffer[i] = 10;
+		}
+
+		bcounter++;
+
+		if (bcounter > new_bperiod)
+		{
+			bcounter = 1;
+		}
+	}
+
+	/////////// INDICATE THAT THIS IS A BEAT //////////////////
+
+	// beat is now
+	beat = 0;
+
+	// offbeat is half of new beat period away
+	m0 = (int) round(((float) new_bperiod)/2);
+}
+
+
+//-------------------------------------------------------------------------------
+// fix tempo to roughly around some value
+void BTrack :: fixtempo(float tempo)
+{
+	// firstly make sure tempo is between 80 and 160 bpm..
+	while (tempo > 160)
+	{
+		tempo = tempo/2;
+	}
+
+	while (tempo < 80)
+	{
+		tempo = tempo * 2;
+	}
+
+	// convert tempo from bpm value to integer index of tempo probability
+	int tempo_index = (int) round((tempo - 80)/2);
+
+	// now set previous fixed previous tempo observation values to zero
+	for (int i=0;i < 41;i++)
+	{
+		prev_delta_fix[i] = 0;
+	}
+
+	// set desired tempo index to 1
+	prev_delta_fix[tempo_index] = 1;
+
+	// set the tempo fix flag
+	tempofix = 1;
+}
+
+//-------------------------------------------------------------------------------
+// do not fix the tempo anymore
+void BTrack :: unfixtempo()
+{
+	// set the tempo fix flag
+	tempofix = 0;
+}
+
+//-------------------------------------------------------------------------------
+// Convert detection function from N samples to 512
+void BTrack :: dfconvert()
+{
+	float output[512];
+
+	double src_ratio = 512.0/((double) dfbuffer_size);
+	int BUFFER_LEN = dfbuffer_size;
+	int output_len;
+	SRC_DATA	src_data ;
+
+	//output_len = (int) floor (((double) BUFFER_LEN) * src_ratio) ;
+	output_len = 512;
+
+	src_data.data_in = dfbuffer;
+	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++)
+	{
+		df512[i] = src_data.data_out[i];
+	}
+}
+
+//-------------------------------------------------------------------------------
+// To calculate the current tempo expressed as the beat period in detection function samples
+void BTrack :: calcTempo()
+{
+	// adaptive threshold on input
+	adapt_thresh(df512,512);
+
+	// calculate auto-correlation function of detection function
+	acf_bal(df512);
+
+	// calculate output of comb filterbank
+	getrcfoutput();
+
+
+	// adaptive threshold on rcf
+	adapt_thresh(rcf,128);
+
+
+	int t_index;
+	int t_index2;
+	// calculate tempo observation vector from bperiod observation vector
+	for (int i = 0;i < 41;i++)
+	{
+		t_index = (int) round(p_fact / ((float) ((2*i)+80)));
+		t_index2 = (int) round(p_fact / ((float) ((4*i)+160)));
+
+
+		t_obs[i] = rcf[t_index-1] + rcf[t_index2-1];
+	}
+
+
+	float maxval;
+	float maxind;
+	float curval;
+
+	// if tempo is fixed then always use a fixed set of tempi as the previous observation probability function
+	if (tempofix == 1)
+	{
+		for (int k = 0;k < 41;k++)
+		{
+			prev_delta[k] = prev_delta_fix[k];
+		}
+	}
+
+	for (int j=0;j < 41;j++)
+	{
+		maxval = -1;
+		for (int i = 0;i < 41;i++)
+		{
+			curval = prev_delta[i]*t_tmat[i][j];
+
+			if (curval > maxval)
+			{
+				maxval = curval;
+			}
+		}
+
+		delta[j] = maxval*t_obs[j];
+	}
+
+
+	normalise(delta,41);
+
+	maxind = -1;
+	maxval = -1;
+
+	for (int j=0;j < 41;j++)
+	{
+		if (delta[j] > maxval)
+		{
+			maxval = delta[j];
+			maxind = j;
+		}
+
+		prev_delta[j] = delta[j];
+	}
+
+	bperiod = round((60.0*44100.0)/(((2*maxind)+80)*((float) framesize)));
+
+	if (bperiod > 0)
+	{
+		est_tempo = 60.0/((((float) framesize) / 44100.0)*bperiod);
+	}
+
+	//cout << bperiod << endl;
+}
+
+//-------------------------------------------------------------------------------
+// calculates an adaptive threshold which is used to remove low level energy from detection function and emphasise peaks
+void BTrack :: adapt_thresh(float x[],int N)
+{
+	//int N = 512; // length of df
+	int i = 0;
+	int k,t = 0;
+	float x_thresh[N];
+
+	int p_post = 7;
+	int p_pre = 8;
+
+	t = min(N,p_post);	// what is smaller, p_post of df size. This is to avoid accessing outside of arrays
+
+	// find threshold for first 't' samples, where a full average cannot be computed yet
+	for (i = 0;i <= t;i++)
+	{
+		k = min((i+p_pre),N);
+		x_thresh[i] = mean_array(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] = mean_array(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 = max((i-p_post),1);
+		x_thresh[i] = mean_array(x,k,N);
+	}
+
+	// subtract the threshold from the detection function and check that it is not less than 0
+	for (i = 0;i < N;i++)
+	{
+		x[i] = x[i] - x_thresh[i];
+		if (x[i] < 0)
+		{
+			x[i] = 0;
+		}
+	}
+}
+
+//-------------------------------------------------------------------------------
+// returns the output of the comb filter
+void BTrack :: getrcfoutput()
+{
+	int numelem;
+
+	for (int i = 0;i < 128;i++)
+	{
+		rcf[i] = 0;
+	}
+
+	numelem = 4;
+
+	for (int i = 2;i <= 127;i++) // max beat period
+	{
+		for (int a = 1;a <= numelem;a++) // number of comb elements
+		{
+			for (int b = 1-a;b <= a-1;b++) // general state using normalisation of comb elements
+			{
+				rcf[i-1] = rcf[i-1] + (acf[(a*i+b)-1]*wv[i-1])/(2*a-1);	// calculate value for comb filter row
+			}
+		}
+	}
+}
+
+//-------------------------------------------------------------------------------
+// calculates the balanced autocorrelation of the smoothed detection function
+void BTrack :: acf_bal(float df_thresh[])
+{
+	int l, n = 0;
+	float 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 = df_thresh[n] * df_thresh[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
+	}
+}
+
+
+//-------------------------------------------------------------------------------
+// calculates the mean of values in an array from index locations [start,end]
+float BTrack :: mean_array(float array[],int start,int end)
+{
+	int i;
+	float sum = 0;
+	int length = end - start + 1;
+
+	// find sum
+	for (i = start;i < end+1;i++)
+	{
+		sum = sum + array[i];
+	}
+
+	return sum / length;	// average and return
+}
+
+//-------------------------------------------------------------------------------
+// normalise the array
+void BTrack :: normalise(float array[],int N)
+{
+	double sum = 0;
+
+	for (int i = 0;i < N;i++)
+	{
+		if (array[i] > 0)
+		{
+			sum = sum + array[i];
+		}
+	}
+
+	if (sum > 0)
+	{
+		for (int i = 0;i < N;i++)
+		{
+			array[i] = array[i] / sum;
+		}
+	}
+}
+
+//-------------------------------------------------------------------------------
+// plot contents of detection function buffer
+void BTrack :: plotdfbuffer()
+{
+	for (int i=0;i < dfbuffer_size;i++)
+	{
+		cout << dfbuffer[i] << endl;
+	}
+
+	cout << "--------------------------------" << endl;
+}
+
+//-------------------------------------------------------------------------------
+// update the cumulative score
+void BTrack :: updatecumscore(float df_sample)
+{
+	int start, end, winsize;
+	float max;
+
+	start = dfbuffer_size - round(2*bperiod);
+	end = dfbuffer_size - round(bperiod/2);
+	winsize = end-start+1;
+
+	float w1[winsize];
+	float v = -2*bperiod;
+	float wcumscore;
+
+
+	// create window
+	for (int i = 0;i < winsize;i++)
+	{
+		w1[i] = exp((-1*pow(tightness*log(-v/bperiod),2))/2);
+		v = v+1;
+	}
+
+	// calculate new cumulative score value
+	max = 0;
+	int n = 0;
+	for (int i=start;i <= end;i++)
+	{
+			wcumscore = cumscore[i]*w1[n];
+
+			if (wcumscore > max)
+			{
+				max = wcumscore;
+			}
+		n++;
+	}
+
+
+	// shift cumulative score back one
+	for (int i = 0;i < (dfbuffer_size-1);i++)
+	{
+		cumscore[i] = cumscore[i+1];
+	}
+
+	// add new value to cumulative score
+	cumscore[dfbuffer_size-1] = ((1-alpha)*df_sample) + (alpha*max);
+
+	cscoreval = cumscore[dfbuffer_size-1];
+
+	//cout << cumscore[dfbuffer_size-1] << endl;
+
+}
+
+//-------------------------------------------------------------------------------
+// plot contents of detection function buffer
+void BTrack :: predictbeat()
+{
+	int winsize = (int) bperiod;
+	float fcumscore[dfbuffer_size + winsize];
+	float w2[winsize];
+	// copy cumscore to first part of fcumscore
+	for (int i = 0;i < dfbuffer_size;i++)
+	{
+		fcumscore[i] = cumscore[i];
+	}
+
+	// create future window
+	float v = 1;
+	for (int i = 0;i < winsize;i++)
+	{
+		w2[i] = exp((-1*pow((v - (bperiod/2)),2))   /  (2*pow((bperiod/2) ,2)));
+		v++;
+	}
+
+	// create past window
+	v = -2*bperiod;
+	int start = dfbuffer_size - round(2*bperiod);
+	int end = dfbuffer_size - round(bperiod/2);
+	int pastwinsize = end-start+1;
+	float w1[pastwinsize];
+
+	for (int i = 0;i < pastwinsize;i++)
+	{
+		w1[i] = exp((-1*pow(tightness*log(-v/bperiod),2))/2);
+		v = v+1;
+	}
+
+
+
+	// calculate future cumulative score
+	float max;
+	int n;
+	float wcumscore;
+	for (int i = dfbuffer_size;i < (dfbuffer_size+winsize);i++)
+	{
+		start = i - round(2*bperiod);
+		end = i - round(bperiod/2);
+
+		max = 0;
+		n = 0;
+		for (int k=start;k <= end;k++)
+		{
+			wcumscore = fcumscore[k]*w1[n];
+
+			if (wcumscore > max)
+			{
+				max = wcumscore;
+			}
+			n++;
+		}
+
+		fcumscore[i] = max;
+	}
+
+
+	// predict beat
+	max = 0;
+	n = 0;
+
+	for (int i = dfbuffer_size;i < (dfbuffer_size+winsize);i++)
+	{
+		wcumscore = fcumscore[i]*w2[n];
+
+		if (wcumscore > max)
+		{
+			max = wcumscore;
+			beat = n;
+		}
+
+		n++;
+	}
+
+
+	// set beat
+	beat = beat;
+
+	// set next prediction time
+	m0 = beat+round(bperiod/2);
+
+
+}
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/BTrack.h	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,99 @@
+//=======================================================================
+/** @file BTrack.h
+ *  @brief Header file for the BTrack beat tracker
+ *  @author Adam Stark
+ *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+//=======================================================================
+
+#ifndef __BTRACK_H
+#define __BTRACK_H
+
+//#include "fftw3.h"
+
+class BTrack {
+
+public:
+	BTrack();				// constructor
+	~BTrack();				// destructor
+
+	void initialise(int fsize);
+	void process(float df_sample);
+	void plotdfbuffer();
+	void updatecumscore(float df_sample);
+	void predictbeat();
+	void dfconvert();
+	void calcTempo();
+	void adapt_thresh(float x[],int N);
+	float mean_array(float array[],int start,int end);
+	void normalise(float array[],int N);
+	void acf_bal(float df_thresh[]);
+	void getrcfoutput();
+	void settempo(float tempo);
+	void fixtempo(float tempo);
+	void unfixtempo();
+
+	int playbeat;
+	float cscoreval;
+	float est_tempo;
+
+private:
+
+	// buffers
+	float *dfbuffer;			// to hold detection function
+	float df512[512];			// to hold resampled detection function
+	float *cumscore;			// to hold cumulative score
+
+	float acf[512];				// to hold autocorrelation function
+
+	float wv[128];				// to hold weighting vector
+
+	float rcf[128];				// to hold comb filter output
+	float t_obs[41];			// to hold tempo version of comb filter output
+
+	float delta[41];			// to hold final tempo candidate array
+	float prev_delta[41];		// previous delta
+	float prev_delta_fix[41];	// fixed tempo version of previous delta
+
+	float t_tmat[41][41];		// transition matrix
+
+
+	// parameters
+	float tightness;
+	float alpha;
+	float bperiod;
+	float tempo;
+
+
+	float p_fact;
+
+
+	//
+	int m0;				// indicates when the next point to predict the next beat is
+	int beat;
+
+	int dfbuffer_size;
+
+
+	int framesize;
+
+
+	int tempofix;
+
+
+};
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/OnsetDetectionFunction.cpp	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,805 @@
+//=======================================================================
+/** @file OnsetDetectionFunction.cpp
+ *  @brief A class for calculating onset detection functions
+ *  @author Adam Stark
+ *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+//=======================================================================
+
+#include <math.h>
+#include <iostream>
+#include "OnsetDetectionFunction.h"
+using namespace std;
+
+//-------------------------------------------------------------------------------
+// Constructor
+OnsetDetectionFunction :: OnsetDetectionFunction()
+{
+	// indicate that we have not initialised yet
+	initialised = 0;
+
+	// set pi
+	pi = 3.14159265358979;
+
+	// initialise with hopsize = 512, framesize = 1024, complex spectral difference DF and hanning window
+	initialise(512,1024,6,1);
+}
+
+
+//-------------------------------------------------------------------------------
+// Constructor (with arguments)
+OnsetDetectionFunction :: OnsetDetectionFunction(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
+{
+	// indicate that we have not initialised yet
+	initialised = 0;
+
+	// set pi
+	pi = 3.14159265358979;
+
+	// initialise with arguments to constructor
+	initialise(arg_hsize,arg_fsize,arg_df_type,arg_win_type);
+}
+
+
+//--------------------------------------------------------------------------------------
+// Destructor
+OnsetDetectionFunction :: ~OnsetDetectionFunction()
+{
+	// destroy fft plan
+    fftw_destroy_plan(p);
+	fftw_free(in);
+	fftw_free(out);
+
+	// deallocate memory
+	delete [] frame;
+	frame = NULL;
+	delete [] window;
+	window = NULL;
+	delete [] wframe;
+	wframe = NULL;
+	delete [] mag;
+	mag = NULL;
+	delete [] mag_old;
+	mag_old = NULL;
+	delete [] phase;
+	phase = NULL;
+	delete [] phase_old;
+	phase_old = NULL;
+	delete [] phase_old_2;
+	phase_old_2 = NULL;
+}
+
+//-------------------------------------------------------------------------------
+// Initialisation
+void OnsetDetectionFunction :: initialise(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type)
+{
+	if (initialised == 1) // if we have already initialised some buffers and an FFT plan
+	{
+		//////////////////////////////////
+		// TIDY UP FIRST - If initialise is called after the class has been initialised
+		// then we want to free up memory and cancel existing FFT plans
+
+		// destroy fft plan
+		fftw_destroy_plan(p);
+		fftw_free(in);
+		fftw_free(out);
+
+
+		// deallocate memory
+		delete [] frame;
+		frame = NULL;
+		delete [] window;
+		window = NULL;
+		delete [] wframe;
+		wframe = NULL;
+		delete [] mag;
+		mag = NULL;
+		delete [] mag_old;
+		mag_old = NULL;
+		delete [] phase;
+		phase = NULL;
+		delete [] phase_old;
+		phase_old = NULL;
+		delete [] phase_old_2;
+		phase_old_2 = NULL;
+
+		////// END TIDY UP ///////////////
+		//////////////////////////////////
+	}
+
+	hopsize = arg_hsize; // set hopsize
+	framesize = arg_fsize; // set framesize
+
+	df_type = arg_df_type; // set detection function type
+
+	// initialise buffers
+	frame = new double[framesize];
+	window = new double[framesize];
+	wframe = new double[framesize];
+
+	mag = new double[framesize];
+	mag_old = new double[framesize];
+
+	phase = new double[framesize];
+	phase_old = new double[framesize];
+	phase_old_2 = new double[framesize];
+
+
+	// set the window to the specified type
+	switch (arg_win_type){
+		case 0:
+			set_win_rectangular();		// Rectangular window
+			break;
+		case 1:
+			set_win_hanning();			// Hanning Window
+			break;
+		case 2:
+			set_win_hamming();			// Hamming Window
+			break;
+		case 3:
+			set_win_blackman();			// Blackman Window
+			break;
+		case 4:
+			set_win_tukey();			// Tukey Window
+			break;
+		default:
+			set_win_hanning();			// DEFAULT: Hanning Window
+	}
+
+
+
+
+	// initialise previous magnitude spectrum to zero
+	for (int i = 0;i < framesize;i++)
+	{
+		mag_old[i] = 0.0;
+		phase_old[i] = 0.0;
+		phase_old_2[i] = 0.0;
+		frame[i] = 0.0;
+	}
+
+	energy_sum_old = 0.0;	// initialise previous energy sum value to zero
+
+	/*  Init fft */
+	in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);		// complex array to hold fft data
+	out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * framesize);	// complex array to hold fft data
+	p = fftw_plan_dft_1d(framesize, in, out, FFTW_FORWARD, FFTW_ESTIMATE);	// FFT plan initialisation
+
+	initialised = 1;
+}
+
+//--------------------------------------------------------------------------------------
+// set the detection function type to that specified by the argument
+void OnsetDetectionFunction :: set_df_type(int arg_df_type)
+{
+	df_type = arg_df_type; // set detection function type
+}
+
+
+//--------------------------------------------------------------------------------------
+// calculates a single detection function sample from a single audio frame.
+double OnsetDetectionFunction :: getDFsample(double inputbuffer[])
+{
+	double df_sample;
+
+	// shift audio samples back in frame by hop size
+	for (int i = 0; i < (framesize-hopsize);i++)
+	{
+		frame[i] = frame[i+hopsize];
+	}
+
+	// add new samples to frame from input buffer
+	int j = 0;
+	for (int i = (framesize-hopsize);i < framesize;i++)
+	{
+		frame[i] = inputbuffer[j];
+		j++;
+	}
+
+	switch (df_type){
+		case 0:
+			df_sample = energy_envelope();	// calculate energy envelope detection function sample
+			break;
+		case 1:
+			df_sample = energy_difference();	// calculate half-wave rectified energy difference detection function sample
+			break;
+		case 2:
+			df_sample = spectral_difference();	// calculate spectral difference detection function sample
+			break;
+		case 3:
+			df_sample = spectral_difference_hwr(); // calculate spectral difference detection function sample (half wave rectified)
+			break;
+		case 4:
+			df_sample = phase_deviation();		// calculate phase deviation detection function sample (half wave rectified)
+			break;
+		case 5:
+			df_sample = complex_spectral_difference(); // calcualte complex spectral difference detection function sample
+			break;
+		case 6:
+			df_sample = complex_spectral_difference_hwr();  // calcualte complex spectral difference detection function sample (half-wave rectified)
+			break;
+		case 7:
+			df_sample = high_frequency_content(); // calculate high frequency content detection function sample
+			break;
+		case 8:
+			df_sample = high_frequency_spectral_difference(); // calculate high frequency spectral difference detection function sample
+			break;
+		case 9:
+			df_sample = high_frequency_spectral_difference_hwr(); // calculate high frequency spectral difference detection function (half-wave rectified)
+			break;
+		default:
+			df_sample = 1.0;
+	}
+
+	return df_sample;
+}
+
+
+//--------------------------------------------------------------------------------------
+// performs the fft, storing the complex result in 'out'
+void OnsetDetectionFunction :: perform_FFT()
+{
+	int fsize2 = (framesize/2);
+
+	// window frame and copy to complex array, swapping the first and second half of the signal
+	for (int i = 0;i < fsize2;i++)
+	{
+		in[i][0] = frame[i+fsize2] * window[i+fsize2];
+		in[i][1] = 0.0;
+		in[i+fsize2][0] = frame[i] * window[i];
+		in[i+fsize2][1] = 0.0;
+	}
+
+	// perform the fft
+	fftw_execute(p);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////// Methods for Detection Functions /////////////////////////////////
+
+//--------------------------------------------------------------------------------------
+// calculates an energy envelope detection function sample
+double OnsetDetectionFunction :: energy_envelope()
+{
+	double sum;
+
+	sum = 0;	// initialise sum
+
+	// sum the squares of the samples
+	for (int i = 0;i < framesize;i++)
+	{
+		sum = sum + (frame[i]*frame[i]);
+	}
+
+	return sum;		// return sum
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a half-wave rectified energy difference detection function sample
+double OnsetDetectionFunction :: energy_difference()
+{
+	double sum;
+	double sample;
+
+	sum = 0;	// initialise sum
+
+	// sum the squares of the samples
+	for (int i = 0;i < framesize;i++)
+	{
+		sum = sum + (frame[i]*frame[i]);
+	}
+
+	sample = sum - energy_sum_old;	// sample is first order difference in energy
+
+	energy_sum_old = sum;	// store energy value for next calculation
+
+	if (sample > 0)
+	{
+		return sample;		// return difference
+	}
+	else
+	{
+		return 0;
+	}
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a spectral difference detection function sample
+double OnsetDetectionFunction :: spectral_difference()
+{
+	double diff;
+	double sum;
+
+	// perform the FFT
+	perform_FFT();
+
+	// compute first (N/2)+1 mag values
+	for (int i = 0;i < (framesize/2)+1;i++)
+	{
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+	}
+	// mag spec symmetric above (N/2)+1 so copy previous values
+	for (int i = (framesize/2)+1;i < framesize;i++)
+	{
+		mag[i] = mag[framesize-i];
+	}
+
+	sum = 0;	// initialise sum to zero
+
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate difference
+		diff = mag[i] - mag_old[i];
+
+		// ensure all difference values are positive
+		if (diff < 0)
+		{
+			diff = diff*-1;
+		}
+
+		// add difference to sum
+		sum = sum+diff;
+
+		// store magnitude spectrum bin for next detection function sample calculation
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a spectral difference detection function sample
+double OnsetDetectionFunction :: spectral_difference_hwr()
+{
+	double diff;
+	double sum;
+
+	// perform the FFT
+	perform_FFT();
+
+	// compute first (N/2)+1 mag values
+	for (int i = 0;i < (framesize/2)+1;i++)
+	{
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+	}
+	// mag spec symmetric above (N/2)+1 so copy previous values
+	for (int i = (framesize/2)+1;i < framesize;i++)
+	{
+		mag[i] = mag[framesize-i];
+	}
+
+	sum = 0;	// initialise sum to zero
+
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate difference
+		diff = mag[i] - mag_old[i];
+
+		// only add up positive differences
+		if (diff > 0)
+		{
+			// add difference to sum
+			sum = sum+diff;
+		}
+
+
+
+		// store magnitude spectrum bin for next detection function sample calculation
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+
+//--------------------------------------------------------------------------------------
+// calculates a phase deviation detection function sample
+double OnsetDetectionFunction :: phase_deviation()
+{
+	double dev,pdev;
+	double sum;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate phase value
+		phase[i] = atan2(out[i][1],out[i][0]);
+
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+
+		// if bin is not just a low energy bin then examine phase deviation
+		if (mag[i] > 0.1)
+		{
+			dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];	// phase deviation
+			pdev = princarg(dev);	// wrap into [-pi,pi] range
+
+			// make all values positive
+			if (pdev < 0)
+			{
+				pdev = pdev*-1;
+			}
+
+			// add to sum
+			sum = sum + pdev;
+		}
+
+		// store values for next calculation
+		phase_old_2[i] = phase_old[i];
+		phase_old[i] = phase[i];
+	}
+
+	return sum;
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a complex spectral difference detection function sample
+double OnsetDetectionFunction :: complex_spectral_difference()
+{
+	double dev,pdev;
+	double sum;
+	double mag_diff,phase_diff;
+	double value;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate phase value
+		phase[i] = atan2(out[i][1],out[i][0]);
+
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+
+		// phase deviation
+		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];
+
+		// wrap into [-pi,pi] range
+		pdev = princarg(dev);
+
+
+		// calculate magnitude difference (real part of Euclidean distance between complex frames)
+		mag_diff = mag[i] - mag_old[i];
+
+		// calculate phase difference (imaginary part of Euclidean distance between complex frames)
+		phase_diff = -mag[i]*sin(pdev);
+
+
+
+		// square real and imaginary parts, sum and take square root
+		value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
+
+
+		// add to sum
+		sum = sum + value;
+
+
+		// store values for next calculation
+		phase_old_2[i] = phase_old[i];
+		phase_old[i] = phase[i];
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a complex spectral difference detection function sample (half-wave rectified)
+double OnsetDetectionFunction :: complex_spectral_difference_hwr()
+{
+	double dev,pdev;
+	double sum;
+	double mag_diff,phase_diff;
+	double value;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate phase value
+		phase[i] = atan2(out[i][1],out[i][0]);
+
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+
+		// phase deviation
+		dev = phase[i] - (2*phase_old[i]) + phase_old_2[i];
+
+		// wrap into [-pi,pi] range
+		pdev = princarg(dev);
+
+
+		// calculate magnitude difference (real part of Euclidean distance between complex frames)
+		mag_diff = mag[i] - mag_old[i];
+
+		// if we have a positive change in magnitude, then include in sum, otherwise ignore (half-wave rectification)
+		if (mag_diff > 0)
+		{
+			// calculate phase difference (imaginary part of Euclidean distance between complex frames)
+			phase_diff = -mag[i]*sin(pdev);
+
+			// square real and imaginary parts, sum and take square root
+			value = sqrt(pow(mag_diff,2) + pow(phase_diff,2));
+
+			// add to sum
+			sum = sum + value;
+		}
+
+		// store values for next calculation
+		phase_old_2[i] = phase_old[i];
+		phase_old[i] = phase[i];
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+
+//--------------------------------------------------------------------------------------
+// calculates a high frequency content detection function sample
+double OnsetDetectionFunction :: high_frequency_content()
+{
+	double sum;
+	double mag_diff;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+
+		sum = sum + (mag[i]*((double) (i+1)));
+
+		// store values for next calculation
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a high frequency spectral difference detection function sample
+double OnsetDetectionFunction :: high_frequency_spectral_difference()
+{
+	double sum;
+	double mag_diff;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+		// calculate difference
+		mag_diff = mag[i] - mag_old[i];
+
+		if (mag_diff < 0)
+		{
+			mag_diff = -mag_diff;
+		}
+
+		sum = sum + (mag_diff*((double) (i+1)));
+
+		// store values for next calculation
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+//--------------------------------------------------------------------------------------
+// calculates a high frequency spectral difference detection function sample (half-wave rectified)
+double OnsetDetectionFunction :: high_frequency_spectral_difference_hwr()
+{
+	double sum;
+	double mag_diff;
+
+	// perform the FFT
+	perform_FFT();
+
+	sum = 0; // initialise sum to zero
+
+	// compute phase values from fft output and sum deviations
+	for (int i = 0;i < framesize;i++)
+	{
+		// calculate magnitude value
+		mag[i] = sqrt(pow(out[i][0],2) + pow(out[i][1],2));
+
+		// calculate difference
+		mag_diff = mag[i] - mag_old[i];
+
+		if (mag_diff > 0)
+		{
+			sum = sum + (mag_diff*((double) (i+1)));
+		}
+
+		// store values for next calculation
+		mag_old[i] = mag[i];
+	}
+
+	return sum;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////// Methods to Calculate Windows ////////////////////////////////////
+
+//--------------------------------------------------------------------------------------
+// HANNING: set the window in the buffer 'window' to a Hanning window
+void OnsetDetectionFunction :: set_win_hanning()
+{
+	double N;		// variable to store framesize minus 1
+
+	N = (double) (framesize-1);	// framesize minus 1
+
+	// Hanning window calculation
+	for (int n = 0;n < framesize;n++)
+	{
+		window[n] = 0.5*(1-cos(2*pi*(n/N)));
+	}
+}
+
+//--------------------------------------------------------------------------------------
+// HAMMING: set the window in the buffer 'window' to a Hanning window
+void OnsetDetectionFunction :: set_win_hamming()
+{
+	double N;		// variable to store framesize minus 1
+	double n_val;	// double version of index 'n'
+
+	N = (double) (framesize-1);	// framesize minus 1
+	n_val = 0;
+
+	// Hamming window calculation
+	for (int n = 0;n < framesize;n++)
+	{
+		window[n] = 0.54 - (0.46*cos(2*pi*(n_val/N)));
+		n_val = n_val+1;
+	}
+}
+
+//--------------------------------------------------------------------------------------
+// BLACKMAN: set the window in the buffer 'window' to a Blackman window
+void OnsetDetectionFunction :: set_win_blackman()
+{
+	double N;		// variable to store framesize minus 1
+	double n_val;	// double version of index 'n'
+
+	N = (double) (framesize-1);	// framesize minus 1
+	n_val = 0;
+
+	// Blackman window calculation
+	for (int n = 0;n < framesize;n++)
+	{
+		window[n] = 0.42 - (0.5*cos(2*pi*(n_val/N))) + (0.08*cos(4*pi*(n_val/N)));
+		n_val = n_val+1;
+	}
+}
+
+//--------------------------------------------------------------------------------------
+// TUKEY: set the window in the buffer 'window' to a Tukey window
+void OnsetDetectionFunction :: set_win_tukey()
+{
+	double N;		// variable to store framesize minus 1
+	double n_val;	// double version of index 'n'
+	double alpha;	// alpha [default value = 0.5];
+
+	int lim1,lim2;
+
+	alpha = 0.5;
+
+	N = (double) (framesize-1);	// framesize minus 1
+
+	// Tukey window calculation
+
+	n_val = (double) (-1*((framesize/2)))+1;
+
+	for (int n = 0;n < framesize;n++)	// left taper
+	{
+		if ((n_val >= 0) && (n_val <= (alpha*(N/2))))
+		{
+			window[n] = 1.0;
+		}
+		else if ((n_val <= 0) && (n_val >= (-1*alpha*(N/2))))
+		{
+			window[n] = 1.0;
+		}
+		else
+		{
+			window[n] = 0.5*(1+cos(pi*(((2*n_val)/(alpha*N))-1)));
+		}
+
+		n_val = n_val+1;
+	}
+
+}
+
+//--------------------------------------------------------------------------------------
+// RECTANGULAR: set the window in the buffer 'window' to a Rectangular window
+void OnsetDetectionFunction :: set_win_rectangular()
+{
+	// Rectangular window calculation
+	for (int n = 0;n < framesize;n++)
+	{
+		window[n] = 1.0;
+	}
+}
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////
+///////////////////////////////// Other Handy Methods //////////////////////////////////////////
+
+//--------------------------------------------------------------------------------------
+// set phase values to the range [-pi,pi]
+double OnsetDetectionFunction :: princarg(double phaseval)
+{
+	// if phase value is less than or equal to -pi then add 2*pi
+	while (phaseval <= (-pi))
+	{
+		phaseval = phaseval + (2*pi);
+	}
+
+	// if phase value is larger than pi, then subtract 2*pi
+	while (phaseval > pi)
+	{
+		phaseval = phaseval - (2*pi);
+	}
+
+	return phaseval;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/OnsetDetectionFunction.h	Tue Jan 21 00:10:11 2014 +0000
@@ -0,0 +1,92 @@
+//=======================================================================
+/** @file OnsetDetectionFunction.h
+ *  @brief A class for calculating onset detection functions
+ *  @author Adam Stark
+ *  @copyright Copyright (C) 2008-2014  Queen Mary University of London
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+//=======================================================================
+
+#ifndef __RTONSETDF_H
+#define __RTONSETDF_H
+
+#include "fftw3.h"
+
+class OnsetDetectionFunction
+{
+public:
+	OnsetDetectionFunction();																// Constructor
+	OnsetDetectionFunction(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type);	// Constructor (with arguments)
+	~OnsetDetectionFunction();																// Destructor
+	void initialise(int arg_hsize,int arg_fsize,int arg_df_type,int arg_win_type);			// Initialisation Function
+
+	double getDFsample(double inputbuffer[]);												// process input buffer and calculate detection function sample
+	void set_df_type(int arg_df_type);														// set the detection function type
+
+private:
+
+	void perform_FFT();																		// perform the FFT on the data in 'frame'
+
+	double energy_envelope();																// calculate energy envelope detection function sample
+	double energy_difference();																// calculate energy difference detection function sample
+	double spectral_difference();															// calculate spectral difference detection function sample
+	double spectral_difference_hwr();														// calculate spectral difference (half wave rectified) detection function sample
+	double phase_deviation();																// calculate phase deviation detection function sample
+	double complex_spectral_difference();													// calculate complex spectral difference detection function sample
+	double complex_spectral_difference_hwr();												// calculate complex spectral difference detection function sample (half-wave rectified)
+	double high_frequency_content();														// calculate high frequency content detection function sample
+	double high_frequency_spectral_difference();											// calculate high frequency spectral difference detection function sample
+	double high_frequency_spectral_difference_hwr();										// calculate high frequency spectral difference detection function sample (half-wave rectified)
+
+	void set_win_rectangular();																// calculate a Rectangular window
+	void set_win_hanning();																	// calculate a Hanning window
+	void set_win_hamming();																	// calculate a Hamming window
+	void set_win_blackman();																// calculate a Blackman window
+	void set_win_tukey();																	// calculate a Tukey window
+
+
+	double princarg(double phaseval);														// set phase values between [-pi, pi]
+
+
+	double pi;																				// pi, the constant
+
+	int framesize;																			// audio framesize
+	int hopsize;																			// audio hopsize
+	int df_type;																			// type of detection function
+
+	fftw_plan p;																			// create fft plan
+	fftw_complex *in;																		// to hold complex fft values for input
+	fftw_complex *out;																		// to hold complex fft values for output
+
+	int initialised;																		// flag indicating whether buffers and FFT plans have been initialised
+
+
+	double *frame;																			// audio frame
+	double *window;																			// window
+	double *wframe;																			// windowed frame
+
+	double energy_sum_old;																	// to hold the previous energy sum value
+
+	double *mag;																			// magnitude spectrum
+	double *mag_old;																		// previous magnitude spectrum
+
+	double *phase;																			// FFT phase values
+	double *phase_old;																		// previous phase values
+	double *phase_old_2;																	// second order previous phase values
+
+};
+
+
+#endif
\ No newline at end of file