source: src/tempo/beattracking.c @ 1e2c82f

feature/autosinkfeature/constantqfeature/pitchshiftfeature/pydocstringsfeature/timestretchpitchshiftsamplertimestretchyinfft+
Last change on this file since 1e2c82f was 1e2c82f, checked in by Paul Brossier <piem@piem.org>, 10 years ago

rename vec_max and vec_max_elem to fvec_max fvec_max_elem

  • Property mode set to 100644
File size: 12.2 KB
Line 
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
28uint_t fvec_gettimesig (fvec_t * acf, uint_t acflen, uint_t gp);
29void aubio_beattracking_checkstate (aubio_beattracking_t * bt);
30
31struct _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
55aubio_beattracking_t *
56new_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
106void
107del_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
121void
122aubio_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 = fvec_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 = fvec_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
254uint_t
255fvec_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
274void
275aubio_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, fvec_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
407smpl_t
408aubio_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
417smpl_t
418aubio_beattracking_get_confidence (aubio_beattracking_t * bt)
419{
420  if (bt->gp) {
421    return fvec_max (bt->acfout);
422  } else {
423    return 0.;
424  }
425}
Note: See TracBrowser for help on using the repository browser.