src/onset/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 "fvec.h"
  22 #include "cvec.h"
  23 #include "spectral/fft.h"
  24 #include "mathutils.h"
  25 #include "utils/hist.h"
  26 #include "onset/onsetdetection.h"
  27
  28
  29 /** structure to store object state */
  30 struct _aubio_onsetdetection_t {
  31   aubio_onsetdetection_type type; /**< onset detection type */
  32   /** Pointer to aubio_onsetdetection_<type> function */
  33   void (*funcpointer)(aubio_onsetdetection_t *o,
  34       cvec_t * fftgrain, fvec_t * onset);
  35   smpl_t threshold;      /**< minimum norm threshold for phase and specdiff */
  36   fvec_t *oldmag;        /**< previous norm 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       // compute the predicted phase
  77       o->dev1->data[i][j] = 2. * o->theta1->data[i][j] - o->theta2->data[i][j];
  78       // compute the euclidean distance in the complex domain
  79       // sqrt ( r_1^2 + r_2^2 - 2 * r_1 * r_2 * \cos ( \phi_1 - \phi_2 ) )
  80       onset->data[i][0] +=
  81         SQRT (ABS (SQR (o->oldmag->data[i][j]) + SQR (fftgrain->norm[i][j])
  82               - 2. * o->oldmag->data[i][j] * fftgrain->norm[i][j]
  83               * COS (o->dev1->data[i][j] - fftgrain->phas[i][j])));
  84       /* swap old phase data (need to remember 2 frames behind)*/
  85       o->theta2->data[i][j] = o->theta1->data[i][j];
  86       o->theta1->data[i][j] = fftgrain->phas[i][j];
  87       /* swap old magnitude data (1 frame is enough) */
  88       o->oldmag->data[i][j] = fftgrain->norm[i][j];
  89     }
  90   }
  91 }
  92
  93
  94 /* Phase Based Method onset detection function */
  95 void aubio_onsetdetection_phase(aubio_onsetdetection_t *o,
  96     cvec_t * fftgrain, fvec_t * onset){
  97   uint_t i, j;
  98   uint_t nbins = fftgrain->length;
  99   for (i=0;i<fftgrain->channels; i++)  {
 100     onset->data[i][0] = 0.0;
 101     o->dev1->data[i][0]=0.;
 102     for ( j=0;j<nbins; j++ )  {
 103       o->dev1->data[i][j] =
 104         aubio_unwrap2pi(
 105             fftgrain->phas[i][j]
 106             -2.0*o->theta1->data[i][j]
 107             +o->theta2->data[i][j]);
 108       if ( o->threshold < fftgrain->norm[i][j] )
 109         o->dev1->data[i][j] = ABS(o->dev1->data[i][j]);
 110       else
 111         o->dev1->data[i][j] = 0.0;
 112       /* keep a track of the past frames */
 113       o->theta2->data[i][j] = o->theta1->data[i][j];
 114       o->theta1->data[i][j] = fftgrain->phas[i][j];
 115     }
 116     /* apply o->histogram */
 117     aubio_hist_dyn_notnull(o->histog,o->dev1);
 118     /* weight it */
 119     aubio_hist_weight(o->histog);
 120     /* its mean is the result */
 121     onset->data[i][0] = aubio_hist_mean(o->histog);
 122     //onset->data[i][0] = fvec_mean(o->dev1);
 123   }
 124 }
 125
 126 /* Spectral difference method onset detection function */
 127 void aubio_onsetdetection_specdiff(aubio_onsetdetection_t *o,
 128     cvec_t * fftgrain, fvec_t * onset){
 129   uint_t i, j;
 130   uint_t nbins = fftgrain->length;
 131   for (i=0;i<fftgrain->channels; i++)  {
 132     onset->data[i][0] = 0.0;
 133     for (j=0;j<nbins; j++)  {
 134       o->dev1->data[i][j] = SQRT(
 135           ABS(SQR( fftgrain->norm[i][j])
 136             - SQR(o->oldmag->data[i][j])));
 137       if (o->threshold < fftgrain->norm[i][j] )
 138         o->dev1->data[i][j] = ABS(o->dev1->data[i][j]);
 139       else
 140         o->dev1->data[i][j] = 0.0;
 141       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 142     }
 143
 144     /* apply o->histogram (act somewhat as a low pass on the
 145      * overall function)*/
 146     aubio_hist_dyn_notnull(o->histog,o->dev1);
 147     /* weight it */
 148     aubio_hist_weight(o->histog);
 149     /* its mean is the result */
 150     onset->data[i][0] = aubio_hist_mean(o->histog);
 151
 152   }
 153 }
 154
 155 /* Kullback Liebler onset detection function
 156  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 157  * negative (1.+) and infinite values (+1.e-10) */
 158 void aubio_onsetdetection_kl(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 159   uint_t i,j;
 160   for (i=0;i<fftgrain->channels;i++) {
 161     onset->data[i][0] = 0.;
 162     for (j=0;j<fftgrain->length;j++) {
 163       onset->data[i][0] += fftgrain->norm[i][j]
 164         *LOG(1.+fftgrain->norm[i][j]/(o->oldmag->data[i][j]+1.e-10));
 165       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 166     }
 167     if (isnan(onset->data[i][0])) onset->data[i][0] = 0.;
 168   }
 169 }
 170
 171 /* Modified Kullback Liebler onset detection function
 172  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 173  * negative (1.+) and infinite values (+1.e-10) */
 174 void aubio_onsetdetection_mkl(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 175   uint_t i,j;
 176   for (i=0;i<fftgrain->channels;i++) {
 177     onset->data[i][0] = 0.;
 178     for (j=0;j<fftgrain->length;j++) {
 179       onset->data[i][0] += LOG(1.+fftgrain->norm[i][j]/(o->oldmag->data[i][j]+1.e-10));
 180       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 181     }
 182     if (isnan(onset->data[i][0])) onset->data[i][0] = 0.;
 183   }
 184 }
 185
 186 /* Spectral flux */
 187 void aubio_onsetdetection_specflux(aubio_onsetdetection_t *o, cvec_t * fftgrain, fvec_t * onset){
 188   uint_t i, j;
 189   for (i=0;i<fftgrain->channels;i++) {
 190     onset->data[i][0] = 0.;
 191     for (j=0;j<fftgrain->length;j++) {
 192       if (fftgrain->norm[i][j] > o->oldmag->data[i][j])
 193         onset->data[i][0] += fftgrain->norm[i][j] - o->oldmag->data[i][j];
 194       o->oldmag->data[i][j] = fftgrain->norm[i][j];
 195     }
 196   }
 197 }
 198
 199 /* Generic function pointing to the choosen one */
 200 void
 201 aubio_onsetdetection(aubio_onsetdetection_t *o, cvec_t * fftgrain,
 202     fvec_t * onset) {
 203   o->funcpointer(o,fftgrain,onset);
 204 }
 205
 206 /* Allocate memory for an onset detection
 207  * depending on the choosen type, allocate memory as needed
 208  */
 209 aubio_onsetdetection_t *
 210 new_aubio_onsetdetection (aubio_onsetdetection_type type,
 211     uint_t size, uint_t channels){
 212   aubio_onsetdetection_t * o = AUBIO_NEW(aubio_onsetdetection_t);
 213   uint_t rsize = size/2+1;
 214   uint_t i;
 215   switch(type) {
 216     /* for both energy and hfc, only fftgrain->norm is required */
 217     case aubio_onset_energy:
 218       break;
 219     case aubio_onset_hfc:
 220       break;
 221       /* the other approaches will need some more memory spaces */
 222     case aubio_onset_complex:
 223       o->oldmag = new_fvec(rsize,channels);
 224       o->dev1   = new_fvec(rsize,channels);
 225       o->theta1 = new_fvec(rsize,channels);
 226       o->theta2 = new_fvec(rsize,channels);
 227       break;
 228     case aubio_onset_phase:
 229       o->dev1   = new_fvec(rsize,channels);
 230       o->theta1 = new_fvec(rsize,channels);
 231       o->theta2 = new_fvec(rsize,channels);
 232       o->histog = new_aubio_hist(0.0, PI, 10, channels);
 233       o->threshold = 0.1;
 234       break;
 235     case aubio_onset_specdiff:
 236       o->oldmag = new_fvec(rsize,channels);
 237       o->dev1   = new_fvec(rsize,channels);
 238       o->histog = new_aubio_hist(0.0, PI, 10, channels);
 239       o->threshold = 0.1;
 240       break;
 241     case aubio_onset_kl:
 242     case aubio_onset_mkl:
 243     case aubio_onset_specflux:
 244       o->oldmag = new_fvec(rsize,channels);
 245       break;
 246     default:
 247       break;
 248   }
 249
 250   /* this switch could be in its own function to change between
 251    * detections on the fly. this would need getting rid of the switch
 252    * above and always allocate all the structure */
 253
 254   switch(type) {
 255     case aubio_onset_energy:
 256       o->funcpointer = aubio_onsetdetection_energy;
 257       break;
 258     case aubio_onset_hfc:
 259       o->funcpointer = aubio_onsetdetection_hfc;
 260       break;
 261     case aubio_onset_complex:
 262       o->funcpointer = aubio_onsetdetection_complex;
 263       break;
 264     case aubio_onset_phase:
 265       o->funcpointer = aubio_onsetdetection_phase;
 266       break;
 267     case aubio_onset_specdiff:
 268       o->funcpointer = aubio_onsetdetection_specdiff;
 269       break;
 270     case aubio_onset_kl:
 271       o->funcpointer = aubio_onsetdetection_kl;
 272       break;
 273     case aubio_onset_mkl:
 274       o->funcpointer = aubio_onsetdetection_mkl;
 275       break;
 276     case aubio_onset_specflux:
 277       o->funcpointer = aubio_onsetdetection_specflux;
 278       break;
 279     default:
 280       break;
 281   }
 282   o->type = type;
 283   return o;
 284 }
 285
 286 void del_aubio_onsetdetection (aubio_onsetdetection_t *o){
 287   switch(o->type) {
 288     /* for both energy and hfc, only fftgrain->norm is required */
 289     case aubio_onset_energy:
 290       break;
 291     case aubio_onset_hfc:
 292       break;
 293       /* the other approaches will need some more memory spaces */
 294     case aubio_onset_complex:
 295       del_fvec(o->oldmag);
 296       del_fvec(o->dev1);
 297       del_fvec(o->theta1);
 298       del_fvec(o->theta2);
 299       break;
 300     case aubio_onset_phase:
 301       del_fvec(o->dev1);
 302       del_fvec(o->theta1);
 303       del_fvec(o->theta2);
 304       del_aubio_hist(o->histog);
 305       break;
 306     case aubio_onset_specdiff:
 307       del_fvec(o->oldmag);
 308       del_fvec(o->dev1);
 309       del_aubio_hist(o->histog);
 310       break;
 311     case aubio_onset_kl:
 312     case aubio_onset_mkl:
 313     case aubio_onset_specflux:
 314       del_fvec(o->oldmag);
 315       break;
 316     default:
 317       break;
 318   }
 319   AUBIO_FREE(o);
 320 }