source: src/tempo/beattracking.c @ acd97d1

feature/autosinkfeature/cnnfeature/cnn_orgfeature/constantqfeature/crepefeature/crepe_orgfeature/pitchshiftfeature/pydocstringsfeature/timestretchfix/ffmpeg5pitchshiftsamplertimestretchyinfft+
Last change on this file since acd97d1 was acd97d1, checked in by Paul Brossier <piem@piem.org>, 10 years ago

src/mathutils.{c,h}: remove fvec_quadint, use fvec_quadratic_peak_pos

  • Property mode set to 100644
File size: 12.1 KB
RevLine 
[b78805a]1/*
[e6a78ea]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
[b78805a]19*/
20
21#include "aubio_priv.h"
[6c7d49b]22#include "fvec.h"
[b78805a]23#include "mathutils.h"
[32d6958]24#include "tempo/beattracking.h"
[b78805a]25
[17b7d66]26/** define to 1 to print out tracking difficulties */
[1812f49]27#define AUBIO_BEAT_WARNINGS 0
[17b7d66]28
[7bf3dcb]29uint_t fvec_gettimesig (fvec_t * acf, uint_t acflen, uint_t gp);
30void aubio_beattracking_checkstate (aubio_beattracking_t * bt);
31
32struct _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;
[b78805a]54};
55
[7bf3dcb]56aubio_beattracking_t *
[d207300]57new_aubio_beattracking (uint_t winlen)
[7bf3dcb]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;
[d207300]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);
[7bf3dcb]89
90  p->timesig = 0;
91
92  /* exponential weighting, dfwv = 0.5 when i =  43 */
93  for (i = 0; i < winlen; i++) {
[d207300]94    p->dfwv->data[i] = (EXP ((LOG (2.0) / rayparam) * (i + 1)))
[7bf3dcb]95        / dfwvnorm;
96  }
97
98  for (i = 0; i < (laglen); i++) {
[d207300]99    p->rwv->data[i] = ((smpl_t) (i + 1.) / SQR ((smpl_t) rayparam)) *
[7bf3dcb]100        EXP ((-SQR ((smpl_t) (i + 1.)) / (2. * SQR ((smpl_t) rayparam))));
101  }
102
103  return p;
[b78805a]104
105}
106
[7bf3dcb]107void
108del_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);
[b78805a]119}
120
121
[7bf3dcb]122void
123aubio_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
[72db1cf]138  uint_t a, b;                  // used to build shift invariant comb filterbank
[7bf3dcb]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++) {
[72db1cf]162      for (b = 1; b < 2 * a; b++) {
163        bt->acfout->data[i] += bt->acf->data[i * a + b - 1]
[7bf3dcb]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 */
[1e2c82f]172  maxindex = fvec_max_elem (bt->acfout);
[acd97d1]173  bt->rp = maxindex ? fvec_quadratic_peak_pos (bt->acfout, maxindex) : 1;
[7bf3dcb]174  //rp = (maxindex==127) ? 43 : maxindex; //rayparam
175  bt->rp = (maxindex == bt->acfout->length - 1) ? bt->rayparam : maxindex;      //rayparam
176
177  /* activate biased filterbank */
178  aubio_beattracking_checkstate (bt);
179#if 0                           // debug metronome mode
180  bt->bp = 36.9142;
[4e19e5b]181#endif
[7bf3dcb]182  bp = bt->bp;
183  /* end of biased filterbank */
184
185
186  /* deliberate integer operation, could be set to 3 max eventually */
187  kmax = FLOOR (winlen / bp);
188
189  /* initialize output */
190  fvec_zeros (bt->phout);
191  for (i = 0; i < bp; i++) {
192    for (k = 0; k < kmax; k++) {
[d207300]193      bt->phout->data[i] += bt->dfrev->data[i + (uint_t) ROUND (bp * k)];
[7bf3dcb]194    }
195  }
196  fvec_weight (bt->phout, bt->phwv);
197
198  /* find Rayleigh period */
[1e2c82f]199  maxindex = fvec_max_elem (bt->phout);
[3de10bb]200  if (maxindex >= winlen - 1) {
[17b7d66]201#if AUBIO_BEAT_WARNINGS
202    AUBIO_WRN ("no idea what this groove's phase is\n");
203#endif /* AUBIO_BEAT_WARNINGS */
[3de10bb]204    phase = step - bt->lastbeat;
205  } else {
[acd97d1]206    phase = fvec_quadratic_peak_pos (bt->phout, maxindex);
[3de10bb]207  }
[17b7d66]208  /* take back one frame delay */
209  phase += 1.;
[7bf3dcb]210#if 0                           // debug metronome mode
211  phase = step - bt->lastbeat;
[4e19e5b]212#endif
[b78805a]213
[7bf3dcb]214  /* reset output */
215  fvec_zeros (output);
216
217  i = 1;
218  beat = bp - phase;
[17b7d66]219
220  // AUBIO_DBG ("bp: %f, phase: %f, lastbeat: %f, step: %d, winlen: %d\n",
221  //    bp, phase, bt->lastbeat, step, winlen);
222
223  /* the next beat will be earlier than 60% of the tempo period
224    skip this one */
225  if ( ( step - bt->lastbeat - phase ) < -0.40 * bp ) {
226#if AUBIO_BEAT_WARNINGS
227    AUBIO_WRN ("back off-beat error, skipping this beat\n");
228#endif /* AUBIO_BEAT_WARNINGS */
229    beat += bp;
230  }
231
[7bf3dcb]232  /* start counting the beats */
[3de10bb]233  while (beat + bp < 0) {
234    beat += bp;
235  }
236
[7bf3dcb]237  if (beat >= 0) {
[3de10bb]238    //AUBIO_DBG ("beat: %d, %f, %f\n", i, bp, beat);
[d207300]239    output->data[i] = beat;
[7bf3dcb]240    i++;
241  }
242
243  while (beat + bp <= step) {
244    beat += bp;
[3de10bb]245    //AUBIO_DBG ("beat: %d, %f, %f\n", i, bp, beat);
[d207300]246    output->data[i] = beat;
[7bf3dcb]247    i++;
248  }
249
250  bt->lastbeat = beat;
[17b7d66]251  /* store the number of beats in this frame as the first element */
[d207300]252  output->data[0] = i;
[b78805a]253}
254
[7bf3dcb]255uint_t
256fvec_gettimesig (fvec_t * acf, uint_t acflen, uint_t gp)
257{
258  sint_t k = 0;
259  smpl_t three_energy = 0., four_energy = 0.;
260  if (acflen > 6 * gp + 2) {
261    for (k = -2; k < 2; k++) {
[d207300]262      three_energy += acf->data[3 * gp + k];
263      four_energy += acf->data[4 * gp + k];
[7bf3dcb]264    }
265  } else {
266    /*Expanded to be more accurate in time sig estimation */
267    for (k = -2; k < 2; k++) {
[d207300]268      three_energy += acf->data[3 * gp + k] + acf->data[6 * gp + k];
269      four_energy += acf->data[4 * gp + k] + acf->data[2 * gp + k];
[7bf3dcb]270    }
271  }
272  return (three_energy > four_energy) ? 3 : 4;
[b78805a]273}
274
[7bf3dcb]275void
276aubio_beattracking_checkstate (aubio_beattracking_t * bt)
277{
[72db1cf]278  uint_t i, j, a, b;
[7bf3dcb]279  uint_t flagconst = 0;
280  sint_t counter = bt->counter;
281  uint_t flagstep = bt->flagstep;
282  smpl_t gp = bt->gp;
283  smpl_t bp = bt->bp;
284  smpl_t rp = bt->rp;
285  smpl_t rp1 = bt->rp1;
286  smpl_t rp2 = bt->rp2;
287  uint_t laglen = bt->rwv->length;
288  uint_t acflen = bt->acf->length;
289  uint_t step = bt->step;
290  fvec_t *acf = bt->acf;
291  fvec_t *acfout = bt->acfout;
292
293  if (gp) {
294    // doshiftfbank again only if context dependent model is in operation
295    //acfout = doshiftfbank(acf,gwv,timesig,laglen,acfout);
296    //don't need acfout now, so can reuse vector
297    // gwv is, in first loop, definitely all zeros, but will have
298    // proper values when context dependent model is activated
299    fvec_zeros (acfout);
300    for (i = 1; i < laglen - 1; i++) {
301      for (a = 1; a <= bt->timesig; a++) {
[72db1cf]302        for (b = 1; b < 2 * a; b++) {
303          acfout->data[i] += acf->data[i * a + b - 1];
[b78805a]304        }
[7bf3dcb]305      }
306    }
307    fvec_weight (acfout, bt->gwv);
[acd97d1]308    gp = fvec_quadratic_peak_pos (acfout, fvec_max_elem (acfout));
[7bf3dcb]309    /*
310       while(gp<32) gp =gp*2;
311       while(gp>64) gp = gp/2;
312     */
313  } else {
314    //still only using general model
315    gp = 0;
316  }
317
318  //now look for step change - i.e. a difference between gp and rp that
319  // is greater than 2*constthresh - always true in first case, since gp = 0
320  if (counter == 0) {
321    if (ABS (gp - rp) > 2. * bt->g_var) {
322      flagstep = 1;             // have observed  step change.
323      counter = 3;              // setup 3 frame counter
324    } else {
325      flagstep = 0;
326    }
327  }
328  //i.e. 3rd frame after flagstep initially set
329  if (counter == 1 && flagstep == 1) {
330    //check for consistency between previous beatperiod values
331    if (ABS (2. * rp - rp1 - rp2) < bt->g_var) {
332      //if true, can activate context dependent model
333      flagconst = 1;
334      counter = 0;              // reset counter and flagstep
335    } else {
336      //if not consistent, then don't flag consistency!
337      flagconst = 0;
338      counter = 2;              // let it look next time
339    }
340  } else if (counter > 0) {
341    //if counter doesn't = 1,
342    counter = counter - 1;
343  }
344
345  rp2 = rp1;
346  rp1 = rp;
347
348  if (flagconst) {
349    /* first run of new hypothesis */
350    gp = rp;
351    bt->timesig = fvec_gettimesig (acf, acflen, gp);
352    for (j = 0; j < laglen; j++)
[d207300]353      bt->gwv->data[j] =
[7bf3dcb]354          EXP (-.5 * SQR ((smpl_t) (j + 1. - gp)) / SQR (bt->g_var));
355    flagconst = 0;
356    bp = gp;
357    /* flat phase weighting */
358    fvec_ones (bt->phwv);
359  } else if (bt->timesig) {
360    /* context dependant model */
361    bp = gp;
362    /* gaussian phase weighting */
363    if (step > bt->lastbeat) {
364      for (j = 0; j < 2 * laglen; j++) {
[d207300]365        bt->phwv->data[j] =
[7bf3dcb]366            EXP (-.5 * SQR ((smpl_t) (1. + j - step +
367                    bt->lastbeat)) / (bp / 8.));
368      }
369    } else {
370      //AUBIO_DBG("NOT using phase weighting as step is %d and lastbeat %d \n",
371      //                step,bt->lastbeat);
372      fvec_ones (bt->phwv);
373    }
374  } else {
375    /* initial state */
376    bp = rp;
377    /* flat phase weighting */
378    fvec_ones (bt->phwv);
379  }
380
381  /* do some further checks on the final bp value */
382
383  /* if tempo is > 206 bpm, half it */
384  while (bp < 25) {
[17b7d66]385#if AUBIO_BEAT_WARNINGS
386    AUBIO_WRN ("doubling from %f (%f bpm) to %f (%f bpm)\n",
387        bp, 60.*44100./512./bp, bp/2., 60.*44100./512./bp/2. );
[7bf3dcb]388    //AUBIO_DBG("warning, halving the tempo from %f\n", 60.*samplerate/hopsize/bp);
[17b7d66]389#endif /* AUBIO_BEAT_WARNINGS */
[7bf3dcb]390    bp = bp * 2;
391  }
392
393  //AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp);
394
395  /* smoothing */
396  //bp = (uint_t) (0.8 * (smpl_t)bp + 0.2 * (smpl_t)bp2);
397  //AUBIO_DBG("tempo:\t%3.5f bpm smoothed | bp2 %d | bp %d | ", 5168./bp, bp2, bp);
398  //bp2 = bp;
399  //AUBIO_DBG("time signature: %d \n", bt->timesig);
400  bt->counter = counter;
401  bt->flagstep = flagstep;
402  bt->gp = gp;
403  bt->bp = bp;
404  bt->rp1 = rp1;
405  bt->rp2 = rp2;
[b78805a]406}
[416c0b5]407
[7bf3dcb]408smpl_t
409aubio_beattracking_get_bpm (aubio_beattracking_t * bt)
410{
[acd97d1]411  if (bt->bp != 0 && bt->timesig != 0 && bt->counter == 0 && bt->flagstep == 0) {
412    return 5168. / fvec_quadratic_peak_pos (bt->acfout, bt->bp);
[7bf3dcb]413  } else {
414    return 0.;
415  }
[416c0b5]416}
[e34b010]417
[7bf3dcb]418smpl_t
419aubio_beattracking_get_confidence (aubio_beattracking_t * bt)
420{
421  if (bt->gp) {
[1e2c82f]422    return fvec_max (bt->acfout);
[7bf3dcb]423  } else {
424    return 0.;
425  }
[e34b010]426}
Note: See TracBrowser for help on using the repository browser.