source: src/tempo/beattracking.c @ 7bf3dcb

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

src/tempo/beattracking.c: no changes, indent

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