src/mathutils.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 /* see in mathutils.h for doc */
  21
  22 #include "aubio_priv.h"
  23 #include "fvec.h"
  24 #include "mathutils.h"
  25 #include "config.h"
  26
  27 fvec_t * new_aubio_window(uint_t size, aubio_window_type wintype) {
  28   // create fvec of size x 1 channel
  29   fvec_t * win = new_fvec( size, 1);
  30   smpl_t * w = win->data[0];
  31   uint_t i;
  32   switch(wintype) {
  33     case aubio_win_rectangle:
  34       for (i=0;i<size;i++)
  35         w[i] = 0.5;
  36       break;
  37     case aubio_win_hamming:
  38       for (i=0;i<size;i++)
  39         w[i] = 0.54 - 0.46 * COS(TWO_PI * i / (size));
  40       break;
  41     case aubio_win_hanning:
  42       for (i=0;i<size;i++)
  43         w[i] = 0.5 - (0.5 * COS(TWO_PI * i / (size)));
  44       break;
  45     case aubio_win_hanningz:
  46       for (i=0;i<size;i++)
  47         w[i] = 0.5 * (1.0 - COS(TWO_PI * i / (size)));
  48       break;
  49     case aubio_win_blackman:
  50       for (i=0;i<size;i++)
  51         w[i] = 0.42
  52           - 0.50 * COS(    TWO_PI*i/(size-1.0))
  53           + 0.08 * COS(2.0*TWO_PI*i/(size-1.0));
  54       break;
  55     case aubio_win_blackman_harris:
  56       for (i=0;i<size;i++)
  57         w[i] = 0.35875
  58           - 0.48829 * COS(    TWO_PI*i/(size-1.0))
  59           + 0.14128 * COS(2.0*TWO_PI*i/(size-1.0))
  60           - 0.01168 * COS(3.0*TWO_PI*i/(size-1.0));
  61       break;
  62     case aubio_win_gaussian:
  63       for (i=0;i<size;i++)
  64         w[i] = EXP(- 1.0 / SQR(size) * SQR(2.0*i-size));
  65       break;
  66     case aubio_win_welch:
  67       for (i=0;i<size;i++)
  68         w[i] = 1.0 - SQR((2*i-size)/(size+1.0));
  69       break;
  70     case aubio_win_parzen:
  71       for (i=0;i<size;i++)
  72         w[i] = 1.0 - ABS((2*i-size)/(size+1.0));
  73       break;
  74     default:
  75       break;
  76   }
  77   return win;
  78 }
  79
  80 smpl_t aubio_unwrap2pi(smpl_t phase) {
  81   /* mod(phase+pi,-2pi)+pi */
  82   return phase + TWO_PI * (1. + FLOOR(-(phase+PI)/TWO_PI));
  83 }
  84
  85 smpl_t vec_mean(fvec_t *s) {
  86   uint_t i,j;
  87   smpl_t tmp = 0.0f;
  88   for (i=0; i < s->channels; i++)
  89     for (j=0; j < s->length; j++)
  90       tmp += s->data[i][j];
  91   return tmp/(smpl_t)(s->length);
  92 }
  93
  94 smpl_t vec_sum(fvec_t *s) {
  95   uint_t i,j;
  96   smpl_t tmp = 0.0f;
  97   for (i=0; i < s->channels; i++)
  98     for (j=0; j < s->length; j++)
  99       tmp += s->data[i][j];
 100   return tmp;
 101 }
 102
 103 smpl_t vec_max(fvec_t *s) {
 104   uint_t i,j;
 105   smpl_t tmp = 0.0f;
 106   for (i=0; i < s->channels; i++)
 107     for (j=0; j < s->length; j++)
 108       tmp = (tmp > s->data[i][j])? tmp : s->data[i][j];
 109   return tmp;
 110 }
 111
 112 smpl_t vec_min(fvec_t *s) {
 113   uint_t i,j;
 114   smpl_t tmp = s->data[0][0];
 115   for (i=0; i < s->channels; i++)
 116     for (j=0; j < s->length; j++)
 117       tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
 118   return tmp;
 119 }
 120
 121 uint_t vec_min_elem(fvec_t *s) {
 122   uint_t i,j=0, pos=0.;
 123   smpl_t tmp = s->data[0][0];
 124   for (i=0; i < s->channels; i++)
 125     for (j=0; j < s->length; j++) {
 126       pos = (tmp < s->data[i][j])? pos : j;
 127       tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
 128     }
 129   return pos;
 130 }
 131
 132 uint_t vec_max_elem(fvec_t *s) {
 133   uint_t i,j=0, pos=0.;
 134   smpl_t tmp = 0.0f;
 135   for (i=0; i < s->channels; i++)
 136     for (j=0; j < s->length; j++) {
 137       pos = (tmp > s->data[i][j])? pos : j;
 138       tmp = (tmp > s->data[i][j])? tmp : s->data[i][j]  ;
 139     }
 140   return pos;
 141 }
 142
 143 void vec_shift(fvec_t *s) {
 144   uint_t i,j;
 145   //smpl_t tmp = 0.0f;
 146   for (i=0; i < s->channels; i++)
 147     for (j=0; j < s->length / 2 ; j++) {
 148       //tmp = s->data[i][j];
 149       //s->data[i][j] = s->data[i][j+s->length/2];
 150       //s->data[i][j+s->length/2] = tmp;
 151       ELEM_SWAP(s->data[i][j],s->data[i][j+s->length/2]);
 152     }
 153 }
 154
 155 smpl_t vec_local_energy(fvec_t * f) {
 156   smpl_t locE = 0.;
 157   uint_t i,j;
 158   for (i=0;i<f->channels;i++)
 159     for (j=0;j<f->length;j++)
 160       locE+=SQR(f->data[i][j]);
 161   return locE;
 162 }
 163
 164 smpl_t vec_local_hfc(fvec_t * f) {
 165   smpl_t locE = 0.;
 166   uint_t i,j;
 167   for (i=0;i<f->channels;i++)
 168     for (j=0;j<f->length;j++)
 169       locE+=(i+1)*f->data[i][j];
 170   return locE;
 171 }
 172
 173 smpl_t vec_alpha_norm(fvec_t * DF, smpl_t alpha) {
 174   smpl_t tmp = 0.;
 175   uint_t i,j;
 176   for (i=0;i<DF->channels;i++)
 177     for (j=0;j<DF->length;j++)
 178       tmp += POW(ABS(DF->data[i][j]),alpha);
 179   return POW(tmp/DF->length,1./alpha);
 180 }
 181
 182 void vec_dc_removal(fvec_t * mag) {
 183     smpl_t mini = 0.;
 184     uint_t length = mag->length, i=0, j;
 185     mini = vec_min(mag);
 186     for (j=0;j<length;j++) {
 187       mag->data[i][j] -= mini;
 188     }
 189 }
 190
 191 void vec_alpha_normalise(fvec_t * mag, uint_t alpha) {
 192   smpl_t alphan = 1.;
 193   uint_t length = mag->length, i=0, j;
 194   alphan = vec_alpha_norm(mag,alpha);
 195   for (j=0;j<length;j++){
 196     mag->data[i][j] /= alphan;
 197   }
 198 }
 199
 200 void vec_add(fvec_t * mag, smpl_t threshold) {
 201   uint_t length = mag->length, i=0, j;
 202   for (j=0;j<length;j++) {
 203     mag->data[i][j] += threshold;
 204   }
 205 }
 206
 207 void vec_adapt_thres(fvec_t * vec, fvec_t * tmp,
 208     uint_t post, uint_t pre) {
 209   uint_t length = vec->length, i=0, j;
 210   for (j=0;j<length;j++) {
 211     vec->data[i][j] -= vec_moving_thres(vec, tmp, post, pre, j);
 212   }
 213 }
 214
 215 smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmpvec,
 216     uint_t post, uint_t pre, uint_t pos) {
 217   smpl_t * medar = (smpl_t *)tmpvec->data[0];
 218   uint_t k;
 219   uint_t win_length =  post+pre+1;
 220   uint_t length =  vec->length;
 221   /* post part of the buffer does not exist */
 222   if (pos<post+1) {
 223     for (k=0;k<post+1-pos;k++)
 224       medar[k] = 0.; /* 0-padding at the beginning */
 225     for (k=post+1-pos;k<win_length;k++)
 226       medar[k] = vec->data[0][k+pos-post];
 227   /* the buffer is fully defined */
 228   } else if (pos+pre<length) {
 229     for (k=0;k<win_length;k++)
 230       medar[k] = vec->data[0][k+pos-post];
 231   /* pre part of the buffer does not exist */
 232   } else {
 233     for (k=0;k<length-pos+post;k++)
 234       medar[k] = vec->data[0][k+pos-post];
 235     for (k=length-pos+post;k<win_length;k++)
 236       medar[k] = 0.; /* 0-padding at the end */
 237   }
 238   return vec_median(tmpvec);
 239 }
 240
 241 smpl_t vec_median(fvec_t * input) {
 242   uint_t n = input->length;
 243   smpl_t * arr = (smpl_t *) input->data[0];
 244   uint_t low, high ;
 245   uint_t median;
 246   uint_t middle, ll, hh;
 247
 248   low = 0 ; high = n-1 ; median = (low + high) / 2;
 249   for (;;) {
 250     if (high <= low) /* One element only */
 251       return arr[median] ;
 252
 253     if (high == low + 1) {  /* Two elements only */
 254       if (arr[low] > arr[high])
 255         ELEM_SWAP(arr[low], arr[high]) ;
 256       return arr[median] ;
 257     }
 258
 259     /* Find median of low, middle and high items; swap into position low */
 260     middle = (low + high) / 2;
 261     if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]);
 262     if (arr[low]    > arr[high])    ELEM_SWAP(arr[low],    arr[high]);
 263     if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
 264
 265     /* Swap low item (now in position middle) into position (low+1) */
 266     ELEM_SWAP(arr[middle], arr[low+1]) ;
 267
 268     /* Nibble from each end towards middle, swapping items when stuck */
 269     ll = low + 1;
 270     hh = high;
 271     for (;;) {
 272       do ll++; while (arr[low] > arr[ll]) ;
 273       do hh--; while (arr[hh]  > arr[low]) ;
 274
 275       if (hh < ll)
 276         break;
 277
 278       ELEM_SWAP(arr[ll], arr[hh]) ;
 279     }
 280
 281     /* Swap middle item (in position low) back into correct position */
 282     ELEM_SWAP(arr[low], arr[hh]) ;
 283
 284     /* Re-set active partition */
 285     if (hh <= median)
 286       low = ll;
 287     if (hh >= median)
 288       high = hh - 1;
 289   }
 290 }
 291
 292 smpl_t vec_quadint(fvec_t * x,uint_t pos, uint_t span) {
 293   smpl_t s0, s1, s2;
 294   uint_t x0 = (pos < span) ? pos : pos - span;
 295   uint_t x2 = (pos + span < x->length) ? pos + span : pos;
 296   if (x0 == pos) return (x->data[0][pos] <= x->data[0][x2]) ? pos : x2;
 297   if (x2 == pos) return (x->data[0][pos] <= x->data[0][x0]) ? pos : x0;
 298   s0 = x->data[0][x0];
 299   s1 = x->data[0][pos]     ;
 300   s2 = x->data[0][x2];
 301   return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
 302 }
 303
 304 smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
 305   smpl_t tmp = s0 + (pf/2.) * (pf * ( s0 - 2.*s1 + s2 ) - 3.*s0 + 4.*s1 - s2);
 306   return tmp;
 307 }
 308
 309 uint_t vec_peakpick(fvec_t * onset, uint_t pos) {
 310   uint_t i=0, tmp=0;
 311   /*for (i=0;i<onset->channels;i++)*/
 312   tmp = (onset->data[i][pos] > onset->data[i][pos-1]
 313       &&  onset->data[i][pos] > onset->data[i][pos+1]
 314       &&  onset->data[i][pos] > 0.);
 315   return tmp;
 316 }
 317
 318 smpl_t aubio_freqtomidi(smpl_t freq) {
 319   /* log(freq/A-2)/log(2) */
 320   smpl_t midi = freq/6.875;
 321   midi = LOG(midi)/0.69314718055995;
 322   midi *= 12;
 323   midi -= 3;
 324   return midi;
 325 }
 326
 327 smpl_t aubio_miditofreq(smpl_t midi) {
 328   smpl_t freq = (midi+3.)/12.;
 329   freq = EXP(freq*0.69314718055995);
 330   freq *= 6.875;
 331   return freq;
 332 }
 333
 334 smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
 335   smpl_t freq = samplerate/fftsize;
 336   return freq*bin;
 337 }
 338
 339 smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
 340   smpl_t midi = aubio_bintofreq(bin,samplerate,fftsize);
 341   return aubio_freqtomidi(midi);
 342 }
 343
 344 smpl_t aubio_freqtobin(smpl_t freq, smpl_t samplerate, smpl_t fftsize) {
 345   smpl_t bin = fftsize/samplerate;
 346   return freq*bin;
 347 }
 348
 349 smpl_t aubio_miditobin(smpl_t midi, smpl_t samplerate, smpl_t fftsize) {
 350   smpl_t freq = aubio_miditofreq(midi);
 351   return aubio_freqtobin(freq,samplerate,fftsize);
 352 }
 353
 354 /** returns 1 if wassilence is 0 and RMS(ibuf)<threshold
 355  * \bug mono
 356  */
 357 uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold) {
 358   smpl_t loudness = 0;
 359   uint_t i=0,j;
 360   for (j=0;j<ibuf->length;j++) {
 361     loudness += SQR(ibuf->data[i][j]);
 362   }
 363   loudness = SQRT(loudness);
 364   loudness /= (smpl_t)ibuf->length;
 365   loudness = LIN2DB(loudness);
 366
 367   return (loudness < threshold);
 368 }
 369
 370 /** returns level log(RMS(ibuf)) if < threshold, 1 otherwise
 371  * \bug mono
 372  */
 373 smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold) {
 374   smpl_t loudness = 0;
 375   uint_t i=0,j;
 376   for (j=0;j<ibuf->length;j++) {
 377     loudness += SQR(ibuf->data[i][j]);
 378   }
 379   loudness = SQRT(loudness);
 380   loudness /= (smpl_t)ibuf->length;
 381   loudness = LIN2DB(loudness);
 382
 383   if (loudness < threshold)
 384     return 1.;
 385   else
 386     return loudness;
 387 }
 388
 389 smpl_t aubio_zero_crossing_rate(fvec_t * input) {
 390   uint_t i=0,j;
 391   uint_t zcr = 0;
 392   for ( j = 1; j < input->length; j++ ) {
 393     // previous was strictly negative
 394     if( input->data[i][j-1] < 0. ) {
 395       // current is positive or null
 396       if ( input->data[i][j] >= 0. ) {
 397         zcr += 1;
 398       }
 399     // previous was positive or null
 400     } else {
 401       // current is strictly negative
 402       if ( input->data[i][j] < 0. ) {
 403         zcr += 1;
 404       }
 405     }
 406   }
 407   return zcr/(smpl_t)input->length;
 408 }
 409
 410 void aubio_autocorr(fvec_t * input, fvec_t * output) {
 411   uint_t i = 0, j = 0, length = input->length;
 412   smpl_t * data = input->data[0];
 413   smpl_t * acf = output->data[0];
 414   smpl_t tmp =0.;
 415   for(i=0;i<length;i++){
 416     for(j=i;j<length;j++){
 417       tmp += data[j-i]*data[j];
 418     }
 419     acf[i] = tmp /(smpl_t)(length-i);
 420     tmp = 0.0;
 421   }
 422 }
 423
 424 void aubio_cleanup(void) {
 425 #if HAVE_FFTW3
 426   fftw_cleanup();
 427 #else
 428 #if HAVE_FFTW3F
 429   fftwf_cleanup();
 430 #endif
 431 #endif
 432 }