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

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