source: src/tempo/beattracking.c @ 6e4ef27

feature/cnn_orgfeature/crepe_org
Last change on this file since 6e4ef27 was 69f74f0, checked in by Paul Brossier <piem@piem.org>, 9 years ago

src/tempo/: add const qualifiers

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