[b78805a] | 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 "sample.h" |
---|
| 22 | #include "mathutils.h" |
---|
| 23 | #include "beattracking.h" |
---|
| 24 | |
---|
| 25 | /* maximum length for rp */ |
---|
| 26 | static uint_t MAXRP = 512; |
---|
| 27 | static smpl_t constthresh = 3.901; //empirically derived! |
---|
| 28 | static smpl_t stepthresh = 2.*3.901; |
---|
| 29 | |
---|
| 30 | uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp); |
---|
| 31 | void aubio_beattracking_checkstate(aubio_beattracking_t * bt); |
---|
[2cdae81] | 32 | smpl_t fvec_getperiod(smpl_t * acf, uint_t timesig, uint_t rp); |
---|
[b78805a] | 33 | |
---|
| 34 | /* could move to struct */ |
---|
[2cdae81] | 35 | uint_t gp=0, bp = 0, rp1 = 0, rp2 = 0, bp2 = 0; |
---|
[b78805a] | 36 | smpl_t g_mu =0.; |
---|
| 37 | smpl_t g_var = 3.901; |
---|
| 38 | uint_t flagconst = 0; |
---|
| 39 | uint_t flagstep = 0; |
---|
| 40 | // needs to be a signed ? |
---|
| 41 | sint_t counter = 0; |
---|
| 42 | uint_t maxindex = 0; |
---|
| 43 | uint_t timesig = 0; |
---|
| 44 | uint_t rp = 1; |
---|
| 45 | uint_t lastbeat = 0; |
---|
| 46 | //number of harmonics in shift invariant comb filterbank |
---|
| 47 | uint_t numelem = 4; |
---|
[2cdae81] | 48 | smpl_t myperiod = 0.; |
---|
| 49 | |
---|
| 50 | |
---|
| 51 | // FIXME |
---|
| 52 | uint_t maxnumelem = 4; |
---|
| 53 | uint_t len = 512; |
---|
| 54 | smpl_t inds[4]; //vector for max index outputs for each harmonic |
---|
| 55 | smpl_t local_acf[512]; //vector to store harmonics of filterbank of acf |
---|
| 56 | |
---|
[b78805a] | 57 | |
---|
| 58 | struct _aubio_beattracking_t { |
---|
| 59 | fvec_t * rwv; /* rayleigh weight vector - rayleigh distribution function */ |
---|
| 60 | fvec_t * gwv; /* rayleigh weight vector - rayleigh distribution function */ |
---|
| 61 | fvec_t * dfwv; /* detection function weighting - exponential curve */ |
---|
| 62 | fvec_t * dfrev; /* reversed onset detection function */ |
---|
| 63 | fvec_t * acf; /* vector for autocorrelation function (of current detection function frame) */ |
---|
| 64 | fvec_t * acfout; /* store result of passing acf through s.i.c.f.b. */ |
---|
| 65 | fvec_t * phwv; /* beat expectation alignment weighting */ |
---|
| 66 | fvec_t * phout; |
---|
| 67 | //uint_t timesig; /* time signature of input, set to zero until context dependent model activated */ |
---|
| 68 | uint_t step; |
---|
| 69 | }; |
---|
| 70 | |
---|
| 71 | aubio_beattracking_t * new_aubio_beattracking(uint_t winlen, |
---|
| 72 | uint_t step, |
---|
| 73 | uint_t laglen, |
---|
| 74 | uint_t rayparam, |
---|
| 75 | uint_t channels) { |
---|
| 76 | |
---|
| 77 | aubio_beattracking_t * p = AUBIO_NEW(aubio_beattracking_t); |
---|
| 78 | |
---|
| 79 | uint_t i = 0; |
---|
| 80 | smpl_t dfwvnorm = EXP((LOG(2.0)/rayparam)*(winlen+2)); |
---|
| 81 | |
---|
| 82 | p->step = step; |
---|
| 83 | p->rwv = new_fvec(laglen,channels); |
---|
| 84 | p->gwv = new_fvec(laglen,channels); |
---|
| 85 | p->dfwv = new_fvec(winlen,channels); |
---|
| 86 | p->dfrev = new_fvec(winlen,channels); |
---|
| 87 | p->acf = new_fvec(winlen,channels); |
---|
| 88 | p->acfout = new_fvec(laglen,channels); |
---|
| 89 | p->phwv = new_fvec(2*laglen,channels); |
---|
| 90 | p->phout = new_fvec(MAXRP,channels); |
---|
| 91 | //p->timesig = 0; |
---|
| 92 | |
---|
| 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))) |
---|
| 96 | / dfwvnorm; |
---|
| 97 | } |
---|
| 98 | |
---|
| 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)))); |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | return p; |
---|
| 105 | |
---|
| 106 | } |
---|
| 107 | |
---|
| 108 | void del_aubio_beattracking(aubio_beattracking_t * p) { |
---|
| 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 | |
---|
| 121 | void aubio_beattracking_do(aubio_beattracking_t * bt, fvec_t * dfframe, fvec_t * output) { |
---|
| 122 | |
---|
| 123 | uint_t i,k; |
---|
| 124 | /* current beat period value found using gaussian weighting (from context dependent model) */ |
---|
| 125 | uint_t step = bt->step; |
---|
| 126 | uint_t laglen = bt->rwv->length; |
---|
| 127 | uint_t winlen = bt->dfwv->length; |
---|
| 128 | smpl_t * phout = bt->phout->data[0]; |
---|
| 129 | smpl_t * phwv = bt->phwv->data[0]; |
---|
| 130 | smpl_t * dfrev = bt->dfrev->data[0]; |
---|
| 131 | smpl_t * dfwv = bt->dfwv->data[0]; |
---|
| 132 | smpl_t * rwv = bt->rwv->data[0]; |
---|
| 133 | smpl_t * acfout = bt->acfout->data[0]; |
---|
| 134 | smpl_t * acf = bt->acf->data[0]; |
---|
| 135 | //smpl_t * out = output->data[0]; |
---|
| 136 | |
---|
| 137 | //parameters for making s.i.c.f.b. |
---|
| 138 | uint_t a,b; |
---|
| 139 | //beat alignment |
---|
| 140 | uint_t phase; |
---|
| 141 | uint_t kmax; |
---|
| 142 | uint_t beat; |
---|
| 143 | |
---|
| 144 | for (i = 0; i < winlen; i++){ |
---|
| 145 | dfrev[winlen-1-i] = 0.; |
---|
| 146 | dfrev[winlen-1-i] = dfframe->data[0][i]*dfwv[i]; |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | /* find autocorrelation function */ |
---|
| 150 | aubio_autocorr(dfframe,bt->acf); |
---|
| 151 | /* |
---|
| 152 | for (i = 0; i < winlen; i++){ |
---|
| 153 | AUBIO_DBG("%f,",acf[i]); |
---|
| 154 | } |
---|
| 155 | AUBIO_DBG("\n"); |
---|
| 156 | */ |
---|
| 157 | |
---|
| 158 | /* get acfout - assume Rayleigh weightvector only */ |
---|
| 159 | /* if timesig is unknown, use metrically unbiased version of filterbank */ |
---|
| 160 | //if(!timesig) |
---|
| 161 | // AUBIO_DBG("using unbiased filterbank, timesig: %d\n", timesig); |
---|
| 162 | //else |
---|
| 163 | // AUBIO_DBG("using biased filterbank, timesig: %d\n", timesig); |
---|
| 164 | if(!timesig) |
---|
| 165 | numelem = 4; |
---|
| 166 | else |
---|
| 167 | numelem = timesig; |
---|
| 168 | |
---|
| 169 | /* first and last output values are left intentionally as zero */ |
---|
| 170 | for (i=0; i < bt->acfout->length; i++) |
---|
| 171 | acfout[i] = 0.; |
---|
[2cdae81] | 172 | |
---|
[b78805a] | 173 | for(i=1;i<laglen-1;i++){ |
---|
[2cdae81] | 174 | for (a=1; a<=numelem; a++){ |
---|
| 175 | for(b=(1-a); b<a; b++){ |
---|
[b78805a] | 176 | acfout[i] += acf[a*(i+1)+b-1] |
---|
| 177 | * 1./(2.*a-1.)*rwv[i]; |
---|
| 178 | } |
---|
| 179 | } |
---|
| 180 | } |
---|
| 181 | |
---|
| 182 | /* find non-zero Rayleigh period */ |
---|
| 183 | maxindex = vec_max_elem(bt->acfout); |
---|
| 184 | rp = maxindex ? maxindex : 1; |
---|
| 185 | rp = (maxindex==127) ? 43 : maxindex; //rayparam |
---|
| 186 | |
---|
[2cdae81] | 187 | // get float period |
---|
| 188 | myperiod = fvec_getperiod(acf,timesig,rp); |
---|
| 189 | //AUBIO_DBG("\nrp = %d myperiod = %f\n",rp,myperiod); |
---|
| 190 | //AUBIO_DBG("accurate tempo is %f bpm\n",5168./myperiod); |
---|
| 191 | |
---|
| 192 | |
---|
[b78805a] | 193 | /* activate biased filterbank */ |
---|
| 194 | aubio_beattracking_checkstate(bt); |
---|
| 195 | /* end of biased filterbank */ |
---|
| 196 | |
---|
| 197 | /* initialize output */ |
---|
| 198 | for(i=0;i<MAXRP;i++) {phout[i] = 0.;} |
---|
| 199 | |
---|
| 200 | /* deliberate integer operation */ |
---|
| 201 | /* could be set to 3 max eventually */ |
---|
| 202 | kmax = winlen/bp; |
---|
| 203 | |
---|
| 204 | for(i=0;i<bp;i++){ |
---|
| 205 | phout[i] = 0.; |
---|
| 206 | for(k=0;k<kmax;k++){ |
---|
| 207 | phout[i] += dfrev[i+bp*k] * phwv[i]; |
---|
| 208 | } |
---|
| 209 | } |
---|
| 210 | |
---|
| 211 | /* find Rayleigh period */ |
---|
| 212 | maxindex = vec_max_elem(bt->phout); |
---|
| 213 | if (maxindex == winlen-1) maxindex = 0; |
---|
| 214 | phase = 1 + maxindex; |
---|
| 215 | |
---|
| 216 | /* debug */ |
---|
| 217 | //AUBIO_DBG("beat period = %d, rp1 = %d, rp2 = %d\n", bp, rp1, rp2); |
---|
| 218 | //AUBIO_DBG("rp = %d, gp = %d, phase = %d\n", rp, gp, phase); |
---|
| 219 | |
---|
| 220 | /* reset output */ |
---|
| 221 | for (i = 0; i < laglen; i++) |
---|
| 222 | output->data[0][i] = 0.; |
---|
| 223 | |
---|
| 224 | i = 1; |
---|
| 225 | beat = bp - phase; |
---|
| 226 | /* start counting the beats */ |
---|
[2cdae81] | 227 | if(beat >= 0) |
---|
[b78805a] | 228 | { |
---|
| 229 | output->data[0][i] = (smpl_t)beat; |
---|
| 230 | i++; |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | while( beat+bp < step ) |
---|
| 234 | { |
---|
| 235 | beat += bp; |
---|
| 236 | output->data[0][i] = (smpl_t)beat; |
---|
| 237 | i++; |
---|
| 238 | } |
---|
| 239 | |
---|
| 240 | lastbeat = beat; |
---|
| 241 | |
---|
| 242 | output->data[0][0] = i; |
---|
| 243 | |
---|
| 244 | } |
---|
| 245 | |
---|
| 246 | uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp){ |
---|
| 247 | sint_t k = 0; |
---|
| 248 | smpl_t three_energy = 0., four_energy = 0.; |
---|
| 249 | if( acflen > 6 * gp + 2 ){ |
---|
| 250 | for(k=-2;k<2;k++){ |
---|
| 251 | three_energy += acf[3*gp+k]; |
---|
| 252 | four_energy += acf[4*gp+k]; |
---|
| 253 | } |
---|
| 254 | } |
---|
| 255 | else{ /*Expanded to be more accurate in time sig estimation*/ |
---|
| 256 | for(k=-2;k<2;k++){ |
---|
| 257 | three_energy += acf[3*gp+k]+acf[6*gp+k]; |
---|
| 258 | four_energy += acf[4*gp+k]+acf[2*gp+k]; |
---|
| 259 | } |
---|
| 260 | } |
---|
| 261 | return (three_energy > four_energy) ? 3 : 4; |
---|
| 262 | } |
---|
| 263 | |
---|
[2cdae81] | 264 | smpl_t fvec_getperiod(smpl_t * acf, uint_t timesig, uint_t rp){ |
---|
| 265 | /*function to make a more accurate beat period measurement.*/ |
---|
| 266 | |
---|
| 267 | smpl_t period = 0.; |
---|
| 268 | smpl_t maxval = 0.; |
---|
| 269 | |
---|
| 270 | |
---|
| 271 | |
---|
| 272 | int a,b; |
---|
| 273 | uint_t i,j; |
---|
| 274 | uint_t acfmi = rp; //acfout max index |
---|
| 275 | uint_t maxind = 0; |
---|
| 276 | |
---|
| 277 | if(!timesig) |
---|
| 278 | numelem = 4; |
---|
| 279 | else |
---|
| 280 | numelem = timesig; |
---|
| 281 | |
---|
| 282 | for (i=0;i<numelem;i++) // initialize |
---|
| 283 | inds[i] = 0.; |
---|
| 284 | |
---|
| 285 | for (i=0;i<len;i++) // initialize |
---|
| 286 | local_acf[i] = 0.; |
---|
| 287 | |
---|
| 288 | // get appropriate acf elements from acf and store in local_acf |
---|
| 289 | for (a=1;a<=4;a++){ |
---|
| 290 | for(b=(1-a);b<a;b++){ |
---|
| 291 | local_acf[a*(acfmi)+b-1] = acf[a*(acfmi)+b-1]; |
---|
| 292 | } |
---|
| 293 | } |
---|
| 294 | |
---|
| 295 | for(i=0;i<numelem;i++){ |
---|
| 296 | |
---|
| 297 | maxindex = 0; |
---|
| 298 | maxval = 0.0; |
---|
| 299 | |
---|
| 300 | for (j=0;j<(acfmi*(i+1)+(i)); j++){ |
---|
| 301 | if(local_acf[j]>maxval){ |
---|
| 302 | maxval = local_acf[j]; |
---|
| 303 | maxind = j; |
---|
| 304 | } |
---|
| 305 | local_acf[maxind] = 0.; |
---|
| 306 | } |
---|
| 307 | //printf("\n\n"); |
---|
| 308 | //AUBIO_DBG("\n maxind is %d\n",maxind); |
---|
| 309 | inds[i] = maxind; |
---|
| 310 | |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | for (i=0;i<numelem;i++){ |
---|
| 314 | period += inds[i]/(i+1.);} |
---|
| 315 | |
---|
| 316 | period = period/numelem; |
---|
| 317 | |
---|
| 318 | return (period); |
---|
| 319 | } |
---|
| 320 | |
---|
| 321 | |
---|
[b78805a] | 322 | void aubio_beattracking_checkstate(aubio_beattracking_t * bt) { |
---|
| 323 | uint_t i,j,a,b; |
---|
| 324 | uint_t laglen = bt->rwv->length; |
---|
| 325 | uint_t acflen = bt->acf->length; |
---|
| 326 | uint_t step = bt->step; |
---|
| 327 | smpl_t * acf = bt->acf->data[0]; |
---|
| 328 | smpl_t * acfout = bt->acfout->data[0]; |
---|
| 329 | smpl_t * gwv = bt->gwv->data[0]; |
---|
| 330 | smpl_t * phwv = bt->phwv->data[0]; |
---|
| 331 | |
---|
| 332 | if (gp) { |
---|
| 333 | // doshiftfbank again only if context dependent model is in operation |
---|
| 334 | //acfout = doshiftfbank(acf,gwv,timesig,laglen,acfout); |
---|
| 335 | //don't need acfout now, so can reuse vector |
---|
| 336 | // gwv is, in first loop, definitely all zeros, but will have |
---|
| 337 | // proper values when context dependent model is activated |
---|
| 338 | for (i=0; i < bt->acfout->length; i++) |
---|
| 339 | acfout[i] = 0.; |
---|
| 340 | for(i=1;i<laglen-1;i++){ |
---|
| 341 | for (a=1;a<=timesig;a++){ |
---|
| 342 | for(b=(1-a);b<a;b++){ |
---|
| 343 | acfout[i] += acf[a*(i+1)+b-1] |
---|
| 344 | * 1. * gwv[i]; |
---|
| 345 | } |
---|
| 346 | } |
---|
| 347 | } |
---|
| 348 | gp = vec_max_elem(bt->acfout); |
---|
[2cdae81] | 349 | /* |
---|
| 350 | while(gp<32) gp =gp*2; |
---|
| 351 | while(gp>64) gp = gp/2; |
---|
| 352 | */ |
---|
[b78805a] | 353 | } else { |
---|
| 354 | //still only using general model |
---|
| 355 | gp = 0; |
---|
| 356 | } |
---|
| 357 | |
---|
| 358 | //now look for step change - i.e. a difference between gp and rp that |
---|
| 359 | // is greater than stepthresh - always true in first case, since gp = 0 |
---|
| 360 | if(counter == 0){ |
---|
| 361 | if(ABS(gp - rp) > stepthresh) { |
---|
| 362 | flagstep = 1; // have observed step change. |
---|
| 363 | counter = 3; // setup 3 frame counter |
---|
| 364 | } else { |
---|
| 365 | flagstep = 0; |
---|
| 366 | } |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | //i.e. 3rd frame after flagstep initially set |
---|
| 370 | if (counter==1 && flagstep==1) { |
---|
| 371 | //check for consistency between previous beatperiod values |
---|
| 372 | if(ABS(2.*rp - rp1 -rp2) < constthresh) { |
---|
| 373 | //if true, can activate context dependent model |
---|
| 374 | flagconst = 1; |
---|
| 375 | counter = 0; // reset counter and flagstep |
---|
| 376 | } else { |
---|
| 377 | //if not consistent, then don't flag consistency! |
---|
| 378 | flagconst = 0; |
---|
| 379 | counter = 2; // let it look next time |
---|
| 380 | } |
---|
| 381 | } else if (counter > 0) { |
---|
| 382 | //if counter doesn't = 1, |
---|
| 383 | counter = counter-1; |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | rp2 = rp1; rp1 = rp; |
---|
| 387 | |
---|
| 388 | if (flagconst) { |
---|
| 389 | /* first run of new hypothesis */ |
---|
| 390 | gp = rp; |
---|
| 391 | g_mu = gp; |
---|
| 392 | timesig = fvec_gettimesig(acf,acflen, gp); |
---|
| 393 | for(j=0;j<laglen;j++) |
---|
| 394 | gwv[j] = EXP(-.5*SQR((smpl_t)(j+1.-g_mu))/SQR(g_var)); |
---|
| 395 | flagconst = 0; |
---|
| 396 | bp = gp; |
---|
| 397 | /* flat phase weighting */ |
---|
| 398 | for(j=0;j<2*laglen;j++) {phwv[j] = 1.;} |
---|
| 399 | } else if (timesig) { |
---|
| 400 | /* context dependant model */ |
---|
| 401 | bp = gp; |
---|
| 402 | /* gaussian phase weighting */ |
---|
| 403 | if (step > lastbeat) { |
---|
| 404 | for(j=0;j<2*laglen;j++) { |
---|
| 405 | phwv[j] = EXP(-.5*SQR((smpl_t)(1.+j-step+lastbeat))/(bp/8.)); |
---|
| 406 | } |
---|
| 407 | } else { |
---|
| 408 | AUBIO_DBG("NOT using phase weighting as step is %d and lastbeat %d \n", |
---|
| 409 | step,lastbeat); |
---|
| 410 | for(j=0;j<2*laglen;j++) {phwv[j] = 1.;} |
---|
| 411 | } |
---|
| 412 | } else { |
---|
| 413 | /* initial state */ |
---|
| 414 | bp = rp; |
---|
| 415 | /* flat phase weighting */ |
---|
| 416 | for(j=0;j<2*laglen;j++) {phwv[j] = 1.;} |
---|
| 417 | } |
---|
| 418 | |
---|
| 419 | /* do some further checks on the final bp value */ |
---|
[2cdae81] | 420 | |
---|
[b78805a] | 421 | /* if tempo is > 206 bpm, half it */ |
---|
| 422 | while (bp < 25) { |
---|
[2cdae81] | 423 | AUBIO_DBG("warning, halving the tempo from %f\n", 5168./bp); |
---|
[b78805a] | 424 | bp = bp*2; |
---|
| 425 | } |
---|
| 426 | |
---|
[2cdae81] | 427 | AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp); |
---|
| 428 | |
---|
| 429 | /* smoothing */ |
---|
| 430 | //bp = (uint_t) (0.8 * (smpl_t)bp + 0.2 * (smpl_t)bp2); |
---|
| 431 | //AUBIO_DBG("tempo:\t%3.5f bpm smoothed | bp2 %d | bp %d | ", 5168./bp, bp2, bp); |
---|
| 432 | //bp2 = bp; |
---|
| 433 | AUBIO_DBG("time signature: %d \n", timesig); |
---|
[b78805a] | 434 | |
---|
| 435 | } |
---|