update qm-dsp library
[ardour.git] / libs / qm-dsp / dsp / tempotracking / TempoTrack.cpp
index c36385f250a0403e02eba7ce0cc548995a10c516..389403edaa8af9ac7832d6e59e31e3ed91017030 100644 (file)
-/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */\r
-\r
-/*\r
-    QM DSP Library\r
-\r
-    Centre for Digital Music, Queen Mary, University of London.\r
-    This file copyright 2005-2006 Christian Landone.and Matthew Davies.\r
+/* -*- 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.and Matthew Davies.
 
     This program is free software; you can redistribute it and/or
     modify it under the terms of the GNU General Public License as
     published by the Free Software Foundation; either version 2 of the
     License, or (at your option) any later version.  See the file
-    COPYING included with this distribution for more information.\r
-*/\r
-\r
-#include "TempoTrack.h"\r
-\r
-#include "maths/MathAliases.h"\r
-#include "maths/MathUtilities.h"\r
-\r
-#include <iostream>\r
-\r
-#include <cassert>\r
-\r
-//#define DEBUG_TEMPO_TRACK 1\r
-\r
-\r
-#define RAY43VAL\r
-\r
-//////////////////////////////////////////////////////////////////////\r
-// Construction/Destruction\r
-//////////////////////////////////////////////////////////////////////\r
-\r
-TempoTrack::TempoTrack( TTParams Params )\r
-{\r
-    m_tempoScratch = NULL;\r
-    m_rawDFFrame = NULL;\r
-    m_smoothDFFrame = NULL;\r
-    m_frameACF = NULL;\r
-       m_smoothRCF = NULL;\r
-\r
-    m_dataLength = 0;\r
-    m_winLength = 0;\r
-    m_lagLength = 0;\r
-\r
-    m_rayparam = 0;\r
-    m_sigma = 0;\r
-    m_DFWVNnorm = 0;\r
-\r
-    initialise( Params );\r
-}\r
-\r
-TempoTrack::~TempoTrack()\r
-{\r
-    deInitialise();\r
-}\r
-\r
-void TempoTrack::initialise( TTParams Params )\r
-{      \r
-    m_winLength = Params.winLength;\r
-    m_lagLength = Params.lagLength;\r
-\r
-    m_rayparam  = 43.0;\r
-    m_sigma = sqrt(3.9017);\r
-    m_DFWVNnorm = exp( ( log( 2.0 ) / m_rayparam ) * ( m_winLength + 2 ) );\r
-\r
-    m_rawDFFrame = new double[ m_winLength ];\r
-    m_smoothDFFrame = new double[ m_winLength ];\r
-    m_frameACF = new double[ m_winLength ];\r
-    m_tempoScratch = new double[ m_lagLength ];\r
-       m_smoothRCF = new double[ m_lagLength ];\r
-\r
-\r
-    unsigned int winPre = Params.WinT.pre;\r
-    unsigned int winPost = Params.WinT.post;\r
-\r
-    m_DFFramer.configure( m_winLength, m_lagLength );\r
-       \r
-    m_DFPParams.length = m_winLength;\r
-    m_DFPParams.AlphaNormParam = Params.alpha;\r
-    m_DFPParams.LPOrd = Params.LPOrd;\r
-    m_DFPParams.LPACoeffs = Params.LPACoeffs;\r
-    m_DFPParams.LPBCoeffs = Params.LPBCoeffs;\r
-    m_DFPParams.winPre = Params.WinT.pre;\r
-    m_DFPParams.winPost = Params.WinT.post;\r
-    m_DFPParams.isMedianPositive = true;\r
-       \r
-    m_DFConditioning = new DFProcess( m_DFPParams );\r
-\r
-\r
-       // these are parameters for smoothing m_tempoScratch\r
-    m_RCFPParams.length = m_lagLength;\r
-    m_RCFPParams.AlphaNormParam = Params.alpha;\r
-    m_RCFPParams.LPOrd = Params.LPOrd;\r
-    m_RCFPParams.LPACoeffs = Params.LPACoeffs;\r
-    m_RCFPParams.LPBCoeffs = Params.LPBCoeffs;\r
-    m_RCFPParams.winPre = Params.WinT.pre;\r
-    m_RCFPParams.winPost = Params.WinT.post;\r
-    m_RCFPParams.isMedianPositive = true;\r
-\r
-    m_RCFConditioning = new DFProcess( m_RCFPParams );\r
-\r
-}\r
-\r
-void TempoTrack::deInitialise()\r
-{      \r
-    delete [] m_rawDFFrame;\r
-       \r
-    delete [] m_smoothDFFrame;\r
-\r
-       delete [] m_smoothRCF;  \r
-       \r
-    delete [] m_frameACF;\r
-\r
-    delete [] m_tempoScratch;\r
-\r
-    delete m_DFConditioning;\r
-       \r
-       delete m_RCFConditioning;\r
-\r
-}\r
-\r
-void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag)\r
-{\r
-    unsigned int i;\r
-\r
-    if( beatLag == 0 )\r
-    {\r
-       for( i = 0; i < winLength; i++ )\r
-       {    \r
-           Filter[ i ] = ( ( i + 1 ) / pow( m_rayparam, 2.0) ) * exp( ( -pow(( i + 1 ),2.0 ) / ( 2.0 * pow( m_rayparam, 2.0))));\r
-       }\r
-    }\r
-    else\r
-    {  \r
-       m_sigma = beatLag/4;\r
-       for( i = 0; i < winLength; i++ )\r
-       {\r
-           double dlag = (double)(i+1) - beatLag;\r
-           Filter[ i ] =  exp(-0.5 * pow(( dlag / m_sigma), 2.0) ) / (sqrt( 2 * PI) * m_sigma);\r
-       }\r
-    }\r
-}\r
-\r
-double TempoTrack::tempoMM(double* ACF, double* weight, int tsig)\r
-{\r
-\r
-    double period = 0;\r
-    double maxValRCF = 0.0;\r
-    unsigned int maxIndexRCF = 0;\r
-\r
-    double* pdPeaks;\r
-\r
-    unsigned int maxIndexTemp;\r
-    double     maxValTemp;\r
-    unsigned int count; \r
-       \r
-    unsigned int numelem,i,j;\r
-    int a, b;\r
-\r
-    for( i = 0; i < m_lagLength; i++ )\r
-       m_tempoScratch[ i ] = 0.0;\r
-\r
-    if( tsig == 0 ) \r
-    {\r
-       //if time sig is unknown, use metrically unbiased version of Filterbank\r
-       numelem = 4;\r
-    }\r
-    else\r
-    {\r
-       numelem = tsig;\r
-    }\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-    std::cerr << "tempoMM: m_winLength = " << m_winLength << ", m_lagLength = " << m_lagLength << ", numelem = " << numelem << std::endl;\r
-#endif\r
-\r
-    for(i=1;i<m_lagLength-1;i++)\r
-    {\r
-       //first and last output values are left intentionally as zero\r
-       for (a=1;a<=numelem;a++)\r
-       {\r
-           for(b=(1-a);b<a;b++)\r
-           {\r
-               if( tsig == 0 )\r
-               {                                       \r
-                   m_tempoScratch[i] += ACF[a*(i+1)+b-1] * (1.0 / (2.0 * (double)a-1)) * weight[i];\r
-               }\r
-               else\r
-               {\r
-                   m_tempoScratch[i] += ACF[a*(i+1)+b-1] * 1 * weight[i];\r
-               }\r
-           }\r
-       }\r
-    }\r
-\r
-\r
-       //////////////////////////////////////////////////\r
-       // MODIFIED BEAT PERIOD EXTRACTION //////////////\r
-       /////////////////////////////////////////////////\r
-\r
-       // find smoothed version of RCF ( as applied to Detection Function)\r
-       m_RCFConditioning->process( m_tempoScratch, m_smoothRCF);\r
-\r
-       if (tsig != 0) // i.e. in context dependent state\r
-       {       \r
-//     NOW FIND MAX INDEX OF ACFOUT\r
-            for( i = 0; i < m_lagLength; i++)\r
-            {\r
-                if( m_tempoScratch[ i ] > maxValRCF)\r
-                {\r
-                    maxValRCF = m_tempoScratch[ i ];\r
-                    maxIndexRCF = i;\r
-                }\r
-            }\r
-       }\r
-       else // using rayleigh weighting\r
-       {\r
-               vector <vector<double> > rcfMat;\r
-       \r
-               double sumRcf = 0.;\r
-       \r
-               double maxVal = 0.;\r
-               // now find the two values which minimise rcfMat\r
-               double minVal = 0.;\r
-               int p_i = 1; // periodicity for row i;\r
-               int p_j = 1; //periodicity for column j;\r
-       \r
-       \r
-               for ( i=0; i<m_lagLength; i++)\r
-               {\r
-                       m_tempoScratch[i] =m_smoothRCF[i];\r
-               }       \r
-\r
-               // normalise m_tempoScratch so that it sums to zero.\r
-               for ( i=0; i<m_lagLength; i++)\r
-               {\r
-                       sumRcf += m_tempoScratch[i];\r
-               }       \r
-       \r
-               for( i=0; i<m_lagLength; i++)\r
-               {\r
-                       m_tempoScratch[i] /= sumRcf;\r
-               }       \r
-       \r
-               // create a matrix to store m_tempoScratchValues modified by log2 ratio\r
-               for ( i=0; i<m_lagLength; i++)\r
-               {\r
-                       rcfMat.push_back  ( vector<double>() ); // adds a new row...\r
-               }\r
-       \r
-               for (i=0; i<m_lagLength; i++)\r
-               {\r
-                       for (j=0; j<m_lagLength; j++)\r
-                       {\r
-                               rcfMat[i].push_back (0.);\r
-                       }\r
-               }\r
-       \r
-               // the 'i' and 'j' indices deliberately start from '1' and not '0'\r
-               for ( i=1; i<m_lagLength; i++)\r
-               {\r
-                       for (j=1; j<m_lagLength; j++)\r
-                       {\r
-                               double log2PeriodRatio = log( static_cast<double>(i)/static_cast<double>(j) ) / log(2.0);\r
-                               rcfMat[i][j] = ( abs(1.0-abs(log2PeriodRatio)) );\r
-                               rcfMat[i][j] += ( 0.01*( 1./(m_tempoScratch[i]+m_tempoScratch[j]) ) );\r
-                       }\r
-               }\r
-               \r
-               // set diagonal equal to maximum value in rcfMat \r
-               // we don't want to pick one strong middle peak - we need a combination of two peaks.\r
-       \r
-               for ( i=1; i<m_lagLength; i++)\r
-               {\r
-                       for (j=1; j<m_lagLength; j++)\r
-                       {\r
-                               if (rcfMat[i][j] > maxVal)\r
-                               {       \r
-                                       maxVal = rcfMat[i][j];\r
-                               }\r
-                       }\r
-               }\r
-       \r
-               for ( i=1; i<m_lagLength; i++)\r
-               {\r
-                       rcfMat[i][i] = maxVal;\r
-               }\r
-       \r
-               // now find the row and column number which minimise rcfMat\r
-               minVal = maxVal;\r
-               \r
-               for ( i=1; i<m_lagLength; i++)\r
-               {\r
-                       for ( j=1; j<m_lagLength; j++)\r
-                       {\r
-                               if (rcfMat[i][j] < minVal)\r
-                               {       \r
-                                       minVal = rcfMat[i][j];\r
-                                       p_i = i;\r
-                                       p_j = j;\r
-                               }\r
-                       }\r
-               }\r
-       \r
-       \r
-               // initially choose p_j (arbitrary) - saves on an else statement\r
-               int beatPeriod = p_j;\r
-               if (m_tempoScratch[p_i] > m_tempoScratch[p_j])\r
-               {\r
-                       beatPeriod = p_i;\r
-               }\r
-               \r
-               // now write the output\r
-               maxIndexRCF = static_cast<int>(beatPeriod);\r
-       }\r
-\r
-\r
-    double locked = 5168.f / maxIndexRCF;\r
-    if (locked >= 30 && locked <= 180) {\r
-        m_lockedTempo = locked;\r
-    }\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-    std::cerr << "tempoMM: locked tempo = " << m_lockedTempo << std::endl;\r
-#endif\r
-\r
-    if( tsig == 0 )\r
-       tsig = 4;\r
-\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-std::cerr << "tempoMM: maxIndexRCF = " << maxIndexRCF << std::endl;\r
-#endif\r
-       \r
-    if( tsig == 4 )\r
-    {\r
-#ifdef DEBUG_TEMPO_TRACK\r
-        std::cerr << "tsig == 4" << std::endl;\r
-#endif\r
-\r
-       pdPeaks = new double[ 4 ];\r
-       for( i = 0; i < 4; i++ ){ pdPeaks[ i ] = 0.0;}\r
-\r
-       pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;\r
-\r
-       maxIndexTemp = 0;\r
-       maxValTemp = 0.0;\r
-       count = 0;\r
-\r
-       for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )\r
-       {\r
-           if( ACF[ i ] > maxValTemp )\r
-           {\r
-               maxValTemp = ACF[ i ];\r
-               maxIndexTemp = count;\r
-           }\r
-           count++;\r
-       }\r
-       pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;\r
-\r
-       maxIndexTemp = 0;\r
-       maxValTemp = 0.0;\r
-       count = 0;\r
-\r
-       for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )\r
-       {\r
-           if( ACF[ i ] > maxValTemp )\r
-           {\r
-               maxValTemp = ACF[ i ];\r
-               maxIndexTemp = count;\r
-           }\r
-           count++;\r
-       }\r
-       pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;\r
-\r
-       maxIndexTemp = 0;\r
-       maxValTemp = 0.0;\r
-       count = 0;\r
-\r
-       for( i = ( 4 * maxIndexRCF + 3) - 3; i < ( 4 * maxIndexRCF + 3) + 4; i++ )\r
-       {\r
-           if( ACF[ i ] > maxValTemp )\r
-           {\r
-               maxValTemp = ACF[ i ];\r
-               maxIndexTemp = count;\r
-           }\r
-           count++;\r
-       }\r
-       pdPeaks[ 3 ] = (double)( maxIndexTemp + 1 + ( (4 * maxIndexRCF + 3) - 9 ) + 1 )/4 ;\r
-\r
-\r
-       period = MathUtilities::mean( pdPeaks, 4 );\r
-    }\r
-    else\r
-    { \r
-#ifdef DEBUG_TEMPO_TRACK\r
-       std::cerr << "tsig != 4" << std::endl;\r
-#endif\r
-\r
-       pdPeaks = new double[ 3 ];\r
-       for( i = 0; i < 3; i++ ){ pdPeaks[ i ] = 0.0;}\r
-\r
-       pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;\r
-\r
-       maxIndexTemp = 0;\r
-       maxValTemp = 0.0;\r
-       count = 0;\r
-\r
-       for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )\r
-       {\r
-           if( ACF[ i ] > maxValTemp )\r
-           {\r
-               maxValTemp = ACF[ i ];\r
-               maxIndexTemp = count;\r
-           }\r
-           count++;\r
-       }\r
-       pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;\r
-\r
-       maxIndexTemp = 0;\r
-       maxValTemp = 0.0;\r
-       count = 0;\r
-\r
-       for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )\r
-       {\r
-           if( ACF[ i ] > maxValTemp )\r
-           {\r
-               maxValTemp = ACF[ i ];\r
-               maxIndexTemp = count;\r
-           }\r
-           count++;\r
-       }\r
-       pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;\r
-\r
-\r
-       period = MathUtilities::mean( pdPeaks, 3 );\r
-    }\r
-\r
-    delete [] pdPeaks;\r
-\r
-    return period;\r
-}\r
-\r
-void TempoTrack::stepDetect( double* periodP, double* periodG, int currentIdx, int* flag )\r
-{\r
-    double stepthresh = 1 * 3.9017;\r
-\r
-    if( *flag )\r
-    {\r
-       if(abs(periodG[ currentIdx ] - periodP[ currentIdx ]) > stepthresh)\r
-       {\r
-           // do nuffin'\r
-       }\r
-    }\r
-    else\r
-    {\r
-       if(fabs(periodG[ currentIdx ]-periodP[ currentIdx ]) > stepthresh)\r
-       {\r
-           *flag = 3;\r
-       }\r
-    }\r
-}\r
-\r
-void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag )\r
-{\r
-    double constthresh = 2 * 3.9017;\r
-\r
-    if( fabs( 2 * periodP[ currentIdx ] - periodP[ currentIdx - 1] - periodP[ currentIdx - 2] ) < constthresh)\r
-    {\r
-       *flag = 1;\r
-    }\r
-    else\r
-    {\r
-       *flag = 0;\r
-    }\r
-}\r
-\r
-int TempoTrack::findMeter(double *ACF, unsigned int len, double period)\r
-{\r
-    int i;\r
-    int p = (int)MathUtilities::round( period );\r
-    int tsig;\r
-\r
-    double Energy_3 = 0.0;\r
-    double Energy_4 = 0.0;\r
-\r
-    double temp3A = 0.0;\r
-    double temp3B = 0.0;\r
-    double temp4A = 0.0;\r
-    double temp4B = 0.0;\r
-\r
-    double* dbf = new double[ len ]; int t = 0;\r
-    for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; }\r
-\r
-    if( (double)len < 6 * p + 2 )\r
-    {\r
-       for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp3A += ACF[ i ];\r
-           dbf[ t++ ] = ACF[ i ];\r
-       }\r
-       \r
-       for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp4A += ACF[ i ];\r
-       }\r
-\r
-       Energy_3 = temp3A;\r
-       Energy_4 = temp4A;\r
-    }\r
-    else\r
-    {\r
-       for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp3A += ACF[ i ];\r
-       }\r
-       \r
-       for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp4A += ACF[ i ];\r
-       }\r
-\r
-       for( i = ( 6 * p - 2 ); i < ( 6 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp3B += ACF[ i ];\r
-       }\r
-       \r
-       for( i = ( 2 * p - 2 ); i < ( 2 * p + 2 ) + 1; i++ )\r
-       {\r
-           temp4B += ACF[ i ];\r
-       }\r
-\r
-       Energy_3 = temp3A + temp3B;\r
-       Energy_4 = temp4A + temp4B;\r
-    }\r
-\r
-    if (Energy_3 > Energy_4)\r
-    {\r
-       tsig = 3;\r
-    }\r
-    else\r
-    {\r
-       tsig = 4;\r
-    }\r
-\r
-\r
-    return tsig;\r
-}\r
-\r
-void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat)\r
-{      \r
-    int p = (int)MathUtilities::round( period );\r
-    int predictedOffset = 0;\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-    std::cerr << "TempoTrack::createPhaseExtractor: period = " << period << ", p = " << p << std::endl;\r
-#endif\r
-\r
-    if (p > 10000) {\r
-        std::cerr << "TempoTrack::createPhaseExtractor: WARNING! Highly implausible period value " << p << "!" << std::endl;\r
-        period = 5168 / 120;\r
-    }\r
-\r
-    double* phaseScratch = new double[ p*2 + 2 ];\r
-    for (int i = 0; i < p*2 + 2; ++i) phaseScratch[i] = 0.0;\r
-\r
-       \r
-    if( lastBeat != 0 )\r
-    {\r
-       lastBeat = (int)MathUtilities::round((double)lastBeat );///(double)winLength);\r
-\r
-        predictedOffset = lastBeat + p - fsp;\r
-\r
-        if (predictedOffset < 0) \r
-        {\r
-            lastBeat = 0;\r
-        }\r
-    }\r
-\r
-    if( lastBeat != 0 )\r
-    {\r
-       int mu = p;\r
-       double sigma = (double)p/8;\r
-       double PhaseMin = 0.0;\r
-       double PhaseMax = 0.0;\r
-       unsigned int scratchLength = p*2;\r
-       double temp = 0.0;\r
-\r
-       for(  int i = 0; i < scratchLength; i++ )\r
-       {\r
-           phaseScratch[ i ] = exp( -0.5 * pow( ( i - mu ) / sigma, 2 ) ) / ( sqrt( 2*PI ) *sigma );\r
-       }\r
-\r
-       MathUtilities::getFrameMinMax( phaseScratch, scratchLength, &PhaseMin, &PhaseMax );\r
-                       \r
-       for(int i = 0; i < scratchLength; i ++)\r
-       {\r
-           temp = phaseScratch[ i ];\r
-           phaseScratch[ i ] = (temp - PhaseMin)/PhaseMax;\r
-       }\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-        std::cerr << "predictedOffset = " << predictedOffset << std::endl;\r
-#endif\r
-\r
-       unsigned int index = 0;\r
-       for (int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++)\r
-       {\r
-#ifdef DEBUG_TEMPO_TRACK\r
-            std::cerr << "assigning to filter index " << index << " (size = " << p*2 << ")" << " value " << phaseScratch[i] << " from scratch index " << i << std::endl;\r
-#endif\r
-           Filter[ index++ ] = phaseScratch[ i ];\r
-       }\r
-    }\r
-    else\r
-    {\r
-       for( int i = 0; i < p; i ++)\r
-       {\r
-           Filter[ i ] = 1;\r
-       }\r
-    }\r
-       \r
-    delete [] phaseScratch;\r
-}\r
-\r
-int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period)\r
-{\r
-    int alignment = 0;\r
-    int p = (int)MathUtilities::round( period );\r
-\r
-    double temp = 0.0;\r
-\r
-    double* y = new double[ winLength ];\r
-    double* align = new double[ p ];\r
-\r
-    for( int i = 0; i < winLength; i++ )\r
-    {  \r
-       y[ i ] = (double)( -i + winLength  )/(double)winLength;\r
-       y[ i ] = pow(y [i ],2.0); // raise to power 2.\r
-    }\r
-\r
-    for( int o = 0; o < p; o++ )\r
-    { \r
-       temp = 0.0;\r
-       for(int i = 1 + (o - 1); i< winLength; i += (p + 1))\r
-       {\r
-           temp = temp + DF[ i ] * y[ i ]; \r
-       }\r
-       align[ o ] = temp * weighting[ o ];       \r
-    }\r
-\r
-\r
-    double valTemp = 0.0;\r
-    for(int i = 0; i < p; i++)\r
-    {\r
-       if( align[ i ] > valTemp )\r
-       {\r
-           valTemp = align[ i ];\r
-           alignment = i;\r
-       }\r
-    }\r
-\r
-    delete [] y;\r
-    delete [] align;\r
-\r
-    return alignment;\r
-}\r
-\r
-int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step )\r
-{\r
-    int beat = 0;\r
-\r
-    int p = (int)MathUtilities::round( period );\r
-    int align = (int)MathUtilities::round( alignment );\r
-    int FSP = (int)MathUtilities::round( FSP0 );\r
-\r
-    int FEP = FSP + ( step );\r
-\r
-    beat = FSP + align;\r
-\r
-    m_beats.push_back( beat );\r
-\r
-    while( beat + p < FEP )\r
-    {\r
-       beat += p;\r
-               \r
-       m_beats.push_back( beat );\r
-    }\r
-\r
-    return beat;\r
-}\r
-\r
-\r
-\r
-vector<int> TempoTrack::process( vector <double> DF,\r
-                                 vector <double> *tempoReturn )\r
-{\r
-    m_dataLength = DF.size();\r
-       \r
-    m_lockedTempo = 0.0;\r
-\r
-    double     period = 0.0;\r
-    int stepFlag = 0;\r
-    int constFlag = 0;\r
-    int FSP = 0;\r
-    int tsig = 0;\r
-    int lastBeat = 0;\r
-\r
-    vector <double> causalDF;\r
-\r
-    causalDF = DF;\r
-\r
-    //Prepare Causal Extension DFData\r
-    unsigned int DFCLength = m_dataLength + m_winLength;\r
-       \r
-    for( unsigned int j = 0; j < m_winLength; j++ )\r
-    {\r
-       causalDF.push_back( 0 );\r
-    }\r
-       \r
-       \r
-    double* RW = new double[ m_lagLength ];\r
-    for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}\r
-\r
-    double* GW = new double[ m_lagLength ];\r
-    for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}\r
-\r
-    double* PW = new double[ m_lagLength ];\r
-    for(unsigned clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}\r
-\r
-    m_DFFramer.setSource( &causalDF[0], m_dataLength );\r
-\r
-    unsigned int TTFrames = m_DFFramer.getMaxNoFrames();\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-    std::cerr << "TTFrames = " << TTFrames << std::endl;\r
-#endif\r
-       \r
-    double* periodP = new double[ TTFrames ];\r
-    for(unsigned clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}\r
-       \r
-    double* periodG = new double[ TTFrames ];\r
-    for(unsigned clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}\r
-       \r
-    double* alignment = new double[ TTFrames ];\r
-    for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}\r
-\r
-    m_beats.clear();\r
-\r
-    createCombFilter( RW, m_lagLength, 0, 0 );\r
-\r
-    int TTLoopIndex = 0;\r
-\r
-    for( unsigned int i = 0; i < TTFrames; i++ )\r
-    {\r
-       m_DFFramer.getFrame( m_rawDFFrame );\r
-\r
-       m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame );\r
-\r
-       m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength );\r
-               \r
-       periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 );\r
-\r
-       if( GW[ 0 ] != 0 )\r
-       {\r
-           periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig );\r
-       }\r
-       else\r
-       {\r
-           periodG[ TTLoopIndex ] = 0.0;\r
-       }\r
-\r
-       stepDetect( periodP, periodG, TTLoopIndex, &stepFlag );\r
-\r
-       if( stepFlag == 1)\r
-       {\r
-           constDetect( periodP, TTLoopIndex, &constFlag );\r
-           stepFlag = 0;\r
-       }\r
-       else\r
-       {\r
-           stepFlag -= 1;\r
-       }\r
-\r
-       if( stepFlag < 0 )\r
-       {\r
-           stepFlag = 0;\r
-       }\r
-\r
-       if( constFlag != 0)\r
-       {\r
-           tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] );\r
-       \r
-           createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] );\r
-                       \r
-           periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); \r
-\r
-           period = periodG[ TTLoopIndex ];\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-            std::cerr << "TempoTrack::process: constFlag == " << constFlag << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;\r
-#endif\r
-\r
-           createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); \r
-\r
-           constFlag = 0;\r
-\r
-       }\r
-       else\r
-       {\r
-           if( GW[ 0 ] != 0 )\r
-           {\r
-               period = periodG[ TTLoopIndex ];\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-                std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;\r
-#endif\r
-\r
-                if (period > 10000) {\r
-                    std::cerr << "TempoTrack::process: WARNING!  Highly implausible period value " << period << "!" << std::endl;\r
-                    std::cerr << "periodG contains (of " << TTFrames << " frames): " << std::endl;\r
-                    for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {\r
-                        std::cerr << i << " -> " << periodG[i] << std::endl;\r
-                    }\r
-                    std::cerr << "periodP contains (of " << TTFrames << " frames): " << std::endl;\r
-                    for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {\r
-                        std::cerr << i << " -> " << periodP[i] << std::endl;\r
-                    }\r
-                    period = 5168 / 120;\r
-                }\r
-\r
-               createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); \r
-\r
-           }\r
-           else\r
-           {\r
-               period = periodP[ TTLoopIndex ];\r
-\r
-#ifdef DEBUG_TEMPO_TRACK\r
-                std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodP = " << period << std::endl;\r
-#endif\r
-\r
-               createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); \r
-           }\r
-       }\r
-\r
-       alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); \r
-\r
-       lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength );\r
-\r
-       FSP += (m_lagLength);\r
-\r
-        if (tempoReturn) tempoReturn->push_back(m_lockedTempo);\r
-\r
-       TTLoopIndex++;\r
-    }\r
-\r
-\r
-    delete [] periodP;\r
-    delete [] periodG;\r
-    delete [] alignment;\r
-\r
-    delete [] RW;\r
-    delete [] GW;\r
-    delete [] PW;\r
-\r
-    return m_beats;\r
-}\r
+    COPYING included with this distribution for more information.
+*/
+
+#include "TempoTrack.h"
+
+#include "maths/MathAliases.h"
+#include "maths/MathUtilities.h"
+
+#include <iostream>
+
+#include <cassert>
+
+//#define DEBUG_TEMPO_TRACK 1
+
+
+#define RAY43VAL
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+TempoTrack::TempoTrack( TTParams Params )
+{
+    m_tempoScratch = NULL;
+    m_rawDFFrame = NULL;
+    m_smoothDFFrame = NULL;
+    m_frameACF = NULL;
+       m_smoothRCF = 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 ];
+       m_smoothRCF = 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 );
+
+
+       // these are parameters for smoothing m_tempoScratch
+    m_RCFPParams.length = m_lagLength;
+    m_RCFPParams.AlphaNormParam = Params.alpha;
+    m_RCFPParams.LPOrd = Params.LPOrd;
+    m_RCFPParams.LPACoeffs = Params.LPACoeffs;
+    m_RCFPParams.LPBCoeffs = Params.LPBCoeffs;
+    m_RCFPParams.winPre = Params.WinT.pre;
+    m_RCFPParams.winPost = Params.WinT.post;
+    m_RCFPParams.isMedianPositive = true;
+
+    m_RCFConditioning = new DFProcess( m_RCFPParams );
+
+}
+
+void TempoTrack::deInitialise()
+{      
+    delete [] m_rawDFFrame;
+       
+    delete [] m_smoothDFFrame;
+
+       delete [] m_smoothRCF;  
+       
+    delete [] m_frameACF;
+
+    delete [] m_tempoScratch;
+
+    delete m_DFConditioning;
+       
+       delete m_RCFConditioning;
+
+}
+
+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/4;
+       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,i,j;
+    int 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;
+    }
+
+#ifdef DEBUG_TEMPO_TRACK
+    std::cerr << "tempoMM: m_winLength = " << m_winLength << ", m_lagLength = " << m_lagLength << ", numelem = " << numelem << std::endl;
+#endif
+
+    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];
+               }
+           }
+       }
+    }
+
+
+       //////////////////////////////////////////////////
+       // MODIFIED BEAT PERIOD EXTRACTION //////////////
+       /////////////////////////////////////////////////
+
+       // find smoothed version of RCF ( as applied to Detection Function)
+       m_RCFConditioning->process( m_tempoScratch, m_smoothRCF);
+
+       if (tsig != 0) // i.e. in context dependent state
+       {       
+//     NOW FIND MAX INDEX OF ACFOUT
+            for( i = 0; i < m_lagLength; i++)
+            {
+                if( m_tempoScratch[ i ] > maxValRCF)
+                {
+                    maxValRCF = m_tempoScratch[ i ];
+                    maxIndexRCF = i;
+                }
+            }
+       }
+       else // using rayleigh weighting
+       {
+               vector <vector<double> > rcfMat;
+       
+               double sumRcf = 0.;
+       
+               double maxVal = 0.;
+               // now find the two values which minimise rcfMat
+               double minVal = 0.;
+               int p_i = 1; // periodicity for row i;
+               int p_j = 1; //periodicity for column j;
+       
+       
+               for ( i=0; i<m_lagLength; i++)
+               {
+                       m_tempoScratch[i] =m_smoothRCF[i];
+               }       
+
+               // normalise m_tempoScratch so that it sums to zero.
+               for ( i=0; i<m_lagLength; i++)
+               {
+                       sumRcf += m_tempoScratch[i];
+               }       
+       
+               for( i=0; i<m_lagLength; i++)
+               {
+                       m_tempoScratch[i] /= sumRcf;
+               }       
+       
+               // create a matrix to store m_tempoScratchValues modified by log2 ratio
+               for ( i=0; i<m_lagLength; i++)
+               {
+                       rcfMat.push_back  ( vector<double>() ); // adds a new row...
+               }
+       
+               for (i=0; i<m_lagLength; i++)
+               {
+                       for (j=0; j<m_lagLength; j++)
+                       {
+                               rcfMat[i].push_back (0.);
+                       }
+               }
+       
+               // the 'i' and 'j' indices deliberately start from '1' and not '0'
+               for ( i=1; i<m_lagLength; i++)
+               {
+                       for (j=1; j<m_lagLength; j++)
+                       {
+                               double log2PeriodRatio = log( static_cast<double>(i)/static_cast<double>(j) ) / log(2.0);
+                               rcfMat[i][j] = ( abs(1.0-abs(log2PeriodRatio)) );
+                               rcfMat[i][j] += ( 0.01*( 1./(m_tempoScratch[i]+m_tempoScratch[j]) ) );
+                       }
+               }
+               
+               // set diagonal equal to maximum value in rcfMat 
+               // we don't want to pick one strong middle peak - we need a combination of two peaks.
+       
+               for ( i=1; i<m_lagLength; i++)
+               {
+                       for (j=1; j<m_lagLength; j++)
+                       {
+                               if (rcfMat[i][j] > maxVal)
+                               {       
+                                       maxVal = rcfMat[i][j];
+                               }
+                       }
+               }
+       
+               for ( i=1; i<m_lagLength; i++)
+               {
+                       rcfMat[i][i] = maxVal;
+               }
+       
+               // now find the row and column number which minimise rcfMat
+               minVal = maxVal;
+               
+               for ( i=1; i<m_lagLength; i++)
+               {
+                       for ( j=1; j<m_lagLength; j++)
+                       {
+                               if (rcfMat[i][j] < minVal)
+                               {       
+                                       minVal = rcfMat[i][j];
+                                       p_i = i;
+                                       p_j = j;
+                               }
+                       }
+               }
+       
+       
+               // initially choose p_j (arbitrary) - saves on an else statement
+               int beatPeriod = p_j;
+               if (m_tempoScratch[p_i] > m_tempoScratch[p_j])
+               {
+                       beatPeriod = p_i;
+               }
+               
+               // now write the output
+               maxIndexRCF = static_cast<int>(beatPeriod);
+       }
+
+
+    double locked = 5168.f / maxIndexRCF;
+    if (locked >= 30 && locked <= 180) {
+        m_lockedTempo = locked;
+    }
+
+#ifdef DEBUG_TEMPO_TRACK
+    std::cerr << "tempoMM: locked tempo = " << m_lockedTempo << std::endl;
+#endif
+
+    if( tsig == 0 )
+       tsig = 4;
+
+
+#ifdef DEBUG_TEMPO_TRACK
+std::cerr << "tempoMM: maxIndexRCF = " << maxIndexRCF << std::endl;
+#endif
+       
+    if( tsig == 4 )
+    {
+#ifdef DEBUG_TEMPO_TRACK
+        std::cerr << "tsig == 4" << std::endl;
+#endif
+
+       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
+    { 
+#ifdef DEBUG_TEMPO_TRACK
+       std::cerr << "tsig != 4" << std::endl;
+#endif
+
+       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;
+
+#ifdef DEBUG_TEMPO_TRACK
+    std::cerr << "TempoTrack::createPhaseExtractor: period = " << period << ", p = " << p << std::endl;
+#endif
+
+    if (p > 10000) {
+        std::cerr << "TempoTrack::createPhaseExtractor: WARNING! Highly implausible period value " << p << "!" << std::endl;
+        period = 5168 / 120;
+    }
+
+    double* phaseScratch = new double[ p*2 + 2 ];
+    for (int i = 0; i < p*2 + 2; ++i) phaseScratch[i] = 0.0;
+
+       
+    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/8;
+       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;
+       }
+
+#ifdef DEBUG_TEMPO_TRACK
+        std::cerr << "predictedOffset = " << predictedOffset << std::endl;
+#endif
+
+       unsigned int index = 0;
+       for (int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++)
+       {
+#ifdef DEBUG_TEMPO_TRACK
+            std::cerr << "assigning to filter index " << index << " (size = " << p*2 << ")" << " value " << phaseScratch[i] << " from scratch index " << i << std::endl;
+#endif
+           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;
+       y[ i ] = pow(y [i ],2.0); // raise to power 2.
+    }
+
+    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( vector <double> DF,
+                                 vector <double> *tempoReturn )
+{
+    m_dataLength = DF.size();
+       
+    m_lockedTempo = 0.0;
+
+    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();
+
+#ifdef DEBUG_TEMPO_TRACK
+    std::cerr << "TTFrames = " << TTFrames << std::endl;
+#endif
+       
+    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 ];
+
+#ifdef DEBUG_TEMPO_TRACK
+            std::cerr << "TempoTrack::process: constFlag == " << constFlag << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;
+#endif
+
+           createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); 
+
+           constFlag = 0;
+
+       }
+       else
+       {
+           if( GW[ 0 ] != 0 )
+           {
+               period = periodG[ TTLoopIndex ];
+
+#ifdef DEBUG_TEMPO_TRACK
+                std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;
+#endif
+
+                if (period > 10000) {
+                    std::cerr << "TempoTrack::process: WARNING!  Highly implausible period value " << period << "!" << std::endl;
+                    std::cerr << "periodG contains (of " << TTFrames << " frames): " << std::endl;
+                    for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {
+                        std::cerr << i << " -> " << periodG[i] << std::endl;
+                    }
+                    std::cerr << "periodP contains (of " << TTFrames << " frames): " << std::endl;
+                    for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {
+                        std::cerr << i << " -> " << periodP[i] << std::endl;
+                    }
+                    period = 5168 / 120;
+                }
+
+               createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); 
+
+           }
+           else
+           {
+               period = periodP[ TTLoopIndex ];
+
+#ifdef DEBUG_TEMPO_TRACK
+                std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodP = " << period << std::endl;
+#endif
+
+               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);
+
+        if (tempoReturn) tempoReturn->push_back(m_lockedTempo);
+
+       TTLoopIndex++;
+    }
+
+
+    delete [] periodP;
+    delete [] periodG;
+    delete [] alignment;
+
+    delete [] RW;
+    delete [] GW;
+    delete [] PW;
+
+    return m_beats;
+}
+