annotate dsp/tempotracking/TempoTrack.cpp @ 271:61ac228cc67e

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