src/onsetdetection.c

   1 /*
   2    Copyright (C) 2003 Paul Brossier
   3
   4    This program is free software; you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation; either version 2 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program; if not, write to the Free Software
  16    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include "aubio_priv.h"
  21 #include "sample.h"
  22 #include "fft.h"
  23 #include "mathutils.h"
  24 #include "hist.h"
  25 #include "onsetdetection.h"
  26
  27
  28 /** structure to store object state */
  29 struct _aubio_onsetdetection_t {
  30   aubio_onsetdetection_type type; /**< onset detection type */
  31   /** Pointer to aubio_onsetdetection_<type> function */
  32   void (*funcpointer)(aubio_onsetdetection_t *o,
  33       cvec_t * fftgrain, fvec_t * onset);
  34   smpl_t threshold;      /**< minimum norm threshold for phase and specdiff */
  35   fvec_t *oldmag;        /**< previous norm vector */
  36   fft_data_t *meas;      /**< current onset detection measure complex vector */
  37   fvec_t *dev1 ;         /**< current onset detection measure vector */
  38   fvec_t *theta1;        /**< previous phase vector, one frame behind */
  39   fvec_t *theta2;        /**< previous phase vector, two frames behind */
  40   aubio_hist_t * histog; /**< histogram */
  41 };
  42
  43
  44 /* Energy based onset detection function */
  45 void aubio_onsetdetection_energy  (aubio_onsetdetection_t *o UNUSED,
  46     cvec_t * fftgrain, fvec_t * onset) {
  47   uint_t i,j;
  48   for (i=0;i<fftgrain->channels;i++) {
  49     onset->data[i][0] = 0.;
  50     for (j=0;j<fftgrain->length;j++) {
  51       onset->data[i][0] += SQR(fftgrain->norm[i][j]);
  52     }
  53   }
  54 }
  55
  56 /* High Frequency Content onset detection function */
  57 void aubio_onsetdetection_hfc(aubio_onsetdetection_t *o UNUSED,
  58     cvec_t * fftgrain, fvec_t * onset){
  59   uint_t i,j;
  60   for (i=0;i<fftgrain->channels;i++) {
  61     onset->data[i][0] = 0.;
  62     for (j=0;j<fftgrain->length;j++) {
  63       onset->data[i][0] += (j+1)*fftgrain->norm[i][j];
  64     }
  65   }
  66 }
  67
  68
  69 /* Complex Domain Method onset detection function */
  70 void aubio_onsetdetection_complex (aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset) {
  71   uint_t i, j;
  72   uint_t nbins = fftgrain->length;
  73   for (i=0;i<fftgrain->channels; i++)  {
  74     onset->data[i][0] = 0.;
  75     for (j=0;j<nbins; j++)  {
  76       o->dev1->data[i][j]    = aubio_unwrap2pi(
  77           fftgrain->phas[i][j]
  78           -2.0*o->theta1->data[i][j]+
  79           o->theta2->data[i][j]);
  80 #ifdef HAVE_COMPLEX_H
  81       o->meas[j] = fftgrain->norm[i][j]*CEXPC(I*o->dev1->data[i][j]);
  82       /* sum on all bins */
  83       onset->data[i][0]     += //(fftgrain->norm[i][j]);
  84           SQRT(SQR( REAL(o->oldmag->data[i][j]-o->meas[j]) )
  85             +  SQR( IMAG(o->oldmag->data[i][j]-o->meas[j]) )
  86             );
  87 #else
  88       o->meas[j]             = (fftgrain->norm[i][j])*COS(o->dev1->data[i][j]);
  89       o->meas[(nbins-1)*2-j] = (fftgrain->norm[i][j])*SIN(o->dev1->data[i][j]);
  90       /* sum on all bins */
  91       onset->data[i][0]     += //(fftgrain->norm[i][j]);
  92           SQRT(SQR( (o->oldmag->data[i][j]-o->meas[j]) )
  93             +  SQR( (-o->meas[(nbins-1)*2-j]) )
  94             );
  95 #endif
  96       /* swap old phase data (need to remember 2 frames behind)*/
  97       o->theta2->data[i][j] = o->theta1->data[i][j];
  98       o->theta1->data[i][j] = fftgrain->phas[i][j];
  99       /* swap old magnitude data (1 frame is enough) */
 100       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 101     }
 102   }
 103 }
 104
 105
 106 /* Phase Based Method onset detection function */
 107 void aubio_onsetdetection_phase(aubio_onsetdetection_t *o,
 108     cvec_t * fftgrain, fvec_t * onset){
 109   uint_t i, j;
 110   uint_t nbins = fftgrain->length;
 111   for (i=0;i<fftgrain->channels; i++)  {
 112     onset->data[i][0] = 0.0f;
 113     o->dev1->data[i][0]=0.;
 114     for ( j=0;j<nbins; j++ )  {
 115       o->dev1->data[i][j] =
 116         aubio_unwrap2pi(
 117             fftgrain->phas[i][j]
 118             -2.0*o->theta1->data[i][j]
 119             +o->theta2->data[i][j]);
 120       if ( o->threshold < fftgrain->norm[i][j] )
 121         o->dev1->data[i][j] = ABS(o->dev1->data[i][j]);
 122       else
 123         o->dev1->data[i][j] = 0.0f;
 124       /* keep a track of the past frames */
 125       o->theta2->data[i][j] = o->theta1->data[i][j];
 126       o->theta1->data[i][j] = fftgrain->phas[i][j];
 127     }
 128     /* apply o->histogram */
 129     aubio_hist_dyn_notnull(o->histog,o->dev1);
 130     /* weight it */
 131     aubio_hist_weight(o->histog);
 132     /* its mean is the result */
 133     onset->data[i][0] = aubio_hist_mean(o->histog);
 134     //onset->data[i][0] = vec_mean(o->dev1);
 135   }
 136 }
 137
 138 /* Spectral difference method onset detection function */
 139 void aubio_onsetdetection_specdiff(aubio_onsetdetection_t *o,
 140     cvec_t * fftgrain, fvec_t * onset){
 141   uint_t i, j;
 142   uint_t nbins = fftgrain->length;
 143   for (i=0;i<fftgrain->channels; i++)  {
 144     onset->data[i][0] = 0.0f;
 145     for (j=0;j<nbins; j++)  {
 146       o->dev1->data[i][j] = SQRT(
 147           ABS(SQR( fftgrain->norm[i][j])
 148             - SQR(o->oldmag->data[i][j])));
 149       if (o->threshold < fftgrain->norm[i][j] )
 150         o->dev1->data[i][j] = ABS(o->dev1->data[i][j]);
 151       else
 152         o->dev1->data[i][j] = 0.0f;
 153       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 154     }
 155
 156     /* apply o->histogram (act somewhat as a low pass on the
 157      * overall function)*/
 158     aubio_hist_dyn_notnull(o->histog,o->dev1);
 159     /* weight it */
 160     aubio_hist_weight(o->histog);
 161     /* its mean is the result */
 162     onset->data[i][0] = aubio_hist_mean(o->histog);
 163
 164   }
 165 }
 166
 167 /* Kullback Liebler onset detection function
 168  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 169  * negative (1.+) and infinite values (+1.e-10) */
 170 void aubio_onsetdetection_kl(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 171   uint_t i,j;
 172   for (i=0;i<fftgrain->channels;i++) {
 173     onset->data[i][0] = 0.;
 174     for (j=0;j<fftgrain->length;j++) {
 175       onset->data[i][0] += fftgrain->norm[i][j]
 176         *LOG(1.+fftgrain->norm[i][j]/(o->oldmag->data[i][j]+1.e-10));
 177       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 178     }
 179     if (isnan(onset->data[i][0])) onset->data[i][0] = 0.;
 180   }
 181 }
 182
 183 /* Modified Kullback Liebler onset detection function
 184  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 185  * negative (1.+) and infinite values (+1.e-10) */
 186 void aubio_onsetdetection_mkl(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 187   uint_t i,j;
 188   for (i=0;i<fftgrain->channels;i++) {
 189     onset->data[i][0] = 0.;
 190     for (j=0;j<fftgrain->length;j++) {
 191       onset->data[i][0] += LOG(1.+fftgrain->norm[i][j]/(o->oldmag->data[i][j]+1.e-10));
 192       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 193     }
 194     if (isnan(onset->data[i][0])) onset->data[i][0] = 0.;
 195   }
 196 }
 197
 198 /* Spectral flux */
 199 void aubio_onsetdetection_specflux(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 200   uint_t i, j;
 201   for (i=0;i<fftgrain->channels;i++) {
 202     onset->data[i][0] = 0.;
 203     for (j=0;j<fftgrain->length;j++) {
 204       if (fftgrain->norm[i][j] > o->oldmag->data[i][j])
 205         onset->data[i][0] += fftgrain->norm[i][j] - o->oldmag->data[i][j];
 206       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 207     }
 208   }
 209 }
 210
 211 /* Generic function pointing to the choosen one */
 212 void
 213 aubio_onsetdetection(aubio_onsetdetection_t *o, cvec_t * fftgrain,
 214     fvec_t * onset) {
 215   o->funcpointer(o,fftgrain,onset);
 216 }
 217
 218 /* Allocate memory for an onset detection
 219  * depending on the choosen type, allocate memory as needed
 220  */
 221 aubio_onsetdetection_t *
 222 new_aubio_onsetdetection (aubio_onsetdetection_type type,
 223     uint_t size, uint_t channels){
 224   aubio_onsetdetection_t * o = AUBIO_NEW(aubio_onsetdetection_t);
 225   uint_t rsize = size/2+1;
 226   uint_t i;
 227   switch(type) {
 228     /* for both energy and hfc, only fftgrain->norm is required */
 229     case aubio_onset_energy:
 230       break;
 231     case aubio_onset_hfc:
 232       break;
 233       /* the other approaches will need some more memory spaces */
 234     case aubio_onset_complex:
 235       o->oldmag = new_fvec(rsize,channels);
 236       /** bug: must be complex array */
 237       o->meas = AUBIO_ARRAY(fft_data_t,size+1);
 238       for (i=0; i<size+1; i++) o->meas[i] = 0;
 239       o->dev1   = new_fvec(rsize,channels);
 240       o->theta1 = new_fvec(rsize,channels);
 241       o->theta2 = new_fvec(rsize,channels);
 242       break;
 243     case aubio_onset_phase:
 244       o->dev1   = new_fvec(rsize,channels);
 245       o->theta1 = new_fvec(rsize,channels);
 246       o->theta2 = new_fvec(rsize,channels);
 247       o->histog = new_aubio_hist(0.0f, PI, 10, channels);
 248       o->threshold = 0.1;
 249       break;
 250     case aubio_onset_specdiff:
 251       o->oldmag = new_fvec(rsize,channels);
 252       o->dev1   = new_fvec(rsize,channels);
 253       o->histog = new_aubio_hist(0.0f, PI, 10, channels);
 254       o->threshold = 0.1;
 255       break;
 256     case aubio_onset_kl:
 257     case aubio_onset_mkl:
 258     case aubio_onset_specflux:
 259       o->oldmag = new_fvec(rsize,channels);
 260       break;
 261     default:
 262       break;
 263   }
 264
 265   /* this switch could be in its own function to change between
 266    * detections on the fly. this would need getting rid of the switch
 267    * above and always allocate all the structure */
 268
 269   switch(type) {
 270     case aubio_onset_energy:
 271       o->funcpointer = aubio_onsetdetection_energy;
 272       break;
 273     case aubio_onset_hfc:
 274       o->funcpointer = aubio_onsetdetection_hfc;
 275       break;
 276     case aubio_onset_complex:
 277       o->funcpointer = aubio_onsetdetection_complex;
 278       break;
 279     case aubio_onset_phase:
 280       o->funcpointer = aubio_onsetdetection_phase;
 281       break;
 282     case aubio_onset_specdiff:
 283       o->funcpointer = aubio_onsetdetection_specdiff;
 284       break;
 285     case aubio_onset_kl:
 286       o->funcpointer = aubio_onsetdetection_kl;
 287       break;
 288     case aubio_onset_mkl:
 289       o->funcpointer = aubio_onsetdetection_mkl;
 290       break;
 291     case aubio_onset_specflux:
 292       o->funcpointer = aubio_onsetdetection_specflux;
 293       break;
 294     default:
 295       break;
 296   }
 297   o->type = type;
 298   return o;
 299 }
 300
 301 void del_aubio_onsetdetection (aubio_onsetdetection_t *o){
 302   switch(o->type) {
 303     /* for both energy and hfc, only fftgrain->norm is required */
 304     case aubio_onset_energy:
 305       break;
 306     case aubio_onset_hfc:
 307       break;
 308       /* the other approaches will need some more memory spaces */
 309     case aubio_onset_complex:
 310       AUBIO_FREE(o->meas);
 311       del_fvec(o->oldmag);
 312       del_fvec(o->dev1);
 313       del_fvec(o->theta1);
 314       del_fvec(o->theta2);
 315       break;
 316     case aubio_onset_phase:
 317       del_fvec(o->dev1);
 318       del_fvec(o->theta1);
 319       del_fvec(o->theta2);
 320       del_aubio_hist(o->histog);
 321       break;
 322     case aubio_onset_specdiff:
 323       del_fvec(o->oldmag);
 324       del_fvec(o->dev1);
 325       del_aubio_hist(o->histog);
 326       break;
 327     case aubio_onset_kl:
 328       del_fvec(o->oldmag);
 329       break;
 330     case aubio_onset_mkl:
 331       del_fvec(o->oldmag);
 332       break;
 333     default:
 334       break;
 335   }
 336   AUBIO_FREE(o);
 337 }