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