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