source: src/tempo/beattracking.c @ 6c7d49b

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

remove src/sample.h

  • Property mode set to 100644
File size: 16.0 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(smpl_t * acf, uint_t acflen, uint_t gp);
26void aubio_beattracking_checkstate(aubio_beattracking_t * bt);
27smpl_t fvec_getperiod(aubio_beattracking_t * bt);
28
29struct _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 */
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        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 */
43        uint_t lastbeat;
44        sint_t counter;
45        uint_t flagstep;
46        smpl_t g_var;
47        uint_t gp;
48        uint_t bp;
49        uint_t rp;
50        uint_t rp1;
51        uint_t rp2;
52};
53
54aubio_beattracking_t * new_aubio_beattracking(uint_t winlen,
55                uint_t channels) {
56
57        aubio_beattracking_t * p = AUBIO_NEW(aubio_beattracking_t);
58        uint_t i        = 0;
59        /* parameter for rayleigh weight vector - sets preferred tempo to
60         * 120bpm [43] */
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 */
68
69        uint_t maxnumelem = 4; /* max number of index output */
70        p->lastbeat = 0;
71        p->counter = 0;
72        p->flagstep = 0;
73        p->g_var = 3.901; // constthresh empirically derived!
74        p->rp = 1;
75        p->gp = 0;
76
77        p->rayparam = rayparam;
78        p->step    = step;
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);
87
88        p->timesig = 0;
89
90        p->inds    = new_fvec(maxnumelem,channels);
91        p->locacf  = new_fvec(winlen,channels); 
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
108void 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        del_fvec(p->locacf);
118        del_fvec(p->inds);
119        AUBIO_FREE(p);
120}
121
122
123void aubio_beattracking_do(aubio_beattracking_t * bt, fvec_t * dfframe, fvec_t * output) {
124
125        uint_t i,k;
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];
137        uint_t maxindex = 0;
138        //number of harmonics in shift invariant comb filterbank
139        uint_t numelem  = 4;
140
141        //smpl_t myperiod   = 0.;
142        //smpl_t * out    = output->data[0];
143
144        //parameters for making s.i.c.f.b.
145        uint_t a,b; 
146        //beat alignment
147        uint_t phase; 
148        uint_t kmax;
149        sint_t beat; 
150        uint_t bp;
151
152        for (i = 0; i < winlen; i++){
153                dfrev[winlen-1-i] = 0.;
154                dfrev[winlen-1-i] = dfframe->data[0][i]*dfwv[i];
155        }
156
157        /* find autocorrelation function */
158        aubio_autocorr(dfframe,bt->acf); 
159        /*
160        for (i = 0; i < winlen; i++){
161                AUBIO_DBG("%f,",acf[i]);
162        }
163        AUBIO_DBG("\n");
164        */
165
166        /* get acfout - assume Rayleigh weightvector only */
167        /* if timesig is unknown, use metrically unbiased version of filterbank */
168        if(!bt->timesig) 
169                numelem = 4;
170        //        AUBIO_DBG("using unbiased filterbank, timesig: %d\n", timesig);
171        else
172                numelem = bt->timesig;
173        //        AUBIO_DBG("using biased filterbank, timesig: %d\n", timesig);
174
175        /* first and last output values are left intentionally as zero */
176        for (i=0; i < bt->acfout->length; i++)
177                acfout[i] = 0.;
178
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];
184                        }
185                }
186        }
187
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
193
194        // get float period
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);
198
199        /* activate biased filterbank */
200        aubio_beattracking_checkstate(bt);
201        bp = bt->bp;
202        /* end of biased filterbank */
203
204        /* initialize output */
205        for(i=0;i<bt->phout->length;i++)     {phout[i] = 0.;} 
206
207        /* deliberate integer operation, could be set to 3 max eventually */
208        kmax = winlen/bp;
209
210        for(i=0;i<bp;i++){
211                phout[i] = 0.;
212                for(k=0;k<kmax;k++){
213                        phout[i] += dfrev[i+bp*k] * phwv[i];
214                }
215        }
216
217        /* find Rayleigh period */
218        maxindex = vec_max_elem(bt->phout);
219        if (maxindex == winlen-1) maxindex = 0;
220        phase =  1 + maxindex;
221
222        /* debug */
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);
225
226        /* reset output */
227        for (i = 0; i < laglen; i++)
228                output->data[0][i] = 0.;
229
230        i = 1;
231        beat =  bp - phase;
232        /* start counting the beats */
233        if(beat >= 0)
234        {
235                output->data[0][i] = (smpl_t)beat;
236                i++;
237        }
238
239        while( beat+bp < step )
240        {
241                beat += bp;
242                output->data[0][i] = (smpl_t)beat;
243                i++;
244        }
245
246        bt->lastbeat = beat;
247        /* store the number of beat found in this frame as the first element */
248        output->data[0][0] = i;
249}
250
251uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp){
252        sint_t k = 0;
253        smpl_t three_energy = 0., four_energy = 0.;
254        if( acflen > 6 * gp + 2 ){
255                for(k=-2;k<2;k++){
256                        three_energy += acf[3*gp+k];
257                        four_energy += acf[4*gp+k];
258                }
259        }
260        else{ /*Expanded to be more accurate in time sig estimation*/
261                for(k=-2;k<2;k++){
262                        three_energy += acf[3*gp+k]+acf[6*gp+k];
263                        four_energy += acf[4*gp+k]+acf[2*gp+k];
264                }
265        }
266        return (three_energy > four_energy) ? 3 : 4;
267}
268
269smpl_t fvec_getperiod(aubio_beattracking_t * bt){
270        /*function to make a more accurate beat period measurement.*/
271
272        smpl_t period = 0.;
273        smpl_t maxval = 0.;
274        uint_t numelem    = 4;
275       
276        sint_t a,b;
277        uint_t i,j;     
278        uint_t acfmi = bt->rp; //acfout max index
279        uint_t maxind = 0;
280
281        if(!bt->timesig)
282                numelem = 4;
283        else
284                numelem = bt->timesig;
285
286        for (i=0;i<numelem;i++) // initialize
287        bt->inds->data[0][i] = 0.;
288
289        for (i=0;i<bt->locacf->length;i++) // initialize
290                bt->locacf->data[0][i] = 0.;
291       
292        // get appropriate acf elements from acf and store in locacf
293        for (a=1;a<=4;a++){
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];                                                     
297                }
298        }
299
300        for(i=0;i<numelem;i++){
301
302                maxind = 0;
303                maxval = 0.0;
304       
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];
308                                maxind = j;
309                        }
310                        //bt->locacf->data[0][maxind] = 0.;
311                        bt->locacf->data[0][j] = 0.;
312                }
313                //AUBIO_DBG("\n maxind is  %d\n",maxind);
314                bt->inds->data[0][i] = maxind;
315
316        }
317
318        for (i=0;i<numelem;i++){
319                period += bt->inds->data[0][i]/(i+1.);}
320
321        period = period/numelem;
322
323        return (period);
324}
325
326
327void aubio_beattracking_checkstate(aubio_beattracking_t * bt) {
328        uint_t i,j,a,b;
329        uint_t flagconst  = 0;
330        sint_t counter  = bt->counter;
331        uint_t flagstep = bt->flagstep;
332        uint_t gp       = bt->gp;
333        uint_t bp       = bt->bp;
334        uint_t rp       = bt->rp;
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];
344
345        if (gp) {
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++)
352                       acfout[i] = 0.;
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] 
357                                                * 1. * gwv[i];
358                                }
359                        }
360                }
361                gp = vec_max_elem(bt->acfout);
362                /*
363                while(gp<32) gp =gp*2;
364                while(gp>64) gp = gp/2;
365                */
366        } else {
367                //still only using general model
368                gp = 0; 
369        }
370
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
373        if(counter == 0){
374                if(ABS(gp - rp) > 2.*bt->g_var) {
375                        flagstep = 1; // have observed  step change.
376                        counter  = 3; // setup 3 frame counter
377                } else {
378                        flagstep = 0;
379                }
380        }
381
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
387                        flagconst = 1;
388                        counter   = 0; // reset counter and flagstep
389                } else {
390                        //if not consistent, then don't flag consistency!
391                        flagconst = 0;
392                        counter   = 2; // let it look next time
393                }
394        } else if (counter > 0) {
395                //if counter doesn't = 1,
396                counter = counter-1;
397        }
398
399        rp2 = rp1; rp1 = rp; 
400
401        if (flagconst) {
402                /* first run of new hypothesis */
403                gp = rp;
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));
407                flagconst = 0;
408                bp = gp;
409                /* flat phase weighting */
410                for(j=0;j<2*laglen;j++)  {phwv[j] = 1.;} 
411        } else if (bt->timesig) {
412                /* context dependant model */
413                bp = gp;
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.));
418                        }
419                } else { 
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.;} 
423                }
424        } else {
425                /* initial state */ 
426                bp = rp;
427                /* flat phase weighting */
428                for(j=0;j<2*laglen;j++)  {phwv[j] = 1.;} 
429        }
430
431        /* do some further checks on the final bp value */
432
433        /* if tempo is > 206 bpm, half it */
434        while (bp < 25) {
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);
437                bp = bp*2;
438        }
439       
440        //AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp);
441
442        /* smoothing */
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);
445        //bp2 = bp;
446        //AUBIO_DBG("time signature: %d \n", bt->timesig);
447        bt->counter = counter;
448        bt->flagstep = flagstep;
449        bt->gp = gp;
450        bt->bp = bp;
451        bt->rp1 = rp1;
452        bt->rp2 = rp2;
453
454}
455
456smpl_t aubio_beattracking_get_bpm(aubio_beattracking_t * bt) {
457        if (bt->timesig != 0 && bt->counter == 0 && bt->flagstep == 0) {
458          return 5168. / (smpl_t)bt->gp;
459        } else {
460          return 0.;
461        }
462}
463
464smpl_t aubio_beattracking_get_confidence(aubio_beattracking_t * bt) {
465        if (bt->gp) return vec_max(bt->acfout);
466        else return 0.;
467}
Note: See TracBrowser for help on using the repository browser.