src/spectral/specdesc.c

   1 /*
   2   Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
   3
   4   This file is part of aubio.
   5
   6   aubio 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 3 of the License, or
   9   (at your option) any later version.
  10
  11   aubio 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.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with aubio.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 */
  20
  21 #include "aubio_priv.h"
  22 #include "fvec.h"
  23 #include "cvec.h"
  24 #include "spectral/fft.h"
  25 #include "spectral/specdesc.h"
  26 #include "mathutils.h"
  27 #include "utils/hist.h"
  28
  29 void aubio_specdesc_energy(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  30 void aubio_specdesc_hfc(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  31 void aubio_specdesc_complex(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  32 void aubio_specdesc_phase(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  33 void aubio_specdesc_wphase(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  34 void aubio_specdesc_specdiff(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  35 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  36 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  37 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset);
  38
  39 extern void aubio_specdesc_centroid (aubio_specdesc_t * o, const cvec_t * spec,
  40     fvec_t * desc);
  41 extern void aubio_specdesc_spread (aubio_specdesc_t * o, const cvec_t * spec,
  42     fvec_t * desc);
  43 extern void aubio_specdesc_skewness (aubio_specdesc_t * o, const cvec_t * spec,
  44     fvec_t * desc);
  45 extern void aubio_specdesc_kurtosis (aubio_specdesc_t * o, const cvec_t * spec,
  46     fvec_t * desc);
  47 extern void aubio_specdesc_slope (aubio_specdesc_t * o, const cvec_t * spec,
  48     fvec_t * desc);
  49 extern void aubio_specdesc_decrease (aubio_specdesc_t * o, const cvec_t * spec,
  50     fvec_t * desc);
  51 extern void aubio_specdesc_rolloff (aubio_specdesc_t * o, const cvec_t * spec,
  52     fvec_t * desc);
  53
  54 /** onsetdetection types */
  55 typedef enum {
  56         aubio_onset_energy,         /**< energy based */
  57         aubio_onset_specdiff,       /**< spectral diff */
  58         aubio_onset_hfc,            /**< high frequency content */
  59         aubio_onset_complex,        /**< complex domain */
  60         aubio_onset_phase,          /**< phase fast */
  61         aubio_onset_wphase,         /**< weighted phase */
  62         aubio_onset_kl,             /**< Kullback Liebler */
  63         aubio_onset_mkl,            /**< modified Kullback Liebler */
  64         aubio_onset_specflux,       /**< spectral flux */
  65         aubio_specmethod_centroid,  /**< spectral centroid */
  66         aubio_specmethod_spread,    /**< spectral spread */
  67         aubio_specmethod_skewness,  /**< spectral skewness */
  68         aubio_specmethod_kurtosis,  /**< spectral kurtosis */
  69         aubio_specmethod_slope,     /**< spectral kurtosis */
  70         aubio_specmethod_decrease,  /**< spectral decrease */
  71         aubio_specmethod_rolloff,   /**< spectral rolloff */
  72         aubio_onset_default = aubio_onset_hfc, /**< default mode, set to hfc */
  73 } aubio_specdesc_type;
  74
  75 /** structure to store object state */
  76 struct _aubio_specdesc_t {
  77   aubio_specdesc_type onset_type; /**< onset detection type */
  78   /** Pointer to aubio_specdesc_<type> function */
  79   void (*funcpointer)(aubio_specdesc_t *o,
  80       const cvec_t * fftgrain, fvec_t * onset);
  81   smpl_t threshold;      /**< minimum norm threshold for phase and specdiff */
  82   fvec_t *oldmag;        /**< previous norm vector */
  83   fvec_t *dev1 ;         /**< current onset detection measure vector */
  84   fvec_t *theta1;        /**< previous phase vector, one frame behind */
  85   fvec_t *theta2;        /**< previous phase vector, two frames behind */
  86   aubio_hist_t * histog; /**< histogram */
  87 };
  88
  89
  90 /* Energy based onset detection function */
  91 void aubio_specdesc_energy  (aubio_specdesc_t *o UNUSED,
  92     const cvec_t * fftgrain, fvec_t * onset) {
  93   uint_t j;
  94   onset->data[0] = 0.;
  95   for (j=0;j<fftgrain->length;j++) {
  96     onset->data[0] += SQR(fftgrain->norm[j]);
  97   }
  98 }
  99
 100 /* High Frequency Content onset detection function */
 101 void aubio_specdesc_hfc(aubio_specdesc_t *o UNUSED,
 102     const cvec_t * fftgrain, fvec_t * onset){
 103   uint_t j;
 104   onset->data[0] = 0.;
 105   for (j=0;j<fftgrain->length;j++) {
 106     onset->data[0] += (j+1)*fftgrain->norm[j];
 107   }
 108 }
 109
 110
 111 /* Complex Domain Method onset detection function */
 112 void aubio_specdesc_complex (aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset) {
 113   uint_t j;
 114   uint_t nbins = fftgrain->length;
 115   onset->data[0] = 0.;
 116   for (j=0;j<nbins; j++)  {
 117     // compute the predicted phase
 118     o->dev1->data[j] = 2. * o->theta1->data[j] - o->theta2->data[j];
 119     // compute the euclidean distance in the complex domain
 120     // sqrt ( r_1^2 + r_2^2 - 2 * r_1 * r_2 * \cos ( \phi_1 - \phi_2 ) )
 121     onset->data[0] +=
 122       SQRT (ABS (SQR (o->oldmag->data[j]) + SQR (fftgrain->norm[j])
 123             - 2 * o->oldmag->data[j] * fftgrain->norm[j]
 124             * COS (o->dev1->data[j] - fftgrain->phas[j])));
 125     /* swap old phase data (need to remember 2 frames behind)*/
 126     o->theta2->data[j] = o->theta1->data[j];
 127     o->theta1->data[j] = fftgrain->phas[j];
 128     /* swap old magnitude data (1 frame is enough) */
 129     o->oldmag->data[j] = fftgrain->norm[j];
 130   }
 131 }
 132
 133
 134 /* Phase Based Method onset detection function */
 135 void aubio_specdesc_phase(aubio_specdesc_t *o,
 136     const cvec_t * fftgrain, fvec_t * onset){
 137   uint_t j;
 138   uint_t nbins = fftgrain->length;
 139   onset->data[0] = 0.0;
 140   o->dev1->data[0]=0.;
 141   for ( j=0;j<nbins; j++ )  {
 142     o->dev1->data[j] =
 143       aubio_unwrap2pi(
 144           fftgrain->phas[j]
 145           -2.0*o->theta1->data[j]
 146           +o->theta2->data[j]);
 147     if ( o->threshold < fftgrain->norm[j] )
 148       o->dev1->data[j] = ABS(o->dev1->data[j]);
 149     else
 150       o->dev1->data[j] = 0.0;
 151     /* keep a track of the past frames */
 152     o->theta2->data[j] = o->theta1->data[j];
 153     o->theta1->data[j] = fftgrain->phas[j];
 154   }
 155   /* apply o->histogram */
 156   aubio_hist_dyn_notnull(o->histog,o->dev1);
 157   /* weight it */
 158   aubio_hist_weight(o->histog);
 159   /* its mean is the result */
 160   onset->data[0] = aubio_hist_mean(o->histog);
 161   //onset->data[0] = fvec_mean(o->dev1);
 162 }
 163
 164 /* weighted phase */
 165 void
 166 aubio_specdesc_wphase(aubio_specdesc_t *o,
 167     const cvec_t *fftgrain, fvec_t *onset) {
 168   uint_t i;
 169   aubio_specdesc_phase(o, fftgrain, onset);
 170   for (i = 0; i < fftgrain->length; i++) {
 171     o->dev1->data[i] *= fftgrain->norm[i];
 172   }
 173   /* apply o->histogram */
 174   aubio_hist_dyn_notnull(o->histog,o->dev1);
 175   /* weight it */
 176   aubio_hist_weight(o->histog);
 177   /* its mean is the result */
 178   onset->data[0] = aubio_hist_mean(o->histog);
 179 }
 180
 181 /* Spectral difference method onset detection function */
 182 void aubio_specdesc_specdiff(aubio_specdesc_t *o,
 183     const cvec_t * fftgrain, fvec_t * onset){
 184   uint_t j;
 185   uint_t nbins = fftgrain->length;
 186     onset->data[0] = 0.0;
 187     for (j=0;j<nbins; j++)  {
 188       o->dev1->data[j] = SQRT(
 189           ABS(SQR( fftgrain->norm[j])
 190             - SQR(o->oldmag->data[j])));
 191       if (o->threshold < fftgrain->norm[j] )
 192         o->dev1->data[j] = ABS(o->dev1->data[j]);
 193       else
 194         o->dev1->data[j] = 0.0;
 195       o->oldmag->data[j] = fftgrain->norm[j];
 196     }
 197
 198     /* apply o->histogram (act somewhat as a low pass on the
 199      * overall function)*/
 200     aubio_hist_dyn_notnull(o->histog,o->dev1);
 201     /* weight it */
 202     aubio_hist_weight(o->histog);
 203     /* its mean is the result */
 204     onset->data[0] = aubio_hist_mean(o->histog);
 205 }
 206
 207 /* Kullback Liebler onset detection function
 208  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 209  * negative (1.+) and infinite values (+1.e-10) */
 210 void aubio_specdesc_kl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
 211   uint_t j;
 212     onset->data[0] = 0.;
 213     for (j=0;j<fftgrain->length;j++) {
 214       onset->data[0] += fftgrain->norm[j]
 215         *LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
 216       o->oldmag->data[j] = fftgrain->norm[j];
 217     }
 218     if (isnan(onset->data[0])) onset->data[0] = 0.;
 219 }
 220
 221 /* Modified Kullback Liebler onset detection function
 222  * note we use ln(1+Xn/(Xn-1+0.0001)) to avoid
 223  * negative (1.+) and infinite values (+1.e-10) */
 224 void aubio_specdesc_mkl(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
 225   uint_t j;
 226     onset->data[0] = 0.;
 227     for (j=0;j<fftgrain->length;j++) {
 228       onset->data[0] += LOG(1.+fftgrain->norm[j]/(o->oldmag->data[j]+1.e-1));
 229       o->oldmag->data[j] = fftgrain->norm[j];
 230     }
 231     if (isnan(onset->data[0])) onset->data[0] = 0.;
 232 }
 233
 234 /* Spectral flux */
 235 void aubio_specdesc_specflux(aubio_specdesc_t *o, const cvec_t * fftgrain, fvec_t * onset){
 236   uint_t j;
 237   onset->data[0] = 0.;
 238   for (j=0;j<fftgrain->length;j++) {
 239     if (fftgrain->norm[j] > o->oldmag->data[j])
 240       onset->data[0] += fftgrain->norm[j] - o->oldmag->data[j];
 241     o->oldmag->data[j] = fftgrain->norm[j];
 242   }
 243 }
 244
 245 /* Generic function pointing to the choosen one */
 246 void
 247 aubio_specdesc_do (aubio_specdesc_t *o, const cvec_t * fftgrain,
 248     fvec_t * onset) {
 249   o->funcpointer(o,fftgrain,onset);
 250 }
 251
 252 /* Allocate memory for an onset detection
 253  * depending on the choosen type, allocate memory as needed
 254  */
 255 aubio_specdesc_t *
 256 new_aubio_specdesc (const char_t * onset_mode, uint_t size){
 257   aubio_specdesc_t * o = AUBIO_NEW(aubio_specdesc_t);
 258   uint_t rsize = size/2+1;
 259   aubio_specdesc_type onset_type;
 260   if (strcmp (onset_mode, "energy") == 0)
 261       onset_type = aubio_onset_energy;
 262   else if (strcmp (onset_mode, "specdiff") == 0)
 263       onset_type = aubio_onset_specdiff;
 264   else if (strcmp (onset_mode, "hfc") == 0)
 265       onset_type = aubio_onset_hfc;
 266   else if (strcmp (onset_mode, "complexdomain") == 0)
 267       onset_type = aubio_onset_complex;
 268   else if (strcmp (onset_mode, "complex") == 0)
 269       onset_type = aubio_onset_complex;
 270   else if (strcmp (onset_mode, "phase") == 0)
 271       onset_type = aubio_onset_phase;
 272   else if (strcmp (onset_mode, "wphase") == 0)
 273       onset_type = aubio_onset_wphase;
 274   else if (strcmp (onset_mode, "mkl") == 0)
 275       onset_type = aubio_onset_mkl;
 276   else if (strcmp (onset_mode, "kl") == 0)
 277       onset_type = aubio_onset_kl;
 278   else if (strcmp (onset_mode, "specflux") == 0)
 279       onset_type = aubio_onset_specflux;
 280   else if (strcmp (onset_mode, "centroid") == 0)
 281       onset_type = aubio_specmethod_centroid;
 282   else if (strcmp (onset_mode, "spread") == 0)
 283       onset_type = aubio_specmethod_spread;
 284   else if (strcmp (onset_mode, "skewness") == 0)
 285       onset_type = aubio_specmethod_skewness;
 286   else if (strcmp (onset_mode, "kurtosis") == 0)
 287       onset_type = aubio_specmethod_kurtosis;
 288   else if (strcmp (onset_mode, "slope") == 0)
 289       onset_type = aubio_specmethod_slope;
 290   else if (strcmp (onset_mode, "decrease") == 0)
 291       onset_type = aubio_specmethod_decrease;
 292   else if (strcmp (onset_mode, "rolloff") == 0)
 293       onset_type = aubio_specmethod_rolloff;
 294   else if (strcmp (onset_mode, "old_default") == 0)
 295       onset_type = aubio_onset_default;
 296   else if (strcmp (onset_mode, "default") == 0)
 297       onset_type = aubio_onset_default;
 298   else {
 299       AUBIO_ERR("specdesc: unknown spectral descriptor type '%s'\n",
 300           onset_mode);
 301       AUBIO_FREE(o);
 302       return NULL;
 303   }
 304   switch(onset_type) {
 305     /* for both energy and hfc, only fftgrain->norm is required */
 306     case aubio_onset_energy:
 307       break;
 308     case aubio_onset_hfc:
 309       break;
 310       /* the other approaches will need some more memory spaces */
 311     case aubio_onset_complex:
 312       o->oldmag = new_fvec(rsize);
 313       o->dev1   = new_fvec(rsize);
 314       o->theta1 = new_fvec(rsize);
 315       o->theta2 = new_fvec(rsize);
 316       break;
 317     case aubio_onset_phase:
 318     case aubio_onset_wphase:
 319       o->dev1   = new_fvec(rsize);
 320       o->theta1 = new_fvec(rsize);
 321       o->theta2 = new_fvec(rsize);
 322       o->histog = new_aubio_hist(0.0, PI, 10);
 323       o->threshold = 0.1;
 324       break;
 325     case aubio_onset_specdiff:
 326       o->oldmag = new_fvec(rsize);
 327       o->dev1   = new_fvec(rsize);
 328       o->histog = new_aubio_hist(0.0, PI, 10);
 329       o->threshold = 0.1;
 330       break;
 331     case aubio_onset_kl:
 332     case aubio_onset_mkl:
 333     case aubio_onset_specflux:
 334       o->oldmag = new_fvec(rsize);
 335       break;
 336     default:
 337       break;
 338   }
 339
 340   switch(onset_type) {
 341     case aubio_onset_energy:
 342       o->funcpointer = aubio_specdesc_energy;
 343       break;
 344     case aubio_onset_hfc:
 345       o->funcpointer = aubio_specdesc_hfc;
 346       break;
 347     case aubio_onset_complex:
 348       o->funcpointer = aubio_specdesc_complex;
 349       break;
 350     case aubio_onset_phase:
 351       o->funcpointer = aubio_specdesc_phase;
 352       break;
 353     case aubio_onset_wphase:
 354       o->funcpointer = aubio_specdesc_wphase;
 355       break;
 356     case aubio_onset_specdiff:
 357       o->funcpointer = aubio_specdesc_specdiff;
 358       break;
 359     case aubio_onset_kl:
 360       o->funcpointer = aubio_specdesc_kl;
 361       break;
 362     case aubio_onset_mkl:
 363       o->funcpointer = aubio_specdesc_mkl;
 364       break;
 365     case aubio_onset_specflux:
 366       o->funcpointer = aubio_specdesc_specflux;
 367       break;
 368     case aubio_specmethod_centroid:
 369       o->funcpointer = aubio_specdesc_centroid;
 370       break;
 371     case aubio_specmethod_spread:
 372       o->funcpointer = aubio_specdesc_spread;
 373       break;
 374     case aubio_specmethod_skewness:
 375       o->funcpointer = aubio_specdesc_skewness;
 376       break;
 377     case aubio_specmethod_kurtosis:
 378       o->funcpointer = aubio_specdesc_kurtosis;
 379       break;
 380     case aubio_specmethod_slope:
 381       o->funcpointer = aubio_specdesc_slope;
 382       break;
 383     case aubio_specmethod_decrease:
 384       o->funcpointer = aubio_specdesc_decrease;
 385       break;
 386     case aubio_specmethod_rolloff:
 387       o->funcpointer = aubio_specdesc_rolloff;
 388       break;
 389     default:
 390       break;
 391   }
 392   o->onset_type = onset_type;
 393   return o;
 394 }
 395
 396 void del_aubio_specdesc (aubio_specdesc_t *o){
 397   switch(o->onset_type) {
 398     case aubio_onset_energy:
 399       break;
 400     case aubio_onset_hfc:
 401       break;
 402     case aubio_onset_complex:
 403       del_fvec(o->oldmag);
 404       del_fvec(o->dev1);
 405       del_fvec(o->theta1);
 406       del_fvec(o->theta2);
 407       break;
 408     case aubio_onset_phase:
 409     case aubio_onset_wphase:
 410       del_fvec(o->dev1);
 411       del_fvec(o->theta1);
 412       del_fvec(o->theta2);
 413       del_aubio_hist(o->histog);
 414       break;
 415     case aubio_onset_specdiff:
 416       del_fvec(o->oldmag);
 417       del_fvec(o->dev1);
 418       del_aubio_hist(o->histog);
 419       break;
 420     case aubio_onset_kl:
 421     case aubio_onset_mkl:
 422     case aubio_onset_specflux:
 423       del_fvec(o->oldmag);
 424       break;
 425     default:
 426       break;
 427   }
 428   AUBIO_FREE(o);
 429 }