source: src/tempo/beattracking.c @ 54f2a70

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

src/tempo/beattracking.c: remove unused fvec_getperiod function

  • Property mode set to 100644
File size: 13.9 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);
27
28struct _aubio_beattracking_t {
29        fvec_t * rwv;    /** rayleigh weight vector - rayleigh distribution function */                   
30        fvec_t * gwv;    /** rayleigh weight vector - rayleigh distribution function */                   
31        fvec_t * dfwv;   /** detection function weighting - exponential curve */
32        fvec_t * dfrev;  /** reversed onset detection function */
33        fvec_t * acf;    /** vector for autocorrelation function (of current detection function frame) */
34        fvec_t * acfout; /** store result of passing acf through s.i.c.f.b. */
35        fvec_t * phwv;   /** beat expectation alignment weighting */
36        fvec_t * phout;
37        uint_t timesig;  /** time signature of input, set to zero until context dependent model activated */
38        uint_t step;
39        uint_t rayparam; /** Rayleigh parameter */
40        uint_t lastbeat;
41        sint_t counter;
42        uint_t flagstep;
43        smpl_t g_var;
44        uint_t gp;
45        uint_t bp;
46        uint_t rp;
47        uint_t rp1;
48        uint_t rp2;
49};
50
51aubio_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;
56        /* parameter for rayleigh weight vector - sets preferred tempo to
57         * 120bpm [43] */
58        smpl_t rayparam = 48./512. * winlen;
59        smpl_t dfwvnorm = EXP((LOG(2.0)/rayparam)*(winlen+2));
60        /** length over which beat period is found [128] */
61        uint_t laglen   = winlen/4;
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 */
65
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
73        p->rayparam = rayparam;
74        p->step    = step;
75        p->rwv     = new_fvec(laglen,channels);
76        p->gwv     = new_fvec(laglen,channels);
77        p->dfwv    = new_fvec(winlen,channels);
78        p->dfrev   = new_fvec(winlen,channels);
79        p->acf     = new_fvec(winlen,channels);
80        p->acfout  = new_fvec(laglen,channels);
81        p->phwv    = new_fvec(2*laglen,channels);
82        p->phout   = new_fvec(winlen,channels);
83
84        p->timesig = 0;
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
101void 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
114void aubio_beattracking_do(aubio_beattracking_t * bt, fvec_t * dfframe, fvec_t * output) {
115
116        uint_t i,k;
117        /* current beat period value found using gaussian weighting (from context dependent model) */
118        uint_t step     = bt->step;
119        uint_t laglen   = bt->rwv->length;
120        uint_t winlen   = bt->dfwv->length;
121        smpl_t * phout  = bt->phout->data[0];
122        smpl_t * phwv   = bt->phwv->data[0];
123        smpl_t * dfrev  = bt->dfrev->data[0];
124        smpl_t * dfwv   = bt->dfwv->data[0];
125        smpl_t * rwv    = bt->rwv->data[0];
126        smpl_t * acfout = bt->acfout->data[0];
127        smpl_t * acf    = bt->acf->data[0];
128        uint_t maxindex = 0;
129        //number of harmonics in shift invariant comb filterbank
130        uint_t numelem  = 4;
131
132        //parameters for making s.i.c.f.b.
133        uint_t a,b; 
134        //beat alignment
135        uint_t phase; 
136        uint_t kmax;
137        sint_t beat; 
138        uint_t bp;
139
140        for (i = 0; i < winlen; i++){
141                dfrev[winlen-1-i] = 0.;
142                dfrev[winlen-1-i] = dfframe->data[0][i]*dfwv[i];
143        }
144
145        /* find autocorrelation function */
146        aubio_autocorr(dfframe,bt->acf); 
147        /*
148        for (i = 0; i < winlen; i++){
149                AUBIO_DBG("%f,",acf[i]);
150        }
151        AUBIO_DBG("\n");
152        */
153
154        /* get acfout - assume Rayleigh weightvector only */
155        /* if timesig is unknown, use metrically unbiased version of filterbank */
156        if(!bt->timesig) 
157                numelem = 4;
158        //        AUBIO_DBG("using unbiased filterbank, timesig: %d\n", timesig);
159        else
160                numelem = bt->timesig;
161        //        AUBIO_DBG("using biased filterbank, timesig: %d\n", timesig);
162
163        /* first and last output values are left intentionally as zero */
164        for (i=0; i < bt->acfout->length; i++)
165                acfout[i] = 0.;
166
167        for(i=1;i<laglen-1;i++){ 
168                for (a=1; a<=numelem; a++){
169                        for(b=(1-a); b<a; b++){
170                                acfout[i] += acf[a*(i+1)+b-1] 
171                                        * 1./(2.*a-1.)*rwv[i];
172                        }
173                }
174        }
175
176        /* find non-zero Rayleigh period */
177        maxindex = vec_max_elem(bt->acfout);
178        bt->rp = maxindex ? maxindex : 1;
179        //rp = (maxindex==127) ? 43 : maxindex; //rayparam
180        bt->rp = (maxindex==bt->acfout->length-1) ? bt->rayparam : maxindex; //rayparam
181
182        /* activate biased filterbank */
183        aubio_beattracking_checkstate(bt);
184        bp = bt->bp;
185        /* end of biased filterbank */
186
187        /* initialize output */
188        for(i=0;i<bt->phout->length;i++)     {phout[i] = 0.;} 
189
190        /* deliberate integer operation, could be set to 3 max eventually */
191        kmax = winlen/bp;
192
193        for(i=0;i<bp;i++){
194                phout[i] = 0.;
195                for(k=0;k<kmax;k++){
196                        phout[i] += dfrev[i+bp*k] * phwv[i];
197                }
198        }
199
200        /* find Rayleigh period */
201        maxindex = vec_max_elem(bt->phout);
202        if (maxindex == winlen-1) maxindex = 0;
203        phase =  1 + maxindex;
204
205        /* debug */
206        //AUBIO_DBG("beat period = %d, rp1 = %d, rp2 = %d\n", bp, rp1, rp2);
207        //AUBIO_DBG("rp = %d, gp = %d, phase = %d\n", bt->rp, bt->gp, phase);
208
209        /* reset output */
210        for (i = 0; i < laglen; i++)
211                output->data[0][i] = 0.;
212
213        i = 1;
214        beat =  bp - phase;
215        /* start counting the beats */
216        if(beat >= 0)
217        {
218                output->data[0][i] = (smpl_t)beat;
219                i++;
220        }
221
222        while( beat+bp < step )
223        {
224                beat += bp;
225                output->data[0][i] = (smpl_t)beat;
226                i++;
227        }
228
229        bt->lastbeat = beat;
230        /* store the number of beat found in this frame as the first element */
231        output->data[0][0] = i;
232}
233
234uint_t fvec_gettimesig(smpl_t * acf, uint_t acflen, uint_t gp){
235        sint_t k = 0;
236        smpl_t three_energy = 0., four_energy = 0.;
237        if( acflen > 6 * gp + 2 ){
238                for(k=-2;k<2;k++){
239                        three_energy += acf[3*gp+k];
240                        four_energy += acf[4*gp+k];
241                }
242        }
243        else{ /*Expanded to be more accurate in time sig estimation*/
244                for(k=-2;k<2;k++){
245                        three_energy += acf[3*gp+k]+acf[6*gp+k];
246                        four_energy += acf[4*gp+k]+acf[2*gp+k];
247                }
248        }
249        return (three_energy > four_energy) ? 3 : 4;
250}
251
252void aubio_beattracking_checkstate(aubio_beattracking_t * bt) {
253        uint_t i,j,a,b;
254        uint_t flagconst  = 0;
255        sint_t counter  = bt->counter;
256        uint_t flagstep = bt->flagstep;
257        uint_t gp       = bt->gp;
258        uint_t bp       = bt->bp;
259        uint_t rp       = bt->rp;
260        uint_t rp1      = bt->rp1;
261        uint_t rp2      = bt->rp2;
262        uint_t laglen   = bt->rwv->length;
263        uint_t acflen   = bt->acf->length;
264        uint_t step     = bt->step;
265        smpl_t * acf    = bt->acf->data[0];
266        smpl_t * acfout = bt->acfout->data[0];
267        smpl_t * gwv    = bt->gwv->data[0];
268        smpl_t * phwv   = bt->phwv->data[0];
269
270        if (gp) {
271                // doshiftfbank again only if context dependent model is in operation
272                //acfout = doshiftfbank(acf,gwv,timesig,laglen,acfout);
273                //don't need acfout now, so can reuse vector
274                // gwv is, in first loop, definitely all zeros, but will have
275                // proper values when context dependent model is activated
276                for (i=0; i < bt->acfout->length; i++)
277                       acfout[i] = 0.;
278                for(i=1;i<laglen-1;i++){ 
279                        for (a=1;a<=bt->timesig;a++){
280                                for(b=(1-a);b<a;b++){
281                                        acfout[i] += acf[a*(i+1)+b-1] 
282                                                * 1. * gwv[i];
283                                }
284                        }
285                }
286                gp = vec_max_elem(bt->acfout);
287                /*
288                while(gp<32) gp =gp*2;
289                while(gp>64) gp = gp/2;
290                */
291        } else {
292                //still only using general model
293                gp = 0; 
294        }
295
296        //now look for step change - i.e. a difference between gp and rp that
297        // is greater than 2*constthresh - always true in first case, since gp = 0
298        if(counter == 0){
299                if(ABS(gp - rp) > 2.*bt->g_var) {
300                        flagstep = 1; // have observed  step change.
301                        counter  = 3; // setup 3 frame counter
302                } else {
303                        flagstep = 0;
304                }
305        }
306
307        //i.e. 3rd frame after flagstep initially set
308        if (counter==1 && flagstep==1) {
309                //check for consistency between previous beatperiod values
310                if(ABS(2.*rp - rp1 -rp2) < bt->g_var) {
311                        //if true, can activate context dependent model
312                        flagconst = 1;
313                        counter   = 0; // reset counter and flagstep
314                } else {
315                        //if not consistent, then don't flag consistency!
316                        flagconst = 0;
317                        counter   = 2; // let it look next time
318                }
319        } else if (counter > 0) {
320                //if counter doesn't = 1,
321                counter = counter-1;
322        }
323
324        rp2 = rp1; rp1 = rp; 
325
326        if (flagconst) {
327                /* first run of new hypothesis */
328                gp = rp;
329                bt->timesig = fvec_gettimesig(acf,acflen, gp);
330                for(j=0;j<laglen;j++)
331                        gwv[j] = EXP(-.5*SQR((smpl_t)(j+1.-gp))/SQR(bt->g_var));
332                flagconst = 0;
333                bp = gp;
334                /* flat phase weighting */
335                for(j=0;j<2*laglen;j++)  {phwv[j] = 1.;} 
336        } else if (bt->timesig) {
337                /* context dependant model */
338                bp = gp;
339                /* gaussian phase weighting */
340                if (step > bt->lastbeat) {
341                        for(j=0;j<2*laglen;j++)  {
342                                phwv[j] = EXP(-.5*SQR((smpl_t)(1.+j-step+bt->lastbeat))/(bp/8.));
343                        }
344                } else { 
345                        //AUBIO_DBG("NOT using phase weighting as step is %d and lastbeat %d \n",
346                        //                step,bt->lastbeat);
347                        for(j=0;j<2*laglen;j++)  {phwv[j] = 1.;} 
348                }
349        } else {
350                /* initial state */ 
351                bp = rp;
352                /* flat phase weighting */
353                for(j=0;j<2*laglen;j++)  {phwv[j] = 1.;} 
354        }
355
356        /* do some further checks on the final bp value */
357
358        /* if tempo is > 206 bpm, half it */
359        while (bp < 25) {
360                //AUBIO_DBG("warning, doubling the beat period from %d\n", bp);
361                //AUBIO_DBG("warning, halving the tempo from %f\n", 60.*samplerate/hopsize/bp);
362                bp = bp*2;
363        }
364       
365        //AUBIO_DBG("tempo:\t%3.5f bpm | ", 5168./bp);
366
367        /* smoothing */
368        //bp = (uint_t) (0.8 * (smpl_t)bp + 0.2 * (smpl_t)bp2);
369        //AUBIO_DBG("tempo:\t%3.5f bpm smoothed | bp2 %d | bp %d | ", 5168./bp, bp2, bp);
370        //bp2 = bp;
371        //AUBIO_DBG("time signature: %d \n", bt->timesig);
372        bt->counter = counter;
373        bt->flagstep = flagstep;
374        bt->gp = gp;
375        bt->bp = bp;
376        bt->rp1 = rp1;
377        bt->rp2 = rp2;
378
379}
380
381smpl_t aubio_beattracking_get_bpm(aubio_beattracking_t * bt) {
382        if (bt->timesig != 0 && bt->counter == 0 && bt->flagstep == 0) {
383          return 5168. / (smpl_t)bt->gp;
384        } else {
385          return 0.;
386        }
387}
388
389smpl_t aubio_beattracking_get_confidence(aubio_beattracking_t * bt) {
390        if (bt->gp) return vec_max(bt->acfout);
391        else return 0.;
392}
Note: See TracBrowser for help on using the repository browser.