X-Git-Url: https://main.carlh.net/gitweb/?a=blobdiff_plain;f=libs%2Fqm-dsp%2Fdsp%2Ftempotracking%2FTempoTrack.cpp;h=389403edaa8af9ac7832d6e59e31e3ed91017030;hb=f68d2e06bcfb81efda107d3b4c3aa7dbc2d73bc2;hp=c36385f250a0403e02eba7ce0cc548995a10c516;hpb=2a27cc475867612afd261e5bf3b2a1a42b9c75cc;p=ardour.git diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp index c36385f250..389403edaa 100644 --- a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp +++ b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp @@ -1,869 +1,870 @@ -/* -*- 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. +/* -*- 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. -*/ - -#include "TempoTrack.h" - -#include "maths/MathAliases.h" -#include "maths/MathUtilities.h" - -#include - -#include - -//#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;iprocess( 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 > 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() ); // adds a new row... - } - - for (i=0; i(i)/static_cast(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 maxVal) - { - maxVal = rcfMat[i][j]; - } - } - } - - for ( i=1; i m_tempoScratch[p_j]) - { - beatPeriod = p_i; - } - - // now write the output - maxIndexRCF = static_cast(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 TempoTrack::process( vector DF, - vector *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 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; -} + COPYING included with this distribution for more information. +*/ + +#include "TempoTrack.h" + +#include "maths/MathAliases.h" +#include "maths/MathUtilities.h" + +#include + +#include + +//#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;iprocess( 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 > 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() ); // adds a new row... + } + + for (i=0; i(i)/static_cast(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 maxVal) + { + maxVal = rcfMat[i][j]; + } + } + } + + for ( i=1; i m_tempoScratch[p_j]) + { + beatPeriod = p_i; + } + + // now write the output + maxIndexRCF = static_cast(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 TempoTrack::process( vector DF, + vector *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 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; +} +