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