1 /*
2 SuperCollider real time audio synthesis system
3 Copyright (c) 2002 James McCartney. All rights reserved.
4 http://www.audiosynth.com
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 // BeatTrack UGen implemented by Nick Collins (http://www.informatics.sussex.ac.uk/users/nc81/)
22 // post FFT UGen version 1 Nov 2007
23
24 // conversion of Matthew Davies autocorrelation beat tracking model, adapted for real-time use
25 // currently using QMUL complex domain onset detection function model
26
27 //#include "SC_PlugIn.h"
28 //#include <vecLib/vecLib.h>
29 //#include <string.h>
30 //#include <math.h>
31 //#include <stdlib.h>
32 //#include <stdio.h>
33
34 #include "ML.h"
35
36 // FFT data
37 //#define N 1024 //FFT size
38 // FFT size over 2
39 #define NOVER2 512
40 //#define NOVER4 256 //FFT size
41 //#define OVERLAP 512
42 //#define OVERLAPINDEX 512
43 //#define HOPSIZE 512
44
45 //#define FS 44100 //assumes fixed sampling rate
46 //#define FRAMESR 86.1328
47 // converted for different sampling rates
48 #define FRAMEPERIOD 0.01161
49 #define SKIP 128
50 //#define TIMEELAPSED 1.48608
51
52 // this data assumes LAGS is 128
53 static float g_m[128] = {
54 0.00054069, 0.00108050, 0.00161855, 0.00215399, 0.00268594, 0.00321356, 0.00373600, 0.00425243, 0.00476204,
55 0.00526404, 0.00575765, 0.00624213, 0.00671675, 0.00718080, 0.00763362, 0.00807455, 0.00850299, 0.00891836,
56 0.00932010, 0.00970771, 0.01008071, 0.01043866, 0.01078115, 0.01110782, 0.01141834, 0.01171242, 0.01198982,
57 0.01225033, 0.01249378, 0.01272003, 0.01292899, 0.01312061, 0.01329488, 0.01345182, 0.01359148, 0.01371396,
58 0.01381939, 0.01390794, 0.01397980, 0.01403520, 0.01407439, 0.01409768, 0.01410536, 0.01409780, 0.01407534,
59 0.01403838, 0.01398734, 0.01392264, 0.01384474, 0.01375410, 0.01365120, 0.01353654, 0.01341062, 0.01327397,
60 0.01312710, 0.01297054, 0.01280484, 0.01263053, 0.01244816, 0.01225827, 0.01206139, 0.01185807, 0.01164884,
61 0.01143424, 0.01121478, 0.01099099, 0.01076337, 0.01053241, 0.01029861, 0.01006244, 0.00982437, 0.00958484,
62 0.00934429, 0.00910314, 0.00886181, 0.00862067, 0.00838011, 0.00814049, 0.00790214, 0.00766540, 0.00743057,
63 0.00719793, 0.00696778, 0.00674036, 0.00651591, 0.00629466, 0.00607682, 0.00586256, 0.00565208, 0.00544551,
64 0.00524301, 0.00504470, 0.00485070, 0.00466109, 0.00447597, 0.00429540, 0.00411944, 0.00394813, 0.00378151,
65 0.00361959, 0.00346238, 0.00330989, 0.00316210, 0.00301899, 0.00288053, 0.00274669, 0.00261741, 0.00249266,
66 0.00237236, 0.00225646, 0.00214488, 0.00203755, 0.00193440, 0.00183532, 0.00174025, 0.00164909, 0.00156174,
67 0.00147811, 0.00139810, 0.00132161, 0.00124854, 0.00117880, 0.00111228, 0.00104887, 0.00098848, 0.00093100,
68 0.00087634, 0.00082438,
69 };
70 static float g_mg[257] = {
71 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
72 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
73 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
74 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
75 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
76 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
77 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
78 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
79 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
80 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
81 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
82 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
83 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
84 0.00000004, 0.00000055, 0.00000627, 0.00005539, 0.00037863, 0.00200318, 0.00820201, 0.02599027, 0.06373712,
85 0.12096648, 0.17767593, 0.20196826, 0.17767593, 0.12096648, 0.06373712, 0.02599027, 0.00820201, 0.00200318,
86 0.00037863, 0.00005539, 0.00000627, 0.00000055, 0.00000004, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
87 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
88 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
89 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
90 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
91 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
92 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
93 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
94 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
95 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
96 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
97 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
98 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
99 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.00000000,
100 };
101
102 // other functions
103 static void BeatTrack_dofft(BeatTrack* unit, uint32);
104 static void complexdf(BeatTrack* unit);
105 static void finaldecision(BeatTrack* unit);
106
107 // amortisation
108 static void autocorr(BeatTrack* unit, int j);
109 static void beatperiod(BeatTrack* unit, int j, int whichm);
110 static float findtor(BeatTrack* unit);
111
112 // as many amortisation steps as tor
113 static void findphase(BeatTrack* unit, int j, int gaussflag, int predicted);
114
115 static int detectperiodchange(BeatTrack* unit);
116 static void findmeter(BeatTrack* unit);
117 static void setupphaseexpectation(BeatTrack* unit); // create Gaussian focussed matrix for phase
118
119
BeatTrack_Ctor(BeatTrack * unit)120 void BeatTrack_Ctor(BeatTrack* unit) {
121 ///////
122 // check sampling rate and establish multipliers on estimates and FFT window size
123 // down sampling by factor of two automatic
124
125 unit->m_srate = unit->mWorld->mFullRate.mSampleRate;
126
127 // if sample rate is 88200 or 96000, assume taking double size FFT to start with
128 if (unit->m_srate > (44100.0 * 1.5))
129 unit->m_srate = unit->m_srate * 0.5;
130
131 unit->m_srateconversion = unit->m_srate / 44100.0;
132
133 // assumes base of 1024 FFT
134 unit->m_frameperiod =
135 (FRAMEPERIOD / unit->m_srateconversion); // in seconds //(int) ((FRAMEPERIOD/unit->m_srateconversion) +0.5);
136
137 printf("srate %f conversion factor %f frame period %f \n", unit->m_srate, unit->m_srateconversion,
138 unit->m_frameperiod);
139
140 unit->m_prevmag = (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
141 unit->m_prevphase = (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
142 unit->m_predict = (float*)RTAlloc(unit->mWorld, NOVER2 * sizeof(float));
143
144 ////////time positions//////////
145 unit->m_frame = 1; // don't decide immediately, wait for maximum period!
146
147 /////////df////////
148 unit->m_dfcounter = DFSTORE - 1;
149 // random uncorrelated noise df store for initialisation
150 // RGen& rgen = *unit->mParent->mRGen;
151
152 // don't want this noise, want consistent starting point!
153 for (int j = 0; j < DFSTORE; ++j) {
154 unit->m_df[j] = 0.0; //(2*rgen.frand() - 1.0);
155 }
156
157 unit->m_dfmemorycounter = 14;
158
159 Clear(15, unit->m_dfmemory);
160
161 /////////tempo assess///////////
162 unit->m_currtempo = 2;
163
164 ////////phase assess///////////
165
166 unit->m_currphase = 0.0;
167
168 unit->m_phase = 0.0;
169
170 // default of 2bps
171 unit->m_phaseperblock = ((float)unit->mWorld->mFullRate.mBufLength * 2) / ((float)unit->mWorld->mSampleRate);
172
173 unit->m_outputphase = unit->m_phase;
174 unit->m_outputtempo = unit->m_currtempo;
175 unit->m_outputphaseperblock = unit->m_phaseperblock;
176
177 unit->halftrig = 0;
178 unit->q1trig = 0;
179 unit->q2trig = 0;
180
181 // amortisation and states
182 unit->m_amortisationstate = 0; // off
183 unit->m_stateflag = 0;
184 unit->m_timesig = 4;
185 unit->m_flagstep = 0;
186
187 unit->mCalcFunc = (UnitCalcFunc)&BeatTrack_next;
188
189 // initialize outputs
190 ZOUT0(0) = 0.0;
191 ZOUT0(1) = 0.0;
192 ZOUT0(2) = 0.0;
193 ZOUT0(3) = unit->m_outputtempo;
194 }
195
196
BeatTrack_Dtor(BeatTrack * unit)197 void BeatTrack_Dtor(BeatTrack* unit) {
198 RTFree(unit->mWorld, unit->m_prevmag);
199 RTFree(unit->mWorld, unit->m_prevphase);
200 RTFree(unit->mWorld, unit->m_predict);
201 }
202
203
BeatTrack_next(BeatTrack * unit,int wrongNumSamples)204 void BeatTrack_next(BeatTrack* unit, int wrongNumSamples) {
205 // float *in = IN(0);
206
207 // printf("%d \n",wrongNumSamples);
208 // int numSamples = unit->mWorld->mFullRate.mBufLength;
209
210 // conditions in reverse order to avoid immediate spillover
211 // printf("state %d \n",unit->m_amortisationstate);
212
213 // keeps incrementing but will be reset with each calculation run
214 unit->m_amortisationsteps = unit->m_amortisationsteps + 1;
215
216 // if state nonzero do something
217 switch (unit->m_amortisationstate) {
218 case 0:
219 break; // do nothing case
220 case 1: // calculate acf
221 autocorr(unit, unit->m_amortcount);
222
223 unit->m_amortcount = unit->m_amortcount + 1;
224
225 if (unit->m_amortcount == unit->m_amortlength) {
226 unit->m_amortisationstate = 2;
227 unit->m_amortlength = 128;
228 unit->m_amortcount = 0;
229
230 unit->m_bestcolumn = 0;
231 unit->m_besttorsum = -1000.0;
232 }
233
234 break;
235 case 2: // periodp
236 beatperiod(unit, unit->m_amortcount, 0);
237
238 unit->m_amortcount = unit->m_amortcount + 1;
239
240 if (unit->m_amortcount == unit->m_amortlength) {
241 unit->m_periodp = findtor(unit);
242
243 if (unit->m_stateflag == 1) {
244 unit->m_amortisationstate = 3;
245 unit->m_amortlength = 128;
246 unit->m_amortcount = 0;
247
248 unit->m_bestcolumn = 0;
249 unit->m_besttorsum = -1000.0;
250
251 } else {
252 unit->m_periodg = -1000; // will always trigger initially
253 unit->m_amortisationstate = 4;
254 }
255 }
256
257 break;
258 case 3: // periodg
259 beatperiod(unit, unit->m_amortcount, 1);
260 unit->m_amortcount = unit->m_amortcount + 1;
261
262 if (unit->m_amortcount == unit->m_amortlength) {
263 unit->m_periodg = findtor(unit);
264
265 unit->m_amortisationstate = 4;
266 }
267
268 break;
269 case 4: // stepdetect/constdetect
270
271 if (detectperiodchange(unit)) {
272 unit->m_amortisationstate = 5;
273 unit->m_amortlength = 128;
274 unit->m_amortcount = 0;
275
276 unit->m_bestcolumn = 0;
277 unit->m_besttorsum = -1000.0;
278
279 unit->m_stateflag = 1;
280 findmeter(unit);
281
282 // set up Gaussian weighting centred on periodp
283 int startindex = 128 - ((int)(unit->m_periodp + 0.5));
284
285 float* mg = unit->m_mg;
286
287 for (int ii = 0; ii < 128; ++ii) {
288 mg[ii] = g_mg[startindex + ii];
289 }
290
291 } else {
292 if (unit->m_stateflag == 1)
293 unit->m_tor = unit->m_periodg;
294 else
295 unit->m_tor = unit->m_periodp;
296
297 unit->m_torround = int(unit->m_tor + 0.5);
298
299 unit->m_amortisationstate = 7;
300 unit->m_amortlength = unit->m_torround;
301 unit->m_amortcount = 0;
302 }
303
304 break;
305
306 case 5: // redo periodg calculation
307 beatperiod(unit, unit->m_amortcount, 1);
308 unit->m_amortcount = unit->m_amortcount + 1;
309
310 if (unit->m_amortcount == unit->m_amortlength) {
311 unit->m_periodg = findtor(unit);
312
313 unit->m_tor = unit->m_periodg;
314 unit->m_torround = int(unit->m_tor + 0.5f);
315
316 unit->m_amortisationstate = 6;
317 unit->m_amortlength = unit->m_torround;
318 unit->m_amortcount = 0;
319
320 setupphaseexpectation(unit);
321
322 // don't need to reset change flag since it isn't stored
323 }
324
325 break;
326 case 6: // flat phase calc after move to context, avoids bias
327 findphase(unit, unit->m_amortcount, 0, 0);
328 unit->m_amortcount = unit->m_amortcount + 1;
329
330 if (unit->m_amortcount == unit->m_amortlength) {
331 unit->m_amortisationstate = 8; // final state
332 }
333
334 break;
335
336 case 7: // phase calc with possible gaussian narrowing of the allowed phases
337
338 findphase(unit, unit->m_amortcount, unit->m_stateflag, (int)(unit->m_currphase * unit->m_torround + 0.5f));
339 unit->m_amortcount = unit->m_amortcount + 1;
340
341 if (unit->m_amortcount == unit->m_amortlength) {
342 unit->m_amortisationstate = 8; // final state
343 }
344
345 break;
346 case 8:
347
348 finaldecision(unit);
349 unit->m_amortisationstate = 0;
350 break;
351
352 default:
353 break;
354 }
355
356
357 // MUST CHECK IF INCIDENT FFT IS >1, if so update buffer with appropriate coefficients
358
359 float fbufnum = ZIN0(0);
360
361 // next FFT bufffer ready, update
362 // assuming at this point that buffer precalculated for any resampling
363 if (!(fbufnum < 0)) {
364 unit->m_frame = unit->m_frame + 1;
365 BeatTrack_dofft(unit, (uint32)fbufnum);
366 }
367
368
369 // test if impulse to output
370 unit->m_phase += unit->m_phaseperblock;
371
372 // if not locked, update output phase from model phase, else keep a separate output phase
373
374 float lock = ZIN0(1);
375 // printf("lock %f \n",lock);
376
377 if (lock < 0.5f) {
378 unit->m_outputphase = unit->m_phase;
379 unit->m_outputtempo = unit->m_currtempo;
380 unit->m_outputphaseperblock = unit->m_phaseperblock;
381 } else
382 unit->m_outputphase += unit->m_outputphaseperblock;
383
384 if (unit->m_phase >= 1.f)
385 unit->m_phase -= 1.f;
386
387 // 0 is beat, 1 is quaver, 2 is semiquaver, 3 is actual current tempo in bps
388 // so no audio accuracy with beats, just asap, may as well be control rate
389 ZOUT0(0) = 0.0;
390 ZOUT0(1) = 0.0;
391 ZOUT0(2) = 0.0;
392 ZOUT0(3) = unit->m_outputtempo; //*0.016666667;
393
394 // output beat
395 if (unit->m_outputphase >= 1.f) {
396 // printf("beat \n");
397
398 unit->m_outputphase -= 1.f;
399 ZOUT0(0) = 1.0;
400 ZOUT0(1) = 1.0;
401 ZOUT0(2) = 1.0;
402 unit->halftrig = 0;
403 unit->q1trig = 0;
404 unit->q2trig = 0;
405 }
406
407 if (unit->m_outputphase >= 0.5f && unit->halftrig == 0) {
408 ZOUT0(1) = 1.0;
409 ZOUT0(2) = 1.0;
410 unit->halftrig = 1;
411 }
412
413 if (unit->m_outputphase >= 0.25f && unit->q1trig == 0) {
414 ZOUT0(2) = 1.0;
415 unit->q1trig = 1;
416 }
417
418 if (unit->m_outputphase >= 0.75f && unit->q2trig == 0) {
419 ZOUT0(2) = 1.0;
420 unit->q2trig = 1;
421 }
422 }
423
424
425 //
426
427 // calculation function once FFT data ready
BeatTrack_dofft(BeatTrack * unit,uint32 ibufnum)428 void BeatTrack_dofft(BeatTrack* unit, uint32 ibufnum) {
429 World* world = unit->mWorld;
430 SndBuf* buf;
431 if (ibufnum >= world->mNumSndBufs) {
432 int localBufNum = ibufnum - world->mNumSndBufs;
433 Graph* parent = unit->mParent;
434 if (localBufNum <= parent->localBufNum) {
435 buf = parent->mLocalSndBufs + localBufNum;
436 } else {
437 buf = world->mSndBufs;
438 }
439 } else {
440 buf = world->mSndBufs + ibufnum;
441 }
442 LOCK_SNDBUF(buf);
443 // int numbins = buf->samples - 2 >> 1;
444
445 unit->m_FFTBuf = buf->data; // just assign it!
446
447 // transfer data to fftbuf in the format expected by this plugin
448
449 // ideally, should do this part separate to plug-in as well, so can compare different detection functions;
450 // also, can run multiple in parallel with own autocorrelations; committee? Committee.ar(period1, phase1, period2,
451 // phase2, period3, phase3)... chooses predominant estimate? feature detection function
452 complexdf(unit);
453
454 if (unit->m_frame % SKIP == 0) {
455 // printf("amortisation time \n");
456
457 // amortisation- 8 control periods in a frame
458 // have 2000 calcs to do, split over 100 control periods = 6400 samples, ie one tempo per control period
459
460 unit->m_bestcolumn = 0;
461 unit->m_besttorsum = -1000.0;
462
463 unit->m_bestphasescore = -1000.0;
464 unit->m_bestphase = 0;
465
466 // state 0 is do nothing
467 unit->m_amortisationstate = 1;
468 unit->m_amortcount = 0;
469 unit->m_amortlength = 128;
470 unit->m_amortisationsteps = 0;
471
472 // fix time reference for calculations, so it doesn't update during the amortisation- this is the beginning of
473 // the df frame
474 unit->m_storedfcounter = unit->m_dfcounter + DFSTORE - DFFRAMELENGTH;
475
476 // ref for phase calculations
477 unit->m_storedfcounterend = unit->m_dfcounter;
478
479 // unit->m_fftstoreposhold= unit->m_fftstorepos;
480
481 unit->m_currphase = unit->m_phase;
482 }
483 }
484
485
autocorr(BeatTrack * unit,int j)486 void autocorr(BeatTrack* unit, int j) {
487 int baseframe = unit->m_storedfcounter + DFSTORE;
488 float* df = unit->m_df;
489 float* acf = unit->m_acf;
490
491 // work out four lags each time
492 for (int k = 0; k < 4; ++k) {
493 int lag = 4 * j + k;
494
495 int correction = abs(lag - DFFRAMELENGTH);
496
497 float sum = 0.0;
498
499 for (int i = lag; i < DFFRAMELENGTH; ++i) {
500 float val1 = df[(i + baseframe) % DFSTORE];
501 float val2 = df[(i + baseframe - lag) % DFSTORE];
502
503 sum += val1 * val2;
504 }
505
506 acf[lag] = sum * correction;
507 }
508 }
509
510
511 // timesig 4 has one more sum term
512 // indices as MATLAB but need to correct maxinds to be in range of tested, not in global range
findtor(BeatTrack * unit)513 float findtor(BeatTrack* unit) {
514 float maxval, val;
515 int ind2, ind3, ind4;
516
517 // put into MATLAB indexing, from 1 to 512
518 int ind = unit->m_bestcolumn + 1;
519
520 float* acf = unit->m_acf - 1;
521
522 ind2 = 0;
523 maxval = -1000;
524
525 for (int i = 2 * ind - 1; i <= (2 * ind + 1); ++i) {
526 val = acf[i];
527
528 if (val > maxval) {
529 maxval = val;
530 ind2 = i - (2 * ind - 1) + 1;
531 }
532 }
533
534 //[val2,ind2] = max(acf(2*ind-1:2*ind+1));
535 ind2 = ind2 + 2 * (ind + 1) - 2;
536
537 ind3 = 0;
538 maxval = -1000;
539
540 for (int i = 3 * ind - 2; i <= (3 * ind + 2); ++i) {
541 val = acf[i];
542
543 if (val > maxval) {
544 maxval = val;
545 ind3 = i - (3 * ind - 2) + 1;
546 }
547 }
548
549 //[val3,ind3] = max(acf(3*ind-2:3*ind+2));
550 ind3 = ind3 + 3 * ind - 4;
551
552 float period;
553
554 if (unit->m_timesig == 4) {
555 ind4 = 0;
556 maxval = -1000;
557
558 for (int i = 4 * ind - 3; i <= 4 * ind + 3; ++i) {
559 val = acf[i];
560
561 if (val > maxval) {
562 maxval = val;
563 ind4 = i - (4 * ind - 3) + 1;
564 }
565 }
566
567 //[val4,ind4] = max(acf(4*ind-3:4*ind+3));
568 ind4 = ind4 + 4 * ind - 9;
569
570 period = (ind + ind2 * 0.5 + ind3 / 3.f + ind4 * 0.25) * 0.25;
571
572 } else
573
574 period = (ind + ind2 * 0.5 + ind3 / 3.f) * 0.3333333;
575
576
577 // printf("period %f ind %d ind2 %d ind3 %d ind4 %d \n",period, ind,ind2,ind3,ind4);
578
579 // unit->m_tor=period;
580 // unit->m_torround= int(period+0.5);
581 //
582
583 return period;
584 }
585
586
587 // 128 calculation calls for multiplying M and acf, calculates M as it goes apart from precalculated Gaussian or Raleigh
588 // distribution
beatperiod(BeatTrack * unit,int j,int whichm)589 void beatperiod(BeatTrack* unit, int j, int whichm) {
590 float* acf = unit->m_acf;
591
592 // int startindex= 512*j;
593 // int endindex=startindex+512;
594
595 float sum = 0.0;
596
597 // unit->m_timesig harmonics
598 for (int i = 1; i <= (unit->m_timesig); ++i) {
599 int num = 2 * i - 1;
600 float wt = 1.0 / (float)num;
601
602 for (int k = 0; k < num; ++k) {
603 int pos = k + (i * j);
604
605 // m[startindex+pos]
606 if (pos < 512)
607 sum += acf[pos] * wt;
608 }
609 }
610
611 // assumes Mg appropriately rotated already
612 float* m;
613
614 if (whichm)
615 m = g_m; // Gaussian weighted context model
616 else
617 m = unit->m_mg; // general model even weighting
618
619 sum = sum * m[j];
620
621 if (sum > unit->m_besttorsum) {
622 unit->m_besttorsum = sum;
623 unit->m_bestcolumn = j;
624 }
625 }
626
627
628 // j out of unit->m_torround
629 // differs to Davies original in that weight the most recent events more- want minimum reaction time
findphase(BeatTrack * unit,int j,int gaussflag,int predicted)630 void findphase(BeatTrack* unit, int j, int gaussflag, int predicted) {
631 float* df = unit->m_df;
632
633 int period = unit->m_torround;
634 int baseframe = unit->m_storedfcounterend + DFSTORE;
635
636 int numfit = -1;
637
638 if (period != 0)
639 // round down
640 numfit = (int)(DFFRAMELENGTH / period) - 1;
641
642 // testing backwards from the baseframe, weighting goes down as 1/k
643 float sum = 0.0;
644
645 for (int k = 0; k < numfit; ++k) {
646 // j is phase to test
647 int location = (baseframe - (period * k) - j) % DFSTORE;
648
649 sum += df[location] / ((float)(k + 1));
650 }
651
652 // Gaussian focus weighting if desired
653 if (gaussflag) {
654 // difference of predicted from j, min distance within period
655 int diff = sc_min(abs(predicted - j), abs(period - predicted + j));
656
657 sum *= unit->m_phaseweights[diff];
658 }
659
660 if (sum > unit->m_bestphasescore) {
661 unit->m_bestphasescore = sum;
662 unit->m_bestphase = j;
663 }
664 }
665
666 //, int predicted
setupphaseexpectation(BeatTrack * unit)667 void setupphaseexpectation(BeatTrack* unit) // create Gaussian focussed matrix for phase
668 {
669 float* wts = unit->m_phaseweights;
670
671 float sigma = unit->m_torround * 0.25f;
672 // float mu=period;
673
674 float mult = 1.0 / (2.5066283 * sigma);
675 float mult2 = 1.0 / (2.0 * sigma * sigma);
676
677 // unit->m_torround
678 for (int i = 0; i < 128; ++i) {
679 wts[i] = mult * (exp(-(i * i) * mult2));
680 }
681 }
682
683
684 // why force a countdown each time? Why not keep a continuous buffer of previous periodp, periodg?
detectperiodchange(BeatTrack * unit)685 int detectperiodchange(BeatTrack* unit) {
686 // stepthresh = 3.9017;
687
688 if (unit->m_flagstep == 0) {
689 if (fabs(unit->m_periodg - unit->m_periodp) > 3.9017f) {
690 unit->m_flagstep = 3;
691 }
692
693 } else {
694 unit->m_flagstep = unit->m_flagstep - 1;
695 }
696
697 if (unit->m_flagstep) {
698 unit->m_prevperiodp[unit->m_flagstep - 1] = unit->m_periodp;
699 }
700
701 if (unit->m_flagstep == 1) {
702 unit->m_flagstep = 0;
703
704 if (fabs(2 * unit->m_prevperiodp[0] - unit->m_prevperiodp[1] - unit->m_prevperiodp[2]) < 7.8034f) //(2*3.9017)
705 return 1;
706 }
707
708 return 0;
709 }
710
711 // add test
findmeter(BeatTrack * unit)712 void findmeter(BeatTrack* unit) {
713 // int i;
714
715 // float * acf= unit->m_acf;
716
717 // float * acf= unit->m_acf-1;
718 //
719 //
720 // int period = ((int)(unit->m_periodp+0.5));
721 //
722 // float three_energy=0.0;
723 // float four_energy=0.0;
724 //
725 // for(i=(3*period-2);i<(3*period+3);++i)
726 // three_energy += acf[i];
727 //
728 // for(i=(4*period-2);i<(4*period+3);++i)
729 // four_energy += acf[i];
730 //
731 // if((6*period+2)<512) {
732 //
733 // for(i=(6*period-2);i<(6*period+3);++i)
734 // three_energy += acf[i];
735 //
736 // for(i=(2*period-2);i<(2*period+3);++i)
737 // four_energy += acf[i];
738 // }
739 //
740 // if (three_energy > four_energy)
741 // unit->m_timesig = 3;
742 // else
743
744 // worked better in evaluation without any 3/4 at all!
745 unit->m_timesig = 4;
746
747 // printf("time sig %d \n",unit->m_timesig);
748 }
749
750
751 // period is unit->m_tor, phase is unit->m_bestphase
752 // float m_tor; int m_torround;
finaldecision(BeatTrack * unit)753 void finaldecision(BeatTrack* unit) {
754 // int i,j;
755
756 unit->m_currtempo = 1.0 / (unit->m_tor * unit->m_frameperiod);
757
758 unit->m_phaseperblock =
759 ((float)unit->mWorld->mFullRate.mBufLength * (unit->m_currtempo)) / ((float)unit->mWorld->mSampleRate);
760
761 // printf("SAMPLErate %f %f %f", unit->mWorld->mSampleRate,unit->m_phaseperblock,unit->m_currtempo);
762
763 // unit->m_amortisationstate control periods worth = 512/64 = 8
764 // float frameselapsed= 0.125*unit->m_amortisationstate;
765 // float timeelapsed= frameselapsed*unit->m_frameperiod;
766
767 float timeelapsed = ((float)(unit->m_amortisationsteps) * (unit->mWorld->mFullRate.mBufLength)
768 / ((float)unit->mWorld->mSampleRate));
769
770 timeelapsed += 7 * unit->m_frameperiod; // compensation for detection function being delayed by 7 frames
771
772 float phaseelapsed = timeelapsed * (unit->m_currtempo);
773
774 float phasebeforeamort = ((float)unit->m_bestphase / unit->m_torround);
775
776 // add phase to compensate for ELAPSEDTIME
777 unit->m_currphase = unit->m_phase = fmod(phasebeforeamort + phaseelapsed, (float)1.0);
778 }
779
780
781 // Now the format is standardised for the SC FFT UGen as
782 // dc, nyquist and then real/imag pairs for each bin going up successively in frequency
783
complexdf(BeatTrack * unit)784 void complexdf(BeatTrack* unit) {
785 float* fftbuf = unit->m_FFTBuf;
786
787 float* prevmag = unit->m_prevmag;
788 float* prevphase = unit->m_prevphase;
789 float* predict = unit->m_predict;
790
791 float sum = 0.0;
792
793 // printf("complex df time \n");
794
795 // sum bins 2 to 256
796 for (int k = 1; k < NOVER2; ++k) {
797 // Change to fftw
798 int index = 2 * k; // k; //2*k;
799
800 float real = fftbuf[index];
801 // N=1024 conventionally here
802 float imag = fftbuf[index + 1]; // fftbuf[N-index];
803
804 float mag = sqrt(
805 real * real
806 + imag * imag); // was 0.5*sqrt(real*real+ imag*imag); reduce by factor of 2 because of altivec side effect
807 float qmag = prevmag[k];
808
809 prevmag[k] = mag;
810
811 float phase = atan2(imag, real);
812 float oldphase = predict[k];
813
814 predict[k] = 2 * phase - prevphase[k];
815 prevphase[k] = phase;
816
817 float phasediff = phase - oldphase;
818
819 // if(k==2) printf("%f %f\n",phase, phasediff);
820
821 // tables for cos/sin/sqrt speeds up? sqrt(1-c*c) slower than sin
822
823 float realpart = (qmag - (mag * cos(phasediff)));
824 float imagpart = (mag * sin(phasediff)); // no need for negative
825
826 float detect = sqrt(realpart * realpart + imagpart * imagpart);
827
828 // detect is always positive
829 // if(k==1)
830 sum += detect; //(fmod(phase+(16*pi),twopi)); //detect;
831
832 // if(k==1) sum+=mag;
833 }
834
835
836 // smoothing and peak picking operation, delay of 8 frames, must be taken account of in final phase correction
837
838 unit->m_dfmemorycounter = (unit->m_dfmemorycounter + 1) % 15;
839 unit->m_dfmemory[unit->m_dfmemorycounter] = sum; // divide by num of bands to get a dB answer
840
841 float rating = 0.0;
842
843 float* dfmemory = unit->m_dfmemory;
844
845 int refpos = unit->m_dfmemorycounter + 15;
846 int centrepos = (refpos - 7) % 15;
847 float centreval = dfmemory[centrepos];
848
849 for (int k = 0; k < 15; ++k) {
850 int pos = (refpos - k) % 15;
851
852 float nextval = centreval - dfmemory[pos];
853
854 if (nextval < 0.0)
855 nextval = nextval * 10;
856
857 rating += nextval;
858 }
859
860 if (rating < 0.0)
861 rating = 0.0;
862
863 // increment first so this frame is unit->m_loudnesscounterdfcounter
864 unit->m_dfcounter = (unit->m_dfcounter + 1) % DFSTORE;
865
866 unit->m_df[unit->m_dfcounter] = rating * 0.1f; // sum //divide by num of bands to get a dB answer
867 }
868