Forgot a help fix: Drag a dock's title bar, not divider, to reposition
[supercollider.git] / server / plugins / FeatureDetection.cpp
blob89b9a035fc6d55a18ef07b9352902961b4dd5f2a
1 /*
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
35 float * m_prevframe;
36 int m_numbins;
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;
44 int m_fourkindex;
48 //FFT onset detector combining 2 advised features from Hainsworth PhD
49 struct PV_HainsworthFoote : public PV_OnsetDetectionBase
51 float m_prevNorm;
52 int m_5kindex,m_30Hzindex;
56 //for time domain onset detection/RMS
57 struct RunningSum : public Unit {
58 int msamp, mcount;
59 float msum,msum2;
60 //float mmeanmult;
61 float* msquares;
64 extern "C"
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; \
84 int bufOK = 1; \
85 World *world = unit->mWorld; \
86 SndBuf *buf; \
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; \
92 } else { \
93 bufOK = 0; \
94 buf = world->mSndBufs; \
95 if(unit->mWorld->mVerbosity > -1){ Print("FFT Ctor error: Buffer number overrun: %i\n", ibufnum); } \
96 } \
97 } else { \
98 buf = world->mSndBufs + ibufnum; \
99 } \
100 int numbins = buf->samples - 2 >> 1; \
101 if (!buf->data) { \
102 if(unit->mWorld->mVerbosity > -1){ Print("FFT Ctor error: Buffer %i not initialised.\n", ibufnum); } \
103 bufOK = 0; \
106 #define PV_FEAT_GET_BUF \
107 PV_FEAT_GET_BUF_UNLOCKED \
108 LOCK_SNDBUF(buf);
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);
120 if(bufOK) {
121 unit->m_prevframe = (float*)RTAlloc(unit->mWorld, insize);
122 memset(unit->m_prevframe, 0, insize);
125 unit->m_waiting=0;
126 unit->m_waitSamp=0;
127 unit->m_waitLen=0;
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);
142 unit->m_hfc= 0.0;
143 unit->m_hfe= 0.0;
144 unit->m_sf= 0.0;
145 unit->m_sc= 0.0;
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)
160 float outval=0.0;
161 float fbufnum = ZIN0(0);
163 if(unit->m_waiting==1) {
164 unit->m_waitSamp+=inNumSamples;
165 if(unit->m_waitSamp>=unit->m_waitLen)
166 unit->m_waiting=0;
169 if (!(fbufnum < 0.f))
170 //if buffer ready to process
172 PV_FEAT_GET_BUF
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);
187 int k= i+1;
188 float qmag= q[i];
189 magsum+=mag;
190 magsumk+=k*mag;
191 magsumkk+=k*k*mag;
192 sfsum+= fabs(mag- (qmag));
193 if(i>k4) hfesum+=mag;
196 float binmult= 1.f/numbins;
197 //normalise
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);
205 //if(sc<0.0) sc=0.0;
206 //if(hfe<0.0) hfe=0.0;
207 //if(hfc<0.0) hfc=0.0;
208 //if(sf<0.0) sf=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
217 unit->m_sc=sc;
218 unit->m_hfe=hfe;
219 unit->m_hfc=hfc;
220 unit->m_sf=sf;
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);
225 //does this trigger?
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");
234 outval=1.0;
235 unit->m_waiting=1;
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)
243 q[i]= p->bin[i].mag;
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)
277 float outval=0.0;
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
289 PV_FEAT_GET_BUF
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);
302 float qmag= q[i];
304 if(i>=h30 && i<k5) {
305 prevmag= qmag;
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;
316 norm+=mag*mag;
317 footesum+=mag*qmag;
321 mkl= mkl/(k5-h30);
322 //Foote measure- footediv will be zero initially
323 float footediv= ((sqrt(norm))*(sqrt(unit->m_prevNorm)));
324 if(footediv<0.0001f)
325 footediv=0.0001f;
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)) {
336 outval=1.0;
337 unit->m_waiting=1;
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)
345 q[i]= p->bin[i].mag;
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);
359 unit->msum=0.0f;
360 unit->msum2=0.0f;
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 )
379 float *in = ZIN(0);
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;
390 int todo=0;
391 int done=0;
392 while(done<inNumSamples) {
393 todo= sc_min(inNumSamples-done,samp-count);
395 for(int j=0;j<todo;++j) {
396 sum -=data[count];
397 float next= ZXP(in);
398 data[count]= next;
399 sum += next;
400 sum2 +=next;
401 ZXP(out) = sum;
402 ++count;
405 done+=todo;
407 if( count == samp ) {
408 count = 0;
409 sum = sum2;
410 sum2 = 0;
415 unit->mcount =count;
416 unit->msum = sum;
417 unit->msum2 = sum2;
420 void initFeatureDetectors(InterfaceTable *it)
422 DefineDtorUnit(PV_JensenAndersen);
423 DefineDtorUnit(PV_HainsworthFoote);
424 DefineDtorUnit(RunningSum);