source: src/beattracking.c @ 0b9e1da5

feature/autosinkfeature/cnnfeature/cnn_orgfeature/constantqfeature/crepefeature/crepe_orgfeature/pitchshiftfeature/pydocstringsfeature/timestretchfix/ffmpeg5pitchshiftsamplertimestretchyinfft+
Last change on this file since 0b9e1da5 was 0b9e1da5, checked in by Paul Brossier <piem@altern.org>, 18 years ago

fix signed/unsigned mismatches in beattracking
fix signed/unsigned mismatches in beattracking

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