# HG changeset patch
# User Chris Cannam <c.cannam@qmul.ac.uk>
# Date 1144258559 0
# Node ID 49844bc8a89583c99cae1d50b980320c96201cbe

* Queen Mary C++ DSP library

diff -r 000000000000 -r 49844bc8a895 base/Pitch.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/Pitch.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,41 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+    This file Copyright 2006 Chris Cannam.
+    All rights reserved.
+*/
+
+#include "Pitch.h"
+
+#include <cmath>
+
+float
+Pitch::getFrequencyForPitch(int midiPitch,
+			    float centsOffset,
+			    float concertA)
+{
+    float p = float(midiPitch) + (centsOffset / 100);
+    return concertA * powf(2.0, (p - 69.0) / 12.0);
+}
+
+int
+Pitch::getPitchForFrequency(float frequency,
+			    float *centsOffsetReturn,
+			    float concertA)
+{
+    float p = 12.0 * (log(frequency / (concertA / 2.0)) / log(2.0)) + 57.0;
+
+    int midiPitch = int(p + 0.00001);
+    float centsOffset = (p - midiPitch) * 100.0;
+
+    if (centsOffset >= 50.0) {
+	midiPitch = midiPitch + 1;
+	centsOffset = -(100.0 - centsOffset);
+    }
+    
+    if (centsOffsetReturn) *centsOffsetReturn = centsOffset;
+    return midiPitch;
+}
+
diff -r 000000000000 -r 49844bc8a895 base/Pitch.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/Pitch.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,26 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+    This file Copyright 2006 Chris Cannam.
+    All rights reserved.
+*/
+
+#ifndef _PITCH_H_
+#define _PITCH_H_
+
+class Pitch
+{
+public:
+    static float getFrequencyForPitch(int midiPitch,
+				      float centsOffset = 0,
+				      float concertA = 440.0);
+
+    static int getPitchForFrequency(float frequency,
+				    float *centsOffsetReturn = 0,
+				    float concertA = 440.0);
+};
+
+
+#endif
diff -r 000000000000 -r 49844bc8a895 base/Window.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/Window.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,120 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP library
+    Centre for Digital Music, Queen Mary, University of London.
+    This file Copyright 2006 Chris Cannam.
+    All rights reserved.
+*/
+
+#ifndef _WINDOW_H_
+#define _WINDOW_H_
+
+#include <cmath>
+#include <iostream>
+#include <map>
+
+enum WindowType {
+    RectangularWindow,
+    BartlettWindow,
+    HammingWindow,
+    HanningWindow,
+    BlackmanWindow,
+    GaussianWindow,
+    ParzenWindow
+};
+
+template <typename T>
+class Window
+{
+public:
+    /**
+     * Construct a windower of the given type.
+     */
+    Window(WindowType type, size_t size) : m_type(type), m_size(size) { encache(); }
+    Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); }
+    Window &operator=(const Window &w) {
+	if (&w == this) return *this;
+	m_type = w.m_type;
+	m_size = w.m_size;
+	encache();
+	return *this;
+    }
+    virtual ~Window() { delete m_cache; }
+    
+    void cut(T *src) const { cut(src, src); }
+    void cut(T *src, T *dst) const {
+	for (size_t i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
+    }
+
+    WindowType getType() const { return m_type; }
+    size_t getSize() const { return m_size; }
+
+protected:
+    WindowType m_type;
+    size_t m_size;
+    T *m_cache;
+    
+    void encache();
+};
+
+template <typename T>
+void Window<T>::encache()
+{
+    size_t n = m_size;
+    T *mult = new T[n];
+    size_t i;
+    for (i = 0; i < n; ++i) mult[i] = 1.0;
+
+    switch (m_type) {
+		
+    case RectangularWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * 0.5;
+	}
+	break;
+	    
+    case BartlettWindow:
+	for (i = 0; i < n/2; ++i) {
+	    mult[i] = mult[i] * (i / T(n/2));
+	    mult[i + n/2] = mult[i + n/2] * (1.0 - (i / T(n/2)));
+	}
+	break;
+	    
+    case HammingWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n));
+	}
+	break;
+	    
+    case HanningWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n));
+	}
+	break;
+	    
+    case BlackmanWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n)
+				 + 0.08 * cos(4 * M_PI * i / n));
+	}
+	break;
+	    
+    case GaussianWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * exp((-1.0 / (n*n)) * ((T(2*i) - n) *
+						      (T(2*i) - n)));
+	}
+	break;
+	    
+    case ParzenWindow:
+	for (i = 0; i < n; ++i) {
+	    mult[i] = mult[i] * (1.0 - fabs((T(2*i) - n) / T(n + 1)));
+	}
+	break;
+    }
+	
+    m_cache = mult;
+}
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/ChromaProcess.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/ChromaProcess.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,213 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+
+#include "ChromaProcess.h"
+#include "dsp/maths/Histogram.h"
+#include <math.h>
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+ChromaProcess::ChromaProcess()
+{
+
+}
+
+ChromaProcess::~ChromaProcess()
+{
+
+}
+
+int ChromaProcess::findChromaBias( vector<double> chromaVector, unsigned int BPO, unsigned int frames )
+{
+    vector<double> newChroma;
+    vector<int> peakIndex;
+    vector<int> modPeakIndex;
+
+    unsigned int chromaLength = chromaVector.size(); 
+
+    unsigned int newLength = chromaLength + (2*BPO);
+
+    newChroma.resize( newLength );
+    newChroma.clear();
+
+    modPeakIndex.resize( newLength );
+    modPeakIndex.clear();
+
+    //adds last row at the top and first row at the bottom to create 
+    //circularity - effectively adds 2 to the bpo-length of the vectors:
+
+    for( unsigned int i = 0; i < BPO; i++ )
+    {
+	newChroma.push_back( chromaVector[ chromaLength - BPO + i ] );
+    }
+
+    for( unsigned i = 0; i < chromaLength; i++ )
+    {
+	newChroma.push_back( chromaVector[ i ] );
+    }
+
+    for( unsigned i = 0; i < BPO; i++ )
+    {
+	newChroma.push_back( chromaVector[ i ] );
+    }
+
+    // pick peaks in the chroma
+    peakIndex = getPeaks( newChroma, BPO );
+
+    // modularises to res = bpo/12 bins:
+    // corrects the mod value for bin 3
+    modPeakIndex = mod( peakIndex, 3 );
+
+    // finds the highest concentration of peaks on the bpo/12 bin resolution
+    THistogram<int> m_hist(3);
+
+    double ave, adev, sdev, var, skew, ccurt;
+
+    m_hist.compute(modPeakIndex);
+
+    m_hist.getMoments( modPeakIndex, ave, adev, sdev, var, skew, ccurt );
+
+    vector <double> histogram = m_hist.geTHistogramD();
+    //////////////////////////////////////////////////////////////////////////////
+
+    ///////////////////////////////////////////////////////////////////////////
+	// Find actual bias from histogram
+	int minIdx, maxIdx;
+	double min, max;
+
+	findHistMaxMin( histogram, &max, &maxIdx, &min, &minIdx );
+
+/*
+  FILE* foutchroma = fopen("../testdata/newchroma.bin","wb");
+  FILE* foutpeaks = fopen("../testdata/peaks.bin","wb");
+
+
+  fwrite( &chromaVector[0], sizeof(double), chromaVector.size(), foutchroma );
+  fwrite( &histogram[0], sizeof(double), histogram.size(), foutpeaks );
+
+  fclose( foutchroma );
+  fclose( foutpeaks );
+*/
+	return maxIdx - 1;
+}
+
+
+vector <int> ChromaProcess::getPeaks(vector <double> chroma, unsigned int BPO)
+{
+    vector <int> peaks;
+	
+    double pre = 0;
+    double post = 0;
+    double current = 0;
+
+    unsigned int BPOCounter = 0;
+    unsigned int mult = 0;
+    unsigned int idx = 0;
+
+    for( unsigned int i = 0; i < chroma.size() - 0; i++ )
+    {
+	BPOCounter++;
+
+	pre = chroma[ i ];
+	current = chroma[ i + 1 ];
+	post = chroma[ i + 2 ];
+		
+	if( (current > 0) && (current > pre) && (current > post) )
+	{
+	    peaks.push_back( BPOCounter + 1);
+	}
+
+		
+	if( BPOCounter == (BPO - 2 ) )
+	{
+	    BPOCounter = 0;
+	    i+=2;
+	}
+		
+    }
+
+    /*
+      for( unsigned int i = 1; i < chroma.size() - 1; i++ )
+      {
+      BPOCounter++ ;
+
+      pre = chroma[ i - 1 ];
+      current = chroma[ i ];
+      post = chroma[ i + 1 ];
+		
+      if( (current > 0) && (current > pre) && (current > post) )
+      {
+      peaks.push_back( BPOCounter + 1 );
+      }
+
+      if( BPOCounter == (PO - 1) )
+      {
+      BPOCounter = 1;
+      i+=2;
+      }
+      }
+    */
+    return peaks;
+}
+
+vector <int> ChromaProcess::mod(vector <int> input, int res)
+{
+    vector <int> result;
+
+    for( unsigned int i = 0; i < input.size(); i++ )
+    {
+	int val = input[ i ];
+	int res = val - res * floor( (double)val / (double)res );
+
+	if( val != 0 )
+	{
+	    if( res == 0 )
+		res = 3;
+
+	    result.push_back( res );
+	}
+	else
+	{
+	    result.push_back( val );
+	}
+    }
+    return result;
+}
+
+void ChromaProcess::findHistMaxMin( vector<double> hist, double* max, int* maxIdx, double* min, int* minIdx )
+{
+    double temp = 0.0;
+    unsigned int vecLength = hist.size();
+
+    *minIdx = 0;
+    *maxIdx = 0;
+
+    *min = hist[0];
+    *max = *min;
+
+    for( unsigned int u = 0; u < vecLength; u++ )
+    {
+	temp = hist[ u ];
+
+	if( temp < *min )
+	{
+	    *min =  temp ;
+	    *minIdx = u;
+	}
+	if( temp > *max )
+	{
+	    *max =  temp ;
+	    *maxIdx = u;
+	}
+
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/ChromaProcess.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/ChromaProcess.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,30 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef CHROMAPROCESS_H
+#define CHROMAPROCESS_H
+
+#include <vector>
+
+using namespace std;
+
+class ChromaProcess  
+{
+public:
+    void findHistMaxMin( vector<double> hist, double* max, int*maxIdx, double* min, int* minIdx );
+    vector <int> mod( vector <int> input, int res );
+    vector <int> getPeaks( vector <double> chroma, unsigned int BPO );
+    int findChromaBias( vector<double> chromaVector, unsigned int BPO, unsigned int frames );
+    ChromaProcess();
+    virtual ~ChromaProcess();
+
+};
+
+#endif // !defined(CHROMAPROCESS_H)
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/Chromagram.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/Chromagram.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,146 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include <iostream>
+#include <cmath>
+#include "dsp/maths/MathUtilities.h"
+#include "Chromagram.h"
+
+//----------------------------------------------------------------------------
+
+Chromagram::Chromagram( ChromaConfig Config ) 
+{
+    initialise( Config );
+}
+
+int Chromagram::initialise( ChromaConfig Config )
+{	
+    m_FMin = Config.min;		// min freq
+    m_FMax = Config.max;		// max freq
+    m_BPO  = Config.BPO;		// bins per octave
+    isNormalised = Config.isNormalised; // true if frame normalisation is required
+
+    // No. of constant Q bins
+    m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0));	
+
+    // Create array for chroma result
+    m_chromadata = new double[ m_BPO ];
+
+    // Initialise FFT object	
+    m_FFT = new FFT;
+
+    // Create Config Structure for ConstantQ operator
+    CQConfig ConstantQConfig;
+
+    // Populate CQ config structure with parameters
+    // inherited from the Chroma config
+    ConstantQConfig.FS	 = Config.FS;
+    ConstantQConfig.min = m_FMin;
+    ConstantQConfig.max = m_FMax;
+    ConstantQConfig.BPO = m_BPO;
+    ConstantQConfig.CQThresh = Config.CQThresh;
+	
+    // Initialise ConstantQ operator
+    m_ConstantQ = new ConstantQ( ConstantQConfig );
+
+    // Initialise working arrays
+    m_frameSize = m_ConstantQ->getfftlength();
+    m_hopSize = m_ConstantQ->gethop();
+
+    m_FFTRe = new double[ m_frameSize ];
+    m_FFTIm = new double[ m_frameSize ];
+    m_CQRe  = new double[ m_uK ];
+    m_CQIm  = new double[ m_uK ];
+
+    // Generate CQ Kernel 
+    m_ConstantQ->sparsekernel();
+    return 1;
+}
+
+Chromagram::~Chromagram()
+{
+    deInitialise();
+}
+
+int Chromagram::deInitialise()
+{
+    delete [] m_chromadata;
+
+    delete m_FFT;
+
+    delete m_ConstantQ;
+
+    delete [] m_FFTRe;
+    delete [] m_FFTIm;
+    delete [] m_CQRe;
+    delete [] m_CQIm;
+    return 1;
+}
+
+//----------------------------------------------------------------------------------
+// returns the absolute value of complex number xx + i*yy
+double Chromagram::kabs(double xx, double yy)
+{
+    double ab = sqrt(xx*xx + yy*yy);
+    return(ab);
+}
+//-----------------------------------------------------------------------------------
+
+
+void Chromagram::unityNormalise(double *src)
+{
+    double min, max;
+
+    double val = 0;
+
+    MathUtilities::getFrameMinMax( src, m_BPO, & min, &max );
+
+    for( unsigned int i = 0; i < m_BPO; i++ )
+    {
+	val = src[ i ] / max;
+
+	src[ i ] = val;
+    }
+}
+
+
+double* Chromagram::process( double *data )
+{
+    //initialise chromadata to 0
+    for (unsigned i=0; i<m_BPO; i++) 
+	m_chromadata[i]=0;
+
+    double cmax = 0.0;
+    double cval = 0;
+
+    // FFT of current frame
+    m_FFT->process( m_frameSize, 0, data, NULL, m_FFTRe, m_FFTIm ); 
+
+    // Calculate ConstantQ frame
+    m_ConstantQ->process( m_FFTRe, m_FFTIm, m_CQRe, m_CQIm );
+	
+    // add each octave of cq data into Chromagram
+    const unsigned octaves = (int)floor(double( m_uK/m_BPO))-1;
+    for (unsigned octave=0; octave<=octaves; octave++) 
+    {
+	unsigned firstBin = octave*m_BPO;
+	for (unsigned i=0; i<m_BPO; i++) 
+	{
+	    m_chromadata[i] += kabs( m_CQRe[ firstBin + i ], m_CQIm[ firstBin + i ]);
+	}
+    }
+
+    if( isNormalised )
+	unityNormalise( m_chromadata );
+
+    return m_chromadata;
+}
+
+
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/Chromagram.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/Chromagram.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,69 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef CHROMAGRAM_H
+#define CHROMAGRAM_H
+
+#include "dsp/transforms/FFT.h"
+#include "ConstantQ.h"
+
+struct ChromaConfig{
+    unsigned int FS;
+    double min;
+    double max;
+    unsigned int BPO;
+    double CQThresh;
+    bool isNormalised;
+};
+
+class Chromagram 
+{
+
+public:	
+    Chromagram( ChromaConfig Config );
+    ~Chromagram();
+	
+    double* process( double *data );
+    void unityNormalise( double* src );
+
+    // Complex arithmetic
+    double kabs( double real, double imag );
+	
+    // Results
+    unsigned int getK() { return m_uK;}
+    unsigned int getFrameSize() { return m_frameSize; }
+    unsigned int getHopSize()   { return m_hopSize; }
+
+private:
+    int initialise( ChromaConfig Config );
+    int deInitialise();
+	
+    double* m_chromadata;
+    double m_FMin;
+    double m_FMax;
+    unsigned int m_BPO;
+    unsigned int m_uK;
+
+    bool isNormalised;
+
+    unsigned int m_frameSize;
+    unsigned int m_hopSize;
+
+    FFT*	 m_FFT;
+    ConstantQ* m_ConstantQ;
+
+    double* m_FFTRe;
+    double* m_FFTIm;
+    double* m_CQRe;
+    double* m_CQIm;
+
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/ConstantQ.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/ConstantQ.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,193 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "ConstantQ.h"
+#include "dsp/transforms/FFT.h"
+
+//---------------------------------------------------------------------------
+// nextpow2 returns the smallest integer n such that 2^n >= x.
+static double nextpow2(double x) {
+    double y = ceil(log(x)/log(2.0));
+    return(y);
+}
+
+static double squaredModule(const double & xx, const double & yy) {
+    return xx*xx + yy*yy;
+}
+
+//----------------------------------------------------------------------------
+
+ConstantQ::ConstantQ( CQConfig Config ) 
+{
+    initialise( Config );
+}
+
+ConstantQ::~ConstantQ()
+{
+    deInitialise();
+}
+
+//----------------------------------------------------------------------------
+void ConstantQ::sparsekernel()
+{
+    //generates spectral kernel matrix (upside down?)
+    // initialise temporal kernel with zeros, twice length to deal w. complex numbers
+
+    double* hammingWindowRe = new double [ m_FFTLength ];
+    double* hammingWindowIm = new double [ m_FFTLength ];
+    double* transfHammingWindowRe = new double [ m_FFTLength ];
+    double* transfHammingWindowIm = new double [ m_FFTLength ];
+
+    for (unsigned u=0; u < m_FFTLength; u++) 
+    {
+	hammingWindowRe[u] = 0;
+	hammingWindowIm[u] = 0;
+    }
+
+
+    // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
+    // The matrix K x fftlength but the non-zero cells are an antialiased
+    // square root function. So mostly is a line, with some grey point.
+    m_sparseKernelIs.reserve( m_FFTLength*2 );
+    m_sparseKernelJs.reserve( m_FFTLength*2 );
+    m_sparseKernelRealValues.reserve( m_FFTLength*2 );
+    m_sparseKernelImagValues.reserve( m_FFTLength*2 );
+	
+    // for each bin value K, calculate temporal kernel, take its fft to
+    //calculate the spectral kernel then threshold it to make it sparse and 
+    //add it to the sparse kernels matrix
+    double squareThreshold = m_CQThresh * m_CQThresh;
+
+    FFT m_FFT;
+	
+    for (unsigned k = m_uK; k--; ) 
+    {
+	// Computing a hamming window
+	const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
+	for (unsigned i=0; i<hammingLength; i++) 
+	{
+	    const double angle = 2*PI*m_dQ*i/hammingLength;
+	    const double real = cos(angle);
+	    const double imag = sin(angle);
+	    const double absol = hamming(hammingLength, i)/hammingLength;
+	    hammingWindowRe[ i ] = absol*real;
+	    hammingWindowIm[ i ] = absol*imag;
+	}
+
+	//do fft of hammingWindow
+	m_FFT.process( m_FFTLength, 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
+
+		
+	for (unsigned j=0; j<( m_FFTLength ); j++) 
+	{
+	    // perform thresholding
+	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
+	    if (squaredBin <= squareThreshold) continue;
+		
+	    // Insert non-zero position indexes, doubled because they are floats
+	    m_sparseKernelIs.push_back(j);
+	    m_sparseKernelJs.push_back(k);
+
+	    // take conjugate, normalise and add to array sparkernel
+	    m_sparseKernelRealValues.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
+	    m_sparseKernelImagValues.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
+	}
+
+    }
+
+    delete [] hammingWindowRe;
+    delete [] hammingWindowIm;
+    delete [] transfHammingWindowRe;
+    delete [] transfHammingWindowIm;
+
+}
+
+//-----------------------------------------------------------------------------
+double* ConstantQ::process( double* fftdata )
+{
+    for (unsigned row=0; row<2*m_uK; row++) 
+    {
+	m_CQdata[ row ] = 0;
+	m_CQdata[ row+1 ] = 0;
+    }
+    const unsigned *fftbin = &(m_sparseKernelIs[0]);
+    const unsigned *cqbin  = &(m_sparseKernelJs[0]);
+    const double   *real   = &(m_sparseKernelRealValues[0]);
+    const double   *imag   = &(m_sparseKernelImagValues[0]);
+    const unsigned int sparseCells = m_sparseKernelRealValues.size();
+	
+    for (unsigned i = 0; i<sparseCells; i++)
+    {
+	const unsigned row = cqbin[i];
+	const unsigned col = fftbin[i];
+	const double & r1  = real[i];
+	const double & i1  = imag[i];
+	const double & r2  = fftdata[ (2*m_FFTLength) - 2*col];
+	const double & i2  = fftdata[ (2*m_FFTLength) - 2*col+1];
+	// add the multiplication
+	m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
+	m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
+    }
+
+    return m_CQdata;
+}
+
+
+void ConstantQ::initialise( CQConfig Config )
+{
+    m_FS = Config.FS;
+    m_FMin = Config.min;		// min freq
+    m_FMax = Config.max;		// max freq
+    m_BPO = Config.BPO;		// bins per octave
+    m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
+
+    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);	// Work out Q value for Filter bank
+    m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));	// No. of constant Q bins
+
+    // work out length of fft required for this constant Q Filter bank
+    m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin )));
+
+    m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32
+
+    // allocate memory for cqdata
+    m_CQdata = new double [2*m_uK];
+}
+
+void ConstantQ::deInitialise()
+{
+    delete [] m_CQdata;
+}
+
+void ConstantQ::process(double *FFTRe, double* FFTIm, double *CQRe, double *CQIm)
+{
+    for (unsigned row=0; row<m_uK; row++) 
+    {
+	CQRe[ row ] = 0;
+	CQIm[ row ] = 0;
+    }
+
+    const unsigned *fftbin = &(m_sparseKernelIs[0]);
+    const unsigned *cqbin  = &(m_sparseKernelJs[0]);
+    const double   *real   = &(m_sparseKernelRealValues[0]);
+    const double   *imag   = &(m_sparseKernelImagValues[0]);
+    const unsigned int sparseCells = m_sparseKernelRealValues.size();
+	
+    for (unsigned i = 0; i<sparseCells; i++)
+    {
+	const unsigned row = cqbin[i];
+	const unsigned col = fftbin[i];
+	const double & r1  = real[i];
+	const double & i1  = imag[i];
+	const double & r2  = FFTRe[ m_FFTLength- col];
+	const double & i2  = FFTIm[ m_FFTLength - col];
+	// add the multiplication
+	CQRe[ row ] += (r1*r2 - i1*i2);
+	CQIm[ row ] += (r1*i2 + i1*r2);
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/chromagram/ConstantQ.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/chromagram/ConstantQ.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,73 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef CONSTANTQ_H
+#define CONSTANTQ_H
+
+#include <vector>
+#include "dsp/maths/MathAliases.h"
+#include "dsp/maths/MathUtilities.h"
+
+struct CQConfig{
+    unsigned int FS;
+    double min;
+    double max;
+    unsigned int BPO;
+    double CQThresh;
+};
+
+class ConstantQ {
+	
+//public functions incl. sparsekernel so can keep out of loop in main
+public:
+    void process( double* FFTRe, double* FFTIm, double* CQRe, double* CQIm );
+
+    ConstantQ( CQConfig Config );
+    ~ConstantQ();
+
+    double* process( double* FFTData );
+
+    void sparsekernel();
+
+    double hamming(int len, int n) {
+	double out = 0.54 - 0.46*cos(2*PI*n/len);
+	return(out);
+    }
+	
+    int getnumwin() { return m_numWin;}
+    double getQ() { return m_dQ;}
+    int getK() {return m_uK ;}
+    int getfftlength() { return m_FFTLength;}
+    int gethop() { return m_hop;}
+
+private:
+    void initialise( CQConfig Config );
+    void deInitialise();
+	
+    double* m_CQdata;
+    unsigned int m_FS;
+    double m_FMin;
+    double m_FMax;
+    double m_dQ;
+    double m_CQThresh;
+    unsigned int m_numWin;
+    unsigned int m_hop;
+    unsigned int m_BPO;
+    unsigned int m_FFTLength;
+    unsigned int m_uK;
+    std::vector<unsigned> m_sparseKernelIs;
+    std::vector<unsigned> m_sparseKernelJs;
+    std::vector<double> m_sparseKernelImagValues;
+    std::vector<double> m_sparseKernelRealValues;
+};
+
+
+#endif//CONSTANTQ_H
+
diff -r 000000000000 -r 49844bc8a895 dsp/maths/Correlation.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/Correlation.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,51 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "Correlation.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Correlation::Correlation()
+{
+
+}
+
+Correlation::~Correlation()
+{
+
+}
+
+void Correlation::doAutoUnBiased(double *src, double *dst, unsigned int length)
+{
+    double tmp = 0.0;
+    double outVal = 0.0;
+
+    unsigned int i,j;
+
+    for( i = 0; i <  length; i++)
+    {
+	for( j = i; j < length; j++)
+	{
+	    tmp += src[ j-i ] * src[ j ]; 
+	}
+
+
+	outVal = tmp / ( length - i );
+
+	if( outVal <= 0 )
+	    dst[ i ] = EPS;
+	else
+	    dst[ i ] = outVal;
+
+	tmp = 0.0;
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/maths/Correlation.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/Correlation.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,25 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef CORRELATION_H
+#define CORRELATION_H
+
+#define  EPS  2.2204e-016
+
+class Correlation  
+{
+public:
+    void doAutoUnBiased( double* src, double* dst, unsigned int length );
+    Correlation();
+    virtual ~Correlation();
+
+};
+
+#endif // 
diff -r 000000000000 -r 49844bc8a895 dsp/maths/Histogram.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/Histogram.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,390 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+// Histogram.h: interface for the THistogram class.
+//
+//////////////////////////////////////////////////////////////////////
+
+
+#ifndef HISTOGRAM_H
+#define HISTOGRAM_H
+
+
+#include <valarray>
+
+/*! \brief A histogram class
+
+  This class computes the histogram of a vector.
+
+\par Template parameters
+	
+	- T type of input data (can be any: float, double, int, UINT, etc...)
+	- TOut type of output data: float or double. (default is double)
+
+\par Moments:
+
+	The moments (average, standard deviation, skewness, etc.) are computed using 
+the algorithm of the Numerical recipies (see Numerical recipies in C, Chapter 14.1, pg 613).
+
+\par Example:
+
+	This example shows the typical use of the class:
+\code
+	// a vector containing the data
+	vector<float> data;
+	// Creating histogram using float data and with 101 containers,
+	THistogram<float> histo(101);
+	// computing the histogram
+	histo.compute(data);
+\endcode
+
+Once this is done, you can get a vector with the histogram or the normalized histogram (such that it's area is 1):
+\code
+	// getting normalized histogram
+	vector<float> v=histo.getNormalizedHistogram();
+\endcode
+
+\par Reference
+	
+	Equally spaced acsissa function integration (used in #GetArea): Numerical Recipies in C, Chapter 4.1, pg 130.
+
+\author Jonathan de Halleux, dehalleux@auto.ucl.ac.be, 2002
+*/
+
+template<class T, class TOut = double>
+class THistogram  
+{
+public:
+    //! \name Constructors
+    //@{
+    /*! Default constructor
+      \param counters the number of histogram containers (default value is 10)
+    */
+    THistogram(unsigned int counters = 10);
+    virtual ~THistogram()				{ clear();};
+    //@}
+
+    //! \name Histogram computation, update
+    //@{
+    /*! Computes histogram of vector v
+      \param v vector to compute histogram
+      \param computeMinMax set to true if min/max of v have to be used to get the histogram limits	
+	
+      This function computes the histogram of v and stores it internally.
+      \sa Update, GeTHistogram
+    */
+    void compute( const std::vector<T>& v, bool computeMinMax = true);
+    //! Update histogram with the vector v
+    void update( const std::vector<T>& v);
+    //! Update histogram with t
+    void update( const T& t);
+    //@}
+
+    //! \name Resetting functions
+    //@{
+    //! Resize the histogram. Warning this function clear the histogram.
+    void resize( unsigned int counters );
+    //! Clears the histogram
+    void clear()						{ m_counters.clear();};
+    //@}
+
+    //! \name Setters
+    //@{
+    /*! This function sets the minimum of the histogram spectrum. 
+      The spectrum is not recomputed, use it with care
+    */
+    void setMinSpectrum( const T& min )					{	m_min = min; computeStep();};
+    /*! This function sets the minimum of the histogram spectrum. 
+      The spectrum is not recomputed, use it with care
+    */
+    void setMaxSpectrum( const T& max )					{	m_max = max; computeStep();};
+    //@}
+    //! \name Getters
+    //@{
+    //! return minimum of histogram spectrum
+    const T& getMinSpectrum() const							{	return m_min;};
+    //! return maximum of histogram spectrum
+    const T& getMaxSpectrum() const							{	return m_max;};
+    //! return step size of histogram containers
+    TOut getStep() const									{	return m_step;};
+    //! return number of points in histogram
+    unsigned int getSum() const;
+    /*! \brief returns area under the histogram 
+
+    The Simpson rule is used to integrate the histogram.
+    */
+    TOut getArea() const;
+
+    /*! \brief Computes the moments of the histogram
+
+    \param data dataset
+    \param ave mean
+    \f[ \bar x = \frac{1}{N} \sum_{j=1}^N x_j\f]
+    \param adev mean absolute deviation
+    \f[ adev(X) = \frac{1}{N} \sum_{j=1}^N | x_j - \bar x |\f]
+    \param var average deviation:
+    \f[ \mbox{Var}(X) = \frac{1}{N-1} \sum_{j=1}^N (x_j - \bar x)^2\f]
+    \param sdev standard deviation:
+    \f[ \sigma(X) = \sqrt{var(\bar x) }\f]
+    \param skew skewness
+    \f[ \mbox{Skew}(X) = \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^3\f]
+    \param kurt kurtosis
+    \f[ \mbox{Kurt}(X) = \left\{ \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^4 \right\} - 3\f]
+
+    */
+    static void getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt);
+
+    //! return number of containers
+    unsigned int getSize() const							{	return m_counters.size();};
+    //! returns i-th counter
+    unsigned int operator [] (unsigned int i) const					{	ASSERT( i < m_counters.size() ); return m_counters[i];};
+    //! return the computed histogram
+    const std::vector<unsigned int>& geTHistogram() const	{	return m_counters;};
+    //! return the computed histogram, in TOuts
+    std::vector<TOut> geTHistogramD() const;
+    /*! return the normalized computed histogram 
+
+    \return the histogram such that the area is equal to 1
+    */
+    std::vector<TOut> getNormalizedHistogram() const;
+    //! returns left containers position
+    std::vector<TOut> getLeftContainers() const;
+    //! returns center containers position
+    std::vector<TOut> getCenterContainers() const;
+    //@}
+protected:
+    //! Computes the step
+    void computeStep()								{	m_step = (TOut)(((TOut)(m_max-m_min)) / (m_counters.size()-1));};
+    //! Data accumulators
+    std::vector<unsigned int> m_counters;
+    //! minimum of dataset
+    T m_min;
+    //! maximum of dataset
+    T m_max;
+    //! width of container
+    TOut m_step;
+};
+
+template<class T, class TOut>
+THistogram<T,TOut>::THistogram(unsigned int counters)
+    : m_counters(counters,0), m_min(0), m_max(0), m_step(0)
+{
+
+}
+
+template<class T, class TOut>
+void THistogram<T,TOut>::resize( unsigned int counters )
+{
+    clear();
+
+    m_counters.resize(counters,0);
+
+    computeStep();
+}
+
+template<class T, class TOut>
+void THistogram<T,TOut>::compute( const std::vector<T>& v, bool computeMinMax)
+{
+    using namespace std;
+    unsigned int i;
+    int index;
+
+    if (m_counters.empty())
+	return;
+
+    if (computeMinMax)
+    {
+	m_max = m_min = v[0];
+	for (i=1;i<v.size();i++)
+	{
+	    m_max = std::max( m_max, v[i]);
+	    m_min = std::min( m_min, v[i]);
+	}
+    }
+
+    computeStep();
+
+    for (i = 0;i < v.size() ; i++)
+    {
+	index=(int) floor( ((TOut)(v[i]-m_min))/m_step ) ;
+
+	if (index >= m_counters.size() || index < 0)
+	    return;
+
+	m_counters[index]++;
+    }
+}
+
+template<class T, class TOut>
+void THistogram<T,TOut>::update( const std::vector<T>& v)
+{
+    if (m_counters.empty())
+	return;
+
+    computeStep();
+
+    TOut size = m_counters.size();
+
+    int index;
+    for (unsigned int i = 0;i < size ; i++)
+    {
+	index = (int)floor(((TOut)(v[i]-m_min))/m_step);
+
+	if (index >= m_counters.size() || index < 0)
+	    return;
+
+	m_counters[index]++;
+    }
+}
+
+template<class T, class TOut>
+void THistogram<T,TOut>::update( const T& t)
+{	
+    int index=(int) floor( ((TOut)(t-m_min))/m_step ) ;
+
+    if (index >= m_counters.size() || index < 0)
+	return;
+
+    m_counters[index]++;
+};
+
+template<class T, class TOut>
+std::vector<TOut> THistogram<T,TOut>::geTHistogramD() const
+{
+    std::vector<TOut> v(m_counters.size());
+    for (unsigned int i = 0;i<m_counters.size(); i++)
+	v[i]=(TOut)m_counters[i];
+
+    return v;
+}
+
+template <class T, class TOut>
+std::vector<TOut> THistogram<T,TOut>::getLeftContainers() const
+{
+    std::vector<TOut> x( m_counters.size());
+
+    for (unsigned int i = 0;i<m_counters.size(); i++)
+	x[i]= m_min + i*m_step;
+
+    return x;
+}
+
+template <class T, class TOut>
+std::vector<TOut> THistogram<T,TOut>::getCenterContainers() const
+{
+    std::vector<TOut> x( m_counters.size());
+
+    for (unsigned int i = 0;i<m_counters.size(); i++)
+	x[i]= m_min + (i+0.5)*m_step;
+
+    return x;
+}
+
+template <class T, class TOut>
+unsigned int THistogram<T,TOut>::getSum() const
+{
+    unsigned int sum = 0;
+    for (unsigned int i = 0;i<m_counters.size(); i++)
+	sum+=m_counters[i];
+
+    return sum;
+}
+
+template <class T, class TOut>
+TOut THistogram<T,TOut>::getArea() const
+{
+    const size_t n=m_counters.size();
+    TOut area=0;
+
+    if (n>6)
+    {
+	area=3.0/8.0*(m_counters[0]+m_counters[n-1])
+	    +7.0/6.0*(m_counters[1]+m_counters[n-2])
+	    +23.0/24.0*(m_counters[2]+m_counters[n-3]);
+	for (unsigned int i=3;i<n-3;i++)
+	{
+	    area+=m_counters[i];
+	}
+    }
+    else if (n>4)
+    {
+	area=5.0/12.0*(m_counters[0]+m_counters[n-1])
+	    +13.0/12.0*(m_counters[1]+m_counters[n-2]);
+	for (unsigned int i=2;i<n-2;i++)
+	{
+	    area+=m_counters[i];
+	}
+    }
+    else if (n>1)
+    {
+	area=1/2.0*(m_counters[0]+m_counters[n-1]);
+	for (unsigned int i=1;i<n-1;i++)
+	{
+	    area+=m_counters[i];
+	}
+    }
+    else 
+	area=0;
+
+    return area*m_step;
+}
+
+template <class T, class TOut>
+std::vector<TOut> THistogram<T,TOut>::getNormalizedHistogram() const
+{
+    std::vector<TOut> normCounters( m_counters.size());
+    TOut area = (TOut)getArea();
+
+    for (unsigned int i = 0;i<m_counters.size(); i++)
+    {
+	normCounters[i]= (TOut)m_counters[i]/area;
+    }
+
+    return normCounters;
+};
+
+template <class T, class TOut>
+void THistogram<T,TOut>::getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt)
+{
+    int j;
+    double ep=0.0,s,p;
+    const size_t n = data.size();
+
+    if (n <= 1)
+	// nrerror("n must be at least 2 in moment");
+	return;
+
+    s=0.0; // First pass to get the mean.
+    for (j=0;j<n;j++)
+	s += data[j];
+	
+    ave=s/(n);
+    adev=var=skew=kurt=0.0; 
+    /* Second pass to get the first (absolute), second,
+       third, and fourth moments of the
+       deviation from the mean. */
+
+    for (j=0;j<n;j++) 
+    {
+	adev += fabs(s=data[j]-(ave));
+	ep += s;
+	var += (p=s*s);
+	skew += (p *= s);
+	kurt += (p *= s);
+    }
+
+
+    adev /= n;
+    var=(var-ep*ep/n)/(n-1); // Corrected two-pass formula.
+    sdev=sqrt(var); // Put the pieces together according to the conventional definitions. 
+    if (var) 
+    {
+	skew /= (n*(var)*(sdev));
+	kurt=(kurt)/(n*(var)*(var))-3.0;
+    } 
+    else
+	//nrerror("No skew/kurtosis when variance = 0 (in moment)");
+	return;
+}
+
+#endif
+
diff -r 000000000000 -r 49844bc8a895 dsp/maths/MathAliases.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/MathAliases.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,55 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef MATHALIASES_H
+#define MATHALIASES_H
+
+#include <cmath>
+#include <complex>
+
+using namespace std;
+typedef complex<double> ComplexData;
+
+
+#ifndef PI
+#define PI (3.14159265358979232846)
+#endif
+
+#define TWO_PI 		(*2.PI)
+
+#define EPS 2.2204e-016
+
+/* aliases to math.h functions */
+#define EXP				exp
+#define COS				cos
+#define SIN				sin
+#define ABS				fabs
+#define POW				powf
+#define SQRT			sqrtf
+#define LOG10			log10f
+#define LOG				logf
+#define FLOOR			floorf
+#define TRUNC			truncf
+
+/* aliases to complex.h functions */
+/** sample = EXPC(complex) */
+#define EXPC			cexpf
+/** complex = CEXPC(complex) */
+#define CEXPC			cexp
+/** sample = ARGC(complex) */
+#define ARGC			cargf
+/** sample = ABSC(complex) norm */
+#define ABSC			cabsf
+/** sample = REAL(complex) */
+#define REAL			crealf
+/** sample = IMAG(complex) */
+#define IMAG			cimagf
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/maths/MathUtilities.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/MathUtilities.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,170 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "MathUtilities.h"
+
+#include <iostream>
+#include <cmath>
+
+
+double MathUtilities::mod(double x, double y)
+{
+    double a = floor( x / y );
+
+    double b = x - ( y * a );
+    return b;
+}
+
+double MathUtilities::princarg(double ang)
+{
+    double ValOut;
+
+    ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
+
+    return ValOut;
+}
+
+void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
+{
+    unsigned int i;
+    double temp = 0.0;
+    double a=0.0;
+	
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+		
+	a  += ::pow( fabs(temp), alpha );
+    }
+    a /= ( double )len;
+    a = ::pow( a, ( 1.0 / (double) alpha ) );
+
+    *ANorm = a;
+}
+
+double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
+{
+    unsigned int i;
+    unsigned int len = data.size();
+    double temp = 0.0;
+    double a=0.0;
+	
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+		
+	a  += ::pow( fabs(temp), alpha );
+    }
+    a /= ( double )len;
+    a = ::pow( a, ( 1.0 / (double) alpha ) );
+
+    return a;
+}
+
+double MathUtilities::round(double x)
+{
+    double val = (double)floor(x + 0.5);
+  
+    return val;
+}
+
+double MathUtilities::median(const double *src, unsigned int len)
+{
+    unsigned int i, j;
+    double tmp = 0.0;
+    double tempMedian;
+    double medianVal;
+ 
+    double* scratch = new double[ len ];//Vector < double > sortedX = Vector < double > ( size );
+
+    for ( i = 0; i < len; i++ )
+    {
+	scratch[i] = src[i];
+    }
+
+    for ( i = 0; i < len - 1; i++ )
+    {
+	for ( j = 0; j < len - 1 - i; j++ )
+	{
+	    if ( scratch[j + 1] < scratch[j] )
+	    {
+		// compare the two neighbors
+		tmp = scratch[j]; // swap a[j] and a[j+1]
+		scratch[j] = scratch[j + 1];
+		scratch[j + 1] = tmp;
+	    }
+	}
+    }
+    int middle;
+    if ( len % 2 == 0 )
+    {
+	middle = len / 2;
+	tempMedian = ( scratch[middle] + scratch[middle - 1] ) / 2;
+    }
+    else
+    {
+	middle = ( int )floor( len / 2.0 );
+	tempMedian = scratch[middle];
+    }
+
+    medianVal = tempMedian;
+
+    delete [] scratch;
+    return medianVal;
+}
+
+double MathUtilities::sum(const double *src, unsigned int len)
+{
+    unsigned int i ;
+    double retVal =0.0;
+
+    for(  i = 0; i < len; i++)
+    {
+	retVal += src[ i ];
+    }
+
+    return retVal;
+}
+
+double MathUtilities::mean(const double *src, unsigned int len)
+{
+    double retVal =0.0;
+
+    double s = sum( src, len );
+	
+    retVal =  s  / (double)len;
+
+    return retVal;
+}
+
+void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
+{
+    unsigned int i;
+    double temp = 0.0;
+    double a=0.0;
+	
+    *min = data[0];
+    *max = data[0];
+
+    for( i = 0; i < len; i++)
+    {
+	temp = data[ i ];
+
+	if( temp < *min )
+	{
+	    *min =  temp ;
+	}
+	if( temp > *max )
+	{
+	    *max =  temp ;
+	}
+		
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/maths/MathUtilities.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/MathUtilities.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,30 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef MATHUTILITIES_H
+#define MATHUTILITIES_H
+
+#include <vector>
+
+class MathUtilities  
+{
+public:	
+    static double round( double x );
+    static void	  getFrameMinMax( const double* data, unsigned int len,  double* min, double* max );
+    static double mean( const double* src, unsigned int len );
+    static double sum( const double* src, unsigned int len );
+    static double princarg( double ang );
+    static double median( const double* src, unsigned int len );
+    static double mod( double x, double y);
+    static void	  getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm);
+    static double getAlphaNorm(const std::vector <double> &data, unsigned int alpha );
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/maths/Polyfit.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/maths/Polyfit.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,398 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+//---------------------------------------------------------------------------
+
+#ifndef PolyfitHPP
+#define PolyfitHPP
+//---------------------------------------------------------------------------
+// Least-squares curve fitting class for arbitrary data types
+/*
+
+{ ******************************************
+  ****   Scientific Subroutine Library  ****
+  ****         for C++ Builder          ****
+  ******************************************
+
+  The following programs were written by Allen Miller and appear in the
+  book entitled "Pascal Programs For Scientists And Engineers" which is
+  published by Sybex, 1981.
+  They were originally typed and submitted to MTPUG in Oct. 1982
+    Juergen Loewner
+    Hoher Heckenweg 3
+    D-4400 Muenster
+  They have had minor corrections and adaptations for Turbo Pascal by
+    Jeff Weiss
+    1572 Peacock Ave.
+    Sunnyvale, CA 94087.
+
+
+2000 Oct 28  Updated for Delphi 4, open array parameters.
+             This allows the routine to be generalised so that it is no longer
+             hard-coded to make a specific order of best fit or work with a
+             specific number of points.
+2001 Jan 07  Update Web site address
+
+Copyright © David J Taylor, Edinburgh and others listed above
+Web site:  www.satsignal.net
+E-mail:    davidtaylor@writeme.com
+}*/
+
+ ///////////////////////////////////////////////////////////////////////////////
+ // Modified by CLandone for VC6 Aug 2004
+ ///////////////////////////////////////////////////////////////////////////////
+
+
+using std::vector;
+
+
+class TPolyFit
+{
+    typedef vector<vector<double> > Matrix;
+public:
+
+    static double PolyFit2 (const vector<double> &x,  // does the work
+			    const vector<double> &y,
+			    vector<double> &coef);
+
+                   
+private:
+    TPolyFit &operator = (const TPolyFit &);   // disable assignment
+    TPolyFit();                                // and instantiation
+    TPolyFit(const TPolyFit&);                 // and copying
+
+  
+    static void Square (const Matrix &x,              // Matrix multiplication routine
+			const vector<double> &y,
+			Matrix &a,                    // A = transpose X times X
+			vector<double> &g,         // G = Y times X
+			const int nrow, const int ncol);
+    // Forms square coefficient matrix
+
+    static bool GaussJordan (Matrix &b,                  // square matrix of coefficients
+			     const vector<double> &y, // constant vector
+			     vector<double> &coef);   // solution vector
+    // returns false if matrix singular
+
+    static bool GaussJordan2(Matrix &b,
+			     const vector<double> &y,
+			     Matrix &w,
+			     vector<vector<int> > &index);
+};
+
+// some utility functions
+
+namespace NSUtility
+{
+    inline void swap(double &a, double &b) {double t = a; a = b; b = t;}
+    void zeroise(vector<double> &array, int n);
+    void zeroise(vector<int> &array, int n);
+    void zeroise(vector<vector<double> > &matrix, int m, int n);
+    void zeroise(vector<vector<int> > &matrix, int m, int n);
+    inline double sqr(const double &x) {return x * x;}
+};
+
+//---------------------------------------------------------------------------
+// Implementation
+//---------------------------------------------------------------------------
+using namespace NSUtility;
+//------------------------------------------------------------------------------------------
+
+
+// main PolyFit routine
+
+double TPolyFit::PolyFit2 (const vector<double> &x,
+			   const vector<double> &y,
+			   vector<double> &coefs)
+// nterms = coefs.size()
+// npoints = x.size()
+{
+    int i, j;
+    double xi, yi, yc, srs, sum_y, sum_y2;
+    Matrix xmatr;        // Data matrix
+    Matrix a;
+    vector<double> g;      // Constant vector
+    const int npoints(x.size());
+    const int nterms(coefs.size());
+    double correl_coef;
+    zeroise(g, nterms);
+    zeroise(a, nterms, nterms);
+    zeroise(xmatr, npoints, nterms);
+    if(nterms < 1)
+	throw "PolyFit called with less than one term";
+    if(npoints < 2)
+	throw "PolyFit called with less than two points";
+    if(npoints != y.size())
+	throw "PolyFit called with x and y of unequal size";
+    for(i = 0; i < npoints; ++i)
+    {
+	//      { setup x matrix }
+	xi = x[i];
+	xmatr[i][0] = 1.0;	   //     { first column }
+	for(j = 1; j < nterms; ++j)
+	    xmatr[i][j] = xmatr [i][j - 1] * xi;
+    }
+    Square (xmatr, y, a, g, npoints, nterms);
+    if(!GaussJordan (a, g, coefs))
+	return -1;
+    sum_y = 0.0;
+    sum_y2 = 0.0;
+    srs = 0.0;
+    for(i = 0; i < npoints; ++i)
+    {
+	yi = y[i];
+	yc = 0.0;
+	for(j = 0; j < nterms; ++j)
+	    yc += coefs [j] * xmatr [i][j];
+	srs += sqr (yc - yi);
+	sum_y += yi;
+	sum_y2 += yi * yi;
+    }
+
+    // If all Y values are the same, avoid dividing by zero
+    correl_coef = sum_y2 - sqr (sum_y) / npoints;
+    // Either return 0 or the correct value of correlation coefficient
+    if (correl_coef != 0)
+	correl_coef = srs / correl_coef;
+    if (correl_coef >= 1)
+	correl_coef = 0.0;
+    else
+	correl_coef = sqrt (1.0 - correl_coef);
+    return correl_coef;
+}
+
+
+//------------------------------------------------------------------------
+
+// Matrix multiplication routine
+// A = transpose X times X
+// G = Y times X
+
+// Form square coefficient matrix
+
+void TPolyFit::Square (const Matrix &x,
+		       const vector<double> &y,
+		       Matrix &a,
+		       vector<double> &g,
+		       const int nrow,
+		       const int ncol)
+{
+    int i, k, l;
+    for(k = 0; k < ncol; ++k)
+    {
+	for(l = 0; l < k + 1; ++l)
+	{
+	    a [k][l] = 0.0;
+	    for(i = 0; i < nrow; ++i)
+	    {
+		a[k][l] += x[i][l] * x [i][k];
+		if(k != l)
+		    a[l][k] = a[k][l];
+	    }
+	}
+	g[k] = 0.0;
+	for(i = 0; i < nrow; ++i)
+	    g[k] += y[i] * x[i][k];
+    }
+}
+//---------------------------------------------------------------------------------
+
+
+bool TPolyFit::GaussJordan (Matrix &b,
+			    const vector<double> &y,
+			    vector<double> &coef)
+//b square matrix of coefficients
+//y constant vector
+//coef solution vector
+//ncol order of matrix got from b.size()
+
+
+{
+/*
+  { Gauss Jordan matrix inversion and solution }
+
+  { B (n, n) coefficient matrix becomes inverse }
+  { Y (n) original constant vector }
+  { W (n, m) constant vector(s) become solution vector }
+  { DETERM is the determinant }
+  { ERROR = 1 if singular }
+  { INDEX (n, 3) }
+  { NV is number of constant vectors }
+*/
+
+    int ncol(b.size());
+    int irow, icol;
+    vector<vector<int> >index;
+    Matrix w;
+
+    zeroise(w, ncol, ncol);
+    zeroise(index, ncol, 3);
+
+    if(!GaussJordan2(b, y, w, index))
+	return false;
+
+    // Interchange columns
+    int m;
+    for (int i = 0; i <  ncol; ++i)
+    {
+	m = ncol - i - 1;
+	if(index [m][0] != index [m][1])
+	{
+	    irow = index [m][0];
+	    icol = index [m][1];
+	    for(int k = 0; k < ncol; ++k)
+		swap (b[k][irow], b[k][icol]);
+	}
+    }
+
+    for(int k = 0; k < ncol; ++k)
+    {
+	if(index [k][2] != 0)
+	{
+	    throw "Error in GaussJordan: matrix is singular";
+	}
+    }
+
+    for( int i = 0; i < ncol; ++i)
+	coef[i] = w[i][0];
+ 
+ 
+    return true;
+}   // end;	{ procedure GaussJordan }
+//----------------------------------------------------------------------------------------------
+
+
+bool TPolyFit::GaussJordan2(Matrix &b,
+			    const vector<double> &y,
+			    Matrix &w,
+			    vector<vector<int> > &index)
+{
+    //GaussJordan2;         // first half of GaussJordan
+    // actual start of gaussj
+ 
+    double big, t;
+    double pivot;
+    double determ;
+    int irow, icol;
+    int ncol(b.size());
+    int nv = 1;                  // single constant vector
+    for(int i = 0; i < ncol; ++i)
+    {
+	w[i][0] = y[i];      // copy constant vector
+	index[i][2] = -1;
+    }
+    determ = 1.0;
+    for(int i = 0; i < ncol; ++i)
+    {
+	// Search for largest element
+	big = 0.0;
+	for(int j = 0; j < ncol; ++j)
+	{
+	    if(index[j][2] != 0)
+	    {
+		for(int k = 0; k < ncol; ++k)
+		{
+		    if(index[k][2] > 0)
+			throw "Error in GaussJordan2: matrix is singular";
+
+		    if(index[k][2] < 0 && fabs(b[j][k]) > big)
+		    {
+			irow = j;
+			icol = k;
+			big = fabs(b[j][k]);
+		    }
+		} //   { k-loop }
+	    }
+	}  // { j-loop }
+	index [icol][2] = index [icol][2] + 1;
+	index [i][0] = irow;
+	index [i][1] = icol;
+
+	// Interchange rows to put pivot on diagonal
+	// GJ3
+	if(irow != icol)
+	{
+	    determ = -determ;
+	    for(int m = 0; m < ncol; ++m)
+		swap (b [irow][m], b[icol][m]);
+	    if (nv > 0)
+		for (int m = 0; m < nv; ++m)
+		    swap (w[irow][m], w[icol][m]);
+	} // end GJ3
+
+	// divide pivot row by pivot column
+	pivot = b[icol][icol];
+	determ *= pivot;
+	b[icol][icol] = 1.0;
+
+	for(int m = 0; m < ncol; ++m)
+	    b[icol][m] /= pivot;
+	if(nv > 0)
+	    for(int m = 0; m < nv; ++m)
+		w[icol][m] /= pivot;
+
+	// Reduce nonpivot rows
+	for(int n = 0; n < ncol; ++n)
+	{
+	    if(n != icol)
+	    {
+		t = b[n][icol];
+		b[n][icol] = 0.0;
+		for(int m = 0; m < ncol; ++m)
+		    b[n][m] -= b[icol][m] * t;
+		if(nv > 0)
+		    for(int m = 0; m < nv; ++m)
+			w[n][m] -= w[icol][m] * t;
+	    }
+	}
+    } // { i-loop }
+    return true;
+}
+//----------------------------------------------------------------------------------------------
+
+//------------------------------------------------------------------------------------
+
+// Utility functions
+//--------------------------------------------------------------------------
+
+// fills a vector with zeros.
+void NSUtility::zeroise(vector<double> &array, int n)
+{
+    array.clear();
+    for(int j = 0; j < n; ++j)
+	array.push_back(0);
+}
+//--------------------------------------------------------------------------
+
+// fills a vector with zeros.
+void NSUtility::zeroise(vector<int> &array, int n)
+{
+    array.clear();
+    for(int j = 0; j < n; ++j)
+	array.push_back(0);
+}
+//--------------------------------------------------------------------------
+
+// fills a (m by n) matrix with zeros.
+void NSUtility::zeroise(vector<vector<double> > &matrix, int m, int n)
+{
+    vector<double> zero;
+    zeroise(zero, n);
+    matrix.clear();
+    for(int j = 0; j < m; ++j)
+	matrix.push_back(zero);
+}
+//--------------------------------------------------------------------------
+
+// fills a (m by n) matrix with zeros.
+void NSUtility::zeroise(vector<vector<int> > &matrix, int m, int n)
+{
+    vector<int> zero;
+    zeroise(zero, n);
+    matrix.clear();
+    for(int j = 0; j < m; ++j)
+	matrix.push_back(zero);
+}
+//--------------------------------------------------------------------------
+
+
+#endif
+ 
diff -r 000000000000 -r 49844bc8a895 dsp/onsets/DetectionFunction.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/onsets/DetectionFunction.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,204 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "DetectionFunction.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+DetectionFunction::DetectionFunction( DFConfig Config ) :
+    m_window(0)
+{
+    magHistory = NULL;
+    phaseHistory = NULL;
+    phaseHistoryOld = NULL;
+    j = ComplexData( 0, 1 );
+
+    initialise( Config );
+}
+
+DetectionFunction::~DetectionFunction()
+{
+    deInitialise();
+}
+
+
+void DetectionFunction::initialise( DFConfig Config )
+{
+    m_dataLength = Config.frameLength;
+    m_halfLength = m_dataLength/2;
+    m_DFType = Config.DFType;
+
+    magHistory = new double[ m_halfLength ];
+    memset(magHistory,0, m_halfLength*sizeof(double));
+		
+    phaseHistory = new double[ m_halfLength ];
+    memset(phaseHistory,0, m_halfLength*sizeof(double));
+
+    phaseHistoryOld = new double[ m_halfLength ];
+    memset(phaseHistoryOld,0, m_halfLength*sizeof(double));
+
+    m_phaseVoc = new PhaseVocoder;
+
+    m_DFWindowedFrame = new double[ m_dataLength ];
+    m_magnitude = new double[ m_halfLength ];
+    m_thetaAngle = new double[ m_halfLength ];
+
+    m_window = new Window<double>(HanningWindow, m_dataLength);
+}
+
+void DetectionFunction::deInitialise()
+{
+    delete [] magHistory ;
+    delete [] phaseHistory ;
+    delete [] phaseHistoryOld ;
+
+    delete m_phaseVoc;
+
+    delete [] m_DFWindowedFrame;
+    delete [] m_magnitude;
+    delete [] m_thetaAngle;
+
+    delete m_window;
+}
+
+double DetectionFunction::process( double *TDomain )
+{
+    double retVal = 0;
+
+    m_window->cut( TDomain, m_DFWindowedFrame );
+	
+    m_phaseVoc->process( m_dataLength, m_DFWindowedFrame, m_magnitude, m_thetaAngle );
+
+    switch( m_DFType )
+    {
+    case DF_HFC:
+	retVal = HFC( m_halfLength, m_magnitude);
+	break;
+	
+    case  DF_SPECDIFF:
+	retVal = specDiff( m_halfLength, m_magnitude);
+	break;
+	
+    case DF_PHASEDEV:
+	retVal = phaseDev( m_halfLength, m_magnitude, m_thetaAngle);
+	break;
+	
+    case DF_COMPLEXSD:
+	retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
+	break;
+    }
+	
+    return retVal;
+}
+
+double DetectionFunction::HFC(unsigned int length, double *src)
+{
+    unsigned int i;
+    double val = 0;
+
+    for( i = 0; i < length; i++)
+    {
+	val += src[ i ] * ( i + 1);
+    }
+    return val;
+}
+
+double DetectionFunction::specDiff(unsigned int length, double *src)
+{
+    unsigned int i;
+    double val = 0.0;
+    double temp = 0.0;
+    double diff = 0.0;
+
+    for( i = 0; i < length; i++)
+    {
+	temp = fabs( (src[ i ] * src[ i ]) - (magHistory[ i ] * magHistory[ i ]) );
+		
+	diff= sqrt(temp);
+
+	if( src[ i ] > 0.1)
+	{
+	    val += diff;
+	}
+
+	magHistory[ i ] = src[ i ];
+    }
+
+    return val;
+}
+
+
+double DetectionFunction::phaseDev(unsigned int length, double *srcMagnitude, double *srcPhase)
+{
+    unsigned int i;
+    double tmpPhase = 0;
+    double tmpVal = 0;
+    double val = 0;
+
+    double dev = 0;
+
+    for( i = 0; i < length; i++)
+    {
+	tmpPhase = (srcPhase[ i ]- 2*phaseHistory[ i ]+phaseHistoryOld[ i ]);
+	dev = MathUtilities::princarg( tmpPhase );
+		
+	if( srcMagnitude[ i  ] > 0.1)
+	{
+	    tmpVal  = fabs( dev);
+	    val += tmpVal ;
+	}
+
+	phaseHistoryOld[ i ] = phaseHistory[ i ] ;
+	phaseHistory[ i ] = srcPhase[ i ];
+    }
+	
+	
+    return val;
+}
+
+
+double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
+{
+    unsigned int i;
+    double val = 0;
+    double tmpPhase = 0;
+    double tmpReal = 0;
+    double tmpImag = 0;
+   
+    double dev = 0;
+    ComplexData meas = ComplexData( 0, 0 );
+
+    for( i = 0; i < length; i++)
+    {
+	tmpPhase = (srcPhase[ i ]- 2*phaseHistory[ i ]+phaseHistoryOld[ i ]);
+	dev= MathUtilities::princarg( tmpPhase );
+		
+	meas = magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
+
+	tmpReal = real( meas );
+	tmpImag = imag( meas );
+
+	val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
+		
+	phaseHistoryOld[ i ] = phaseHistory[ i ] ;
+	phaseHistory[ i ] = srcPhase[ i ];
+	magHistory[ i ] = srcMagnitude[ i ];
+    }
+
+    return val;
+}
+
+double* DetectionFunction::getSpectrumMagnitude()
+{
+    return m_magnitude;
+}
+
diff -r 000000000000 -r 49844bc8a895 dsp/onsets/DetectionFunction.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/onsets/DetectionFunction.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,72 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef DETECTIONFUNCTION_H
+#define DETECTIONFUNCTION_H
+
+#include "dsp/maths/MathUtilities.h"
+#include "dsp/maths/MathAliases.h"
+#include "dsp/phasevocoder/PhaseVocoder.h"
+#include "base/Window.h"
+
+#define DF_HFC (1)
+#define DF_SPECDIFF (2)
+#define DF_PHASEDEV (3)
+#define DF_COMPLEXSD (4)
+
+struct DFConfig{
+    double stepSecs; // DF step in seconds
+    unsigned int stepSize; // DF step in samples
+    unsigned int frameLength; // DF analysis window - usually 2*step
+    int DFType; // type of detection function ( see defines )
+};
+
+class DetectionFunction  
+{
+public:
+    double* getSpectrumMagnitude();
+    DetectionFunction( DFConfig Config );
+    virtual ~DetectionFunction();
+    double process( double* TDomain );
+
+private:
+    double HFC( unsigned int length, double* src);
+    double specDiff( unsigned int length, double* src);
+    double phaseDev(unsigned int length, double *srcMagnitude, double *srcPhase);
+    double complexSD(unsigned int length, double *srcMagnitude, double *srcPhase);
+	
+private:
+    void initialise( DFConfig Config );
+    void deInitialise();
+
+    int m_DFType;
+    unsigned int m_dataLength;
+    unsigned int m_halfLength;
+
+    double* magHistory;
+    double* phaseHistory;
+    double* phaseHistoryOld;
+
+    double* m_DFWindowedFrame; // Array for windowed analysis frame
+    double* m_magnitude; // Magnitude of analysis frame ( frequency domain )
+    double* m_thetaAngle;// Phase of analysis frame ( frequency domain )
+
+
+    vector < ComplexData > meas ;
+	
+    ComplexData j;
+
+    Window<double> *m_window;
+
+    PhaseVocoder*	m_phaseVoc;	// Phase Vocoder
+
+};
+
+#endif 
diff -r 000000000000 -r 49844bc8a895 dsp/onsets/PeakPicking.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/onsets/PeakPicking.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,137 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "PeakPicking.h"
+#include "dsp/maths/Polyfit.h"
+
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+PeakPicking::PeakPicking( PPickParams Config )
+{
+    m_workBuffer = NULL;
+    initialise( Config );
+}
+
+PeakPicking::~PeakPicking()
+{
+    deInitialise();
+}
+
+void PeakPicking::initialise( PPickParams Config )
+{
+    m_DFLength = Config.length ;
+    Qfilta = Config.QuadThresh.a ;
+    Qfiltb = Config.QuadThresh.b ;
+    Qfiltc = Config.QuadThresh.c ;
+	
+    m_DFProcessingParams.length = m_DFLength; 
+    m_DFProcessingParams.LPOrd = Config.LPOrd; 
+    m_DFProcessingParams.LPACoeffs = Config.LPACoeffs; 
+    m_DFProcessingParams.LPBCoeffs = Config.LPBCoeffs; 
+    m_DFProcessingParams.winPre  = Config.WinT.pre;
+    m_DFProcessingParams.winPost = Config.WinT.post; 
+    m_DFProcessingParams.AlphaNormParam = Config.alpha;
+    m_DFProcessingParams.isMedianPositive = false;
+	
+    m_DFSmoothing = new DFProcess( m_DFProcessingParams );
+
+    m_workBuffer = new double[ m_DFLength ];
+    memset( m_workBuffer, 0, sizeof(double)*m_DFLength);
+}
+
+void PeakPicking::deInitialise()
+{
+    delete [] m_workBuffer;
+    delete m_DFSmoothing;
+    m_workBuffer = NULL;
+}
+
+void PeakPicking::process( double* src, unsigned int len, vector<int> &onsets )
+{
+    vector <double> m_maxima;	
+
+    // Signal conditioning 
+    m_DFSmoothing->process( src, m_workBuffer );
+	
+    for( unsigned int u = 0; u < len; u++)
+    {
+	m_maxima.push_back( m_workBuffer[ u ] );		
+    }
+	
+    quadEval( m_maxima, onsets );
+
+    for( int b = 0; b <  m_maxima.size(); b++)
+    {
+	src[ b ] = m_maxima[ b ];
+    }
+}
+
+int PeakPicking::quadEval( vector<double> &src, vector<int> &idx )
+{
+    unsigned int maxLength;
+
+    vector <int> m_maxIndex;
+    vector <int> m_onsetPosition;
+	
+    vector <double> m_maxFit;
+    vector <double> m_poly;
+    vector <double> m_err;
+
+    double p;
+
+    m_poly.push_back(0);
+    m_poly.push_back(0);
+    m_poly.push_back(0);
+
+    for(  int t = -2; t < 3; t++)
+    {
+	m_err.push_back( (double)t );
+    }
+    for( unsigned int i = 2; i < src.size() - 2; i++)
+    {
+	if( (src[ i ] > src[ i  - 1 ]) && (src[ i ] > src[ i  + 1 ]) && (src[ i ] > 0) )
+	{
+	    m_maxIndex.push_back(  i + 1 );
+	}
+    }
+
+    maxLength = m_maxIndex.size();
+
+    double selMax = 0;
+
+    for( unsigned int j = 0; j < maxLength ; j++)
+    {
+	for(  int k = -3; k < 2; k++)
+	{
+	    selMax = src[ m_maxIndex[j] + k ] ;
+	    m_maxFit.push_back(selMax);			
+	}
+
+	p = TPolyFit::PolyFit2( m_err, m_maxFit, m_poly);
+
+	double f = m_poly[0];
+	double g = m_poly[1];
+	double h = m_poly[2];
+
+	int kk = m_poly.size();
+		
+	if( h < -Qfilta || f  >  Qfiltc)
+	{
+	    idx.push_back( m_maxIndex[j] );
+	}
+
+	m_maxFit.erase( m_maxFit.begin(), m_maxFit.end() );
+    }
+
+    return 1;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/onsets/PeakPicking.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/onsets/PeakPicking.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,77 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+// PeakPicking.h: interface for the PeakPicking class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#ifndef PEAKPICKING_H
+#define PEAKPICKING_H
+
+#include "dsp/maths/MathUtilities.h"
+#include "dsp/maths/MathAliases.h"
+#include "dsp/signalconditioning/DFProcess.h"
+
+
+struct WinThresh
+{
+    unsigned int pre;
+    unsigned int  post;
+};
+
+struct QFitThresh
+{
+    double a;
+    double b;
+    double c;
+};
+
+struct PPickParams
+{
+    unsigned int length; //Detection FunctionLength
+    double tau; // time resolution of the detection function:
+    unsigned int alpha; //alpha-norm parameter
+    double cutoff;//low-pass Filter cutoff freq
+    unsigned int LPOrd; // low-pass Filter order
+    double* LPACoeffs; //low pass Filter den coefficients
+    double* LPBCoeffs; //low pass Filter num coefficients
+    WinThresh WinT;//window size in frames for adaptive thresholding [pre post]:
+    QFitThresh QuadThresh;
+};
+
+class PeakPicking  
+{
+public:
+    PeakPicking( PPickParams Config );
+    virtual ~PeakPicking();
+	
+    void process( double* src, unsigned int len, vector<int> &onsets  );
+
+
+private:
+    void initialise( PPickParams Config  );
+    void deInitialise();
+    int  quadEval( vector<double> &src, vector<int> &idx );
+	
+    DFProcConfig m_DFProcessingParams;
+
+    unsigned int m_DFLength ;
+    double m_alphaNormParam ;
+    double Qfilta ;
+    double Qfiltb;
+    double Qfiltc;
+
+
+    double* m_workBuffer;
+	
+    DFProcess*	m_DFSmoothing;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/phasevocoder/PhaseVocoder.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/phasevocoder/PhaseVocoder.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,97 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "PhaseVocoder.h"
+#include "dsp/transforms/FFT.h"
+#include <math.h>
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+PhaseVocoder::PhaseVocoder()
+{
+
+}
+
+PhaseVocoder::~PhaseVocoder()
+{
+
+}
+
+void PhaseVocoder::FFTShift(unsigned int size, double *src)
+{
+    // IN-place Rotation of FFT arrays
+    unsigned int i;
+
+    shiftBuffer = new double[size/2];
+
+    for( i = 0; i < size/2; i++)
+    {
+	shiftBuffer[ i ] = src[ i ];
+	src[ i ] = src[ i + size/2];
+    }
+	
+    for( i =size/2; i < size;  i++)
+    {
+	src[ i ] = shiftBuffer[ i -(size/2)];
+    }
+
+    delete [] shiftBuffer;
+
+}
+
+void PhaseVocoder::process(unsigned int size, double *src, double *mag, double *theta)
+{
+
+    // Primary Interface to Phase Vocoder
+    realOut = new double[ size ];
+    imagOut = new double[ size ];
+
+    FFTShift( size, src);
+	
+    coreFFT( size, src, 0, realOut, imagOut);
+
+    getMagnitude( size/2, mag, realOut, imagOut);
+    getPhase( size/2, theta, realOut, imagOut);
+
+    delete [] realOut;
+    delete [] imagOut;
+}
+
+
+void PhaseVocoder::coreFFT( unsigned int NumSamples, double *RealIn, double* ImagIn, double *RealOut, double *ImagOut)
+{
+    // This function is taken from a standard freeware implementation defined in FFT.h
+    // TODO: Use  FFTW
+    FFT::process( NumSamples,0, RealIn, ImagIn, RealOut, ImagOut );
+}
+
+void PhaseVocoder::getMagnitude(unsigned int size, double *mag, double *real, double *imag)
+{	
+    unsigned int j;
+
+    for( j = 0; j < size; j++)
+    {
+	mag[ j ] = sqrt( real[ j ] * real[ j ] + imag[ j ] * imag[ j ]);
+    }
+}
+
+void PhaseVocoder::getPhase(unsigned int size, double *theta, double *real, double *imag)
+{
+    unsigned int k;
+
+    // Phase Angle "matlab" style 
+    //Watch out for quadrant mapping  !!!
+    for( k = 0; k < size; k++)
+    {
+	theta[ k ] = atan2( -imag[ k ], real[ k ]);
+    }	
+}
diff -r 000000000000 -r 49844bc8a895 dsp/phasevocoder/PhaseVocoder.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/phasevocoder/PhaseVocoder.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,35 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef PHASEVOCODER_H
+#define PHASEVOCODER_H
+
+
+class PhaseVocoder  
+{
+public:
+    PhaseVocoder();
+    virtual ~PhaseVocoder();
+
+    void process( unsigned int size, double* src, double* mag, double* theta);
+    void FFTShift( unsigned int size, double* src);
+
+protected:
+    void getPhase(unsigned int size, double *theta, double *real, double *imag);
+    void coreFFT( unsigned int NumSamples, double *RealIn, double* ImagIn, double *RealOut, double *ImagOut);
+    void getMagnitude( unsigned int size, double* mag, double* real, double* imag);
+
+    double* shiftBuffer;
+    double* imagOut;
+    double* realOut;
+
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/rateconversion/Decimator.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/rateconversion/Decimator.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,154 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "Decimator.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Decimator::Decimator( unsigned int inLength, unsigned int decFactor )
+{
+
+    m_inputLength = 0;
+    m_outputLength = 0;
+    m_decFactor = 1;
+
+    initialise( inLength, decFactor );
+}
+
+Decimator::~Decimator()
+{
+    deInitialise();
+}
+
+void Decimator::initialise( unsigned int inLength, unsigned int decFactor)
+{
+    m_inputLength = inLength;
+    m_decFactor = decFactor;
+    m_outputLength = m_inputLength / m_decFactor;
+
+    decBuffer = new double[ m_inputLength ];
+
+    if( m_decFactor == 4 )
+    {
+	//////////////////////////////////////////////////
+	b[ 0 ] = 0.10133306904918619;
+	b[ 1 ] = -0.2447523353702363;
+	b[ 2 ] = 0.33622528590120965;
+	b[ 3 ] = -0.13936581560633518;
+	b[ 4 ] = -0.13936581560633382;
+	b[ 5 ] = 0.3362252859012087;
+	b[ 6 ] = -0.2447523353702358;
+	b[ 7 ] = 0.10133306904918594;
+
+	a[ 0 ] = 1;
+	a[ 1 ] = -3.9035590278139427;
+	a[ 2 ] = 7.5299379980621133;
+	a[ 3 ] = -8.6890803793177511;
+	a[ 4 ] = 6.4578667096099176;
+	a[ 5 ] = -3.0242979431223631;
+	a[ 6 ] = 0.83043385136748382;
+	a[ 7 ] = -0.094420800837809335;
+	//////////////////////////////////////////////////
+    }
+    else if( m_decFactor == 2 )
+    {
+	//////////////////////////////////////////////////
+	b[ 0 ] = 0.20898944260075727;
+	b[ 1 ] = 0.40011234879814367;
+	b[ 2 ] = 0.819741973072733;
+	b[ 3 ] = 1.0087419911682323;
+	b[ 4 ] = 1.0087419911682325;
+	b[ 5 ] = 0.81974197307273156;
+	b[ 6 ] = 0.40011234879814295;
+	b[ 7 ] = 0.20898944260075661;
+
+	a[ 0 ] = 1;
+	a[ 1 ] = 0.0077331184208358217;
+	a[ 2 ] = 1.9853971155964376;
+	a[ 3 ] = 0.19296739275341004;
+	a[ 4 ] = 1.2330748872852182;
+	a[ 5 ] = 0.18705341389316466;
+	a[ 6 ] = 0.23659265908013868;
+	a[ 7 ] = 0.032352924250533946;
+    }
+    else
+    {
+	//////////////////////////////////////////////////
+	b[ 0 ] = 1;
+	b[ 1 ] = 0;
+	b[ 2 ] = 0;
+	b[ 3 ] = 0;
+	b[ 4 ] = 0;
+	b[ 5 ] = 0;
+	b[ 6 ] = 0;
+	b[ 7 ] = 0;
+
+	a[ 0 ] = 1;
+	a[ 1 ] = 0;
+	a[ 2 ] = 0;
+	a[ 3 ] = 0;
+	a[ 4 ] = 0;
+	a[ 5 ] = 0;
+	a[ 6 ] = 0;
+	a[ 7 ] = 0;
+    }
+
+    resetFilter();
+}
+
+void Decimator::deInitialise()
+{
+    delete [] decBuffer;
+}
+
+void Decimator::resetFilter()
+{
+    Input = Output = 0;
+
+    o1=o2=o3=o4=o5=o6=o7=0;
+}
+
+void Decimator::doAntiAlias(double *src, double *dst, unsigned int length)
+{
+
+    for( unsigned int i = 0; i < length; i++ )
+    {
+	Input = (double)src[ i ];
+
+	Output = Input * b[ 0 ] + o1;
+
+	o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
+	o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
+	o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
+	o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
+	o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
+	o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
+	o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
+
+	dst[ i ] = Output;
+    }
+
+}
+
+void Decimator::process(double *src, double *dst)
+{
+    if( m_decFactor != 1 )
+    {
+	doAntiAlias( src, decBuffer, m_inputLength );
+    }
+    unsigned idx = 0;
+
+    for( unsigned int i = 0; i < m_outputLength; i++ )
+    {
+	dst[ idx++ ] = decBuffer[ m_decFactor * i ];
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/rateconversion/Decimator.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/rateconversion/Decimator.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,42 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef DECIMATOR_H
+#define DECIMATOR_H
+
+class Decimator  
+{
+public:
+    void process( double* src, double* dst );
+    void doAntiAlias( double* src, double* dst, unsigned int length );
+
+    Decimator( unsigned int inLength, unsigned int decFactor );
+    virtual ~Decimator();
+
+private:
+    void resetFilter();
+    void deInitialise();
+    void initialise( unsigned int inLength, unsigned int decFactor );
+
+    unsigned int m_inputLength;
+    unsigned int m_outputLength;
+    unsigned int m_decFactor;
+
+    double Input;
+    double Output ;
+
+    double o1,o2,o3,o4,o5,o6,o7;
+
+    double a[ 9 ];
+    double b[ 9 ];
+	
+    double* decBuffer;
+};
+
+#endif // 
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/DFProcess.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/DFProcess.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,172 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "DFProcess.h"
+#include "dsp/maths/MathUtilities.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+DFProcess::DFProcess( DFProcConfig Config )
+{
+    filtSrc = NULL;
+    filtDst = NULL;	
+    m_filtScratchIn = NULL;
+    m_filtScratchOut = NULL;
+
+    m_FFOrd = 0;
+
+    initialise( Config );
+}
+
+DFProcess::~DFProcess()
+{
+    deInitialise();
+}
+
+void DFProcess::initialise( DFProcConfig Config )
+{
+    m_length = Config.length;
+    m_winPre = Config.winPre;
+    m_winPost = Config.winPost;
+    m_alphaNormParam = Config.AlphaNormParam;
+
+    m_isMedianPositive = Config.isMedianPositive;
+
+    filtSrc = new double[ m_length ];
+    filtDst = new double[ m_length ];
+
+	
+    //Low Pass Smoothing Filter Config
+    m_FilterConfigParams.ord = Config.LPOrd;
+    m_FilterConfigParams.ACoeffs = Config.LPACoeffs;
+    m_FilterConfigParams.BCoeffs = Config.LPBCoeffs;
+	
+    m_FiltFilt = new FiltFilt( m_FilterConfigParams );	
+}
+
+void DFProcess::deInitialise()
+{
+    delete [] filtSrc;
+
+    delete [] filtDst;
+
+    delete [] m_filtScratchIn;
+
+    delete [] m_filtScratchOut;
+
+    delete m_FiltFilt;
+}
+
+void DFProcess::process(double *src, double* dst)
+{
+    removeDCNormalize( src, filtSrc );
+
+    m_FiltFilt->process( filtSrc, filtDst, m_length );
+
+    medianFilter( filtDst, dst );
+}
+
+
+void DFProcess::medianFilter(double *src, double *dst)
+{
+    unsigned int i,k,j,l;
+    unsigned int index = 0;
+
+    double val = 0;
+
+    double* y = new double[ m_winPost + m_winPre + 1];
+    memset( y, 0, sizeof( double ) * ( m_winPost + m_winPre + 1) );
+
+    double* scratch = new double[ m_length ];
+
+    for( i = 0; i < m_winPre; i++)
+    {
+	k = i + m_winPost + 1;
+
+	for( j = 0; j < k; j++)
+	{
+	    y[ j ] = src[ j ];
+	}
+	scratch[ index ] = MathUtilities::median( y, k );
+	index++;
+    }
+
+    for(  i = 0; i < ( m_length - ( m_winPost + m_winPre ) ); i ++)
+    {
+			 
+	l = 0;
+	for(  j  = i; j < ( i + m_winPost + m_winPre + 1); j++)
+	{
+	    y[ l ] = src[ j ];
+	    l++;
+	}
+
+	scratch[ index++ ] = MathUtilities::median( y, (m_winPost + m_winPre + 1 ));
+    }
+
+    for( i = std::max( m_length - m_winPost, (unsigned)1); i < m_length; i++)
+    {
+	k = std::max( i - m_winPre, (unsigned)1);
+
+	l = 0;
+	for( j = k; j < m_length; j++)
+	{
+	    y[ l ] = src[ j ];
+
+	    l++;
+	}
+		
+	scratch[ index++ ] = MathUtilities::median( y, l); 
+    }
+
+
+    for( i = 0; i < m_length; i++ )
+    {
+	val = src[ i ] - scratch[ i ];// - 0.033;
+		
+	if( m_isMedianPositive )
+	{
+	    if( val > 0 )
+	    {
+		dst[ i ]  = val;
+	    }
+	    else
+	    {
+		dst[ i ]  = 0;
+	    }
+	}
+	else
+	{
+	    dst[ i ]  = val;
+	}
+    }
+	
+    delete [] y;
+    delete [] scratch;
+}
+
+
+void DFProcess::removeDCNormalize( double *src, double*dst )
+{
+    double DFmax = 0;
+    double DFMin = 0;
+    double DFAlphaNorm = 0;
+
+    MathUtilities::getFrameMinMax( src, m_length, &DFMin, &DFmax );
+
+    MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm );
+
+    for( unsigned int i = 0; i< m_length; i++)
+    {
+	dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; 
+    }
+}
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/DFProcess.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/DFProcess.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,63 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef CDFPROCESS_H
+#define CDFPROCESS_H
+
+#include <stdio.h>
+#include "FiltFilt.h"
+
+struct DFProcConfig{
+    unsigned int length; 
+    unsigned int LPOrd; 
+    double *LPACoeffs; 
+    double *LPBCoeffs; 
+    unsigned int winPre;
+    unsigned int winPost; 
+    double AlphaNormParam;
+    bool isMedianPositive;};
+
+class DFProcess  
+{
+public:
+    DFProcess( DFProcConfig Config );
+    virtual ~DFProcess();
+
+    void process( double* src, double* dst );
+
+	
+private:
+    void initialise( DFProcConfig Config );
+    void deInitialise();
+    void removeDCNormalize( double *src, double*dst );
+    void medianFilter( double* src, double* dst );
+
+    unsigned int m_length;
+    unsigned int m_FFOrd;
+
+    unsigned int m_winPre;
+    unsigned int m_winPost;
+
+    double m_alphaNormParam;
+
+    double* filtSrc;
+    double* filtDst;
+
+    double* m_filtScratchIn;
+    double* m_filtScratchOut;
+
+    FiltFiltConfig m_FilterConfigParams;
+
+    FiltFilt* m_FiltFilt;
+
+    bool m_isMedianPositive;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/FiltFilt.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/FiltFilt.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,124 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "FiltFilt.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+FiltFilt::FiltFilt( FiltFiltConfig Config )
+{
+    m_filtScratchIn = NULL;
+    m_filtScratchOut = NULL;
+    m_ord = 0;
+	
+    initialise( Config );
+}
+
+FiltFilt::~FiltFilt()
+{
+    deInitialise();
+}
+
+void FiltFilt::initialise( FiltFiltConfig Config )
+{
+    m_ord = Config.ord;
+    m_filterConfig.ord = Config.ord;
+    m_filterConfig.ACoeffs = Config.ACoeffs;
+    m_filterConfig.BCoeffs = Config.BCoeffs;
+	
+    m_filter = new Filter( m_filterConfig );
+}
+
+void FiltFilt::deInitialise()
+{
+    delete m_filter;
+}
+
+
+void FiltFilt::process(double *src, double *dst, unsigned int length)
+{	
+    unsigned int i;
+
+    unsigned int nFilt = m_ord + 1;
+    unsigned int nFact = 3 * ( nFilt - 1);
+    unsigned int nExt	= length + 2 * nFact;
+
+
+    m_filtScratchIn = new double[ nExt ];
+    m_filtScratchOut = new double[ nExt ];
+
+	
+    for( i = 0; i< nExt; i++ ) 
+    {
+	m_filtScratchIn[ i ] = 0.0;
+	m_filtScratchOut[ i ] = 0.0;
+    }
+
+    // Edge transients reflection
+    double sample0 = 2 * src[ 0 ];
+    double sampleN = 2 * src[ length - 1 ];
+
+    unsigned int index = 0;
+    for( i = nFact; i > 0; i-- )
+    {
+	m_filtScratchIn[ index++ ] = sample0 - src[ i ];
+    }
+    index = 0;
+    for( i = 0; i < nFact; i++ )
+    {
+	m_filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ];
+    }
+
+    index = 0;
+    for( i = 0; i < length; i++ )
+    {
+	m_filtScratchIn[ i + nFact ] = src[ i ];
+    }
+	
+    ////////////////////////////////
+    // Do  0Ph filtering
+    m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt);
+	
+    // reverse the series for FILTFILT 
+    for ( i = 0; i < nExt; i++)
+    { 
+	m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1];
+    }
+
+    // do FILTER again 
+    m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt);
+	
+    // reverse the series back 
+    for ( i = 0; i < nExt; i++)
+    {
+	m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1 ];
+    }
+    for ( i = 0;i < nExt; i++)
+    {
+	m_filtScratchOut[ i ] = m_filtScratchIn[ i ];
+    }
+
+    index = 0;
+    for( i = 0; i < length; i++ )
+    {
+	dst[ index++ ] = m_filtScratchOut[ i + nFact ];
+    }	
+
+    delete [] m_filtScratchIn;
+    delete [] m_filtScratchOut;
+
+}
+
+void FiltFilt::reset()
+{
+
+}
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/FiltFilt.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/FiltFilt.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,45 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef FILTFILT_H
+#define FILTFILT_H
+
+#include "Filter.h"
+
+struct FiltFiltConfig{
+    unsigned int ord;
+    double* ACoeffs;
+    double* BCoeffs;
+};
+
+class FiltFilt  
+{
+public:
+    FiltFilt( FiltFiltConfig Config );
+    virtual ~FiltFilt();
+
+    void reset();
+    void process( double* src, double* dst, unsigned int length );
+
+private:
+    void initialise( FiltFiltConfig Config );
+    void deInitialise();
+
+    unsigned int m_ord;
+
+    Filter* m_filter;
+
+    double* m_filtScratchIn;
+    double* m_filtScratchOut;
+
+    FilterConfig m_filterConfig;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/Filter.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/Filter.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,82 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "Filter.h"
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Filter::Filter( FilterConfig Config )
+{
+    m_ord = 0;
+    m_outBuffer = NULL;
+    m_inBuffer = NULL;
+
+    initialise( Config );
+}
+
+Filter::~Filter()
+{
+    deInitialise();
+}
+
+void Filter::initialise( FilterConfig Config )
+{
+    m_ord = Config.ord;
+    m_ACoeffs = Config.ACoeffs;
+    m_BCoeffs = Config.BCoeffs;
+
+    m_inBuffer = new double[ m_ord + 1 ];
+    m_outBuffer = new double[ m_ord + 1 ];
+
+    reset();
+}
+
+void Filter::deInitialise()
+{
+    delete[] m_inBuffer;
+    delete[] m_outBuffer;
+}
+
+void Filter::reset()
+{
+    for( unsigned int i = 0; i < m_ord+1; i++ ){ m_inBuffer[ i ] = 0.0; }
+    for(unsigned int  i = 0; i < m_ord+1; i++ ){ m_outBuffer[ i ] = 0.0; }
+}
+
+void Filter::process( double *src, double *dst, unsigned int length )
+{
+    unsigned int SP,i,j;
+
+    double xin,xout;
+
+    for (SP=0;SP<length;SP++)
+    {
+        xin=src[SP];
+        /* move buffer */
+        for ( i = 0; i < m_ord; i++) {m_inBuffer[ m_ord - i ]=m_inBuffer[ m_ord - i - 1 ];}
+        m_inBuffer[0]=xin;
+
+        xout=0.0;
+        for (j=0;j< m_ord + 1; j++)
+	    xout = xout + m_BCoeffs[ j ] * m_inBuffer[ j ];
+        for (j = 0; j < m_ord; j++)
+	    xout= xout - m_ACoeffs[ j + 1 ] * m_outBuffer[ j ];
+
+        dst[ SP ] = xout;
+        for ( i = 0; i < m_ord - 1; i++ ) { m_outBuffer[ m_ord - i - 1 ] = m_outBuffer[ m_ord - i - 2 ];}
+        m_outBuffer[0]=xout;
+
+    } /* end of SP loop */
+}
+
+
+
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/Filter.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/Filter.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,48 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef FILTER_H
+#define FILTER_H
+
+#ifndef NULL
+#define NULL 0
+#endif
+
+struct FilterConfig{
+    unsigned int ord;
+    double* ACoeffs;
+    double* BCoeffs;
+};
+
+class Filter  
+{
+public:
+    Filter( FilterConfig Config );
+    virtual ~Filter();
+
+    void reset();
+
+    void process( double *src, double *dst, unsigned int length );
+	
+
+private:
+    void initialise( FilterConfig Config );
+    void deInitialise();
+
+    unsigned int m_ord;
+
+    double* m_inBuffer;
+    double* m_outBuffer;
+
+    double* m_ACoeffs;
+    double* m_BCoeffs;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/Framer.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/Framer.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,104 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "Framer.h"
+#include <math.h>
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Framer::Framer()
+{
+    m_dataFrame = NULL;
+    m_strideFrame = NULL;
+}
+
+Framer::~Framer()
+{
+    if( m_dataFrame != NULL )
+	delete [] m_dataFrame;
+
+    if( m_strideFrame != NULL )
+	delete [] m_strideFrame;
+}
+
+void Framer::configure( unsigned int frameLength, unsigned int hop )
+{
+    m_frameLength = frameLength;
+    m_stepSize = hop;
+
+    resetCounters();
+
+    if( m_dataFrame != NULL )
+    {
+	delete [] m_dataFrame;	
+	m_dataFrame = NULL;
+    }
+    m_dataFrame = new double[ m_frameLength ];
+
+    if( m_strideFrame != NULL )
+    {
+	delete [] m_strideFrame;	
+	m_strideFrame = NULL;
+    }
+    m_strideFrame = new double[ m_stepSize ];
+}
+
+void Framer::getFrame(double *dst)
+{
+
+    if( (m_ulSrcIndex + ( m_frameLength) ) < m_ulSampleLen )
+    {
+	for( unsigned int u = 0; u < m_frameLength; u++)
+	{
+	    dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; 
+	}	
+	m_ulSrcIndex -= ( m_frameLength - m_stepSize );
+    }
+    else
+    {
+	unsigned int rem = (m_ulSampleLen - m_ulSrcIndex );
+	unsigned int zero = m_frameLength - rem;
+
+	for( unsigned int u = 0; u < rem; u++ )
+	{
+	    dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ];
+	}
+		
+	for( unsigned int u = 0; u < zero; u++ )
+	{
+	    dst[ rem + u ] = 0;
+	}
+
+	m_ulSrcIndex -= (( rem - m_stepSize ) );
+    }
+
+    m_framesRead++;
+}
+
+void Framer::resetCounters()
+{
+    m_framesRead = 0;
+    m_ulSrcIndex = 0;
+}
+
+unsigned int Framer::getMaxNoFrames()
+{
+    return m_maxFrames;
+}
+
+void Framer::setSource(double *src, unsigned int length)
+{
+    m_srcBuffer = src;
+    m_ulSampleLen = length;
+
+    m_maxFrames = (unsigned int)ceil( (double)m_ulSampleLen/(double)m_stepSize ) ;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/signalconditioning/Framer.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/signalconditioning/Framer.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,47 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef FRAMER_H
+#define FRAMER_H
+
+//#include <io.h>
+#include <fcntl.h>
+#include <stdio.h>
+
+
+class Framer  
+{
+public:
+    void setSource( double* src, unsigned int length );
+    unsigned int getMaxNoFrames();
+    void getFrame( double* dst );
+    void configure( unsigned int frameLength, unsigned int hop );
+    Framer();
+    virtual ~Framer();
+
+    void resetCounters();
+
+private:
+
+    unsigned long	m_ulSampleLen;		// DataLength (samples)
+    unsigned int	m_framesRead;		// Read Frames Index
+
+    double*			m_srcBuffer;
+    double*			m_dataFrame;		// Analysis Frame Buffer
+    double*			m_strideFrame;		// Stride Frame Buffer
+    unsigned int	m_frameLength;		// Analysis Frame Length
+    unsigned int	m_stepSize;		// Analysis Frame Stride
+
+    unsigned int	m_maxFrames;
+
+    unsigned long	m_ulSrcIndex;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/tempotracking/TempoTrack.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tempotracking/TempoTrack.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,774 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#include "TempoTrack.h"
+
+#include "dsp/maths/MathAliases.h"
+#include "dsp/maths/MathUtilities.h"
+
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+TempoTrack::TempoTrack( TTParams Params )
+{
+    m_tempoScratch = NULL;
+    m_rawDFFrame = NULL;
+    m_smoothDFFrame = NULL;
+    m_frameACF = NULL;
+
+    m_dataLength = 0;
+    m_winLength = 0;
+    m_lagLength = 0;
+
+    m_rayparam = 0;
+    m_sigma = 0;
+    m_DFWVNnorm = 0;
+
+    initialise( Params );
+}
+
+TempoTrack::~TempoTrack()
+{
+    deInitialise();
+}
+
+void TempoTrack::initialise( TTParams Params )
+{	
+    m_winLength = Params.winLength;
+    m_lagLength = Params.lagLength;
+
+    m_rayparam	 = 43.0;
+    m_sigma = sqrt(3.9017);
+    m_DFWVNnorm = exp( ( log( 2.0 ) / m_rayparam ) * ( m_winLength + 2 ) );
+
+    m_rawDFFrame = new double[ m_winLength ];
+    m_smoothDFFrame = new double[ m_winLength ];
+    m_frameACF = new double[ m_winLength ];
+    m_tempoScratch = new double[ m_lagLength ];
+
+    unsigned int winPre = Params.WinT.pre;
+    unsigned int winPost = Params.WinT.post;
+
+    m_DFFramer.configure( m_winLength, m_lagLength );
+	
+    m_DFPParams.length = m_winLength;
+    m_DFPParams.AlphaNormParam = Params.alpha;
+    m_DFPParams.LPOrd = Params.LPOrd;
+    m_DFPParams.LPACoeffs = Params.LPACoeffs;
+    m_DFPParams.LPBCoeffs = Params.LPBCoeffs;
+    m_DFPParams.winPre = Params.WinT.pre;
+    m_DFPParams.winPost = Params.WinT.post;
+    m_DFPParams.isMedianPositive = true;
+	
+    m_DFConditioning = new DFProcess( m_DFPParams );
+
+}
+
+void TempoTrack::deInitialise()
+{	
+    delete [] m_rawDFFrame;
+	
+    delete [] m_smoothDFFrame;
+	
+    delete [] m_frameACF;
+
+    delete [] m_tempoScratch;
+
+    delete m_DFConditioning;
+}
+
+void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag)
+{
+    unsigned int i;
+
+    if( beatLag == 0 )
+    {
+	for( i = 0; i < winLength; i++ )
+	{    
+	    Filter[ i ] = ( ( i + 1 ) / pow( m_rayparam, 2.0) ) * exp( ( -pow(( i + 1 ),2.0 ) / ( 2.0 * pow( m_rayparam, 2.0))));
+	}
+    }
+    else
+    {	
+	m_sigma = beatLag/8;
+	for( i = 0; i < winLength; i++ )
+	{
+	    double dlag = (double)(i+1) - beatLag;
+	    Filter[ i ] =  exp(-0.5 * pow(( dlag / m_sigma), 2.0) ) / (sqrt( 2 * PI) * m_sigma);
+	}
+    }
+}
+
+double TempoTrack::tempoMM(double* ACF, double* weight, int tsig)
+{
+
+    double period = 0;
+    double maxValRCF = 0.0;
+    unsigned int maxIndexRCF = 0;
+
+    double* pdPeaks;
+
+    unsigned int maxIndexTemp;
+    double	maxValTemp;
+    unsigned int count; 
+	
+    unsigned int numelem;
+    int i, a, b;
+
+    for( i = 0; i < m_lagLength; i++ )
+	m_tempoScratch[ i ] = 0.0;
+
+    if( tsig == 0 ) 
+    {
+	//if time sig is unknown, use metrically unbiased version of Filterbank
+	numelem = 4;
+    }
+    else
+    {
+	numelem = tsig;
+    }
+
+    for(i=1;i<m_lagLength-1;i++)
+    {
+	//first and last output values are left intentionally as zero
+	for (a=1;a<=numelem;a++)
+	{
+	    for(b=(1-a);b<a;b++)
+	    {
+		if( tsig == 0 )
+		{					
+		    m_tempoScratch[i] += ACF[a*(i+1)+b-1] * (1.0 / (2.0 * (double)a-1)) * weight[i];
+		}
+		else
+		{
+		    m_tempoScratch[i] += ACF[a*(i+1)+b-1] * 1 * weight[i];
+		}
+	    }
+	}
+    }
+
+
+    //NOW FIND MAX INDEX OF ACFOUT
+    for( i = 0; i < m_lagLength; i++)
+    {
+	if( m_tempoScratch[ i ] > maxValRCF)
+	{
+	    maxValRCF = m_tempoScratch[ i ];
+	    maxIndexRCF = i;
+	}
+    }
+	
+    if( tsig == 0 )
+	tsig = 4;
+
+	
+    if( tsig == 4 )
+    {
+	pdPeaks = new double[ 4 ];
+	for( i = 0; i < 4; i++ ){ pdPeaks[ i ] = 0.0;}
+
+	pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;
+
+	maxIndexTemp = 0;
+	maxValTemp = 0.0;
+	count = 0;
+
+	for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )
+	{
+	    if( ACF[ i ] > maxValTemp )
+	    {
+		maxValTemp = ACF[ i ];
+		maxIndexTemp = count;
+	    }
+	    count++;
+	}
+	pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;
+
+	maxIndexTemp = 0;
+	maxValTemp = 0.0;
+	count = 0;
+
+	for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )
+	{
+	    if( ACF[ i ] > maxValTemp )
+	    {
+		maxValTemp = ACF[ i ];
+		maxIndexTemp = count;
+	    }
+	    count++;
+	}
+	pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;
+
+	maxIndexTemp = 0;
+	maxValTemp = 0.0;
+	count = 0;
+
+	for( i = ( 4 * maxIndexRCF + 3) - 3; i < ( 4 * maxIndexRCF + 3) + 4; i++ )
+	{
+	    if( ACF[ i ] > maxValTemp )
+	    {
+		maxValTemp = ACF[ i ];
+		maxIndexTemp = count;
+	    }
+	    count++;
+	}
+	pdPeaks[ 3 ] = (double)( maxIndexTemp + 1 + ( (4 * maxIndexRCF + 3) - 9 ) + 1 )/4 ;
+
+
+	period = MathUtilities::mean( pdPeaks, 4 );
+    }
+    else
+    {
+	pdPeaks = new double[ 3 ];
+	for( i = 0; i < 3; i++ ){ pdPeaks[ i ] = 0.0;}
+
+	pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;
+
+	maxIndexTemp = 0;
+	maxValTemp = 0.0;
+	count = 0;
+
+	for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )
+	{
+	    if( ACF[ i ] > maxValTemp )
+	    {
+		maxValTemp = ACF[ i ];
+		maxIndexTemp = count;
+	    }
+	    count++;
+	}
+	pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;
+
+	maxIndexTemp = 0;
+	maxValTemp = 0.0;
+	count = 0;
+
+	for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )
+	{
+	    if( ACF[ i ] > maxValTemp )
+	    {
+		maxValTemp = ACF[ i ];
+		maxIndexTemp = count;
+	    }
+	    count++;
+	}
+	pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;
+
+
+	period = MathUtilities::mean( pdPeaks, 3 );
+    }
+
+    delete [] pdPeaks;
+
+    return period;
+}
+
+void TempoTrack::stepDetect( double* periodP, double* periodG, int currentIdx, int* flag )
+{
+    double stepthresh = 1 * 3.9017;
+
+    if( *flag )
+    {
+	if(abs(periodG[ currentIdx ] - periodP[ currentIdx ]) > stepthresh)
+	{
+	    // do nuffin'
+	}
+    }
+    else
+    {
+	if(fabs(periodG[ currentIdx ]-periodP[ currentIdx ]) > stepthresh)
+	{
+	    *flag = 3;
+	}
+    }
+}
+
+void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag )
+{
+    double constthresh = 2 * 3.9017;
+
+    if( fabs( 2 * periodP[ currentIdx ] - periodP[ currentIdx - 1] - periodP[ currentIdx - 2] ) < constthresh)
+    {
+	*flag = 1;
+    }
+    else
+    {
+	*flag = 0;
+    }
+}
+
+int TempoTrack::findMeter(double *ACF, unsigned int len, double period)
+{
+    int i;
+    int p = (int)MathUtilities::round( period );
+    int tsig;
+
+    double Energy_3 = 0.0;
+    double Energy_4 = 0.0;
+
+    double temp3A = 0.0;
+    double temp3B = 0.0;
+    double temp4A = 0.0;
+    double temp4B = 0.0;
+
+    double* dbf = new double[ len ]; int t = 0;
+    for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; }
+
+    if( (double)len < 6 * p + 2 )
+    {
+	for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )
+	{
+	    temp3A += ACF[ i ];
+	    dbf[ t++ ] = ACF[ i ];
+	}
+	
+	for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )
+	{
+	    temp4A += ACF[ i ];
+	}
+
+	Energy_3 = temp3A;
+	Energy_4 = temp4A;
+    }
+    else
+    {
+	for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )
+	{
+	    temp3A += ACF[ i ];
+	}
+	
+	for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )
+	{
+	    temp4A += ACF[ i ];
+	}
+
+	for( i = ( 6 * p - 2 ); i < ( 6 * p + 2 ) + 1; i++ )
+	{
+	    temp3B += ACF[ i ];
+	}
+	
+	for( i = ( 2 * p - 2 ); i < ( 2 * p + 2 ) + 1; i++ )
+	{
+	    temp4B += ACF[ i ];
+	}
+
+	Energy_3 = temp3A + temp3B;
+	Energy_4 = temp4A + temp4B;
+    }
+
+    if (Energy_3 > Energy_4)
+    {
+	tsig = 3;
+    }
+    else
+    {
+	tsig = 4;
+    }
+
+
+    return tsig;
+}
+
+void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat)
+{	
+    int p = (int)MathUtilities::round( period );
+    int predictedOffset = 0;
+
+    double* phaseScratch = new double[ p*2 ];
+
+	
+    if( lastBeat != 0 )
+    {
+	lastBeat = (int)MathUtilities::round((double)lastBeat );///(double)winLength);
+
+	    predictedOffset = lastBeat + p - fsp;
+
+	    if (predictedOffset < 0) 
+	    {
+		lastBeat = 0;
+	    }
+    }
+
+    if( lastBeat != 0 )
+    {
+	int mu = p;
+	double sigma = (double)p/4;
+	double PhaseMin = 0.0;
+	double PhaseMax = 0.0;
+	unsigned int scratchLength = p*2;
+	double temp = 0.0;
+
+	for(  int i = 0; i < scratchLength; i++ )
+	{
+	    phaseScratch[ i ] = exp( -0.5 * pow( ( i - mu ) / sigma, 2 ) ) / ( sqrt( 2*PI ) *sigma );
+	}
+
+	MathUtilities::getFrameMinMax( phaseScratch, scratchLength, &PhaseMin, &PhaseMax );
+			
+	for(int i = 0; i < scratchLength; i ++)
+	{
+	    temp = phaseScratch[ i ];
+	    phaseScratch[ i ] = (temp - PhaseMin)/PhaseMax;
+	}
+
+	unsigned int index = 0;
+	for(int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++)
+	{
+	    Filter[ index++ ] = phaseScratch[ i ];
+	}
+    }
+    else
+    {
+	for( int i = 0; i < p; i ++)
+	{
+	    Filter[ i ] = 1;
+	}
+    }
+	
+    delete [] phaseScratch;
+}
+
+int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period)
+{
+    int alignment = 0;
+    int p = (int)MathUtilities::round( period );
+
+    double temp = 0.0;
+
+    double* y = new double[ winLength ];
+    double* align = new double[ p ];
+
+    for( int i = 0; i < winLength; i++ )
+    {	
+	y[ i ] = (double)( -i + winLength  )/(double)winLength;
+    }
+
+    for( int o = 0; o < p; o++ )
+    { 
+	temp = 0.0;
+	for(int i = 1 + (o - 1); i< winLength; i += (p + 1))
+	{
+	    temp = temp + DF[ i ] * y[ i ]; 
+	}
+	align[ o ] = temp * weighting[ o ];       
+    }
+
+
+    double valTemp = 0.0;
+    for(int i = 0; i < p; i++)
+    {
+	if( align[ i ] > valTemp )
+	{
+	    valTemp = align[ i ];
+	    alignment = i;
+	}
+    }
+
+    delete [] y;
+    delete [] align;
+
+    return alignment;
+}
+
+int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step )
+{
+    int beat = 0;
+
+    int p = (int)MathUtilities::round( period );
+    int align = (int)MathUtilities::round( alignment );
+    int FSP = (int)MathUtilities::round( FSP0 );
+
+    int FEP = FSP + ( step );
+
+    beat = FSP + align;
+
+    m_beats.push_back( beat );
+
+    while( beat + p < FEP )
+    {
+	beat += p;
+		
+	m_beats.push_back( beat );
+    }
+
+    return beat;
+}
+
+vector<int> TempoTrack::process(double *DF, unsigned int length)
+{
+    m_dataLength = length;
+	
+    double	period = 0.0;
+    int stepFlag = 0;
+    int constFlag = 0;
+    int FSP = 0;
+    int tsig = 0;
+    int lastBeat = 0;
+
+	
+    double* RW = new double[ m_lagLength ];
+    for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}
+
+    double* GW = new double[ m_lagLength ];
+    for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}
+
+    double* PW = new double[ m_lagLength ];
+    for(unsigned int clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}
+
+    m_DFFramer.setSource( DF, m_dataLength );
+
+    unsigned int TTFrames = m_DFFramer.getMaxNoFrames();
+	
+    double* periodP = new double[ TTFrames ];
+    for(unsigned int  clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}
+	
+    double* periodG = new double[ TTFrames ];
+    for(unsigned int  clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}
+	
+    double* alignment = new double[ TTFrames ];
+    for(unsigned int  clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}
+
+    m_beats.clear();
+
+    createCombFilter( RW, m_lagLength, 0, 0 );
+
+    int TTLoopIndex = 0;
+
+    for( unsigned int i = 0; i < TTFrames; i++ )
+    {
+	m_DFFramer.getFrame( m_rawDFFrame );
+
+	m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame );
+
+	m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength );
+		
+	periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 );
+
+	if( GW[ 0 ] != 0 )
+	{
+	    periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig );
+	}
+	else
+	{
+	    periodG[ TTLoopIndex ] = 0.0;
+	}
+
+	stepDetect( periodP, periodG, TTLoopIndex, &stepFlag );
+
+	if( stepFlag == 1)
+	{
+	    constDetect( periodP, TTLoopIndex, &constFlag );
+	    stepFlag = 0;
+	}
+	else
+	{
+	    stepFlag -= 1;
+	}
+
+	if( stepFlag < 0 )
+	{
+	    stepFlag = 0;
+	}
+
+	if( constFlag != 0)
+	{
+	    tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] );
+	
+	    createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] );
+			
+	    periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); 
+
+	    period = periodG[ TTLoopIndex ];
+
+	    createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); 
+
+	    constFlag = 0;
+
+	}
+	else
+	{
+	    if( GW[ 0 ] != 0 )
+	    {
+		period = periodG[ TTLoopIndex ];
+		createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); 
+
+	    }
+	    else
+	    {
+		period = periodP[ TTLoopIndex ];
+		createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); 
+	    }
+	}
+
+	alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); 
+
+	lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength );
+
+	FSP += (m_lagLength);
+
+	TTLoopIndex++;
+    }
+
+
+    delete [] periodP;
+    delete [] periodG;
+    delete [] alignment;
+
+    delete [] RW;
+    delete [] GW;
+    delete [] PW;
+
+    return m_beats;
+}
+
+
+
+
+
+vector<int> TempoTrack::process( vector <double> DF )
+{
+    m_dataLength = DF.size();
+	
+    double	period = 0.0;
+    int stepFlag = 0;
+    int constFlag = 0;
+    int FSP = 0;
+    int tsig = 0;
+    int lastBeat = 0;
+
+    vector <double> causalDF;
+
+    causalDF = DF;
+
+    //Prepare Causal Extension DFData
+    unsigned int DFCLength = m_dataLength + m_winLength;
+	
+    for( unsigned int j = 0; j < m_winLength; j++ )
+    {
+	causalDF.push_back( 0 );
+    }
+	
+	
+    double* RW = new double[ m_lagLength ];
+    for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}
+
+    double* GW = new double[ m_lagLength ];
+    for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}
+
+    double* PW = new double[ m_lagLength ];
+    for(unsigned clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}
+
+    m_DFFramer.setSource( &causalDF[0], m_dataLength );
+
+    unsigned int TTFrames = m_DFFramer.getMaxNoFrames();
+	
+    double* periodP = new double[ TTFrames ];
+    for(unsigned clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}
+	
+    double* periodG = new double[ TTFrames ];
+    for(unsigned clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}
+	
+    double* alignment = new double[ TTFrames ];
+    for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}
+
+    m_beats.clear();
+
+    createCombFilter( RW, m_lagLength, 0, 0 );
+
+    int TTLoopIndex = 0;
+
+    for( unsigned int i = 0; i < TTFrames; i++ )
+    {
+	m_DFFramer.getFrame( m_rawDFFrame );
+
+	m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame );
+
+	m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength );
+		
+	periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 );
+
+	if( GW[ 0 ] != 0 )
+	{
+	    periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig );
+	}
+	else
+	{
+	    periodG[ TTLoopIndex ] = 0.0;
+	}
+
+	stepDetect( periodP, periodG, TTLoopIndex, &stepFlag );
+
+	if( stepFlag == 1)
+	{
+	    constDetect( periodP, TTLoopIndex, &constFlag );
+	    stepFlag = 0;
+	}
+	else
+	{
+	    stepFlag -= 1;
+	}
+
+	if( stepFlag < 0 )
+	{
+	    stepFlag = 0;
+	}
+
+	if( constFlag != 0)
+	{
+	    tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] );
+	
+	    createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] );
+			
+	    periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); 
+
+	    period = periodG[ TTLoopIndex ];
+
+	    createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); 
+
+	    constFlag = 0;
+
+	}
+	else
+	{
+	    if( GW[ 0 ] != 0 )
+	    {
+		period = periodG[ TTLoopIndex ];
+		createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); 
+
+	    }
+	    else
+	    {
+		period = periodP[ TTLoopIndex ];
+		createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); 
+	    }
+	}
+
+	alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); 
+
+	lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength );
+
+	FSP += (m_lagLength);
+
+	TTLoopIndex++;
+    }
+
+
+    delete [] periodP;
+    delete [] periodG;
+    delete [] alignment;
+
+    delete [] RW;
+    delete [] GW;
+    delete [] PW;
+
+    return m_beats;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/tempotracking/TempoTrack.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tempotracking/TempoTrack.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,97 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2005-2006 Christian Landone.
+    All rights reserved.
+*/
+
+#ifndef TEMPOTRACK_H
+#define TEMPOTRACK_H
+
+
+#include <stdio.h>
+#include <vector>
+
+#include "dsp/signalconditioning/DFProcess.h"
+#include "dsp/maths/Correlation.h"
+#include "dsp/signalconditioning/Framer.h"
+
+
+
+using std::vector;
+
+struct WinThresh
+{
+    unsigned int pre;
+    unsigned int  post;
+};
+
+struct TTParams
+{
+    unsigned int winLength; //Analysis window length
+    unsigned int lagLength; //Lag & Stride size
+    unsigned int alpha; //alpha-norm parameter
+    unsigned int LPOrd; // low-pass Filter order
+    double* LPACoeffs; //low pass Filter den coefficients
+    double* LPBCoeffs; //low pass Filter num coefficients
+    WinThresh WinT;//window size in frames for adaptive thresholding [pre post]:
+};
+
+
+class TempoTrack  
+{
+public:
+    TempoTrack( TTParams Params );
+    virtual ~TempoTrack();
+
+    vector<int> process( vector <double> DF );
+    vector<int> process( double* DF, unsigned int length );
+
+	
+private:
+    void initialise( TTParams Params );
+    void deInitialise();
+
+    int beatPredict( unsigned int FSP, double alignment, double period, unsigned int step);
+    int phaseMM( double* DF, double* weighting, unsigned int winLength, double period );
+    void createPhaseExtractor( double* Filter, unsigned int winLength,  double period,  unsigned int fsp, unsigned int lastBeat );
+    int findMeter( double* ACF,  unsigned int len, double period );
+    void constDetect( double* periodP, int currentIdx, int* flag );
+    void stepDetect( double* periodP, double* periodG, int currentIdx, int* flag );
+    void createCombFilter( double* Filter, unsigned int winLength, unsigned int TSig, double beatLag );
+    double tempoMM( double* ACF, double* weight, int sig );
+	
+    unsigned int m_dataLength;
+    unsigned int m_winLength;
+    unsigned int m_lagLength;
+
+    double		 m_rayparam;
+    double		 m_sigma;
+    double		 m_DFWVNnorm;
+
+    vector<int>	 m_beats; // Vector of detected beats
+
+    double* m_tempoScratch;
+	
+    // Processing Buffers 
+    double* m_rawDFFrame; // Original Detection Function Analysis Frame
+    double* m_smoothDFFrame; // Smoothed Detection Function Analysis Frame
+    double* m_frameACF; // AutoCorrelation of Smoothed Detection Function 
+
+    //Low Pass Coefficients for DF Smoothing
+    double* m_ACoeffs;
+    double* m_BCoeffs;
+	
+    // Objetcs/operators declaration
+    Framer m_DFFramer;
+    DFProcess* m_DFConditioning;
+    Correlation m_correlator;
+
+    // Config structure for DFProcess
+    DFProcConfig m_DFPParams;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/ChangeDetectionFunction.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/ChangeDetectionFunction.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,139 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#include "ChangeDetectionFunction.h"
+
+#ifndef PI
+#define PI (3.14159265358979232846)
+#endif
+
+
+
+ChangeDetectionFunction::ChangeDetectionFunction(ChangeDFConfig config) :
+	m_dFilterSigma(0.0), m_iFilterWidth(0)
+{
+	setFilterWidth(config.smoothingWidth);
+}
+
+ChangeDetectionFunction::~ChangeDetectionFunction()
+{
+}
+
+void ChangeDetectionFunction::setFilterWidth(const int iWidth)
+{
+	m_iFilterWidth = iWidth*2+1;
+	
+	// it is assumed that the gaussian is 0 outside of +/- FWHM
+	// => filter width = 2*FWHM = 2*2.3548*sigma
+	m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548);
+	m_vaGaussian.resize(m_iFilterWidth);
+	
+	double dScale = 1.0 / (m_dFilterSigma*sqrt(2*PI));
+	
+	for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++)
+	{
+		double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) );
+		m_vaGaussian[x + (m_iFilterWidth-1)/2] = w;
+	}
+	
+#ifdef DEBUG_CHANGE_DETECTION_FUNCTION
+	std::cout << "Filter sigma: " << m_dFilterSigma << std::endl;
+	std::cout << "Filter width: " << m_iFilterWidth << std::endl;
+#endif
+}
+
+
+ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram)
+{
+	ChangeDistance retVal;
+	retVal.resize(rTCSGram.getSize(), 0.0);
+	
+	TCSGram smoothedTCSGram;
+
+	for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
+	{
+		int iSkipLower = 0;
+	
+		int iLowerPos = iPosition - (m_iFilterWidth-1)/2;
+		int iUpperPos = iPosition + (m_iFilterWidth-1)/2;
+	
+		if (iLowerPos < 0)
+		{
+			iSkipLower = -iLowerPos;
+			iLowerPos = 0;
+		}
+	
+		if (iUpperPos >= rTCSGram.getSize())
+		{
+			int iMaxIndex = rTCSGram.getSize() - 1;
+			iUpperPos = iMaxIndex;
+		}
+	
+		TCSVector smoothedVector;
+
+		// for every bin of the vector, calculate the smoothed value
+		for (int iPC = 0; iPC < 6; iPC++)
+		{	
+			size_t j = 0;
+			double dSmoothedValue = 0.0;
+			TCSVector rCV;
+		
+			for (int i = iLowerPos; i <= iUpperPos; i++)
+			{
+				rTCSGram.getTCSVector(i, rCV);
+				dSmoothedValue += m_vaGaussian[iSkipLower + j++] * rCV[iPC];
+			}
+
+			smoothedVector[iPC] = dSmoothedValue;
+		}
+		
+		smoothedTCSGram.addTCSVector(smoothedVector);
+	}
+
+	for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
+	{
+		/*
+			TODO: calculate a confidence measure for the current estimation
+			if the current estimate is not confident enough, look further into the future/the past
+			e.g., High frequency content, zero crossing rate, spectral flatness
+		*/
+		
+		TCSVector nextTCS;
+		TCSVector previousTCS;
+		
+		int iWindow = 1;
+
+		// while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0)
+		{
+			smoothedTCSGram.getTCSVector(iPosition-iWindow, previousTCS);
+			// std::cout << previousTCS.magnitude() << std::endl;
+			iWindow++;
+		}
+		
+		iWindow = 1;
+		
+		// while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) )
+		{
+			smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS);
+			iWindow++;
+		}
+
+		double distance = 0.0;
+		// Euclidean distance
+		for (size_t j = 0; j < 6; j++)
+		{
+			distance += std::pow(nextTCS[j] - previousTCS[j], 2.0);
+		}
+	
+		retVal[iPosition] = std::pow(distance, 0.5);
+	}
+
+	return retVal;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/ChangeDetectionFunction.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/ChangeDetectionFunction.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,43 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#ifndef _CHANGEDETECTIONFUNCTION_
+#define _CHANGEDETECTIONFUNCTION_
+
+#define DEBUG_CHANGE_DETECTION_FUNCTION 1
+
+#include "TCSgram.h"
+
+#include <valarray>
+using std::valarray;
+
+typedef	valarray<double> ChangeDistance;
+
+struct ChangeDFConfig
+{
+	int smoothingWidth;
+};
+
+class ChangeDetectionFunction
+{
+public:
+	ChangeDetectionFunction(ChangeDFConfig);
+	~ChangeDetectionFunction();
+	ChangeDistance process(const TCSGram& rTCSGram);
+private:
+	void setFilterWidth(const int iWidth);
+	
+private:
+	valarray<double> m_vaGaussian;
+	double m_dFilterSigma;
+	int m_iFilterWidth;
+};
+
+#endif // _CHANGDETECTIONFUNCTION_
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/TCSgram.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/TCSgram.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,72 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#include "TCSgram.h"
+
+#include <valarray>
+#include <cmath>
+#include <iostream>
+#include <limits>
+
+#include "dsp/maths/MathUtilities.h"
+
+TCSGram::TCSGram() :
+	m_uNumBins(6)
+{
+}
+
+TCSGram::~TCSGram()
+{
+}
+
+
+void TCSGram::getTCSVector(int iPosition, TCSVector& rTCSVector) const
+{
+	if (iPosition < 0) 
+		rTCSVector = TCSVector();
+	else if (iPosition >= m_VectorList.size())
+		rTCSVector = TCSVector();
+	else
+		rTCSVector = m_VectorList[iPosition].second;
+}
+
+long TCSGram::getTime(size_t uPosition) const
+{
+	return m_VectorList[uPosition].first;
+}
+
+
+void TCSGram::addTCSVector(const TCSVector& rTCSVector)
+{
+	size_t uSize = m_VectorList.size();
+	long lMilliSeconds = static_cast<long>(uSize*m_dFrameDurationMS);
+	std::pair<long, TCSVector> p; 
+	p.first = lMilliSeconds;
+	p.second = rTCSVector;
+	
+	m_VectorList.push_back(p);
+}
+
+long TCSGram::getDuration() const
+{
+	size_t uSize = m_VectorList.size();
+	return static_cast<long>(uSize*m_dFrameDurationMS);
+}
+
+void TCSGram::printDebug()
+{
+	vectorlist_t::iterator vectorIterator = m_VectorList.begin();
+	
+	while (vectorIterator != m_VectorList.end())
+	{
+		vectorIterator->second.printDebug();
+		vectorIterator++;
+	}
+}
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/TCSgram.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/TCSgram.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,44 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#ifndef _TCSGram_
+#define _TCSGram_
+
+#include <vector>
+#include <valarray>
+#include <utility>
+
+#include "TonalEstimator.h"
+
+typedef std::vector<std::pair<long, TCSVector> > vectorlist_t;
+
+class TCSGram
+{
+public:	
+	TCSGram();
+	~TCSGram();
+	void getTCSVector(int, TCSVector&) const;
+	void addTCSVector(const TCSVector&);
+	long getTime(size_t) const;
+	long getDuration() const;
+	void printDebug();
+	int getSize() const { return m_VectorList.size(); }
+	void reserve(size_t uSize) { m_VectorList.reserve(uSize); }
+	void clear() { m_VectorList.clear(); }
+	void setFrameDuration(const double dFrameDurationMS) { m_dFrameDurationMS = dFrameDurationMS; }
+	void setNumBins(const unsigned int uNumBins) { m_uNumBins = uNumBins; }
+	void normalize();
+protected:
+	vectorlist_t m_VectorList;
+	unsigned int m_uNumBins;
+	double m_dFrameDurationMS;
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/TonalEstimator.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/TonalEstimator.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,98 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#include "TonalEstimator.h"
+
+#include <cmath>
+#include <iostream>
+
+#ifndef PI
+#define PI (3.14159265358979232846)
+#endif
+
+TonalEstimator::TonalEstimator()
+{
+	m_Basis.resize(6);
+
+	int i = 0;
+	
+	
+	// circle of fifths
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = std::sin( (7.0 / 6.0) * iP * PI);
+	}
+	
+	i++;
+
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = std::cos( (7.0 / 6.0) * iP * PI);
+	}
+	
+	i++;
+	
+	
+	// circle of major thirds
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = 0.6 * std::sin( (2.0 / 3.0) * iP * PI);
+	}
+	
+	i++;
+
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = 0.6 * std::cos( (2.0 / 3.0) * iP * PI);
+	}
+
+	i++;
+
+
+	// circle of minor thirds
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = 1.1 * std::sin( (3.0 / 2.0) * iP * PI);
+	}
+	
+	i++;
+
+	m_Basis[i].resize(12);
+	for (int iP = 0; iP < 12; iP++)
+	{
+		m_Basis[i][iP] = 1.1 * std::cos( (3.0 / 2.0) * iP * PI);
+	}
+
+}
+
+TonalEstimator::~TonalEstimator()
+{
+}
+
+TCSVector TonalEstimator::transform2TCS(const ChromaVector& rVector)
+{
+	TCSVector vaRetVal;
+	vaRetVal.resize(6, 0.0);
+		
+	for (int i = 0; i < 6; i++)
+	{
+		for (int iP = 0; iP < 12; iP++)
+		{
+			vaRetVal[i] += m_Basis[i][iP] * rVector[iP];
+		}
+	}
+	
+	return vaRetVal;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/tonal/TonalEstimator.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/tonal/TonalEstimator.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,94 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2006 Martin Gasser.
+    All rights reserved.
+*/
+
+#ifndef _TONALESTIMATOR_
+#define _TONALESTIMATOR_
+
+
+#include <valarray>
+#include <numeric>
+#include <algorithm>
+#include <iostream>
+
+class ChromaVector : public std::valarray<double>
+{
+public:
+	ChromaVector(size_t uSize = 12) : std::valarray<double>()
+	{ resize(uSize, 0.0f); }
+	
+	virtual ~ChromaVector() {};
+	
+	void printDebug()
+	{
+		for (int i = 0; i < size(); i++)
+		{
+			std::cout <<  (*this)[i] << ";";
+		}
+		
+		std::cout << std::endl;
+	}
+	
+	void normalizeL1()
+	{
+		// normalize the chroma vector (L1 norm)
+		double dSum = 0.0;
+	
+		for (size_t i = 0; i < 12; (dSum += std::abs((*this)[i++])));
+		for (size_t i = 0; i < 12; dSum > 0.0000001?((*this)[i] /= dSum):(*this)[i]=0.0, i++);
+
+	}
+	
+};
+
+class TCSVector : public std::valarray<double>
+{
+public:
+	TCSVector() : std::valarray<double>()
+	{ resize(6, 0.0f); }
+	
+	virtual ~TCSVector() {};
+
+	void printDebug()
+	{
+		for (int i = 0; i < size(); i++)
+		{
+			std::cout <<  (*this)[i] << ";";
+		}
+		
+		std::cout << std::endl;
+	}
+	
+	double magnitude() const
+	{
+		double dMag = 0.0;
+		
+		for (size_t i = 0; i < 6; i++)
+		{
+			dMag += std::pow((*this)[i], 2.0);
+		}
+		
+		return std::sqrt(dMag);
+	}
+
+};
+
+
+
+class TonalEstimator
+{
+public:
+	TonalEstimator();
+	virtual ~TonalEstimator();
+	TCSVector transform2TCS(const ChromaVector& rVector);
+protected:
+	std::valarray< std::valarray<double> > m_Basis;
+};
+
+#endif // _TONALESTIMATOR_
diff -r 000000000000 -r 49844bc8a895 dsp/transforms/FFT.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/transforms/FFT.cpp	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,154 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file is based on Don Cross's public domain FFT implementation.
+*/
+
+#include "FFT.h"
+#include <cmath>
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+FFT::FFT()
+{
+
+}
+
+FFT::~FFT()
+{
+
+}
+
+void FFT::process(unsigned int p_nSamples, bool p_bInverseTransform, double *p_lpRealIn, double *p_lpImagIn, double *p_lpRealOut, double *p_lpImagOut)
+{
+
+    if(!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
+
+
+    unsigned int NumBits;
+    unsigned int i, j, k, n;
+    unsigned int BlockSize, BlockEnd;
+
+    double angle_numerator = 2.0 * M_PI;
+    double tr, ti;
+
+    if( !isPowerOfTwo(p_nSamples) )
+    {
+	return;
+    }
+
+    if( p_bInverseTransform ) angle_numerator = -angle_numerator;
+
+    NumBits = numberOfBitsNeeded ( p_nSamples );
+
+
+    for( i=0; i < p_nSamples; i++ )
+    {
+	j = reverseBits ( i, NumBits );
+	p_lpRealOut[j] = p_lpRealIn[i];
+	p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
+    }
+
+
+    BlockEnd = 1;
+    for( BlockSize = 2; BlockSize <= p_nSamples; BlockSize <<= 1 )
+    {
+	double delta_angle = angle_numerator / (double)BlockSize;
+	double sm2 = -sin ( -2 * delta_angle );
+	double sm1 = -sin ( -delta_angle );
+	double cm2 = cos ( -2 * delta_angle );
+	double cm1 = cos ( -delta_angle );
+	double w = 2 * cm1;
+	double ar[3], ai[3];
+
+	for( i=0; i < p_nSamples; i += BlockSize )
+	{
+
+	    ar[2] = cm2;
+	    ar[1] = cm1;
+
+	    ai[2] = sm2;
+	    ai[1] = sm1;
+
+	    for ( j=i, n=0; n < BlockEnd; j++, n++ )
+	    {
+
+		ar[0] = w*ar[1] - ar[2];
+		ar[2] = ar[1];
+		ar[1] = ar[0];
+
+		ai[0] = w*ai[1] - ai[2];
+		ai[2] = ai[1];
+		ai[1] = ai[0];
+
+		k = j + BlockEnd;
+		tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
+		ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
+
+		p_lpRealOut[k] = p_lpRealOut[j] - tr;
+		p_lpImagOut[k] = p_lpImagOut[j] - ti;
+
+		p_lpRealOut[j] += tr;
+		p_lpImagOut[j] += ti;
+
+	    }
+	}
+
+	BlockEnd = BlockSize;
+
+    }
+
+
+    if( p_bInverseTransform )
+    {
+	double denom = (double)p_nSamples;
+
+	for ( i=0; i < p_nSamples; i++ )
+	{
+	    p_lpRealOut[i] /= denom;
+	    p_lpImagOut[i] /= denom;
+	}
+    }
+}
+
+bool FFT::isPowerOfTwo(unsigned int p_nX)
+{
+    if( p_nX < 2 ) return false;
+
+    if( p_nX & (p_nX-1) ) return false;
+
+    return true;
+}
+
+unsigned int FFT::numberOfBitsNeeded(unsigned int p_nSamples)
+{	
+    int i;
+
+    if( p_nSamples < 2 )
+    {
+	return 0;
+    }
+
+    for ( i=0; ; i++ )
+    {
+	if( p_nSamples & (1 << i) ) return i;
+    }
+}
+
+unsigned int FFT::reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
+{
+    unsigned int i, rev;
+
+    for(i=rev=0; i < p_nBits; i++)
+    {
+	rev = (rev << 1) | (p_nIndex & 1);
+	p_nIndex >>= 1;
+    }
+
+    return rev;
+}
diff -r 000000000000 -r 49844bc8a895 dsp/transforms/FFT.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dsp/transforms/FFT.h	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,28 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file is based on Don Cross's public domain FFT implementation.
+*/
+
+#ifndef FFT_H
+#define FFT_H
+
+class FFT  
+{
+public:
+    static void process(unsigned int p_nSamples, bool p_bInverseTransform,
+                        double *p_lpRealIn, double *p_lpImagIn,
+                        double *p_lpRealOut, double *p_lpImagOut);
+    FFT();
+    virtual ~FFT();
+
+protected:
+    static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits);
+    static unsigned int numberOfBitsNeeded( unsigned int p_nSamples );
+    static bool isPowerOfTwo( unsigned int p_nX );
+};
+
+#endif
diff -r 000000000000 -r 49844bc8a895 qm-dsp.pro
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/qm-dsp.pro	Wed Apr 05 17:35:59 2006 +0000
@@ -0,0 +1,76 @@
+TEMPLATE = lib
+CONFIG += release warn_on staticlib
+OBJECTS_DIR = tmp_obj
+MOC_DIR = tmp_moc
+DEPENDPATH += base \
+              dsp/chromagram \
+              dsp/maths \
+              dsp/onsets \
+              dsp/phasevocoder \
+              dsp/rateconversion \
+              dsp/signalconditioning \
+              dsp/tempotracking \
+              dsp/tonal \
+              dsp/transforms
+INCLUDEPATH += . \
+               base \
+               dsp/maths \
+               dsp/chromagram \
+               dsp/transforms \
+               dsp/onsets \
+               dsp/phasevocoder \
+               dsp/signalconditioning \
+               dsp/rateconversion \
+               dsp/tempotracking \
+               dsp/tonal
+
+# Input
+HEADERS += base/Pitch.h \
+           base/Window.h \
+           dsp/chromagram/Chromagram.h \
+           dsp/chromagram/ChromaProcess.h \
+           dsp/chromagram/ConstantQ.h \
+           dsp/maths/Correlation.h \
+           dsp/maths/Histogram.h \
+           dsp/maths/MathAliases.h \
+           dsp/maths/MathUtilities.h \
+           dsp/maths/Polyfit.h \
+           dsp/maths/SimpleMatrix.h \
+           dsp/onsets/DetectionFunction.h \
+           dsp/onsets/PeakPicking.h \
+           dsp/phasevocoder/PhaseVocoder.h \
+           dsp/rateconversion/Decimator.h \
+           dsp/signalconditioning/DFProcess.h \
+           dsp/signalconditioning/Filter.h \
+           dsp/signalconditioning/FiltFilt.h \
+           dsp/signalconditioning/Framer.h \
+           dsp/signalconditioning/Windowing.h \
+           dsp/tempotracking/TempoTrack.h \
+           dsp/tonal/ChangeDetectionFunction.h \
+           dsp/tonal/ChromaMatrix.h \
+           dsp/tonal/TCSgram.h \
+           dsp/tonal/TonalEstimator.h \
+           dsp/transforms/FFT.h \
+           dsp/transforms/Fourier.h
+SOURCES += base/Pitch.cpp \
+           dsp/chromagram/Chromagram.cpp \
+           dsp/chromagram/ChromaProcess.cpp \
+           dsp/chromagram/ConstantQ.cpp \
+           dsp/maths/Correlation.cpp \
+           dsp/maths/MathUtilities.cpp \
+           dsp/onsets/DetectionFunction.cpp \
+           dsp/onsets/PeakPicking.cpp \
+           dsp/phasevocoder/PhaseVocoder.cpp \
+           dsp/rateconversion/Decimator.cpp \
+           dsp/signalconditioning/DFProcess.cpp \
+           dsp/signalconditioning/Filter.cpp \
+           dsp/signalconditioning/FiltFilt.cpp \
+           dsp/signalconditioning/Framer.cpp \
+           dsp/signalconditioning/Windowing.cpp \
+           dsp/tempotracking/TempoTrack.cpp \
+           dsp/tonal/ChangeDetectionFunction.cpp \
+           dsp/tonal/ChromaMatrix.cpp \
+           dsp/tonal/TCSgram.cpp \
+           dsp/tonal/TonalEstimator.cpp \
+           dsp/transforms/FFT.cpp \
+           dsp/transforms/Fourier.cpp