2 SuperCollider real time audio synthesis system
3 Copyright (c) 2002 James McCartney. All rights reserved.
4 http://www.audiosynth.com
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.
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.
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
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 Verlag.
24 //Hainsworth, S. (2003) Techniques for the Automated Analysis of Musical Audio. PhD, university of cambridge engineering dept.
26 //possible to make a Goto style Detector for a given band and with history of two samples-
27 //should do separately as PV_GotoBandTrack
29 //next perhaps Duxbury et al/ Mauri et al different conception of high frequency content with ratio of changes
31 #include "FFT_UGens.h"
33 struct PV_OnsetDetectionBase
: public Unit
37 int m_waiting
, m_waitSamp
, m_waitLen
;
40 //FFT onset detector combining 4 advised features from Jensen/Andersen
41 struct PV_JensenAndersen
: public PV_OnsetDetectionBase
43 float m_hfc
,m_hfe
,m_sc
,m_sf
;
48 //FFT onset detector combining 2 advised features from Hainsworth PhD
49 struct PV_HainsworthFoote
: public PV_OnsetDetectionBase
52 int m_5kindex
,m_30Hzindex
;
56 //for time domain onset detection/RMS
57 struct RunningSum
: public Unit
{
66 void PV_OnsetDetectionBase_Ctor(PV_OnsetDetectionBase
*unit
);
67 void PV_OnsetDetectionBase_Dtor(PV_OnsetDetectionBase
*unit
);
69 void PV_JensenAndersen_Ctor(PV_JensenAndersen
*unit
);
70 void PV_JensenAndersen_Dtor(PV_JensenAndersen
*unit
);
71 void PV_JensenAndersen_next(PV_JensenAndersen
*unit
, int inNumSamples
);
73 void PV_HainsworthFoote_Ctor(PV_HainsworthFoote
*unit
);
74 void PV_HainsworthFoote_Dtor(PV_HainsworthFoote
*unit
);
75 void PV_HainsworthFoote_next(PV_HainsworthFoote
*unit
, int inNumSamples
);
77 void RunningSum_next_k(RunningSum
*unit
, int inNumSamples
);
78 void RunningSum_Ctor(RunningSum
* unit
);
79 void RunningSum_Dtor(RunningSum
* unit
);
82 #define PV_FEAT_GET_BUF_UNLOCKED \
83 uint32 ibufnum = (uint32)fbufnum; \
85 World *world = unit->mWorld; \
87 if (ibufnum >= world->mNumSndBufs) { \
88 int localBufNum = ibufnum - world->mNumSndBufs; \
89 Graph *parent = unit->mParent; \
90 if(localBufNum <= parent->localBufNum) { \
91 buf = parent->mLocalSndBufs + localBufNum; \
94 buf = world->mSndBufs; \
95 if(unit->mWorld->mVerbosity > -1){ Print("FFT Ctor error: Buffer number overrun: %i\n", ibufnum); } \
98 buf = world->mSndBufs + ibufnum; \
100 int numbins = buf->samples - 2 >> 1; \
102 if(unit->mWorld->mVerbosity > -1){ Print("FFT Ctor error: Buffer %i not initialised.\n", ibufnum); } \
106 #define PV_FEAT_GET_BUF \
107 PV_FEAT_GET_BUF_UNLOCKED \
111 void PV_OnsetDetectionBase_Ctor(PV_OnsetDetectionBase
*unit
)
113 float fbufnum
= ZIN0(0);
115 PV_FEAT_GET_BUF_UNLOCKED
117 unit
->m_numbins
= numbins
;
118 int insize
= unit
->m_numbins
* sizeof(float);
121 unit
->m_prevframe
= (float*)RTAlloc(unit
->mWorld
, insize
);
122 memset(unit
->m_prevframe
, 0, insize
);
130 void PV_OnsetDetectionBase_Dtor(PV_OnsetDetectionBase
*unit
)
132 if(unit
->m_prevframe
)
133 RTFree(unit
->mWorld
, unit
->m_prevframe
);
138 void PV_JensenAndersen_Ctor(PV_JensenAndersen
*unit
)
140 PV_OnsetDetectionBase_Ctor(unit
);
147 unit
->m_fourkindex
= (int)(4000.0/(unit
->mWorld
->mSampleRate
))*(unit
->m_numbins
);
149 SETCALC(PV_JensenAndersen_next
);
152 void PV_JensenAndersen_Dtor(PV_JensenAndersen
*unit
)
154 PV_OnsetDetectionBase_Dtor(unit
);
158 void PV_JensenAndersen_next(PV_JensenAndersen
*unit
, int inNumSamples
)
161 float fbufnum
= ZIN0(0);
163 if(unit
->m_waiting
==1) {
164 unit
->m_waitSamp
+=inNumSamples
;
165 if(unit
->m_waitSamp
>=unit
->m_waitLen
)
169 if (!(fbufnum
< 0.f
))
170 //if buffer ready to process
174 SCPolarBuf
*p
= ToPolarApx(buf
);
176 //four spectral features useful for onset detection according to Jensen/Andersen
178 float magsum
=0.0, magsumk
=0.0, magsumkk
=0.0, sfsum
=0.0, hfesum
=0.0;
180 float * q
= unit
->m_prevframe
;
182 int k4
= unit
->m_fourkindex
;
184 //ignores dc, nyquist
185 for (int i
=0; i
<numbins
; ++i
) {
186 float mag
= ((p
->bin
[i
]).mag
);
192 sfsum
+= fabs(mag
- (qmag
));
193 if(i
>k4
) hfesum
+=mag
;
196 float binmult
= 1.f
/numbins
;
198 float sc
= (magsumk
/magsum
)*binmult
;
199 float hfe
= hfesum
*binmult
;
200 float hfc
= magsumkk
*binmult
*binmult
*binmult
;
201 float sf
= sfsum
*binmult
;
203 //printf("sc %f hfe %f hfc %f sf %f \n",sc, hfe, hfc, sf);
206 //if(hfe<0.0) hfe=0.0;
207 //if(hfc<0.0) hfc=0.0;
210 //ratio of current to previous frame perhaps better indicator than first derivative difference
211 float scdiff
= sc
-(unit
->m_sc
);
212 float hfediff
= hfe
-(unit
->m_hfe
);
213 float hfcdiff
= hfc
-(unit
->m_hfc
);
214 float sfdiff
= sf
-(unit
->m_sf
);
216 //store as old frame values for taking difference
222 //printf("sc %f hfe %f hfc %f sf %f \n",sc, hfe, hfc, sf);
223 //printf("sc %f hfe %f hfc %f sf %f \n",scdiff, hfediff, hfcdiff, sfdiff);
226 //may need to take derivatives across previous frames by storing old values
228 float sum
= (ZIN0(1)*scdiff
)+(ZIN0(2)*hfediff
)+(ZIN0(3)*hfcdiff
)+(ZIN0(4)*sfdiff
);
230 //printf("sum %f thresh %f \n",sum, ZIN0(7));
232 //if over threshold, may also impose a wait here
233 if(sum
>ZIN0(5) && (unit
->m_waiting
==0)) {//printf("bang! \n");
236 unit
->m_waitSamp
=inNumSamples
;
237 unit
->m_waitLen
=(int)(ZIN0(6)*(world
->mSampleRate
));
240 //take copy of this frame's magnitudes as prevframe
242 for (int i
=0; i
<numbins
; ++i
)
246 Fill(inNumSamples
, &ZOUT0(0), outval
);
250 void PV_HainsworthFoote_Ctor(PV_HainsworthFoote
*unit
)
252 PV_OnsetDetectionBase_Ctor(unit
);
254 World
*world
= unit
->mWorld
;
256 unit
->m_5kindex
= (int)((5000.0/(world
->mSampleRate
))*(unit
->m_numbins
));
257 unit
->m_30Hzindex
= (int)((30.0/(world
->mSampleRate
))*(unit
->m_numbins
));
259 unit
->m_prevNorm
= 1.0;
261 //unit->m_5kindex, unit->m_30Hzindex,
262 //printf("numbins %d sr %d \n", unit->m_numbins, world->mSampleRate);
263 //printf("test %d sr %f 5k %d 30Hz %d\n", unit->m_numbins, world->mSampleRate, unit->m_5kindex, unit->m_30Hzindex);
265 SETCALC(PV_HainsworthFoote_next
);
268 void PV_HainsworthFoote_Dtor(PV_HainsworthFoote
*unit
)
270 PV_OnsetDetectionBase_Dtor(unit
);
273 static const float lmult
= 1.442695040889; //loge(2) reciprocal
275 void PV_HainsworthFoote_next(PV_HainsworthFoote
*unit
, int inNumSamples
)
278 float fbufnum
= ZIN0(0);
280 if(unit
->m_waiting
==1)
282 unit
->m_waitSamp
+=inNumSamples
;
283 if(unit
->m_waitSamp
>=unit
->m_waitLen
) {unit
->m_waiting
=0;}
286 if (!(fbufnum
< 0.f
))
287 //if buffer ready to process
291 SCPolarBuf
*p
= ToPolarApx(buf
);
293 float dnk
, prevmag
, mkl
=0.0, footesum
=0.0, norm
=0.0;
295 float * q
= unit
->m_prevframe
;
297 int k5
= unit
->m_5kindex
;
298 int h30
= unit
->m_30Hzindex
;
300 for (int i
=0; i
<numbins
; ++i
) {
301 float mag
= ((p
->bin
[i
]).mag
);
306 //avoid divide by zero
307 if(prevmag
<0.0001) prevmag
=0.0001;
309 //no log2 in maths library, so use log2(x)= log(x)/log(2) where log is to base e
310 //could just use log and ignore scale factor but hey let's stay accurate to the source for now
311 dnk
= log(mag
/prevmag
)*lmult
;
313 if(dnk
>0.0) mkl
+=dnk
;
322 //Foote measure- footediv will be zero initially
323 float footediv
= ((sqrt(norm
))*(sqrt(unit
->m_prevNorm
)));
326 float foote
= 1.0- (footesum
/footediv
); //1.0 - similarity
327 //printf("mkl %f foote %f \n",mkl, foote);
329 unit
->m_prevNorm
= norm
;
330 float sum
= (ZIN0(1)*mkl
)+(ZIN0(2)*foote
);
332 //printf("sum %f thresh %f \n",sum, ZIN0(7));
334 //if over threshold, may also impose a 50mS wait here
335 if(sum
>ZIN0(3) && (unit
->m_waiting
==0)) {
338 unit
->m_waitSamp
=inNumSamples
;
339 unit
->m_waitLen
=(int)(ZIN0(4)*(unit
->mWorld
->mSampleRate
));
342 //take copy of this frame's magnitudes as prevframe
344 for (int i
=0; i
<numbins
; ++i
)
348 Fill(inNumSamples
, &ZOUT0(0), outval
);
352 void RunningSum_Ctor( RunningSum
* unit
)
354 SETCALC(RunningSum_next_k
);
356 unit
->msamp
= (int) ZIN0(1);
358 //unit->mmeanmult= 1.0f/(unit->msamp);
361 unit
->mcount
=0; //unit->msamp-1;
363 unit
->msquares
= (float*)RTAlloc(unit
->mWorld
, unit
->msamp
* sizeof(float));
364 //initialise to zeroes
365 for(int i
=0; i
<unit
->msamp
; ++i
)
366 unit
->msquares
[i
]=0.f
;
370 void RunningSum_Dtor(RunningSum
*unit
)
372 RTFree(unit
->mWorld
, unit
->msquares
);
375 //RMS is easy because convolution kernel can be updated just by deleting oldest sample and adding newest-
376 //half hanning window convolution etc requires updating values for all samples in memory on each iteration
377 void RunningSum_next_k( RunningSum
*unit
, int inNumSamples
)
380 float *out
= ZOUT(0);
382 int count
= unit
->mcount
;
383 int samp
= unit
->msamp
;
385 float * data
= unit
->msquares
;
386 float sum
= unit
->msum
;
387 //avoids floating point error accumulation over time- thanks to Ross Bencina
388 float sum2
= unit
->msum2
;
392 while(done
<inNumSamples
) {
393 todo
= sc_min(inNumSamples
-done
,samp
-count
);
395 for(int j
=0;j
<todo
;++j
) {
407 if( count
== samp
) {
420 void initFeatureDetectors(InterfaceTable
*it
)
422 DefineDtorUnit(PV_JensenAndersen
);
423 DefineDtorUnit(PV_HainsworthFoote
);
424 DefineDtorUnit(RunningSum
);