2 Copyright (C) 2005 Matthew Davies and Paul Brossier
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.
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.
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.
20 #include "aubio_priv.h"
22 #include "mathutils.h"
23 #include "beattracking.h"
25 uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp);
26 void aubio_beattracking_checkstate(aubio_beattracking_t * bt);
27 smpl_t fvec_getperiod(aubio_beattracking_t * bt);
29 struct _aubio_beattracking_t {
30 fvec_t * rwv; /** rayleigh weight vector - rayleigh distribution function */
31 fvec_t * gwv; /** rayleigh weight vector - rayleigh distribution function */
32 fvec_t * dfwv; /** detection function weighting - exponential curve */
33 fvec_t * dfrev; /** reversed onset detection function */
34 fvec_t * acf; /** vector for autocorrelation function (of current detection function frame) */
35 fvec_t * acfout; /** store result of passing acf through s.i.c.f.b. */
36 fvec_t * phwv; /** beat expectation alignment weighting */
38 uint_t timesig; /** time signature of input, set to zero until context dependent model activated */
40 fvec_t * locacf; /** vector to store harmonics of filterbank of acf */
41 fvec_t * inds; /** vector for max index outputs for each harmonic */
42 uint_t rayparam; /** Rayleigh parameter */
54 aubio_beattracking_t * new_aubio_beattracking(uint_t winlen,
57 aubio_beattracking_t * p = AUBIO_NEW(aubio_beattracking_t);
59 /* parameter for rayleigh weight vector - sets preferred tempo to
61 smpl_t rayparam = 48./512. * winlen;
62 smpl_t dfwvnorm = EXP((LOG(2.0)/rayparam)*(winlen+2));
63 /** length over which beat period is found [128] */
64 uint_t laglen = winlen/4;
65 /** step increment - both in detection function samples -i.e. 11.6ms or
66 * 1 onset frame [128] */
67 uint_t step = winlen/4; /* 1.5 seconds */
69 uint_t maxnumelem = 4; /* max number of index output */
73 p->g_var = 3.901; // constthresh empirically derived!
77 p->rayparam = rayparam;
79 p->rwv = new_fvec(laglen,channels);
80 p->gwv = new_fvec(laglen,channels);
81 p->dfwv = new_fvec(winlen,channels);
82 p->dfrev = new_fvec(winlen,channels);
83 p->acf = new_fvec(winlen,channels);
84 p->acfout = new_fvec(laglen,channels);
85 p->phwv = new_fvec(2*laglen,channels);
86 p->phout = new_fvec(winlen,channels);
90 p->inds = new_fvec(maxnumelem,channels);
91 p->locacf = new_fvec(winlen,channels);
93 /* exponential weighting, dfwv = 0.5 when i = 43 */
94 for (i=0;i<winlen;i++) {
95 p->dfwv->data[0][i] = (EXP((LOG(2.0)/rayparam)*(i+1)))
99 for (i=0;i<(laglen);i++){
100 p->rwv->data[0][i] = ((smpl_t)(i+1.) / SQR((smpl_t)rayparam)) *
101 EXP((-SQR((smpl_t)(i+1.)) / (2.*SQR((smpl_t)rayparam))));
108 void del_aubio_beattracking(aubio_beattracking_t * p) {
123 void aubio_beattracking_do(aubio_beattracking_t * bt, fvec_t * dfframe, fvec_t * output) {
126 /* current beat period value found using gaussian weighting (from context dependent model) */
127 uint_t step = bt->step;
128 uint_t laglen = bt->rwv->length;
129 uint_t winlen = bt->dfwv->length;
130 smpl_t * phout = bt->phout->data[0];
131 smpl_t * phwv = bt->phwv->data[0];
132 smpl_t * dfrev = bt->dfrev->data[0];
133 smpl_t * dfwv = bt->dfwv->data[0];
134 smpl_t * rwv = bt->rwv->data[0];
135 smpl_t * acfout = bt->acfout->data[0];
136 smpl_t * acf = bt->acf->data[0];
138 //number of harmonics in shift invariant comb filterbank
141 //smpl_t myperiod = 0.;
142 //smpl_t * out = output->data[0];
144 //parameters for making s.i.c.f.b.
152 for (i = 0; i < winlen; i++){
153 dfrev[winlen-1-i] = 0.;
154 dfrev[winlen-1-i] = dfframe->data[0][i]*dfwv[i];
157 /* find autocorrelation function */
158 aubio_autocorr(dfframe,bt->acf);
160 for (i = 0; i < winlen; i++){
161 AUBIO_DBG("%f,",acf[i]);
166 /* get acfout - assume Rayleigh weightvector only */
167 /* if timesig is unknown, use metrically unbiased version of filterbank */
170 // AUBIO_DBG("using unbiased filterbank, timesig: %d\n", timesig);
172 numelem = bt->timesig;
173 // AUBIO_DBG("using biased filterbank, timesig: %d\n", timesig);
175 /* first and last output values are left intentionally as zero */
176 for (i=0; i < bt->acfout->length; i++)
179 for(i=1;i<laglen-1;i++){
180 for (a=1; a<=numelem; a++){
181 for(b=(1-a); b<a; b++){
182 acfout[i] += acf[a*(i+1)+b-1]
183 * 1./(2.*a-1.)*rwv[i];
188 /* find non-zero Rayleigh period */
189 maxindex = vec_max_elem(bt->acfout);
190 bt->rp = maxindex ? maxindex : 1;
191 //rp = (maxindex==127) ? 43 : maxindex; //rayparam
192 bt->rp = (maxindex==bt->acfout->length-1) ? bt->rayparam : maxindex; //rayparam
195 //myperiod = fvec_getperiod(bt);
196 //AUBIO_DBG("\nrp = %d myperiod = %f\n",bt->rp,myperiod);
197 //AUBIO_DBG("accurate tempo is %f bpm\n",5168./myperiod);
199 /* activate biased filterbank */
200 aubio_beattracking_checkstate(bt);
202 /* end of biased filterbank */
204 /* initialize output */
205 for(i=0;i<bt->phout->length;i++) {phout[i] = 0.;}
207 /* deliberate integer operation, could be set to 3 max eventually */
213 phout[i] += dfrev[i+bp*k] * phwv[i];
217 /* find Rayleigh period */
218 maxindex = vec_max_elem(bt->phout);
219 if (maxindex == winlen-1) maxindex = 0;
220 phase = 1 + maxindex;
223 //AUBIO_DBG("beat period = %d, rp1 = %d, rp2 = %d\n", bp, rp1, rp2);
224 //AUBIO_DBG("rp = %d, gp = %d, phase = %d\n", bt->rp, bt->gp, phase);
227 for (i = 0; i < laglen; i++)
228 output->data[0][i] = 0.;
232 /* start counting the beats */
235 output->data[0][i] = (smpl_t)beat;
239 while( beat+bp < step )
242 output->data[0][i] = (smpl_t)beat;
247 /* store the number of beat found in this frame as the first element */
248 output->data[0][0] = i;
251 uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp){
253 smpl_t three_energy = 0., four_energy = 0.;
254 if( acflen > 6 * gp + 2 ){
256 three_energy += acf[3*gp+k];
257 four_energy += acf[4*gp+k];
260 else{ /*Expanded to be more accurate in time sig estimation*/
262 three_energy += acf[3*gp+k]+acf[6*gp+k];
263 four_energy += acf[4*gp+k]+acf[2*gp+k];
266 return (three_energy > four_energy) ? 3 : 4;
269 smpl_t fvec_getperiod(aubio_beattracking_t * bt){
270 /*function to make a more accurate beat period measurement.*/
278 uint_t acfmi = bt->rp; //acfout max index
284 numelem = bt->timesig;
286 for (i=0;i<numelem;i++) // initialize
287 bt->inds->data[0][i] = 0.;
289 for (i=0;i<bt->locacf->length;i++) // initialize
290 bt->locacf->data[0][i] = 0.;
292 // get appropriate acf elements from acf and store in locacf
294 for(b=(1-a);b<a;b++){
295 bt->locacf->data[0][a*(acfmi)+b-1] =
296 bt->acf->data[0][a*(acfmi)+b-1];
300 for(i=0;i<numelem;i++){
305 for (j=0;j<(acfmi*(i+1)+(i)); j++){
306 if(bt->locacf->data[0][j]>maxval){
307 maxval = bt->locacf->data[0][j];
310 //bt->locacf->data[0][maxind] = 0.;
311 bt->locacf->data[0][j] = 0.;
313 //AUBIO_DBG("\n maxind is %d\n",maxind);
314 bt->inds->data[0][i] = maxind;
318 for (i=0;i<numelem;i++){
319 period += bt->inds->data[0][i]/(i+1.);}
321 period = period/numelem;
327 void aubio_beattracking_checkstate(aubio_beattracking_t * bt) {
329 uint_t flagconst = 0;
330 sint_t counter = bt->counter;
331 uint_t flagstep = bt->flagstep;
335 uint_t rp1 = bt->rp1;
336 uint_t rp2 = bt->rp2;
337 uint_t laglen = bt->rwv->length;
338 uint_t acflen = bt->acf->length;
339 uint_t step = bt->step;
340 smpl_t * acf = bt->acf->data[0];
341 smpl_t * acfout = bt->acfout->data[0];
342 smpl_t * gwv = bt->gwv->data[0];
343 smpl_t * phwv = bt->phwv->data[0];
346 // doshiftfbank again only if context dependent model is in operation
347 //acfout = doshiftfbank(acf,gwv,timesig,laglen,acfout);
348 //don't need acfout now, so can reuse vector
349 // gwv is, in first loop, definitely all zeros, but will have
350 // proper values when context dependent model is activated
351 for (i=0; i < bt->acfout->length; i++)
353 for(i=1;i<laglen-1;i++){
354 for (a=1;a<=bt->timesig;a++){
355 for(b=(1-a);b<a;b++){
356 acfout[i] += acf[a*(i+1)+b-1]
361 gp = vec_max_elem(bt->acfout);
363 while(gp<32) gp =gp*2;
364 while(gp>64) gp = gp/2;
367 //still only using general model
371 //now look for step change - i.e. a difference between gp and rp that
372 // is greater than 2*constthresh - always true in first case, since gp = 0
374 if(ABS(gp - rp) > 2.*bt->g_var) {
375 flagstep = 1; // have observed step change.
376 counter = 3; // setup 3 frame counter
382 //i.e. 3rd frame after flagstep initially set
383 if (counter==1 && flagstep==1) {
384 //check for consistency between previous beatperiod values
385 if(ABS(2.*rp - rp1 -rp2) < bt->g_var) {
386 //if true, can activate context dependent model
388 counter = 0; // reset counter and flagstep
390 //if not consistent, then don't flag consistency!
392 counter = 2; // let it look next time
394 } else if (counter > 0) {
395 //if counter doesn't = 1,
402 /* first run of new hypothesis */
404 bt->timesig = fvec_gettimesig(acf,acflen, gp);
405 for(j=0;j<laglen;j++)
406 gwv[j] = EXP(-.5*SQR((smpl_t)(j+1.-gp))/SQR(bt->g_var));
409 /* flat phase weighting */
410 for(j=0;j<2*laglen;j++) {phwv[j] = 1.;}
411 } else if (bt->timesig) {
412 /* context dependant model */
414 /* gaussian phase weighting */
415 if (step > bt->lastbeat) {
416 for(j=0;j<2*laglen;j++) {
417 phwv[j] = EXP(-.5*SQR((smpl_t)(1.+j-step+bt->lastbeat))/(bp/8.));
420 //AUBIO_DBG("NOT using phase weighting as step is %d and lastbeat %d \n",
421 // step,bt->lastbeat);
422 for(j=0;j<2*laglen;j++) {phwv[j] = 1.;}
427 /* flat phase weighting */
428 for(j=0;j<2*laglen;j++) {phwv[j] = 1.;}
431 /* do some further checks on the final bp value */
433 /* if tempo is > 206 bpm, half it */
435 //AUBIO_DBG("warning, doubling the beat period from %d\n", bp);
436 //AUBIO_DBG("warning, halving the tempo from %f\n", 60.*samplerate/hopsize/bp);
440 //AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp);
443 //bp = (uint_t) (0.8 * (smpl_t)bp + 0.2 * (smpl_t)bp2);
444 //AUBIO_DBG("tempo:\t%3.5f bpm smoothed | bp2 %d | bp %d | ", 5168./bp, bp2, bp);
446 //AUBIO_DBG("time signature: %d \n", bt->timesig);
447 bt->counter = counter;
448 bt->flagstep = flagstep;