libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */\r
   2 \r
   3 /*\r
   4     QM DSP Library\r
   5 \r
   6     Centre for Digital Music, Queen Mary, University of London.\r
   7     This file copyright 2005-2006 Christian Landone.and Matthew Davies.\r
   8
   9     This program is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU General Public License as
  11     published by the Free Software Foundation; either version 2 of the
  12     License, or (at your option) any later version.  See the file
  13     COPYING included with this distribution for more information.\r
  14 */\r
  15 \r
  16 #include "TempoTrack.h"\r
  17 \r
  18 #include "maths/MathAliases.h"\r
  19 #include "maths/MathUtilities.h"\r
  20 \r
  21 #include <iostream>\r
  22 \r
  23 #include <cassert>\r
  24 \r
  25 //#define DEBUG_TEMPO_TRACK 1\r
  26 \r
  27 \r
  28 #define RAY43VAL\r
  29 \r
  30 //////////////////////////////////////////////////////////////////////\r
  31 // Construction/Destruction\r
  32 //////////////////////////////////////////////////////////////////////\r
  33 \r
  34 TempoTrack::TempoTrack( TTParams Params )\r
  35 {\r
  36     m_tempoScratch = NULL;\r
  37     m_rawDFFrame = NULL;\r
  38     m_smoothDFFrame = NULL;\r
  39     m_frameACF = NULL;\r
  40         m_smoothRCF = NULL;\r
  41 \r
  42     m_dataLength = 0;\r
  43     m_winLength = 0;\r
  44     m_lagLength = 0;\r
  45 \r
  46     m_rayparam = 0;\r
  47     m_sigma = 0;\r
  48     m_DFWVNnorm = 0;\r
  49 \r
  50     initialise( Params );\r
  51 }\r
  52 \r
  53 TempoTrack::~TempoTrack()\r
  54 {\r
  55     deInitialise();\r
  56 }\r
  57 \r
  58 void TempoTrack::initialise( TTParams Params )\r
  59 {       \r
  60     m_winLength = Params.winLength;\r
  61     m_lagLength = Params.lagLength;\r
  62 \r
  63     m_rayparam   = 43.0;\r
  64     m_sigma = sqrt(3.9017);\r
  65     m_DFWVNnorm = exp( ( log( 2.0 ) / m_rayparam ) * ( m_winLength + 2 ) );\r
  66 \r
  67     m_rawDFFrame = new double[ m_winLength ];\r
  68     m_smoothDFFrame = new double[ m_winLength ];\r
  69     m_frameACF = new double[ m_winLength ];\r
  70     m_tempoScratch = new double[ m_lagLength ];\r
  71         m_smoothRCF = new double[ m_lagLength ];\r
  72 \r
  73 \r
  74     unsigned int winPre = Params.WinT.pre;\r
  75     unsigned int winPost = Params.WinT.post;\r
  76 \r
  77     m_DFFramer.configure( m_winLength, m_lagLength );\r
  78         \r
  79     m_DFPParams.length = m_winLength;\r
  80     m_DFPParams.AlphaNormParam = Params.alpha;\r
  81     m_DFPParams.LPOrd = Params.LPOrd;\r
  82     m_DFPParams.LPACoeffs = Params.LPACoeffs;\r
  83     m_DFPParams.LPBCoeffs = Params.LPBCoeffs;\r
  84     m_DFPParams.winPre = Params.WinT.pre;\r
  85     m_DFPParams.winPost = Params.WinT.post;\r
  86     m_DFPParams.isMedianPositive = true;\r
  87         \r
  88     m_DFConditioning = new DFProcess( m_DFPParams );\r
  89 \r
  90 \r
  91         // these are parameters for smoothing m_tempoScratch\r
  92     m_RCFPParams.length = m_lagLength;\r
  93     m_RCFPParams.AlphaNormParam = Params.alpha;\r
  94     m_RCFPParams.LPOrd = Params.LPOrd;\r
  95     m_RCFPParams.LPACoeffs = Params.LPACoeffs;\r
  96     m_RCFPParams.LPBCoeffs = Params.LPBCoeffs;\r
  97     m_RCFPParams.winPre = Params.WinT.pre;\r
  98     m_RCFPParams.winPost = Params.WinT.post;\r
  99     m_RCFPParams.isMedianPositive = true;\r
 100 \r
 101     m_RCFConditioning = new DFProcess( m_RCFPParams );\r
 102 \r
 103 }\r
 104 \r
 105 void TempoTrack::deInitialise()\r
 106 {       \r
 107     delete [] m_rawDFFrame;\r
 108         \r
 109     delete [] m_smoothDFFrame;\r
 110 \r
 111         delete [] m_smoothRCF;  \r
 112         \r
 113     delete [] m_frameACF;\r
 114 \r
 115     delete [] m_tempoScratch;\r
 116 \r
 117     delete m_DFConditioning;\r
 118         \r
 119         delete m_RCFConditioning;\r
 120 \r
 121 }\r
 122 \r
 123 void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag)\r
 124 {\r
 125     unsigned int i;\r
 126 \r
 127     if( beatLag == 0 )\r
 128     {\r
 129         for( i = 0; i < winLength; i++ )\r
 130         {    \r
 131             Filter[ i ] = ( ( i + 1 ) / pow( m_rayparam, 2.0) ) * exp( ( -pow(( i + 1 ),2.0 ) / ( 2.0 * pow( m_rayparam, 2.0))));\r
 132         }\r
 133     }\r
 134     else\r
 135     {   \r
 136         m_sigma = beatLag/4;\r
 137         for( i = 0; i < winLength; i++ )\r
 138         {\r
 139             double dlag = (double)(i+1) - beatLag;\r
 140             Filter[ i ] =  exp(-0.5 * pow(( dlag / m_sigma), 2.0) ) / (sqrt( 2 * PI) * m_sigma);\r
 141         }\r
 142     }\r
 143 }\r
 144 \r
 145 double TempoTrack::tempoMM(double* ACF, double* weight, int tsig)\r
 146 {\r
 147 \r
 148     double period = 0;\r
 149     double maxValRCF = 0.0;\r
 150     unsigned int maxIndexRCF = 0;\r
 151 \r
 152     double* pdPeaks;\r
 153 \r
 154     unsigned int maxIndexTemp;\r
 155     double      maxValTemp;\r
 156     unsigned int count; \r
 157         \r
 158     unsigned int numelem,i,j;\r
 159     int a, b;\r
 160 \r
 161     for( i = 0; i < m_lagLength; i++ )\r
 162         m_tempoScratch[ i ] = 0.0;\r
 163 \r
 164     if( tsig == 0 ) \r
 165     {\r
 166         //if time sig is unknown, use metrically unbiased version of Filterbank\r
 167         numelem = 4;\r
 168     }\r
 169     else\r
 170     {\r
 171         numelem = tsig;\r
 172     }\r
 173 \r
 174 #ifdef DEBUG_TEMPO_TRACK\r
 175     std::cerr << "tempoMM: m_winLength = " << m_winLength << ", m_lagLength = " << m_lagLength << ", numelem = " << numelem << std::endl;\r
 176 #endif\r
 177 \r
 178     for(i=1;i<m_lagLength-1;i++)\r
 179     {\r
 180         //first and last output values are left intentionally as zero\r
 181         for (a=1;a<=numelem;a++)\r
 182         {\r
 183             for(b=(1-a);b<a;b++)\r
 184             {\r
 185                 if( tsig == 0 )\r
 186                 {                                       \r
 187                     m_tempoScratch[i] += ACF[a*(i+1)+b-1] * (1.0 / (2.0 * (double)a-1)) * weight[i];\r
 188                 }\r
 189                 else\r
 190                 {\r
 191                     m_tempoScratch[i] += ACF[a*(i+1)+b-1] * 1 * weight[i];\r
 192                 }\r
 193             }\r
 194         }\r
 195     }\r
 196 \r
 197 \r
 198         //////////////////////////////////////////////////\r
 199         // MODIFIED BEAT PERIOD EXTRACTION //////////////\r
 200         /////////////////////////////////////////////////\r
 201 \r
 202         // find smoothed version of RCF ( as applied to Detection Function)\r
 203         m_RCFConditioning->process( m_tempoScratch, m_smoothRCF);\r
 204 \r
 205         if (tsig != 0) // i.e. in context dependent state\r
 206         {       \r
 207 //     NOW FIND MAX INDEX OF ACFOUT\r
 208             for( i = 0; i < m_lagLength; i++)\r
 209             {\r
 210                 if( m_tempoScratch[ i ] > maxValRCF)\r
 211                 {\r
 212                     maxValRCF = m_tempoScratch[ i ];\r
 213                     maxIndexRCF = i;\r
 214                 }\r
 215             }\r
 216         }\r
 217         else // using rayleigh weighting\r
 218         {\r
 219                 vector <vector<double> > rcfMat;\r
 220         \r
 221                 double sumRcf = 0.;\r
 222         \r
 223                 double maxVal = 0.;\r
 224                 // now find the two values which minimise rcfMat\r
 225                 double minVal = 0.;\r
 226                 int p_i = 1; // periodicity for row i;\r
 227                 int p_j = 1; //periodicity for column j;\r
 228         \r
 229         \r
 230                 for ( i=0; i<m_lagLength; i++)\r
 231                 {\r
 232                         m_tempoScratch[i] =m_smoothRCF[i];\r
 233                 }       \r
 234 \r
 235                 // normalise m_tempoScratch so that it sums to zero.\r
 236                 for ( i=0; i<m_lagLength; i++)\r
 237                 {\r
 238                         sumRcf += m_tempoScratch[i];\r
 239                 }       \r
 240         \r
 241                 for( i=0; i<m_lagLength; i++)\r
 242                 {\r
 243                         m_tempoScratch[i] /= sumRcf;\r
 244                 }       \r
 245         \r
 246                 // create a matrix to store m_tempoScratchValues modified by log2 ratio\r
 247                 for ( i=0; i<m_lagLength; i++)\r
 248                 {\r
 249                         rcfMat.push_back  ( vector<double>() ); // adds a new row...\r
 250                 }\r
 251         \r
 252                 for (i=0; i<m_lagLength; i++)\r
 253                 {\r
 254                         for (j=0; j<m_lagLength; j++)\r
 255                         {\r
 256                                 rcfMat[i].push_back (0.);\r
 257                         }\r
 258                 }\r
 259         \r
 260                 // the 'i' and 'j' indices deliberately start from '1' and not '0'\r
 261                 for ( i=1; i<m_lagLength; i++)\r
 262                 {\r
 263                         for (j=1; j<m_lagLength; j++)\r
 264                         {\r
 265                                 double log2PeriodRatio = log( static_cast<double>(i)/static_cast<double>(j) ) / log(2.0);\r
 266                                 rcfMat[i][j] = ( abs(1.0-abs(log2PeriodRatio)) );\r
 267                                 rcfMat[i][j] += ( 0.01*( 1./(m_tempoScratch[i]+m_tempoScratch[j]) ) );\r
 268                         }\r
 269                 }\r
 270                 \r
 271                 // set diagonal equal to maximum value in rcfMat \r
 272                 // we don't want to pick one strong middle peak - we need a combination of two peaks.\r
 273         \r
 274                 for ( i=1; i<m_lagLength; i++)\r
 275                 {\r
 276                         for (j=1; j<m_lagLength; j++)\r
 277                         {\r
 278                                 if (rcfMat[i][j] > maxVal)\r
 279                                 {       \r
 280                                         maxVal = rcfMat[i][j];\r
 281                                 }\r
 282                         }\r
 283                 }\r
 284         \r
 285                 for ( i=1; i<m_lagLength; i++)\r
 286                 {\r
 287                         rcfMat[i][i] = maxVal;\r
 288                 }\r
 289         \r
 290                 // now find the row and column number which minimise rcfMat\r
 291                 minVal = maxVal;\r
 292                 \r
 293                 for ( i=1; i<m_lagLength; i++)\r
 294                 {\r
 295                         for ( j=1; j<m_lagLength; j++)\r
 296                         {\r
 297                                 if (rcfMat[i][j] < minVal)\r
 298                                 {       \r
 299                                         minVal = rcfMat[i][j];\r
 300                                         p_i = i;\r
 301                                         p_j = j;\r
 302                                 }\r
 303                         }\r
 304                 }\r
 305         \r
 306         \r
 307                 // initially choose p_j (arbitrary) - saves on an else statement\r
 308                 int beatPeriod = p_j;\r
 309                 if (m_tempoScratch[p_i] > m_tempoScratch[p_j])\r
 310                 {\r
 311                         beatPeriod = p_i;\r
 312                 }\r
 313                 \r
 314                 // now write the output\r
 315                 maxIndexRCF = static_cast<int>(beatPeriod);\r
 316         }\r
 317 \r
 318 \r
 319     double locked = 5168.f / maxIndexRCF;\r
 320     if (locked >= 30 && locked <= 180) {\r
 321         m_lockedTempo = locked;\r
 322     }\r
 323 \r
 324 #ifdef DEBUG_TEMPO_TRACK\r
 325     std::cerr << "tempoMM: locked tempo = " << m_lockedTempo << std::endl;\r
 326 #endif\r
 327 \r
 328     if( tsig == 0 )\r
 329         tsig = 4;\r
 330 \r
 331 \r
 332 #ifdef DEBUG_TEMPO_TRACK\r
 333 std::cerr << "tempoMM: maxIndexRCF = " << maxIndexRCF << std::endl;\r
 334 #endif\r
 335         \r
 336     if( tsig == 4 )\r
 337     {\r
 338 #ifdef DEBUG_TEMPO_TRACK\r
 339         std::cerr << "tsig == 4" << std::endl;\r
 340 #endif\r
 341 \r
 342         pdPeaks = new double[ 4 ];\r
 343         for( i = 0; i < 4; i++ ){ pdPeaks[ i ] = 0.0;}\r
 344 \r
 345         pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;\r
 346 \r
 347         maxIndexTemp = 0;\r
 348         maxValTemp = 0.0;\r
 349         count = 0;\r
 350 \r
 351         for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )\r
 352         {\r
 353             if( ACF[ i ] > maxValTemp )\r
 354             {\r
 355                 maxValTemp = ACF[ i ];\r
 356                 maxIndexTemp = count;\r
 357             }\r
 358             count++;\r
 359         }\r
 360         pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;\r
 361 \r
 362         maxIndexTemp = 0;\r
 363         maxValTemp = 0.0;\r
 364         count = 0;\r
 365 \r
 366         for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )\r
 367         {\r
 368             if( ACF[ i ] > maxValTemp )\r
 369             {\r
 370                 maxValTemp = ACF[ i ];\r
 371                 maxIndexTemp = count;\r
 372             }\r
 373             count++;\r
 374         }\r
 375         pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;\r
 376 \r
 377         maxIndexTemp = 0;\r
 378         maxValTemp = 0.0;\r
 379         count = 0;\r
 380 \r
 381         for( i = ( 4 * maxIndexRCF + 3) - 3; i < ( 4 * maxIndexRCF + 3) + 4; i++ )\r
 382         {\r
 383             if( ACF[ i ] > maxValTemp )\r
 384             {\r
 385                 maxValTemp = ACF[ i ];\r
 386                 maxIndexTemp = count;\r
 387             }\r
 388             count++;\r
 389         }\r
 390         pdPeaks[ 3 ] = (double)( maxIndexTemp + 1 + ( (4 * maxIndexRCF + 3) - 9 ) + 1 )/4 ;\r
 391 \r
 392 \r
 393         period = MathUtilities::mean( pdPeaks, 4 );\r
 394     }\r
 395     else\r
 396     { \r
 397 #ifdef DEBUG_TEMPO_TRACK\r
 398        std::cerr << "tsig != 4" << std::endl;\r
 399 #endif\r
 400 \r
 401         pdPeaks = new double[ 3 ];\r
 402         for( i = 0; i < 3; i++ ){ pdPeaks[ i ] = 0.0;}\r
 403 \r
 404         pdPeaks[ 0 ] = ( double )maxIndexRCF + 1;\r
 405 \r
 406         maxIndexTemp = 0;\r
 407         maxValTemp = 0.0;\r
 408         count = 0;\r
 409 \r
 410         for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ )\r
 411         {\r
 412             if( ACF[ i ] > maxValTemp )\r
 413             {\r
 414                 maxValTemp = ACF[ i ];\r
 415                 maxIndexTemp = count;\r
 416             }\r
 417             count++;\r
 418         }\r
 419         pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2;\r
 420 \r
 421         maxIndexTemp = 0;\r
 422         maxValTemp = 0.0;\r
 423         count = 0;\r
 424 \r
 425         for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ )\r
 426         {\r
 427             if( ACF[ i ] > maxValTemp )\r
 428             {\r
 429                 maxValTemp = ACF[ i ];\r
 430                 maxIndexTemp = count;\r
 431             }\r
 432             count++;\r
 433         }\r
 434         pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3;\r
 435 \r
 436 \r
 437         period = MathUtilities::mean( pdPeaks, 3 );\r
 438     }\r
 439 \r
 440     delete [] pdPeaks;\r
 441 \r
 442     return period;\r
 443 }\r
 444 \r
 445 void TempoTrack::stepDetect( double* periodP, double* periodG, int currentIdx, int* flag )\r
 446 {\r
 447     double stepthresh = 1 * 3.9017;\r
 448 \r
 449     if( *flag )\r
 450     {\r
 451         if(abs(periodG[ currentIdx ] - periodP[ currentIdx ]) > stepthresh)\r
 452         {\r
 453             // do nuffin'\r
 454         }\r
 455     }\r
 456     else\r
 457     {\r
 458         if(fabs(periodG[ currentIdx ]-periodP[ currentIdx ]) > stepthresh)\r
 459         {\r
 460             *flag = 3;\r
 461         }\r
 462     }\r
 463 }\r
 464 \r
 465 void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag )\r
 466 {\r
 467     double constthresh = 2 * 3.9017;\r
 468 \r
 469     if( fabs( 2 * periodP[ currentIdx ] - periodP[ currentIdx - 1] - periodP[ currentIdx - 2] ) < constthresh)\r
 470     {\r
 471         *flag = 1;\r
 472     }\r
 473     else\r
 474     {\r
 475         *flag = 0;\r
 476     }\r
 477 }\r
 478 \r
 479 int TempoTrack::findMeter(double *ACF, unsigned int len, double period)\r
 480 {\r
 481     int i;\r
 482     int p = (int)MathUtilities::round( period );\r
 483     int tsig;\r
 484 \r
 485     double Energy_3 = 0.0;\r
 486     double Energy_4 = 0.0;\r
 487 \r
 488     double temp3A = 0.0;\r
 489     double temp3B = 0.0;\r
 490     double temp4A = 0.0;\r
 491     double temp4B = 0.0;\r
 492 \r
 493     double* dbf = new double[ len ]; int t = 0;\r
 494     for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; }\r
 495 \r
 496     if( (double)len < 6 * p + 2 )\r
 497     {\r
 498         for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )\r
 499         {\r
 500             temp3A += ACF[ i ];\r
 501             dbf[ t++ ] = ACF[ i ];\r
 502         }\r
 503         \r
 504         for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )\r
 505         {\r
 506             temp4A += ACF[ i ];\r
 507         }\r
 508 \r
 509         Energy_3 = temp3A;\r
 510         Energy_4 = temp4A;\r
 511     }\r
 512     else\r
 513     {\r
 514         for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ )\r
 515         {\r
 516             temp3A += ACF[ i ];\r
 517         }\r
 518         \r
 519         for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ )\r
 520         {\r
 521             temp4A += ACF[ i ];\r
 522         }\r
 523 \r
 524         for( i = ( 6 * p - 2 ); i < ( 6 * p + 2 ) + 1; i++ )\r
 525         {\r
 526             temp3B += ACF[ i ];\r
 527         }\r
 528         \r
 529         for( i = ( 2 * p - 2 ); i < ( 2 * p + 2 ) + 1; i++ )\r
 530         {\r
 531             temp4B += ACF[ i ];\r
 532         }\r
 533 \r
 534         Energy_3 = temp3A + temp3B;\r
 535         Energy_4 = temp4A + temp4B;\r
 536     }\r
 537 \r
 538     if (Energy_3 > Energy_4)\r
 539     {\r
 540         tsig = 3;\r
 541     }\r
 542     else\r
 543     {\r
 544         tsig = 4;\r
 545     }\r
 546 \r
 547 \r
 548     return tsig;\r
 549 }\r
 550 \r
 551 void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat)\r
 552 {       \r
 553     int p = (int)MathUtilities::round( period );\r
 554     int predictedOffset = 0;\r
 555 \r
 556 #ifdef DEBUG_TEMPO_TRACK\r
 557     std::cerr << "TempoTrack::createPhaseExtractor: period = " << period << ", p = " << p << std::endl;\r
 558 #endif\r
 559 \r
 560     if (p > 10000) {\r
 561         std::cerr << "TempoTrack::createPhaseExtractor: WARNING! Highly implausible period value " << p << "!" << std::endl;\r
 562         period = 5168 / 120;\r
 563     }\r
 564 \r
 565     double* phaseScratch = new double[ p*2 + 2 ];\r
 566     for (int i = 0; i < p*2 + 2; ++i) phaseScratch[i] = 0.0;\r
 567 \r
 568         \r
 569     if( lastBeat != 0 )\r
 570     {\r
 571         lastBeat = (int)MathUtilities::round((double)lastBeat );///(double)winLength);\r
 572 \r
 573         predictedOffset = lastBeat + p - fsp;\r
 574 \r
 575         if (predictedOffset < 0) \r
 576         {\r
 577             lastBeat = 0;\r
 578         }\r
 579     }\r
 580 \r
 581     if( lastBeat != 0 )\r
 582     {\r
 583         int mu = p;\r
 584         double sigma = (double)p/8;\r
 585         double PhaseMin = 0.0;\r
 586         double PhaseMax = 0.0;\r
 587         unsigned int scratchLength = p*2;\r
 588         double temp = 0.0;\r
 589 \r
 590         for(  int i = 0; i < scratchLength; i++ )\r
 591         {\r
 592             phaseScratch[ i ] = exp( -0.5 * pow( ( i - mu ) / sigma, 2 ) ) / ( sqrt( 2*PI ) *sigma );\r
 593         }\r
 594 \r
 595         MathUtilities::getFrameMinMax( phaseScratch, scratchLength, &PhaseMin, &PhaseMax );\r
 596                         \r
 597         for(int i = 0; i < scratchLength; i ++)\r
 598         {\r
 599             temp = phaseScratch[ i ];\r
 600             phaseScratch[ i ] = (temp - PhaseMin)/PhaseMax;\r
 601         }\r
 602 \r
 603 #ifdef DEBUG_TEMPO_TRACK\r
 604         std::cerr << "predictedOffset = " << predictedOffset << std::endl;\r
 605 #endif\r
 606 \r
 607         unsigned int index = 0;\r
 608         for (int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++)\r
 609         {\r
 610 #ifdef DEBUG_TEMPO_TRACK\r
 611             std::cerr << "assigning to filter index " << index << " (size = " << p*2 << ")" << " value " << phaseScratch[i] << " from scratch index " << i << std::endl;\r
 612 #endif\r
 613             Filter[ index++ ] = phaseScratch[ i ];\r
 614         }\r
 615     }\r
 616     else\r
 617     {\r
 618         for( int i = 0; i < p; i ++)\r
 619         {\r
 620             Filter[ i ] = 1;\r
 621         }\r
 622     }\r
 623         \r
 624     delete [] phaseScratch;\r
 625 }\r
 626 \r
 627 int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period)\r
 628 {\r
 629     int alignment = 0;\r
 630     int p = (int)MathUtilities::round( period );\r
 631 \r
 632     double temp = 0.0;\r
 633 \r
 634     double* y = new double[ winLength ];\r
 635     double* align = new double[ p ];\r
 636 \r
 637     for( int i = 0; i < winLength; i++ )\r
 638     {   \r
 639         y[ i ] = (double)( -i + winLength  )/(double)winLength;\r
 640         y[ i ] = pow(y [i ],2.0); // raise to power 2.\r
 641     }\r
 642 \r
 643     for( int o = 0; o < p; o++ )\r
 644     { \r
 645         temp = 0.0;\r
 646         for(int i = 1 + (o - 1); i< winLength; i += (p + 1))\r
 647         {\r
 648             temp = temp + DF[ i ] * y[ i ]; \r
 649         }\r
 650         align[ o ] = temp * weighting[ o ];       \r
 651     }\r
 652 \r
 653 \r
 654     double valTemp = 0.0;\r
 655     for(int i = 0; i < p; i++)\r
 656     {\r
 657         if( align[ i ] > valTemp )\r
 658         {\r
 659             valTemp = align[ i ];\r
 660             alignment = i;\r
 661         }\r
 662     }\r
 663 \r
 664     delete [] y;\r
 665     delete [] align;\r
 666 \r
 667     return alignment;\r
 668 }\r
 669 \r
 670 int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step )\r
 671 {\r
 672     int beat = 0;\r
 673 \r
 674     int p = (int)MathUtilities::round( period );\r
 675     int align = (int)MathUtilities::round( alignment );\r
 676     int FSP = (int)MathUtilities::round( FSP0 );\r
 677 \r
 678     int FEP = FSP + ( step );\r
 679 \r
 680     beat = FSP + align;\r
 681 \r
 682     m_beats.push_back( beat );\r
 683 \r
 684     while( beat + p < FEP )\r
 685     {\r
 686         beat += p;\r
 687                 \r
 688         m_beats.push_back( beat );\r
 689     }\r
 690 \r
 691     return beat;\r
 692 }\r
 693 \r
 694 \r
 695 \r
 696 vector<int> TempoTrack::process( vector <double> DF,\r
 697                                  vector <double> *tempoReturn )\r
 698 {\r
 699     m_dataLength = DF.size();\r
 700         \r
 701     m_lockedTempo = 0.0;\r
 702 \r
 703     double      period = 0.0;\r
 704     int stepFlag = 0;\r
 705     int constFlag = 0;\r
 706     int FSP = 0;\r
 707     int tsig = 0;\r
 708     int lastBeat = 0;\r
 709 \r
 710     vector <double> causalDF;\r
 711 \r
 712     causalDF = DF;\r
 713 \r
 714     //Prepare Causal Extension DFData\r
 715     unsigned int DFCLength = m_dataLength + m_winLength;\r
 716         \r
 717     for( unsigned int j = 0; j < m_winLength; j++ )\r
 718     {\r
 719         causalDF.push_back( 0 );\r
 720     }\r
 721         \r
 722         \r
 723     double* RW = new double[ m_lagLength ];\r
 724     for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}\r
 725 \r
 726     double* GW = new double[ m_lagLength ];\r
 727     for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}\r
 728 \r
 729     double* PW = new double[ m_lagLength ];\r
 730     for(unsigned clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}\r
 731 \r
 732     m_DFFramer.setSource( &causalDF[0], m_dataLength );\r
 733 \r
 734     unsigned int TTFrames = m_DFFramer.getMaxNoFrames();\r
 735 \r
 736 #ifdef DEBUG_TEMPO_TRACK\r
 737     std::cerr << "TTFrames = " << TTFrames << std::endl;\r
 738 #endif\r
 739         \r
 740     double* periodP = new double[ TTFrames ];\r
 741     for(unsigned clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}\r
 742         \r
 743     double* periodG = new double[ TTFrames ];\r
 744     for(unsigned clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}\r
 745         \r
 746     double* alignment = new double[ TTFrames ];\r
 747     for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}\r
 748 \r
 749     m_beats.clear();\r
 750 \r
 751     createCombFilter( RW, m_lagLength, 0, 0 );\r
 752 \r
 753     int TTLoopIndex = 0;\r
 754 \r
 755     for( unsigned int i = 0; i < TTFrames; i++ )\r
 756     {\r
 757         m_DFFramer.getFrame( m_rawDFFrame );\r
 758 \r
 759         m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame );\r
 760 \r
 761         m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength );\r
 762                 \r
 763         periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 );\r
 764 \r
 765         if( GW[ 0 ] != 0 )\r
 766         {\r
 767             periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig );\r
 768         }\r
 769         else\r
 770         {\r
 771             periodG[ TTLoopIndex ] = 0.0;\r
 772         }\r
 773 \r
 774         stepDetect( periodP, periodG, TTLoopIndex, &stepFlag );\r
 775 \r
 776         if( stepFlag == 1)\r
 777         {\r
 778             constDetect( periodP, TTLoopIndex, &constFlag );\r
 779             stepFlag = 0;\r
 780         }\r
 781         else\r
 782         {\r
 783             stepFlag -= 1;\r
 784         }\r
 785 \r
 786         if( stepFlag < 0 )\r
 787         {\r
 788             stepFlag = 0;\r
 789         }\r
 790 \r
 791         if( constFlag != 0)\r
 792         {\r
 793             tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] );\r
 794         \r
 795             createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] );\r
 796                         \r
 797             periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); \r
 798 \r
 799             period = periodG[ TTLoopIndex ];\r
 800 \r
 801 #ifdef DEBUG_TEMPO_TRACK\r
 802             std::cerr << "TempoTrack::process: constFlag == " << constFlag << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;\r
 803 #endif\r
 804 \r
 805             createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); \r
 806 \r
 807             constFlag = 0;\r
 808 \r
 809         }\r
 810         else\r
 811         {\r
 812             if( GW[ 0 ] != 0 )\r
 813             {\r
 814                 period = periodG[ TTLoopIndex ];\r
 815 \r
 816 #ifdef DEBUG_TEMPO_TRACK\r
 817                 std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodG = " << period << std::endl;\r
 818 #endif\r
 819 \r
 820                 if (period > 10000) {\r
 821                     std::cerr << "TempoTrack::process: WARNING!  Highly implausible period value " << period << "!" << std::endl;\r
 822                     std::cerr << "periodG contains (of " << TTFrames << " frames): " << std::endl;\r
 823                     for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {\r
 824                         std::cerr << i << " -> " << periodG[i] << std::endl;\r
 825                     }\r
 826                     std::cerr << "periodP contains (of " << TTFrames << " frames): " << std::endl;\r
 827                     for (int i = 0; i < TTLoopIndex + 3 && i < TTFrames; ++i) {\r
 828                         std::cerr << i << " -> " << periodP[i] << std::endl;\r
 829                     }\r
 830                     period = 5168 / 120;\r
 831                 }\r
 832 \r
 833                 createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); \r
 834 \r
 835             }\r
 836             else\r
 837             {\r
 838                 period = periodP[ TTLoopIndex ];\r
 839 \r
 840 #ifdef DEBUG_TEMPO_TRACK\r
 841                 std::cerr << "TempoTrack::process: GW[0] == " << GW[0] << ", TTLoopIndex = " << TTLoopIndex << ", period from periodP = " << period << std::endl;\r
 842 #endif\r
 843 \r
 844                 createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); \r
 845             }\r
 846         }\r
 847 \r
 848         alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); \r
 849 \r
 850         lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength );\r
 851 \r
 852         FSP += (m_lagLength);\r
 853 \r
 854         if (tempoReturn) tempoReturn->push_back(m_lockedTempo);\r
 855 \r
 856         TTLoopIndex++;\r
 857     }\r
 858 \r
 859 \r
 860     delete [] periodP;\r
 861     delete [] periodG;\r
 862     delete [] alignment;\r
 863 \r
 864     delete [] RW;\r
 865     delete [] GW;\r
 866     delete [] PW;\r
 867 \r
 868     return m_beats;\r
 869 }\r