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 // Feature (Onset) Detection implemented by sick lincoln for sc3
22 // Jensen,K. & Andersen, T. H. (2003). Real-time Beat Estimation Using Feature Extraction.
23 // In Proceedings of the Computer Music Modeling and RetrievalSymposium, Lecture Notes in Computer Science. Springer
24 // Verlag. Hainsworth, S. (2003) Techniques for the Automated Analysis of Musical Audio. PhD, university of cambridge
25 // engineering dept.
26
27 // possible to make a Goto style Detector for a given band and with history of two samples-
28 // should do separately as PV_GotoBandTrack
29
30 // next perhaps Duxbury et al/ Mauri et al different conception of high frequency content with ratio of changes
31
32 #include "FFT_UGens.h"
33
34 struct PV_OnsetDetectionBase : public Unit {
35 float* m_prevframe;
36 int m_numbins;
37 int m_waiting, m_waitSamp, m_waitLen;
38 };
39
40 // FFT onset detector combining 4 advised features from Jensen/Andersen
41 struct PV_JensenAndersen : public PV_OnsetDetectionBase {
42 float m_hfc, m_hfe, m_sc, m_sf;
43 int m_fourkindex;
44 };
45
46
47 // FFT onset detector combining 2 advised features from Hainsworth PhD
48 struct PV_HainsworthFoote : public PV_OnsetDetectionBase {
49 float m_prevNorm;
50 int m_5kindex, m_30Hzindex;
51 };
52
53 // for time domain onset detection/RMS
54 struct RunningSum : public Unit {
55 int msamp, mcount;
56 float msum, msum2;
57 // float mmeanmult;
58 float* msquares;
59 };
60
61 extern "C" {
62 void PV_OnsetDetectionBase_Ctor(PV_OnsetDetectionBase* unit);
63 void PV_OnsetDetectionBase_Dtor(PV_OnsetDetectionBase* unit);
64
65 void PV_JensenAndersen_Ctor(PV_JensenAndersen* unit);
66 void PV_JensenAndersen_Dtor(PV_JensenAndersen* unit);
67 void PV_JensenAndersen_next(PV_JensenAndersen* unit, int inNumSamples);
68
69 void PV_HainsworthFoote_Ctor(PV_HainsworthFoote* unit);
70 void PV_HainsworthFoote_Dtor(PV_HainsworthFoote* unit);
71 void PV_HainsworthFoote_next(PV_HainsworthFoote* unit, int inNumSamples);
72
73 void RunningSum_next_k(RunningSum* unit, int inNumSamples);
74 void RunningSum_Ctor(RunningSum* unit);
75 void RunningSum_Dtor(RunningSum* unit);
76 }
77
78 #define PV_FEAT_GET_BUF_UNLOCKED \
79 uint32 ibufnum = (uint32)fbufnum; \
80 int bufOK = 1; \
81 World* world = unit->mWorld; \
82 SndBuf* buf; \
83 if (ibufnum >= world->mNumSndBufs) { \
84 int localBufNum = ibufnum - world->mNumSndBufs; \
85 Graph* parent = unit->mParent; \
86 if (localBufNum <= parent->localBufNum) { \
87 buf = parent->mLocalSndBufs + localBufNum; \
88 } else { \
89 bufOK = 0; \
90 buf = world->mSndBufs; \
91 if (unit->mWorld->mVerbosity > -1) { \
92 Print("FFT Ctor error: Buffer number overrun: %i\n", ibufnum); \
93 } \
94 } \
95 } else { \
96 buf = world->mSndBufs + ibufnum; \
97 } \
98 int numbins = (buf->samples - 2) >> 1; \
99 if (!buf->data) { \
100 if (unit->mWorld->mVerbosity > -1) { \
101 Print("FFT Ctor error: Buffer %i not initialised.\n", ibufnum); \
102 } \
103 bufOK = 0; \
104 }
105
106 #define PV_FEAT_GET_BUF \
107 PV_FEAT_GET_BUF_UNLOCKED \
108 LOCK_SNDBUF(buf);
109
110
PV_OnsetDetectionBase_Ctor(PV_OnsetDetectionBase * unit)111 void PV_OnsetDetectionBase_Ctor(PV_OnsetDetectionBase* unit) {
112 float fbufnum = ZIN0(0);
113
114 PV_FEAT_GET_BUF_UNLOCKED
115
116 unit->m_numbins = numbins;
117 int insize = unit->m_numbins * sizeof(float);
118
119 if (bufOK) {
120 unit->m_prevframe = (float*)RTAlloc(unit->mWorld, insize);
121 memset(unit->m_prevframe, 0, insize);
122 }
123
124 unit->m_waiting = 0;
125 unit->m_waitSamp = 0;
126 unit->m_waitLen = 0;
127 ClearUnitOutputs(unit, 1);
128 }
129
PV_OnsetDetectionBase_Dtor(PV_OnsetDetectionBase * unit)130 void PV_OnsetDetectionBase_Dtor(PV_OnsetDetectionBase* unit) {
131 if (unit->m_prevframe)
132 RTFree(unit->mWorld, unit->m_prevframe);
133 }
134
135
PV_JensenAndersen_Ctor(PV_JensenAndersen * unit)136 void PV_JensenAndersen_Ctor(PV_JensenAndersen* unit) {
137 PV_OnsetDetectionBase_Ctor(unit);
138
139 unit->m_hfc = 0.0;
140 unit->m_hfe = 0.0;
141 unit->m_sf = 0.0;
142 unit->m_sc = 0.0;
143
144 unit->m_fourkindex = (int)(4000.0 / (unit->mWorld->mSampleRate)) * (unit->m_numbins);
145
146 SETCALC(PV_JensenAndersen_next);
147 }
148
PV_JensenAndersen_Dtor(PV_JensenAndersen * unit)149 void PV_JensenAndersen_Dtor(PV_JensenAndersen* unit) { PV_OnsetDetectionBase_Dtor(unit); }
150
151
PV_JensenAndersen_next(PV_JensenAndersen * unit,int inNumSamples)152 void PV_JensenAndersen_next(PV_JensenAndersen* unit, int inNumSamples) {
153 float outval = 0.0;
154 float fbufnum = ZIN0(0);
155
156 if (unit->m_waiting == 1) {
157 unit->m_waitSamp += inNumSamples;
158 if (unit->m_waitSamp >= unit->m_waitLen)
159 unit->m_waiting = 0;
160 }
161
162 if (!(fbufnum < 0.f))
163 // if buffer ready to process
164 {
165 PV_FEAT_GET_BUF
166
167 SCPolarBuf* p = ToPolarApx(buf);
168
169 // four spectral features useful for onset detection according to Jensen/Andersen
170
171 float magsum = 0.0, magsumk = 0.0, magsumkk = 0.0, sfsum = 0.0, hfesum = 0.0;
172
173 float* q = unit->m_prevframe;
174
175 int k4 = unit->m_fourkindex;
176
177 // ignores dc, nyquist
178 for (int i = 0; i < numbins; ++i) {
179 float mag = ((p->bin[i]).mag);
180 int k = i + 1;
181 float qmag = q[i];
182 magsum += mag;
183 magsumk += k * mag;
184 magsumkk += k * k * mag;
185 sfsum += fabs(mag - (qmag));
186 if (i > k4)
187 hfesum += mag;
188 }
189
190 float binmult = 1.f / numbins;
191 // normalise
192 float sc = (magsumk / magsum) * binmult;
193 float hfe = hfesum * binmult;
194 float hfc = magsumkk * binmult * binmult * binmult;
195 float sf = sfsum * binmult;
196
197 // printf("sc %f hfe %f hfc %f sf %f \n",sc, hfe, hfc, sf);
198
199 // if(sc<0.0) sc=0.0;
200 // if(hfe<0.0) hfe=0.0;
201 // if(hfc<0.0) hfc=0.0;
202 // if(sf<0.0) sf=0.0;
203
204 // ratio of current to previous frame perhaps better indicator than first derivative difference
205 float scdiff = sc - (unit->m_sc);
206 float hfediff = hfe - (unit->m_hfe);
207 float hfcdiff = hfc - (unit->m_hfc);
208 float sfdiff = sf - (unit->m_sf);
209
210 // store as old frame values for taking difference
211 unit->m_sc = sc;
212 unit->m_hfe = hfe;
213 unit->m_hfc = hfc;
214 unit->m_sf = sf;
215
216 // printf("sc %f hfe %f hfc %f sf %f \n",sc, hfe, hfc, sf);
217 // printf("sc %f hfe %f hfc %f sf %f \n",scdiff, hfediff, hfcdiff, sfdiff);
218
219 // does this trigger?
220 // may need to take derivatives across previous frames by storing old values
221
222 float sum = (ZIN0(1) * scdiff) + (ZIN0(2) * hfediff) + (ZIN0(3) * hfcdiff) + (ZIN0(4) * sfdiff);
223
224 // printf("sum %f thresh %f \n",sum, ZIN0(7));
225
226 // if over threshold, may also impose a wait here
227 if (sum > ZIN0(5) && (unit->m_waiting == 0)) { // printf("bang! \n");
228 outval = 1.0;
229 unit->m_waiting = 1;
230 unit->m_waitSamp = inNumSamples;
231 unit->m_waitLen = (int)(ZIN0(6) * (world->mSampleRate));
232 }
233
234 // take copy of this frame's magnitudes as prevframe
235
236 for (int i = 0; i < numbins; ++i)
237 q[i] = p->bin[i].mag;
238 }
239
240 Fill(inNumSamples, &ZOUT0(0), outval);
241 }
242
243
PV_HainsworthFoote_Ctor(PV_HainsworthFoote * unit)244 void PV_HainsworthFoote_Ctor(PV_HainsworthFoote* unit) {
245 PV_OnsetDetectionBase_Ctor(unit);
246
247 World* world = unit->mWorld;
248
249 unit->m_5kindex = (int)((5000.0 / (world->mSampleRate)) * (unit->m_numbins));
250 unit->m_30Hzindex = (int)((30.0 / (world->mSampleRate)) * (unit->m_numbins));
251
252 unit->m_prevNorm = 1.0;
253
254 // unit->m_5kindex, unit->m_30Hzindex,
255 // printf("numbins %d sr %d \n", unit->m_numbins, world->mSampleRate);
256 // printf("test %d sr %f 5k %d 30Hz %d\n", unit->m_numbins, world->mSampleRate, unit->m_5kindex, unit->m_30Hzindex);
257
258 SETCALC(PV_HainsworthFoote_next);
259 }
260
PV_HainsworthFoote_Dtor(PV_HainsworthFoote * unit)261 void PV_HainsworthFoote_Dtor(PV_HainsworthFoote* unit) { PV_OnsetDetectionBase_Dtor(unit); }
262
263 static const float lmult = 1.442695040889; // loge(2) reciprocal
264
PV_HainsworthFoote_next(PV_HainsworthFoote * unit,int inNumSamples)265 void PV_HainsworthFoote_next(PV_HainsworthFoote* unit, int inNumSamples) {
266 float outval = 0.0;
267 float fbufnum = ZIN0(0);
268
269 if (unit->m_waiting == 1) {
270 unit->m_waitSamp += inNumSamples;
271 if (unit->m_waitSamp >= unit->m_waitLen) {
272 unit->m_waiting = 0;
273 }
274 }
275
276 if (!(fbufnum < 0.f))
277 // if buffer ready to process
278 {
279 PV_FEAT_GET_BUF
280
281 SCPolarBuf* p = ToPolarApx(buf);
282
283 float dnk, prevmag, mkl = 0.0, footesum = 0.0, norm = 0.0;
284
285 float* q = unit->m_prevframe;
286
287 int k5 = unit->m_5kindex;
288 int h30 = unit->m_30Hzindex;
289
290 for (int i = 0; i < numbins; ++i) {
291 float mag = ((p->bin[i]).mag);
292 float qmag = q[i];
293
294 if (i >= h30 && i < k5) {
295 prevmag = qmag;
296 // avoid divide by zero
297 if (prevmag < 0.0001)
298 prevmag = 0.0001;
299
300 // no log2 in maths library, so use log2(x)= log(x)/log(2) where log is to base e
301 // could just use log and ignore scale factor but hey let's stay accurate to the source for now
302 dnk = log(mag / prevmag) * lmult;
303
304 if (dnk > 0.0)
305 mkl += dnk;
306 }
307
308 norm += mag * mag;
309 footesum += mag * qmag;
310 }
311
312
313 mkl = mkl / (k5 - h30);
314 // Foote measure- footediv will be zero initially
315 float footediv = ((sqrt(norm)) * (sqrt(unit->m_prevNorm)));
316 if (footediv < 0.0001f)
317 footediv = 0.0001f;
318 float foote = 1.0 - (footesum / footediv); // 1.0 - similarity
319 // printf("mkl %f foote %f \n",mkl, foote);
320
321 unit->m_prevNorm = norm;
322 float sum = (ZIN0(1) * mkl) + (ZIN0(2) * foote);
323
324 // printf("sum %f thresh %f \n",sum, ZIN0(7));
325
326 // if over threshold, may also impose a 50mS wait here
327 if (sum > ZIN0(3) && (unit->m_waiting == 0)) {
328 outval = 1.0;
329 unit->m_waiting = 1;
330 unit->m_waitSamp = inNumSamples;
331 unit->m_waitLen = (int)(ZIN0(4) * (unit->mWorld->mSampleRate));
332 }
333
334 // take copy of this frame's magnitudes as prevframe
335
336 for (int i = 0; i < numbins; ++i)
337 q[i] = p->bin[i].mag;
338 }
339
340 Fill(inNumSamples, &ZOUT0(0), outval);
341 }
342
343
RunningSum_Ctor(RunningSum * unit)344 void RunningSum_Ctor(RunningSum* unit) {
345 SETCALC(RunningSum_next_k);
346
347 unit->msamp = (int)ZIN0(1);
348
349 // unit->mmeanmult= 1.0f/(unit->msamp);
350 unit->msum = 0.0f;
351 unit->msum2 = 0.0f;
352 unit->mcount = 0; // unit->msamp-1;
353
354 unit->msquares = (float*)RTAlloc(unit->mWorld, unit->msamp * sizeof(float));
355 // initialise to zeroes
356 for (int i = 0; i < unit->msamp; ++i)
357 unit->msquares[i] = 0.f;
358 OUT0(0) = 0.f;
359 }
360
RunningSum_Dtor(RunningSum * unit)361 void RunningSum_Dtor(RunningSum* unit) { RTFree(unit->mWorld, unit->msquares); }
362
363 // RMS is easy because convolution kernel can be updated just by deleting oldest sample and adding newest-
364 // half hanning window convolution etc requires updating values for all samples in memory on each iteration
RunningSum_next_k(RunningSum * unit,int inNumSamples)365 void RunningSum_next_k(RunningSum* unit, int inNumSamples) {
366 float* in = ZIN(0);
367 float* out = ZOUT(0);
368
369 int count = unit->mcount;
370 int samp = unit->msamp;
371
372 float* data = unit->msquares;
373 float sum = unit->msum;
374 // avoids floating point error accumulation over time- thanks to Ross Bencina
375 float sum2 = unit->msum2;
376
377 int todo = 0;
378 int done = 0;
379 while (done < inNumSamples) {
380 todo = sc_min(inNumSamples - done, samp - count);
381
382 for (int j = 0; j < todo; ++j) {
383 sum -= data[count];
384 float next = ZXP(in);
385 data[count] = next;
386 sum += next;
387 sum2 += next;
388 ZXP(out) = sum;
389 ++count;
390 }
391
392 done += todo;
393
394 if (count == samp) {
395 count = 0;
396 sum = sum2;
397 sum2 = 0;
398 }
399 }
400
401 unit->mcount = count;
402 unit->msum = sum;
403 unit->msum2 = sum2;
404 }
405
initFeatureDetectors(InterfaceTable * it)406 void initFeatureDetectors(InterfaceTable* it) {
407 DefineDtorUnit(PV_JensenAndersen);
408 DefineDtorUnit(PV_HainsworthFoote);
409 DefineDtorUnit(RunningSum);
410 }
411