src/tempo/beattracking.c

   1 /*
   2    Copyright (C) 2005 Matthew Davies and Paul Brossier
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include "aubio_priv.h"
  21 #include "fvec.h"
  22 #include "mathutils.h"
  23 #include "tempo/beattracking.h"
  24
  25 /** define to 1 to print out tracking difficulties */
  26 #define AUBIO_BEAT_WARNINGS 1
  27
  28 uint_t fvec_gettimesig (fvec_t * acf, uint_t acflen, uint_t gp);
  29 void aubio_beattracking_checkstate (aubio_beattracking_t * bt);
  30
  31 struct _aubio_beattracking_t
  32 {
  33   fvec_t *rwv;           /** rayleigh weighting for beat period in general model */
  34   fvec_t *dfwv;          /** exponential weighting for beat alignment in general model */
  35   fvec_t *gwv;           /** gaussian weighting for beat period in context dependant model */
  36   fvec_t *phwv;          /** gaussian weighting for beat alignment in context dependant model */
  37   fvec_t *dfrev;         /** reversed onset detection function */
  38   fvec_t *acf;           /** vector for autocorrelation function (of current detection function frame) */
  39   fvec_t *acfout;        /** store result of passing acf through s.i.c.f.b. */
  40   fvec_t *phout;
  41   uint_t timesig;        /** time signature of input, set to zero until context dependent model activated */
  42   uint_t step;
  43   uint_t rayparam;       /** Rayleigh parameter */
  44   smpl_t lastbeat;
  45   sint_t counter;
  46   uint_t flagstep;
  47   smpl_t g_var;
  48   smpl_t gp;
  49   smpl_t bp;
  50   smpl_t rp;
  51   smpl_t rp1;
  52   smpl_t rp2;
  53 };
  54
  55 aubio_beattracking_t *
  56 new_aubio_beattracking (uint_t winlen, uint_t channels)
  57 {
  58
  59   aubio_beattracking_t *p = AUBIO_NEW (aubio_beattracking_t);
  60   uint_t i = 0;
  61   /* parameter for rayleigh weight vector - sets preferred tempo to
  62    * 120bpm [43] */
  63   smpl_t rayparam = 48. / 512. * winlen;
  64   smpl_t dfwvnorm = EXP ((LOG (2.0) / rayparam) * (winlen + 2));
  65   /* length over which beat period is found [128] */
  66   uint_t laglen = winlen / 4;
  67   /* step increment - both in detection function samples -i.e. 11.6ms or
  68    * 1 onset frame [128] */
  69   uint_t step = winlen / 4;     /* 1.5 seconds */
  70
  71   p->lastbeat = 0;
  72   p->counter = 0;
  73   p->flagstep = 0;
  74   p->g_var = 3.901;             // constthresh empirically derived!
  75   p->rp = 1;
  76   p->gp = 0;
  77
  78   p->rayparam = rayparam;
  79   p->step = step;
  80   p->rwv = new_fvec (laglen, 1);
  81   p->gwv = new_fvec (laglen, 1);
  82   p->dfwv = new_fvec (winlen, 1);
  83   p->dfrev = new_fvec (winlen, channels);
  84   p->acf = new_fvec (winlen, channels);
  85   p->acfout = new_fvec (laglen, channels);
  86   p->phwv = new_fvec (2 * laglen, 1);
  87   p->phout = new_fvec (winlen, channels);
  88
  89   p->timesig = 0;
  90
  91   /* exponential weighting, dfwv = 0.5 when i =  43 */
  92   for (i = 0; i < winlen; i++) {
  93     p->dfwv->data[0][i] = (EXP ((LOG (2.0) / rayparam) * (i + 1)))
  94         / dfwvnorm;
  95   }
  96
  97   for (i = 0; i < (laglen); i++) {
  98     p->rwv->data[0][i] = ((smpl_t) (i + 1.) / SQR ((smpl_t) rayparam)) *
  99         EXP ((-SQR ((smpl_t) (i + 1.)) / (2. * SQR ((smpl_t) rayparam))));
 100   }
 101
 102   return p;
 103
 104 }
 105
 106 void
 107 del_aubio_beattracking (aubio_beattracking_t * p)
 108 {
 109   del_fvec (p->rwv);
 110   del_fvec (p->gwv);
 111   del_fvec (p->dfwv);
 112   del_fvec (p->dfrev);
 113   del_fvec (p->acf);
 114   del_fvec (p->acfout);
 115   del_fvec (p->phwv);
 116   del_fvec (p->phout);
 117   AUBIO_FREE (p);
 118 }
 119
 120
 121 void
 122 aubio_beattracking_do (aubio_beattracking_t * bt, fvec_t * dfframe,
 123     fvec_t * output)
 124 {
 125
 126   uint_t i, k;
 127   uint_t step = bt->step;
 128   uint_t laglen = bt->rwv->length;
 129   uint_t winlen = bt->dfwv->length;
 130   uint_t maxindex = 0;
 131   //number of harmonics in shift invariant comb filterbank
 132   uint_t numelem = 4;
 133
 134   smpl_t phase;                 // beat alignment (step - lastbeat)
 135   smpl_t beat;                  // beat position
 136   smpl_t bp;                    // beat period
 137   uint_t a, b;                  // used to build shift invariant comb filterbank
 138   uint_t kmax;                  // number of elements used to find beat phase
 139
 140   /* copy dfframe, apply detection function weighting, and revert */
 141   fvec_copy (dfframe, bt->dfrev);
 142   fvec_weight (bt->dfrev, bt->dfwv);
 143   fvec_rev (bt->dfrev);
 144
 145   /* compute autocorrelation function */
 146   aubio_autocorr (dfframe, bt->acf);
 147
 148   /* if timesig is unknown, use metrically unbiased version of filterbank */
 149   if (!bt->timesig) {
 150     numelem = 4;
 151   } else {
 152     numelem = bt->timesig;
 153   }
 154
 155   /* first and last output values are left intentionally as zero */
 156   fvec_zeros (bt->acfout);
 157
 158   /* compute shift invariant comb filterbank */
 159   for (i = 1; i < laglen - 1; i++) {
 160     for (a = 1; a <= numelem; a++) {
 161       for (b = (1 - a); b < a; b++) {
 162         bt->acfout->data[0][i] += bt->acf->data[0][a * (i + 1) + b - 1]
 163             * 1. / (2. * a - 1.);
 164       }
 165     }
 166   }
 167   /* apply Rayleigh weight */
 168   fvec_weight (bt->acfout, bt->rwv);
 169
 170   /* find non-zero Rayleigh period */
 171   maxindex = vec_max_elem (bt->acfout);
 172   bt->rp = maxindex ? vec_quadint (bt->acfout, maxindex, 1) : 1;
 173   //rp = (maxindex==127) ? 43 : maxindex; //rayparam
 174   bt->rp = (maxindex == bt->acfout->length - 1) ? bt->rayparam : maxindex;      //rayparam
 175
 176   /* activate biased filterbank */
 177   aubio_beattracking_checkstate (bt);
 178 #if 0                           // debug metronome mode
 179   bt->bp = 36.9142;
 180 #endif
 181   bp = bt->bp;
 182   /* end of biased filterbank */
 183
 184
 185   /* deliberate integer operation, could be set to 3 max eventually */
 186   kmax = FLOOR (winlen / bp);
 187
 188   /* initialize output */
 189   fvec_zeros (bt->phout);
 190   for (i = 0; i < bp; i++) {
 191     for (k = 0; k < kmax; k++) {
 192       bt->phout->data[0][i] += bt->dfrev->data[0][i + (uint_t) ROUND (bp * k)];
 193     }
 194   }
 195   fvec_weight (bt->phout, bt->phwv);
 196
 197   /* find Rayleigh period */
 198   maxindex = vec_max_elem (bt->phout);
 199   if (maxindex >= winlen - 1) {
 200 #if AUBIO_BEAT_WARNINGS
 201     AUBIO_WRN ("no idea what this groove's phase is\n");
 202 #endif /* AUBIO_BEAT_WARNINGS */
 203     phase = step - bt->lastbeat;
 204   } else {
 205     phase = vec_quadint (bt->phout, maxindex, 1);
 206   }
 207   /* take back one frame delay */
 208   phase += 1.;
 209 #if 0                           // debug metronome mode
 210   phase = step - bt->lastbeat;
 211 #endif
 212
 213   /* reset output */
 214   fvec_zeros (output);
 215
 216   i = 1;
 217   beat = bp - phase;
 218
 219   // AUBIO_DBG ("bp: %f, phase: %f, lastbeat: %f, step: %d, winlen: %d\n",
 220   //    bp, phase, bt->lastbeat, step, winlen);
 221
 222   /* the next beat will be earlier than 60% of the tempo period
 223     skip this one */
 224   if ( ( step - bt->lastbeat - phase ) < -0.40 * bp ) {
 225 #if AUBIO_BEAT_WARNINGS
 226     AUBIO_WRN ("back off-beat error, skipping this beat\n");
 227 #endif /* AUBIO_BEAT_WARNINGS */
 228     beat += bp;
 229   }
 230
 231   /* start counting the beats */
 232   while (beat + bp < 0) {
 233     beat += bp;
 234   }
 235
 236   if (beat >= 0) {
 237     //AUBIO_DBG ("beat: %d, %f, %f\n", i, bp, beat);
 238     output->data[0][i] = beat;
 239     i++;
 240   }
 241
 242   while (beat + bp <= step) {
 243     beat += bp;
 244     //AUBIO_DBG ("beat: %d, %f, %f\n", i, bp, beat);
 245     output->data[0][i] = beat;
 246     i++;
 247   }
 248
 249   bt->lastbeat = beat;
 250   /* store the number of beats in this frame as the first element */
 251   output->data[0][0] = i;
 252 }
 253
 254 uint_t
 255 fvec_gettimesig (fvec_t * acf, uint_t acflen, uint_t gp)
 256 {
 257   sint_t k = 0;
 258   smpl_t three_energy = 0., four_energy = 0.;
 259   if (acflen > 6 * gp + 2) {
 260     for (k = -2; k < 2; k++) {
 261       three_energy += acf->data[0][3 * gp + k];
 262       four_energy += acf->data[0][4 * gp + k];
 263     }
 264   } else {
 265     /*Expanded to be more accurate in time sig estimation */
 266     for (k = -2; k < 2; k++) {
 267       three_energy += acf->data[0][3 * gp + k] + acf->data[0][6 * gp + k];
 268       four_energy += acf->data[0][4 * gp + k] + acf->data[0][2 * gp + k];
 269     }
 270   }
 271   return (three_energy > four_energy) ? 3 : 4;
 272 }
 273
 274 void
 275 aubio_beattracking_checkstate (aubio_beattracking_t * bt)
 276 {
 277   uint_t i, j, a, b;
 278   uint_t flagconst = 0;
 279   sint_t counter = bt->counter;
 280   uint_t flagstep = bt->flagstep;
 281   smpl_t gp = bt->gp;
 282   smpl_t bp = bt->bp;
 283   smpl_t rp = bt->rp;
 284   smpl_t rp1 = bt->rp1;
 285   smpl_t rp2 = bt->rp2;
 286   uint_t laglen = bt->rwv->length;
 287   uint_t acflen = bt->acf->length;
 288   uint_t step = bt->step;
 289   fvec_t *acf = bt->acf;
 290   fvec_t *acfout = bt->acfout;
 291
 292   if (gp) {
 293     // doshiftfbank again only if context dependent model is in operation
 294     //acfout = doshiftfbank(acf,gwv,timesig,laglen,acfout);
 295     //don't need acfout now, so can reuse vector
 296     // gwv is, in first loop, definitely all zeros, but will have
 297     // proper values when context dependent model is activated
 298     fvec_zeros (acfout);
 299     for (i = 1; i < laglen - 1; i++) {
 300       for (a = 1; a <= bt->timesig; a++) {
 301         for (b = (1 - a); b < a; b++) {
 302           acfout->data[0][i] += acf->data[0][a * (i + 1) + b - 1];
 303         }
 304       }
 305     }
 306     fvec_weight (acfout, bt->gwv);
 307     gp = vec_quadint (acfout, vec_max_elem (acfout), 1);
 308     /*
 309        while(gp<32) gp =gp*2;
 310        while(gp>64) gp = gp/2;
 311      */
 312   } else {
 313     //still only using general model
 314     gp = 0;
 315   }
 316
 317   //now look for step change - i.e. a difference between gp and rp that
 318   // is greater than 2*constthresh - always true in first case, since gp = 0
 319   if (counter == 0) {
 320     if (ABS (gp - rp) > 2. * bt->g_var) {
 321       flagstep = 1;             // have observed  step change.
 322       counter = 3;              // setup 3 frame counter
 323     } else {
 324       flagstep = 0;
 325     }
 326   }
 327   //i.e. 3rd frame after flagstep initially set
 328   if (counter == 1 && flagstep == 1) {
 329     //check for consistency between previous beatperiod values
 330     if (ABS (2. * rp - rp1 - rp2) < bt->g_var) {
 331       //if true, can activate context dependent model
 332       flagconst = 1;
 333       counter = 0;              // reset counter and flagstep
 334     } else {
 335       //if not consistent, then don't flag consistency!
 336       flagconst = 0;
 337       counter = 2;              // let it look next time
 338     }
 339   } else if (counter > 0) {
 340     //if counter doesn't = 1,
 341     counter = counter - 1;
 342   }
 343
 344   rp2 = rp1;
 345   rp1 = rp;
 346
 347   if (flagconst) {
 348     /* first run of new hypothesis */
 349     gp = rp;
 350     bt->timesig = fvec_gettimesig (acf, acflen, gp);
 351     for (j = 0; j < laglen; j++)
 352       bt->gwv->data[0][j] =
 353           EXP (-.5 * SQR ((smpl_t) (j + 1. - gp)) / SQR (bt->g_var));
 354     flagconst = 0;
 355     bp = gp;
 356     /* flat phase weighting */
 357     fvec_ones (bt->phwv);
 358   } else if (bt->timesig) {
 359     /* context dependant model */
 360     bp = gp;
 361     /* gaussian phase weighting */
 362     if (step > bt->lastbeat) {
 363       for (j = 0; j < 2 * laglen; j++) {
 364         bt->phwv->data[0][j] =
 365             EXP (-.5 * SQR ((smpl_t) (1. + j - step +
 366                     bt->lastbeat)) / (bp / 8.));
 367       }
 368     } else {
 369       //AUBIO_DBG("NOT using phase weighting as step is %d and lastbeat %d \n",
 370       //                step,bt->lastbeat);
 371       fvec_ones (bt->phwv);
 372     }
 373   } else {
 374     /* initial state */
 375     bp = rp;
 376     /* flat phase weighting */
 377     fvec_ones (bt->phwv);
 378   }
 379
 380   /* do some further checks on the final bp value */
 381
 382   /* if tempo is > 206 bpm, half it */
 383   while (bp < 25) {
 384 #if AUBIO_BEAT_WARNINGS
 385     AUBIO_WRN ("doubling from %f (%f bpm) to %f (%f bpm)\n",
 386         bp, 60.*44100./512./bp, bp/2., 60.*44100./512./bp/2. );
 387     //AUBIO_DBG("warning, halving the tempo from %f\n", 60.*samplerate/hopsize/bp);
 388 #endif /* AUBIO_BEAT_WARNINGS */
 389     bp = bp * 2;
 390   }
 391
 392   //AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp);
 393
 394   /* smoothing */
 395   //bp = (uint_t) (0.8 * (smpl_t)bp + 0.2 * (smpl_t)bp2);
 396   //AUBIO_DBG("tempo:\t%3.5f bpm smoothed | bp2 %d | bp %d | ", 5168./bp, bp2, bp);
 397   //bp2 = bp;
 398   //AUBIO_DBG("time signature: %d \n", bt->timesig);
 399   bt->counter = counter;
 400   bt->flagstep = flagstep;
 401   bt->gp = gp;
 402   bt->bp = bp;
 403   bt->rp1 = rp1;
 404   bt->rp2 = rp2;
 405 }
 406
 407 smpl_t
 408 aubio_beattracking_get_bpm (aubio_beattracking_t * bt)
 409 {
 410   if (bt->timesig != 0 && bt->counter == 0 && bt->flagstep == 0) {
 411     return 5168. / vec_quadint (bt->acfout, bt->bp, 1);
 412   } else {
 413     return 0.;
 414   }
 415 }
 416
 417 smpl_t
 418 aubio_beattracking_get_confidence (aubio_beattracking_t * bt)
 419 {
 420   if (bt->gp) {
 421     return vec_max (bt->acfout);
 422   } else {
 423     return 0.;
 424   }
 425 }