server/plugins/onsetsds.c

   1 /*
   2         OnsetsDS - real time musical onset detection library.
   3     Copyright (c) 2007 Dan Stowell. All rights reserved.
   4
   5     This program is free software; you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation; either version 2 of the License, or
   8     (at your option) any later version.
   9
  10     This program is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License
  16     along with this program; if not, write to the Free Software
  17     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
  18 */
  19
  20 #include "onsetsds.h"
  21
  22
  23 #define ODS_DEBUG_POST_CSV 0
  24
  25 #ifdef _MSC_VER
  26 // msvc doesn't support c99
  27 #define hypotf _hypotf
  28 #define inline /* inline */
  29 #endif
  30
  31 static inline float onsetsds_phase_rewrap(float phase){
  32         return (phase>MINUSPI && phase<PI) ? phase : phase + TWOPI * (1.f + floorf((MINUSPI - phase) * INV_TWOPI));
  33 }
  34
  35
  36 size_t onsetsds_memneeded (int odftype, size_t fftsize, unsigned int medspan){
  37
  38         /*
  39         Need memory for:
  40         - median calculation (2 * medspan floats)
  41         - storing old values (whether as OdsPolarBuf or as weirder float lists)
  42         - storing the OdsPolarBuf (size is NOT sizeof(OdsPolarBuf) but is fftsize)
  43         - storing the PSP (numbins + 2 values)
  44         All these are floats.
  45         */
  46
  47         int numbins = (fftsize >> 1) - 1; // No of bins, not counting DC/nyq
  48
  49         switch(odftype){
  50                 case ODS_ODF_POWER:
  51                 case ODS_ODF_MAGSUM:
  52
  53                         // No old FFT frames needed, easy:
  54                         return (medspan+medspan + fftsize + numbins + 2) * sizeof(float);
  55
  56                 case ODS_ODF_COMPLEX:
  57                 case ODS_ODF_RCOMPLEX:
  58
  59                         return (medspan+medspan + fftsize + numbins + 2
  60                                         // For each bin (NOT dc/nyq) we store mag, phase and d_phase
  61                                         + numbins + numbins + numbins
  62                                 ) * sizeof(float);
  63
  64                 case ODS_ODF_PHASE:
  65                 case ODS_ODF_WPHASE:
  66
  67                         return (medspan+medspan + fftsize + numbins + 2
  68                                         // For each bin (NOT dc/nyq) we store phase and d_phase
  69                                         + numbins + numbins
  70                                 ) * sizeof(float);
  71
  72                 case ODS_ODF_MKL:
  73
  74                         return (medspan+medspan + fftsize + numbins + 2
  75                                         // For each bin (NOT dc/nyq) we store mag
  76                                         + numbins
  77                                 ) * sizeof(float);
  78
  79
  80                         break;
  81
  82         }
  83         return -1; //bleh
  84 }
  85
  86
  87 void onsetsds_init(OnsetsDS *ods, float *odsdata, int fftformat,
  88                            int odftype, size_t fftsize, unsigned int medspan, float srate){
  89
  90    int numbins, realnumbins;
  91
  92         // The main pointer to the processing area - other pointers will indicate areas within this
  93         ods->data = odsdata;
  94         // Set all vals in processing area to zero
  95         memset(odsdata, 0, onsetsds_memneeded(odftype, fftsize, medspan));
  96
  97         ods->srate = srate;
  98
  99         numbins  = (fftsize >> 1) - 1; // No of bins, not counting DC/nyq
 100         realnumbins = numbins + 2;
 101
 102         // Also point the other pointers to the right places
 103         ods->curr     = (OdsPolarBuf*) odsdata;
 104         ods->psp      = odsdata + fftsize;
 105         ods->odfvals  = odsdata + fftsize + realnumbins;
 106         ods->sortbuf  = odsdata + fftsize + realnumbins + medspan;
 107         ods->other    = odsdata + fftsize + realnumbins + medspan + medspan;
 108
 109         // Default settings for Adaptive Whitening, user can set own values after init
 110         onsetsds_setrelax(ods, 1.f, fftsize>>1);
 111         ods->floor    = 0.1;
 112
 113         switch(odftype){
 114                 case ODS_ODF_POWER:
 115                         ods->odfparam = 0.01; // "powthresh" in SC code
 116                         ods->normfactor = 2560.f / (realnumbins * fftsize);
 117                         break;
 118                 case ODS_ODF_MAGSUM:
 119                         ods->odfparam = 0.01; // "powthresh" in SC code
 120                         ods->normfactor = 113.137085f / (realnumbins * sqrt(fftsize));
 121                         break;
 122                 case ODS_ODF_COMPLEX:
 123                         ods->odfparam = 0.01; // "powthresh" in SC code
 124                         ods->normfactor = 231.70475f / pow(fftsize, 1.5);// / fftsize;
 125                         break;
 126                 case ODS_ODF_RCOMPLEX:
 127                         ods->odfparam = 0.01; // "powthresh" in SC code
 128                         ods->normfactor = 231.70475f / pow(fftsize, 1.5);// / fftsize;
 129                         break;
 130                 case ODS_ODF_PHASE:
 131                         ods->odfparam = 0.01; // "powthresh" in SC code
 132                         ods->normfactor = 5.12f / fftsize;// / fftsize;
 133                         break;
 134                 case ODS_ODF_WPHASE:
 135                         ods->odfparam = 0.0001; // "powthresh" in SC code. For WPHASE it's kind of superfluous.
 136                         ods->normfactor = 115.852375f / pow(fftsize, 1.5);// / fftsize;
 137                         break;
 138                 case ODS_ODF_MKL:
 139                         ods->odfparam = 0.01; // EPSILON parameter. Brossier recommends 1e-6 but I (ICMC 2007) found larger vals (e.g 0.01) to work better
 140                         ods->normfactor = 7.68f * 0.25f / fftsize;
 141                         break;
 142                 default:
 143                         printf("onsetsds_init ERROR: \"odftype\" is not a recognised value\n");
 144         }
 145
 146         ods->odfvalpost = 0.f;
 147         ods->odfvalpostprev = 0.f;
 148         ods->thresh   = 0.5f;
 149         ods->logmags = false;
 150
 151         ods->odftype  = odftype;
 152         ods->whtype   = ODS_WH_ADAPT_MAX1;
 153         ods->fftformat = fftformat;
 154
 155         ods->whiten   = (odftype != ODS_ODF_MKL); // Deactivate whitening for MKL by default
 156         ods->detected = false;
 157         ods->med_odd  = (medspan & 1) != 0;
 158
 159         ods->medspan  = medspan;
 160
 161         ods->mingap   = 0;
 162         ods->gapleft  = 0;
 163
 164         ods->fftsize  = fftsize;
 165         ods->numbins  = numbins;
 166
 167         //printf("End of _init: normfactor is %g\n", ods->normfactor);
 168
 169 }
 170
 171 bool onsetsds_process(OnsetsDS* ods, float* fftbuf){
 172         onsetsds_loadframe(ods, fftbuf);
 173
 174         onsetsds_whiten(ods);
 175         onsetsds_odf(ods);
 176         onsetsds_detect(ods);
 177
 178         return ods->detected;
 179 }
 180
 181
 182 void onsetsds_setrelax(OnsetsDS* ods, float time, size_t hopsize){
 183         ods->relaxtime = time;
 184         ods->relaxcoef = (time == 0.0f) ? 0.0f : exp((ods_log1 * hopsize)/(time * ods->srate));
 185 }
 186
 187
 188
 189 void onsetsds_loadframe(OnsetsDS* ods, float* fftbuf){
 190
 191         float *pos, *pos2, imag, real;
 192         int i;
 193
 194         switch(ods->fftformat){
 195                 case ODS_FFT_SC3_POLAR:
 196                         // The format is the same! dc, nyq, mag[1], phase[1], ...
 197                         memcpy(ods->curr, fftbuf, ods->fftsize * sizeof(float));
 198                         break;
 199
 200                 case ODS_FFT_SC3_COMPLEX:
 201
 202                         ods->curr->dc  = fftbuf[0];
 203                         ods->curr->nyq = fftbuf[1];
 204
 205                         // Then convert cartesian to polar:
 206                         pos = fftbuf + 2;
 207                         for(i=0; i< (ods->numbins << 1); i += 2){
 208                                 real = pos[i];
 209                                 imag = pos[i+1]; // Plus 1 rather than increment; seems to avoid LSU reject on my PPC
 210                                 ods->curr->bin[i].mag   = hypotf(imag, real);
 211                                 ods->curr->bin[i].phase = atan2f(imag, real);
 212                         }
 213                         break;
 214
 215                 case ODS_FFT_FFTW3_HC:
 216
 217                         ods->curr->dc  = fftbuf[0];
 218                         ods->curr->nyq = fftbuf[ods->fftsize>>1];
 219
 220                         // Then convert cartesian to polar:
 221                         // (Starting positions: real and imag for bin 1)
 222                         pos  = fftbuf + 1;
 223                         pos2 = fftbuf + ods->fftsize - 1;
 224                         for(i=0; i<ods->numbins; i++){
 225                                 real = *(pos++);
 226                                 imag = *(pos2--);
 227                                 ods->curr->bin[i].mag   = hypotf(imag, real);
 228                                 ods->curr->bin[i].phase = atan2f(imag, real);
 229                         }
 230                         break;
 231
 232                 case ODS_FFT_FFTW3_R2C:
 233
 234                         ods->curr->dc  = fftbuf[0];
 235                         ods->curr->nyq = fftbuf[ods->fftsize];
 236
 237                         // Then convert cartesian to polar:
 238                         pos = fftbuf + 2;
 239                         for(i=0; i<ods->numbins; i++){
 240                                 real = *(pos++);
 241                                 imag = *(pos++);
 242                                 ods->curr->bin[i].mag   = hypotf(imag, real);
 243                                 ods->curr->bin[i].phase = atan2f(imag, real);
 244                         }
 245                         break;
 246
 247         }
 248
 249         // Conversion to log-domain magnitudes, including re-scaling to aim back at the zero-to-one range.
 250         // Not well tested yet.
 251         if(ods->logmags){
 252                 for(i=0; i<ods->numbins; i++){
 253                         ods->curr->bin[i].mag =
 254                                 (log(ods_max(ods->curr->bin[i].mag, ODS_LOG_LOWER_LIMIT)) - ODS_LOGOF_LOG_LOWER_LIMIT) * ODS_ABSINVOF_LOGOF_LOG_LOWER_LIMIT;
 255                 }
 256                 ods->curr->dc =
 257                         (log(ods_max(ods_abs(ods->curr->dc ), ODS_LOG_LOWER_LIMIT)) - ODS_LOGOF_LOG_LOWER_LIMIT) * ODS_ABSINVOF_LOGOF_LOG_LOWER_LIMIT;
 258                 ods->curr->nyq =
 259                         (log(ods_max(ods_abs(ods->curr->nyq), ODS_LOG_LOWER_LIMIT)) - ODS_LOGOF_LOG_LOWER_LIMIT) * ODS_ABSINVOF_LOGOF_LOG_LOWER_LIMIT;
 260         }
 261
 262 }
 263
 264 void onsetsds_whiten(OnsetsDS* ods){
 265
 266         float val,oldval, relaxcoef, floor;
 267         int numbins, i;
 268         OdsPolarBuf *curr;
 269         float *psp;
 270         float *pspp1; // Offset by 1, avoids quite a lot of "+1"s in the following code
 271
 272         if(ods->whtype == ODS_WH_NONE){
 273                 //printf("onsetsds_whiten(): ODS_WH_NONE, skipping\n");
 274                 return;
 275         }
 276
 277         // NB: Apart from the above, ods->whtype is currently IGNORED and only one mode is used.
 278
 279         relaxcoef = ods->relaxcoef;
 280         numbins = ods->numbins;
 281         curr = ods->curr;
 282         psp = ods->psp;
 283         pspp1 = psp + 1;
 284         floor = ods->floor;
 285
 286         //printf("onsetsds_whiten: relaxcoef=%g, relaxtime=%g, floor=%g\n", relaxcoef, ods->relaxtime, floor);
 287
 288         ////////////////////// For each bin, update the record of the peak value /////////////////////
 289
 290         val = fabs(curr->dc);   // Grab current magnitude
 291         oldval = psp[0];
 292         // If it beats the amplitude stored then that's our new amplitude;
 293         // otherwise our new amplitude is a decayed version of the old one
 294         if(val < oldval) {
 295                 val = val + (oldval - val) * relaxcoef;
 296         }
 297         psp[0] = val; // Store the "amplitude trace" back
 298
 299         val = fabs(curr->nyq);
 300         oldval = pspp1[numbins];
 301         if(val < oldval) {
 302                 val = val + (oldval - val) * relaxcoef;
 303         }
 304         pspp1[numbins] = val;
 305
 306         for(i=0; i<numbins; ++i){
 307                 val = fabs(curr->bin[i].mag);
 308                 oldval = pspp1[i];
 309                 if(val < oldval) {
 310                         val = val + (oldval - val) * relaxcoef;
 311                 }
 312                 pspp1[i] = val;
 313         }
 314
 315         //////////////////////////// Now for each bin, rescale the current magnitude ////////////////////////////
 316         curr->dc  /= ods_max(floor, psp[0]);
 317         curr->nyq /= ods_max(floor, pspp1[numbins]);
 318         for(i=0; i<numbins; ++i){
 319                 curr->bin[i].mag /= ods_max(floor, pspp1[i]);
 320         }
 321 }
 322
 323 void onsetsds_odf(OnsetsDS* ods){
 324
 325         int numbins = ods->numbins;
 326         OdsPolarBuf *curr = ods->curr;
 327         float* val = ods->odfvals;
 328         int i, tbpointer;
 329         float deviation, diff, curmag;
 330         double totdev;
 331         float predmag, predphase, yesterphase, yesterphasediff;
 332         float yestermag;
 333
 334
 335         bool rectify = true;
 336
 337         // Here we shunt the "old" ODF values down one place
 338         memcpy(val + 1, val, (ods->medspan - 1)*sizeof(float));
 339
 340         // Now calculate a new value and store in ods->odfvals[0]
 341         switch(ods->odftype){
 342                 case ODS_ODF_POWER:
 343
 344                         *val = (curr->nyq  *  curr->nyq)  +  (curr->dc  *  curr->dc);
 345                         for(i=0; i<numbins; i++){
 346                                 *val += curr->bin[i].mag  *  curr->bin[i].mag;
 347                         }
 348                         break;
 349
 350                 case ODS_ODF_MAGSUM:
 351
 352                         *val = ods_abs(curr->nyq) + ods_abs(curr->dc);
 353
 354                         for(i=0; i<numbins; i++){
 355                                 *val += ods_abs(curr->bin[i].mag);
 356                         }
 357                         break;
 358
 359                 case ODS_ODF_COMPLEX:
 360                         rectify = false;
 361                         // ...and then drop through to:
 362                 case ODS_ODF_RCOMPLEX:
 363
 364                         // Note: "other" buf is stored in this format: mag[0],phase[0],d_phase[0],mag[1],phase[1],d_phase[1], ...
 365
 366                         // Iterate through, calculating the deviation from expected value.
 367                         totdev = 0.0;
 368                         tbpointer = 0;
 369                         for (i=0; i<numbins; ++i) {
 370                                 curmag = ods_abs(curr->bin[i].mag);
 371
 372                                 // Predict mag as yestermag
 373                                 predmag         = ods->other[tbpointer++];
 374                                 yesterphase     = ods->other[tbpointer++];
 375                                 yesterphasediff = ods->other[tbpointer++];
 376
 377                                 // Thresholding as Brossier did - discard (ignore) bin's deviation if bin's power is minimal
 378                                 if(curmag > ods->odfparam) {
 379                                         // If rectifying, ignore decreasing bins
 380                                         if((!rectify) || !(curmag < predmag)){
 381
 382                                                 // Predict phase as yesterval + yesterfirstdiff
 383                                                 predphase = yesterphase + yesterphasediff;
 384
 385                                                 // Here temporarily using the "deviation" var to store the phase difference
 386                                                 //  so that the rewrap macro can use it more efficiently
 387                                                 deviation = predphase - curr->bin[i].phase;
 388
 389                                                 // Deviation is Euclidean distance between predicted and actual.
 390                                                 // In polar coords: sqrt(r1^2 +  r2^2 - r1r2 cos (theta1 - theta2))
 391                                                 deviation = sqrtf(predmag * predmag + curmag * curmag
 392                                                                                   - predmag * curmag * cosf(onsetsds_phase_rewrap(deviation))
 393                                                                                 );
 394
 395                                                 totdev += deviation;
 396                                         }
 397                                 }
 398                         }
 399
 400                         // totdev will be the output, but first we need to fill tempbuf with today's values, ready for tomorrow.
 401                         tbpointer = 0;
 402                         for (i=0; i<numbins; ++i) {
 403                                 ods->other[tbpointer++] = ods_abs(curr->bin[i].mag); // Storing mag
 404                                 diff = curr->bin[i].phase - ods->other[tbpointer]; // Retrieving yesterphase from buf
 405                                 ods->other[tbpointer++] = curr->bin[i].phase; // Storing phase
 406                                 // Wrap onto +-PI range
 407                                 diff = onsetsds_phase_rewrap(diff);
 408
 409                                 ods->other[tbpointer++] = diff; // Storing first diff to buf
 410
 411                         }
 412                         *val = (float)totdev;
 413
 414                         break;
 415
 416
 417                 case ODS_ODF_PHASE:
 418                         rectify = false; // So, actually, "rectify" means "useweighting" in this context
 419                         // ...and then drop through to:
 420                 case ODS_ODF_WPHASE:
 421
 422                         // Note: "other" buf is stored in this format: phase[0],d_phase[0],phase[1],d_phase[1], ...
 423
 424                         // Iterate through, calculating the deviation from expected value.
 425                         totdev = 0.0;
 426                         tbpointer = 0;
 427                         for (i=0; i<numbins; ++i) {
 428                                 // Thresholding as Brossier did - discard (ignore) bin's phase deviation if bin's power is low
 429                                 if(ods_abs(curr->bin[i].mag) > ods->odfparam) {
 430
 431                                         // Deviation is the *second difference* of the phase, which is calc'ed as curval - yesterval - yesterfirstdiff
 432                                         deviation = curr->bin[i].phase - ods->other[tbpointer++] - ods->other[tbpointer++];
 433                                         // Wrap onto +-PI range
 434                                         deviation = onsetsds_phase_rewrap(deviation);
 435
 436                                         if(rectify){ // "rectify" meaning "useweighting"...
 437                                                 totdev += fabs(deviation * ods_abs(curr->bin[i].mag));
 438                                         } else {
 439                                                 totdev += fabs(deviation);
 440                                         }
 441                                 }
 442                         }
 443
 444                         // totdev will be the output, but first we need to fill tempbuf with today's values, ready for tomorrow.
 445                         tbpointer = 0;
 446                         for (i=0; i<numbins; ++i) {
 447                                 diff = curr->bin[i].phase - ods->other[tbpointer]; // Retrieving yesterphase from buf
 448                                 ods->other[tbpointer++] = curr->bin[i].phase; // Storing phase
 449                                 // Wrap onto +-PI range
 450                                 diff = onsetsds_phase_rewrap(diff);
 451
 452                                 ods->other[tbpointer++] = diff; // Storing first diff to buf
 453
 454                         }
 455                         *val = (float)totdev;
 456                         break;
 457
 458
 459                 case ODS_ODF_MKL:
 460
 461                         // Iterate through, calculating the Modified Kullback-Liebler distance
 462                         totdev = 0.0;
 463                         tbpointer = 0;
 464                         for (i=0; i<numbins; ++i) {
 465                                 curmag = ods_abs(curr->bin[i].mag);
 466                                 yestermag = ods->other[tbpointer];
 467
 468                                 // Here's the main implementation of Brossier's MKL eq'n (eqn 2.9 from his thesis):
 469                                 deviation = ods_abs(curmag) / (ods_abs(yestermag) + ods->odfparam);
 470                                 totdev += log(1.f + deviation);
 471
 472                                 // Store the mag as yestermag
 473                                 ods->other[tbpointer++] = curmag;
 474                         }
 475                         *val = (float)totdev;
 476                         break;
 477
 478         }
 479
 480 #if ODS_DEBUG_POST_CSV
 481         printf("%g,", *val);
 482         printf("%g,", ods->odfvals[0] * ods->normfactor);
 483 #endif
 484
 485         ods->odfvals[0] *= ods->normfactor;
 486 }
 487 // End of ODF function
 488
 489 void SelectionSort(float *array, int length);
 490 void SelectionSort(float *array, int length)
 491 {
 492   // Algo is simply based on http://en.wikibooks.org/wiki/Algorithm_implementation/Sorting/Selection_sort
 493   int max, i;
 494   float temp;
 495   while(length > 0)
 496   {
 497     max = 0;
 498     for(i = 1; i < length; i++)
 499       if(array[i] > array[max])
 500         max = i;
 501     temp = array[length-1];
 502     array[length-1] = array[max];
 503     array[max] = temp;
 504     length--;
 505   }
 506 }
 507
 508
 509 void onsetsds_detect(OnsetsDS* ods){
 510
 511         float* sortbuf = ods->sortbuf;
 512         int medspan = ods->medspan;
 513
 514         // Shift the yesterval to its rightful place
 515         ods->odfvalpostprev = ods->odfvalpost;
 516
 517         ///////// MEDIAN REMOVAL ////////////
 518
 519         // Copy odfvals to sortbuf
 520         memcpy(sortbuf, ods->odfvals, medspan * sizeof(float));
 521
 522         // Sort sortbuf
 523         SelectionSort(sortbuf, medspan);
 524
 525         // Subtract the middlest value === the median
 526         if(ods->med_odd){
 527                 ods->odfvalpost = ods->odfvals[0]
 528                            - sortbuf[(medspan - 1) >> 1];
 529         }else{
 530                 ods->odfvalpost = ods->odfvals[0]
 531                            - ((sortbuf[medspan >> 1]
 532                                  + sortbuf[(medspan >> 1) - 1]) * 0.5f);
 533
 534         }
 535
 536         // Detection not allowed if we're too close to a previous detection.
 537         if(ods->gapleft != 0) {
 538                 ods->gapleft--;
 539                 ods->detected = false;
 540         } else {
 541                 // Now do the detection.
 542                 ods->detected = (ods->odfvalpost > ods->thresh) && (ods->odfvalpostprev <= ods->thresh);
 543                 if(ods->detected){
 544                         ods->gapleft = ods->mingap;
 545                 }
 546         }
 547 #if ODS_DEBUG_POST_CSV
 548         printf("%g\n", ods->odfvalpost);
 549 #endif
 550 }
 551
 552