src/tempo/beattracking.c

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