server/plugins/BeatTrack.cpp

   1 /*
   2         SuperCollider real time audio synthesis system
   3  Copyright (c) 2002 James McCartney. All rights reserved.
   4         http://www.audiosynth.com
   5
   6  This program is free software; you can redistribute it and/or modify
   7  it under the terms of the GNU General Public License as published by
   8  the Free Software Foundation; either version 2 of the License, or
   9  (at your option) any later version.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15
  16  You should have received a copy of the GNU General Public License
  17  along with this program; if not, write to the Free Software
  18  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
  19  */
  20
  21 //BeatTrack UGen implemented by Nick Collins (http://www.informatics.sussex.ac.uk/users/nc81/)
  22 //post FFT UGen version 1 Nov 2007
  23
  24 //conversion of Matthew Davies autocorrelation beat tracking model, adapted for real-time use
  25 //currently using QMUL complex domain onset detection function model
  26
  27 //#include "SC_PlugIn.h"
  28 //#include <vecLib/vecLib.h>
  29 //#include <string.h>
  30 //#include <math.h>
  31 //#include <stdlib.h>
  32 //#include <stdio.h>
  33
  34 #include "ML.h"
  35
  36 //FFT data
  37 //#define N 1024  //FFT size
  38 //FFT size over 2
  39 #define NOVER2 512
  40 //#define NOVER4 256  //FFT size
  41 //#define OVERLAP 512
  42 //#define OVERLAPINDEX 512
  43 //#define HOPSIZE 512
  44
  45 //#define FS 44100 //assumes fixed sampling rate
  46 //#define FRAMESR 86.1328
  47 //converted for different sampling rates
  48 #define FRAMEPERIOD 0.01161
  49 #define SKIP 128
  50 //#define TIMEELAPSED 1.48608
  51
  52 //this data assumes LAGS is 128
  53 static float g_m[128]= {0.00054069,0.00108050,0.00161855,0.00215399,0.00268594,0.00321356,0.00373600,0.00425243,0.00476204,0.00526404,0.00575765,0.00624213,0.00671675,0.00718080,0.00763362,0.00807455,0.00850299,0.00891836,0.00932010,0.00970771,0.01008071,0.01043866,0.01078115,0.01110782,0.01141834,0.01171242,0.01198982,0.01225033,0.01249378,0.01272003,0.01292899,0.01312061,0.01329488,0.01345182,0.01359148,0.01371396,0.01381939,0.01390794,0.01397980,0.01403520,0.01407439,0.01409768,0.01410536,0.01409780,0.01407534,0.01403838,0.01398734,0.01392264,0.01384474,0.01375410,0.01365120,0.01353654,0.01341062,0.01327397,0.01312710,0.01297054,0.01280484,0.01263053,0.01244816,0.01225827,0.01206139,0.01185807,0.01164884,0.01143424,0.01121478,0.01099099,0.01076337,0.01053241,0.01029861,0.01006244,0.00982437,0.00958484,0.00934429,0.00910314,0.00886181,0.00862067,0.00838011,0.00814049,0.00790214,0.00766540,0.00743057,0.00719793,0.00696778,0.00674036,0.00651591,0.00629466,0.00607682,0.00586256,0.00565208,0.00544551,0.00524301,0.00504470,0.00485070,0.00466109,0.00447597,0.00429540,0.00411944,0.00394813,0.00378151,0.00361959,0.00346238,0.00330989,0.00316210,0.00301899,0.00288053,0.00274669,0.00261741,0.00249266,0.00237236,0.00225646,0.00214488,0.00203755,0.00193440,0.00183532,0.00174025,0.00164909,0.00156174,0.00147811,0.00139810,0.00132161,0.00124854,0.00117880,0.00111228,0.00104887,0.00098848,0.00093100,0.00087634,0.00082438,};
  54 static float g_mg[257]= {0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000004,0.00000055,0.00000627,0.00005539,0.00037863,0.00200318,0.00820201,0.02599027,0.06373712,0.12096648,0.17767593,0.20196826,0.17767593,0.12096648,0.06373712,0.02599027,0.00820201,0.00200318,0.00037863,0.00005539,0.00000627,0.00000055,0.00000004,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,0.00000000,};
  55
  56 //other functions
  57 static void BeatTrack_dofft(BeatTrack *unit, uint32);
  58 static void complexdf(BeatTrack *unit);
  59 static void finaldecision(BeatTrack *unit);
  60
  61 //amortisation
  62 static void autocorr(BeatTrack *unit,int j);
  63 static void beatperiod(BeatTrack *unit,int j, int whichm);
  64 static float findtor(BeatTrack *unit);
  65
  66 //as many amortisation steps as tor
  67 static void findphase(BeatTrack *unit,int j,int gaussflag,int predicted);
  68
  69 static int detectperiodchange(BeatTrack *unit);
  70 static void findmeter(BeatTrack *unit);
  71 static void setupphaseexpectation(BeatTrack *unit); //create Gaussian focussed matrix for phase
  72
  73
  74 void BeatTrack_Ctor(BeatTrack* unit)
  75 {
  76         ///////
  77         //check sampling rate and establish multipliers on estimates and FFT window size
  78         //down sampling by factor of two automatic
  79
  80         unit->m_srate = unit->mWorld->mFullRate.mSampleRate;
  81
  82         //if sample rate is 88200 or 96000, assume taking double size FFT to start with
  83         if(unit->m_srate > (44100.0*1.5)) unit->m_srate = unit->m_srate*0.5;
  84
  85         unit->m_srateconversion = unit->m_srate/44100.0;
  86
  87         //assumes base of 1024 FFT
  88         unit->m_frameperiod= (FRAMEPERIOD/unit->m_srateconversion); //in seconds //(int) ((FRAMEPERIOD/unit->m_srateconversion) +0.5);
  89
  90         printf("srate %f conversion factor %f frame period %f \n", unit->m_srate, unit->m_srateconversion, unit->m_frameperiod);
  91
  92         unit->m_prevmag= (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
  93         unit->m_prevphase= (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
  94         unit->m_predict= (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
  95
  96         ////////time positions//////////
  97         unit->m_frame=1; //don't decide immediately, wait for maximum period!
  98
  99         /////////df////////
 100         unit->m_dfcounter=DFSTORE-1;
 101         //random uncorrelated noise df store for initialisation
 102         //RGen& rgen = *unit->mParent->mRGen;
 103
 104         //don't want this noise, want consistent starting point!
 105         for(int j=0;j<DFSTORE;++j) {
 106                 unit->m_df[j]=0.0; //(2*rgen.frand() - 1.0);
 107         }
 108
 109         unit->m_dfmemorycounter=14;
 110
 111         Clear(15, unit->m_dfmemory);
 112
 113         /////////tempo assess///////////
 114         unit->m_currtempo=2;
 115
 116         ////////phase assess///////////
 117
 118         unit->m_currphase=0.0;
 119
 120         unit->m_phase=0.0;
 121
 122         //default of 2bps
 123         unit->m_phaseperblock= ((float)unit->mWorld->mFullRate.mBufLength*2)/((float)unit->mWorld->mSampleRate);
 124
 125         unit->m_outputphase= unit->m_phase;
 126         unit->m_outputtempo= unit->m_currtempo;
 127         unit->m_outputphaseperblock= unit->m_phaseperblock;
 128
 129         unit->halftrig=0;
 130     unit->q1trig=0;
 131         unit->q2trig=0;
 132
 133         //amortisation and states
 134         unit->m_amortisationstate=0; //off
 135         unit->m_stateflag=0;
 136         unit->m_timesig=4;
 137         unit->m_flagstep=0;
 138
 139         unit->mCalcFunc = (UnitCalcFunc)&BeatTrack_next;
 140 }
 141
 142
 143
 144 void BeatTrack_Dtor(BeatTrack *unit)
 145 {
 146         RTFree(unit->mWorld, unit->m_prevmag);
 147         RTFree(unit->mWorld, unit->m_prevphase);
 148         RTFree(unit->mWorld, unit->m_predict);
 149 }
 150
 151
 152 void BeatTrack_next(BeatTrack *unit, int wrongNumSamples)
 153 {
 154         //float *in = IN(0);
 155
 156         //printf("%d \n",wrongNumSamples);
 157         //int numSamples = unit->mWorld->mFullRate.mBufLength;
 158
 159         //conditions in reverse order to avoid immediate spillover
 160         //printf("state %d \n",unit->m_amortisationstate);
 161
 162         //keeps incrementing but will be reset with each calculation run
 163         unit->m_amortisationsteps=unit->m_amortisationsteps+1;
 164
 165         //if state nonzero do something
 166         switch(unit->m_amortisationstate) {
 167                 case 0:
 168                         break; //do nothing case
 169                 case 1: //calculate acf
 170                         autocorr(unit,unit->m_amortcount);
 171
 172                         unit->m_amortcount=unit->m_amortcount+1;
 173
 174                         if(unit->m_amortcount==unit->m_amortlength) {
 175                                 unit->m_amortisationstate=2;
 176                                 unit->m_amortlength=128;
 177                                 unit->m_amortcount=0;
 178
 179                                 unit->m_bestcolumn=0;
 180                                 unit->m_besttorsum= -1000.0;
 181
 182                         }
 183
 184                                 break;
 185                 case 2: //periodp
 186                         beatperiod(unit,unit->m_amortcount,0);
 187
 188                         unit->m_amortcount=unit->m_amortcount+1;
 189
 190                         if(unit->m_amortcount==unit->m_amortlength) {
 191
 192                                 unit->m_periodp=findtor(unit);
 193
 194                                 if(unit->m_stateflag==1) {
 195                                         unit->m_amortisationstate=3;
 196                                         unit->m_amortlength=128;
 197                                         unit->m_amortcount=0;
 198
 199                                         unit->m_bestcolumn=0;
 200                                         unit->m_besttorsum= -1000.0;
 201
 202                                 } else {
 203                                         unit->m_periodg= -1000; //will always trigger initially
 204                                         unit->m_amortisationstate=4;
 205                                 }
 206                         }
 207
 208                                 break;
 209                 case 3: //periodg
 210                         beatperiod(unit,unit->m_amortcount,1);
 211                         unit->m_amortcount=unit->m_amortcount+1;
 212
 213                         if(unit->m_amortcount==unit->m_amortlength) {
 214
 215                                 unit->m_periodg=findtor(unit);
 216
 217                                 unit->m_amortisationstate=4;
 218                         }
 219
 220                                 break;
 221                 case 4: //stepdetect/constdetect
 222
 223                         if(detectperiodchange(unit)) {
 224
 225                                 unit->m_amortisationstate=5;
 226                                 unit->m_amortlength=128;
 227                                 unit->m_amortcount=0;
 228
 229                                 unit->m_bestcolumn=0;
 230                                 unit->m_besttorsum= -1000.0;
 231
 232                                 unit->m_stateflag=1;
 233                                 findmeter(unit);
 234
 235                                 //set up Gaussian weighting centred on periodp
 236                                 int startindex= 128- ((int)(unit->m_periodp+0.5));
 237
 238                                 float * mg=unit->m_mg;
 239
 240                                 for (int ii=0; ii<128;++ii){
 241                                         mg[ii]= g_mg[startindex+ii];
 242                                 }
 243
 244                         } else {
 245
 246                                 if(unit->m_stateflag==1)
 247                                         unit->m_tor= unit->m_periodg;
 248                                 else
 249                                         unit->m_tor= unit->m_periodp;
 250
 251                                 unit->m_torround= int(unit->m_tor+0.5);
 252
 253                                 unit->m_amortisationstate=7;
 254                                 unit->m_amortlength=unit->m_torround;
 255                                 unit->m_amortcount=0;
 256                         }
 257
 258                         break;
 259
 260                 case 5: //redo periodg calculation
 261                         beatperiod(unit,unit->m_amortcount,1);
 262                         unit->m_amortcount=unit->m_amortcount+1;
 263
 264                         if(unit->m_amortcount==unit->m_amortlength) {
 265
 266                                 unit->m_periodg=findtor(unit);
 267
 268                                 unit->m_tor= unit->m_periodg;
 269                                 unit->m_torround= int(unit->m_tor+0.5f);
 270
 271                                 unit->m_amortisationstate=6;
 272                                 unit->m_amortlength=unit->m_torround;
 273                                 unit->m_amortcount=0;
 274
 275                                 setupphaseexpectation(unit);
 276
 277                                 //don't need to reset change flag since it isn't stored
 278                         }
 279
 280                                 break;
 281                 case 6: //flat phase calc after move to context, avoids bias
 282                         findphase(unit,unit->m_amortcount,0,0);
 283                         unit->m_amortcount=unit->m_amortcount+1;
 284
 285                         if(unit->m_amortcount==unit->m_amortlength) {
 286
 287                                 unit->m_amortisationstate=8; //final state
 288                         }
 289
 290                                 break;
 291
 292                 case 7: //phase calc with possible gaussian narrowing of the allowed phases
 293
 294                         findphase(unit,unit->m_amortcount,unit->m_stateflag,(int)(unit->m_currphase*unit->m_torround+0.5f));
 295                         unit->m_amortcount=unit->m_amortcount+1;
 296
 297                         if(unit->m_amortcount==unit->m_amortlength) {
 298
 299                                 unit->m_amortisationstate=8; //final state
 300                         }
 301
 302                                 break;
 303                 case 8:
 304
 305                         finaldecision(unit);
 306                         unit->m_amortisationstate=0;
 307                         break;
 308
 309                 default:
 310                         break;
 311         }
 312
 313
 314         //MUST CHECK IF INCIDENT FFT IS >1, if so update buffer with appropriate coefficients
 315
 316         float fbufnum = ZIN0(0);
 317
 318         //next FFT bufffer ready, update
 319         //assuming at this point that buffer precalculated for any resampling
 320         if (!(fbufnum<0)) {
 321
 322                 unit->m_frame= unit->m_frame+1;
 323                 BeatTrack_dofft(unit, (uint32)fbufnum);
 324
 325         }
 326
 327
 328         //test if impulse to output
 329         unit->m_phase+=unit->m_phaseperblock;
 330
 331         //if not locked, update output phase from model phase, else keep a separate output phase
 332
 333         float lock= ZIN0(1);
 334         //printf("lock %f \n",lock);
 335
 336         if (lock<0.5f) {
 337                 unit->m_outputphase= unit->m_phase;
 338                 unit->m_outputtempo= unit->m_currtempo;
 339                 unit->m_outputphaseperblock= unit->m_phaseperblock;
 340         } else
 341                 unit->m_outputphase+=unit->m_outputphaseperblock;
 342
 343         if (unit->m_phase >= 1.f) unit->m_phase-= 1.f;
 344
 345         //0 is beat, 1 is quaver, 2 is semiquaver, 3 is actual current tempo in bps
 346         //so no audio accuracy with beats, just asap, may as well be control rate
 347         ZOUT0(0)=0.0;
 348         ZOUT0(1)=0.0;
 349         ZOUT0(2)=0.0;
 350         ZOUT0(3)=unit->m_outputtempo; //*0.016666667;
 351
 352         //output beat
 353         if (unit->m_outputphase >= 1.f) {
 354
 355                 //printf("beat \n");
 356
 357                 unit->m_outputphase -= 1.f;
 358                 ZOUT0(0)=1.0;
 359                 ZOUT0(1)=1.0;
 360                 ZOUT0(2)=1.0;
 361                 unit->halftrig=0;
 362                 unit->q1trig=0;
 363                 unit->q2trig=0;
 364         }
 365
 366         if (unit->m_outputphase>=0.5f && unit->halftrig==0) {
 367                 ZOUT0(1)=1.0;
 368                 ZOUT0(2)=1.0;
 369                 unit->halftrig=1;
 370         }
 371
 372         if (unit->m_outputphase>=0.25f && unit->q1trig==0) {
 373                 ZOUT0(2)=1.0;
 374                 unit->q1trig=1;
 375         }
 376
 377         if (unit->m_outputphase>=0.75f && unit->q2trig==0) {
 378                 ZOUT0(2)=1.0;
 379                 unit->q2trig=1;
 380         }
 381
 382 }
 383
 384
 385
 386 //
 387
 388 //calculation function once FFT data ready
 389 void BeatTrack_dofft(BeatTrack *unit, uint32 ibufnum)
 390 {
 391         World *world = unit->mWorld;
 392         SndBuf *buf;
 393         if (ibufnum >= world->mNumSndBufs) {
 394                 int localBufNum = ibufnum - world->mNumSndBufs;
 395                 Graph *parent = unit->mParent;
 396                 if(localBufNum <= parent->localBufNum) {
 397                         buf = parent->mLocalSndBufs + localBufNum;
 398                 } else {
 399                         buf = world->mSndBufs;
 400                 }
 401         } else {
 402                 buf = world->mSndBufs + ibufnum;
 403         }
 404         LOCK_SNDBUF(buf);
 405         //int numbins = buf->samples - 2 >> 1;
 406
 407         unit->m_FFTBuf = buf->data; //just assign it!
 408
 409         //transfer data to fftbuf in the format expected by this plugin
 410
 411         //ideally, should do this part separate to plug-in as well, so can compare different detection functions;
 412         //also, can run multiple in parallel with own autocorrelations; committee? Committee.ar(period1, phase1, period2, phase2, period3, phase3)...
 413         //chooses predominant estimate?
 414         //feature detection function
 415         complexdf(unit);
 416
 417         if (unit->m_frame%SKIP==0) {
 418
 419                 //printf("amortisation time \n");
 420
 421                 //amortisation- 8 control periods in a frame
 422                 //have 2000 calcs to do, split over 100 control periods = 6400 samples, ie one tempo per control period
 423
 424                 unit->m_bestcolumn=0;
 425                 unit->m_besttorsum= -1000.0;
 426
 427                 unit->m_bestphasescore = -1000.0;
 428                 unit->m_bestphase = 0;
 429
 430                 //state 0 is do nothing
 431                 unit->m_amortisationstate=1;
 432                 unit->m_amortcount=0;
 433                 unit->m_amortlength=128;
 434                 unit->m_amortisationsteps=0;
 435
 436                 //fix time reference for calculations, so it doesn't update during the amortisation- this is the beginning of the df frame
 437                 unit->m_storedfcounter= unit->m_dfcounter+DFSTORE-DFFRAMELENGTH;
 438
 439                 //ref for phase calculations
 440                 unit->m_storedfcounterend= unit->m_dfcounter;
 441
 442                 //unit->m_fftstoreposhold= unit->m_fftstorepos;
 443
 444                 unit->m_currphase=unit->m_phase;
 445
 446         }
 447
 448 }
 449
 450
 451 void autocorr(BeatTrack *unit,int j)
 452 {
 453         int baseframe=unit->m_storedfcounter+DFSTORE;
 454         float * df= unit->m_df;
 455         float * acf= unit->m_acf;
 456
 457         //work out four lags each time
 458         for (int k=0;k<4;++k) {
 459
 460                 int lag=4*j+k;
 461
 462                 int correction= abs(lag-DFFRAMELENGTH);
 463
 464                 float sum=0.0;
 465
 466                 for (int i=lag;i<DFFRAMELENGTH; ++i) {
 467
 468                         float val1= df[(i+baseframe)%DFSTORE];
 469                         float val2= df[(i+baseframe-lag)%DFSTORE];
 470
 471                         sum+= val1*val2;
 472                 }
 473
 474                 acf[lag]=sum*correction;
 475
 476         }
 477
 478 }
 479
 480
 481
 482 //timesig 4 has one more sum term
 483 //indices as MATLAB but need to correct maxinds to be in range of tested, not in global range
 484 float findtor(BeatTrack *unit)
 485 {
 486         float maxval, val;
 487         int ind2,ind3,ind4;
 488
 489         //put into MATLAB indexing, from 1 to 512
 490         int ind= unit->m_bestcolumn+1;
 491
 492         float * acf= unit->m_acf-1;
 493
 494         ind2=0;
 495         maxval=-1000;
 496
 497         for(int i=2*ind-1;i<=(2*ind+1);++i){
 498
 499                 val=acf[i];
 500
 501                 if(val>maxval) {
 502
 503                         maxval=val;
 504                         ind2=i-(2*ind-1)+1;
 505                 }
 506
 507         }
 508
 509         //[val2,ind2] = max(acf(2*ind-1:2*ind+1));
 510         ind2 = ind2 + 2*(ind+1)-2;
 511
 512         ind3=0;
 513         maxval=-1000;
 514
 515         for(int i=3*ind-2;i<=(3*ind+2);++i){
 516
 517                 val=acf[i];
 518
 519                 if(val>maxval) {
 520
 521                         maxval=val;
 522                         ind3=i-(3*ind-2)+1;
 523                 }
 524
 525         }
 526
 527         //[val3,ind3] = max(acf(3*ind-2:3*ind+2));
 528         ind3 = ind3 + 3*ind-4;
 529
 530         float period;
 531
 532         if (unit->m_timesig==4) {
 533
 534                 ind4=0;
 535                 maxval=-1000;
 536
 537                 for(int i=4*ind-3;i<=4*ind+3;++i){
 538
 539                         val=acf[i];
 540
 541                         if(val>maxval) {
 542
 543                                 maxval=val;
 544                                 ind4=i-(4*ind-3)+1;
 545                         }
 546
 547                 }
 548
 549                 //[val4,ind4] = max(acf(4*ind-3:4*ind+3));
 550                 ind4 = ind4 + 4*ind-9;
 551
 552                 period= (ind+ ind2*0.5+ind3/3.f +ind4*0.25)*0.25;
 553
 554         } else
 555
 556                 period= (ind+ ind2*0.5+ind3/3.f)*0.3333333;
 557
 558
 559         //printf("period %f ind %d ind2 %d ind3 %d ind4 %d \n",period, ind,ind2,ind3,ind4);
 560
 561         //unit->m_tor=period;
 562         //unit->m_torround= int(period+0.5);
 563         //
 564
 565         return period;
 566 }
 567
 568
 569
 570
 571 //128 calculation calls for multiplying M and acf, calculates M as it goes apart from precalculated Gaussian or Raleigh distribution
 572 void beatperiod(BeatTrack *unit,int j, int whichm)
 573 {
 574         int baseframe = unit->m_storedfcounter+DFSTORE;
 575         float * acf = unit->m_acf;
 576
 577         //int startindex= 512*j;
 578         //int endindex=startindex+512;
 579
 580         float sum=0.0;
 581
 582         //unit->m_timesig harmonics
 583         for (int i=1;i<=(unit->m_timesig); ++i) {
 584
 585                 int num = 2*i-1;
 586                 float wt= 1.0/(float)num;
 587
 588                 for (int k=0;k<num; ++k) {
 589
 590                         int pos= k+(i*j);
 591
 592                         //m[startindex+pos]
 593                         if(pos<512)
 594                                 sum+= acf[pos]*wt;
 595                 }
 596
 597         }
 598
 599         //assumes Mg appropriately rotated already
 600         float * m;
 601
 602         if(whichm)
 603                 m=g_m; //Gaussian weighted context model
 604         else
 605                 m=unit->m_mg; //general model even weighting
 606
 607         sum=sum*m[j];
 608
 609         if (sum>unit->m_besttorsum) {
 610                 unit->m_besttorsum=sum;
 611                 unit->m_bestcolumn=j;
 612         }
 613 }
 614
 615
 616 //j out of unit->m_torround
 617 //differs to Davies original in that weight the most recent events more- want minimum reaction time
 618 void findphase(BeatTrack *unit,int j,int gaussflag, int predicted)
 619 {
 620         float * df= unit->m_df;
 621
 622         int period= unit->m_torround;
 623         int baseframe=unit->m_storedfcounterend+DFSTORE;
 624
 625         int numfit= -1;
 626
 627         if(period != 0)
 628         //round down
 629         numfit= (int)(DFFRAMELENGTH/period)-1;
 630
 631         //testing backwards from the baseframe, weighting goes down as 1/k
 632         float sum=0.0;
 633
 634         for (int k=0;k<numfit;++k) {
 635
 636                 //j is phase to test
 637                 int location= (baseframe-(period*k)-j)%DFSTORE;
 638
 639                 sum+= df[location]/((float)(k+1));
 640
 641         }
 642
 643         //Gaussian focus weighting if desired
 644         if (gaussflag) {
 645
 646                 //difference of predicted from j, min distance within period
 647                 int diff= sc_min(abs(predicted-j),abs(period-predicted+j));
 648
 649                 sum *= unit->m_phaseweights[diff];
 650
 651         }
 652
 653         if (sum>unit->m_bestphasescore) {
 654
 655                 unit->m_bestphasescore = sum;
 656                 unit->m_bestphase = j;
 657
 658         }
 659
 660 }
 661
 662 //, int predicted
 663 void setupphaseexpectation(BeatTrack *unit)  //create Gaussian focussed matrix for phase
 664 {
 665         float * wts= unit->m_phaseweights;
 666
 667         float sigma= unit->m_torround * 0.25f;
 668         //float mu=period;
 669
 670         float mult= 1.0/(2.5066283*sigma);
 671         float mult2= 1.0/(2.0*sigma*sigma);
 672
 673         //unit->m_torround
 674         for (int i=0; i<128;++i) {
 675                 wts[i]= mult*(exp(-(i*i)*mult2));
 676         }
 677
 678 }
 679
 680
 681 //why force a countdown each time? Why not keep a continuous buffer of previous periodp, periodg?
 682 int detectperiodchange(BeatTrack *unit)
 683 {
 684         //stepthresh = 3.9017;
 685
 686         if(unit->m_flagstep==0) {
 687
 688                 if(fabs(unit->m_periodg-unit->m_periodp) > 3.9017f) {
 689                         unit->m_flagstep= 3;
 690                 }
 691
 692         } else {
 693                 unit->m_flagstep= unit->m_flagstep-1;
 694         }
 695
 696         if(unit->m_flagstep) {
 697
 698                 unit->m_prevperiodp[unit->m_flagstep-1]=unit->m_periodp;
 699         }
 700
 701         if(unit->m_flagstep==1) {
 702
 703                 unit->m_flagstep= 0;
 704
 705                 if(fabs(2*unit->m_prevperiodp[0] - unit->m_prevperiodp[1] - unit->m_prevperiodp[2]) < 7.8034f) //(2*3.9017)
 706                         return 1;
 707
 708         }
 709
 710         return 0;
 711
 712 }
 713
 714 //add test
 715 void findmeter(BeatTrack *unit)
 716 {
 717
 718         //int i;
 719
 720         //float * acf= unit->m_acf;
 721
 722 //      float * acf= unit->m_acf-1;
 723 //
 724 //
 725 //      int period = ((int)(unit->m_periodp+0.5));
 726 //
 727 //      float three_energy=0.0;
 728 //      float four_energy=0.0;
 729 //
 730 //      for(i=(3*period-2);i<(3*period+3);++i)
 731 //              three_energy += acf[i];
 732 //
 733 //      for(i=(4*period-2);i<(4*period+3);++i)
 734 //              four_energy += acf[i];
 735 //
 736 //      if((6*period+2)<512) {
 737 //
 738 //              for(i=(6*period-2);i<(6*period+3);++i)
 739 //                      three_energy += acf[i];
 740 //
 741 //              for(i=(2*period-2);i<(2*period+3);++i)
 742 //                      four_energy += acf[i];
 743 //      }
 744 //
 745 //      if (three_energy > four_energy)
 746 //              unit->m_timesig = 3;
 747 //      else
 748
 749 //worked better in evaluation without any 3/4 at all!
 750                 unit->m_timesig = 4;
 751
 752         //printf("time sig %d \n",unit->m_timesig);
 753
 754 }
 755
 756
 757 //period is unit->m_tor, phase is unit->m_bestphase
 758 //      float m_tor; int m_torround;
 759 void finaldecision(BeatTrack *unit)
 760 {
 761         //int i,j;
 762
 763         unit->m_currtempo= 1.0/(unit->m_tor*unit->m_frameperiod);
 764
 765         unit->m_phaseperblock= ((float)unit->mWorld->mFullRate.mBufLength*(unit->m_currtempo))/((float)unit->mWorld->mSampleRate);
 766
 767         //printf("SAMPLErate %f %f %f", unit->mWorld->mSampleRate,unit->m_phaseperblock,unit->m_currtempo);
 768
 769         //unit->m_amortisationstate control periods worth = 512/64 = 8
 770         //float frameselapsed= 0.125*unit->m_amortisationstate;
 771         //float timeelapsed= frameselapsed*unit->m_frameperiod;
 772
 773         float timeelapsed= ((float)(unit->m_amortisationsteps)*(unit->mWorld->mFullRate.mBufLength)/((float)unit->mWorld->mSampleRate));
 774
 775         timeelapsed += 7*unit->m_frameperiod; //compensation for detection function being delayed by 7 frames
 776
 777         float phaseelapsed= timeelapsed*(unit->m_currtempo);
 778
 779         float phasebeforeamort= ((float)unit->m_bestphase/unit->m_torround);
 780
 781         //add phase to compensate for ELAPSEDTIME
 782         unit->m_currphase= unit->m_phase = fmod(phasebeforeamort+phaseelapsed,(float)1.0);
 783
 784 }
 785
 786
 787 //Now the format is standardised for the SC FFT UGen as
 788 //dc, nyquist and then real/imag pairs for each bin going up successively in frequency
 789
 790 void complexdf(BeatTrack *unit)
 791 {
 792         float * fftbuf= unit->m_FFTBuf;
 793
 794         float * prevmag= unit->m_prevmag;
 795         float * prevphase= unit->m_prevphase;
 796         float * predict= unit->m_predict;
 797
 798         float sum=0.0;
 799
 800         //printf("complex df time \n");
 801
 802         //sum bins 2 to 256
 803         for (int k=1; k<NOVER2; ++k){
 804
 805                         //Change to fftw
 806                         int index= 2*k; //k; //2*k;
 807
 808                         float real=fftbuf[index];
 809                         //N=1024 conventionally here
 810                         float imag=fftbuf[index+1]; //fftbuf[N-index];
 811
 812                         float mag= sqrt(real*real+ imag*imag); // was  0.5*sqrt(real*real+ imag*imag); reduce by factor of 2 because of altivec side effect
 813                         float qmag= prevmag[k];
 814
 815                         prevmag[k]=mag;
 816
 817                         float phase= atan2(imag,real);
 818                         float oldphase = predict[k];
 819
 820                         predict[k]= 2*phase- prevphase[k];
 821                         prevphase[k]= phase;
 822
 823                         float phasediff= phase-oldphase;
 824
 825                         //if(k==2) printf("%f %f\n",phase, phasediff);
 826
 827                         //tables for cos/sin/sqrt speeds up? sqrt(1-c*c) slower than sin
 828
 829                         float realpart= (qmag-(mag*cos(phasediff)));
 830                         float imagpart= (mag*sin(phasediff)); //no need for negative
 831
 832                         float detect= sqrt(realpart*realpart + imagpart*imagpart);
 833
 834                         //detect is always positive
 835                         //if(k==1)
 836                         sum+=detect; //(fmod(phase+(16*pi),twopi)); //detect;
 837
 838                         //if(k==1) sum+=mag;
 839         }
 840
 841
 842         //smoothing and peak picking operation, delay of 8 frames, must be taken account of in final phase correction
 843
 844         unit->m_dfmemorycounter=(unit->m_dfmemorycounter+1)%15;
 845         unit->m_dfmemory[unit->m_dfmemorycounter]=sum; //divide by num of bands to get a dB answer
 846
 847         float rating=0.0;
 848
 849         float * dfmemory=unit->m_dfmemory;
 850
 851         int refpos=unit->m_dfmemorycounter+15;
 852         int centrepos=(refpos-7)%15;
 853         float centreval=dfmemory[centrepos];
 854
 855         for (int k=0;k<15; ++k) {
 856
 857                 int pos=(refpos-k)%15;
 858
 859                 float nextval= centreval-dfmemory[pos];
 860
 861                 if (nextval<0.0)
 862                         nextval=nextval*10;
 863
 864                 rating+=nextval;
 865         }
 866
 867         if(rating<0.0) rating=0.0;
 868
 869         //increment first so this frame is unit->m_loudnesscounterdfcounter
 870         unit->m_dfcounter=(unit->m_dfcounter+1)%DFSTORE;
 871
 872         unit->m_df[unit->m_dfcounter]=rating*0.1f; //sum //divide by num of bands to get a dB answer
 873
 874 }