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 fvec_mean(fvec_t *s) {
  87   uint_t i,j;
  88   smpl_t tmp = 0.0;
  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 fvec_sum(fvec_t *s) {
  96   uint_t i,j;
  97   smpl_t tmp = 0.0;
  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 fvec_max(fvec_t *s) {
 105   uint_t i,j;
 106   smpl_t tmp = 0.0;
 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 fvec_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 fvec_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 fvec_max_elem(fvec_t *s) {
 134   uint_t i,j=0, pos=0.;
 135   smpl_t tmp = 0.0;
 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 fvec_shift(fvec_t *s) {
 145   uint_t i,j;
 146   //smpl_t tmp = 0.0;
 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 fvec_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 fvec_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 fvec_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
 184 fvec_min_removal (fvec_t * o)
 185 {
 186   uint_t i, j;
 187   smpl_t mini = fvec_min (mag);
 188   for (i = 0; i < o->channels; i++) {
 189     for (j = 0; j < o->length; j++) {
 190       mag->data[i][j] -= mini;
 191     }
 192   }
 193 }
 194
 195 void fvec_alpha_normalise(fvec_t * mag, uint_t alpha) {
 196   smpl_t alphan = 1.;
 197   uint_t length = mag->length, i=0, j;
 198   alphan = fvec_alpha_norm(mag,alpha);
 199   for (j=0;j<length;j++){
 200     mag->data[i][j] /= alphan;
 201   }
 202 }
 203
 204 void fvec_add(fvec_t * mag, smpl_t threshold) {
 205   uint_t length = mag->length, i=0, j;
 206   for (j=0;j<length;j++) {
 207     mag->data[i][j] += threshold;
 208   }
 209 }
 210
 211 void fvec_adapt_thres(fvec_t * vec, fvec_t * tmp,
 212     uint_t post, uint_t pre) {
 213   uint_t length = vec->length, i=0, j;
 214   for (j=0;j<length;j++) {
 215     vec->data[i][j] -= fvec_moving_thres(vec, tmp, post, pre, j);
 216   }
 217 }
 218
 219 smpl_t fvec_moving_thres(fvec_t * vec, fvec_t * tmpvec,
 220     uint_t post, uint_t pre, uint_t pos) {
 221   smpl_t * medar = (smpl_t *)tmpvec->data[0];
 222   uint_t k;
 223   uint_t win_length =  post+pre+1;
 224   uint_t length =  vec->length;
 225   /* post part of the buffer does not exist */
 226   if (pos<post+1) {
 227     for (k=0;k<post+1-pos;k++)
 228       medar[k] = 0.; /* 0-padding at the beginning */
 229     for (k=post+1-pos;k<win_length;k++)
 230       medar[k] = vec->data[0][k+pos-post];
 231   /* the buffer is fully defined */
 232   } else if (pos+pre<length) {
 233     for (k=0;k<win_length;k++)
 234       medar[k] = vec->data[0][k+pos-post];
 235   /* pre part of the buffer does not exist */
 236   } else {
 237     for (k=0;k<length-pos+post;k++)
 238       medar[k] = vec->data[0][k+pos-post];
 239     for (k=length-pos+post;k<win_length;k++)
 240       medar[k] = 0.; /* 0-padding at the end */
 241   }
 242   return fvec_median(tmpvec);
 243 }
 244
 245 smpl_t fvec_median(fvec_t * input) {
 246   uint_t n = input->length;
 247   smpl_t * arr = (smpl_t *) input->data[0];
 248   uint_t low, high ;
 249   uint_t median;
 250   uint_t middle, ll, hh;
 251
 252   low = 0 ; high = n-1 ; median = (low + high) / 2;
 253   for (;;) {
 254     if (high <= low) /* One element only */
 255       return arr[median] ;
 256
 257     if (high == low + 1) {  /* Two elements only */
 258       if (arr[low] > arr[high])
 259         ELEM_SWAP(arr[low], arr[high]) ;
 260       return arr[median] ;
 261     }
 262
 263     /* Find median of low, middle and high items; swap into position low */
 264     middle = (low + high) / 2;
 265     if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]);
 266     if (arr[low]    > arr[high])    ELEM_SWAP(arr[low],    arr[high]);
 267     if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
 268
 269     /* Swap low item (now in position middle) into position (low+1) */
 270     ELEM_SWAP(arr[middle], arr[low+1]) ;
 271
 272     /* Nibble from each end towards middle, swapping items when stuck */
 273     ll = low + 1;
 274     hh = high;
 275     for (;;) {
 276       do ll++; while (arr[low] > arr[ll]) ;
 277       do hh--; while (arr[hh]  > arr[low]) ;
 278
 279       if (hh < ll)
 280         break;
 281
 282       ELEM_SWAP(arr[ll], arr[hh]) ;
 283     }
 284
 285     /* Swap middle item (in position low) back into correct position */
 286     ELEM_SWAP(arr[low], arr[hh]) ;
 287
 288     /* Re-set active partition */
 289     if (hh <= median)
 290       low = ll;
 291     if (hh >= median)
 292       high = hh - 1;
 293   }
 294 }
 295
 296 smpl_t fvec_quadint(fvec_t * x,uint_t pos, uint_t span) {
 297   smpl_t s0, s1, s2;
 298   uint_t x0 = (pos < span) ? pos : pos - span;
 299   uint_t x2 = (pos + span < x->length) ? pos + span : pos;
 300   if (x0 == pos) return (x->data[0][pos] <= x->data[0][x2]) ? pos : x2;
 301   if (x2 == pos) return (x->data[0][pos] <= x->data[0][x0]) ? pos : x0;
 302   s0 = x->data[0][x0];
 303   s1 = x->data[0][pos]     ;
 304   s2 = x->data[0][x2];
 305   return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
 306 }
 307
 308 smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
 309   smpl_t tmp = s0 + (pf/2.) * (pf * ( s0 - 2.*s1 + s2 ) - 3.*s0 + 4.*s1 - s2);
 310   return tmp;
 311 }
 312
 313 uint_t fvec_peakpick(fvec_t * onset, uint_t pos) {
 314   uint_t i=0, tmp=0;
 315   /*for (i=0;i<onset->channels;i++)*/
 316   tmp = (onset->data[i][pos] > onset->data[i][pos-1]
 317       &&  onset->data[i][pos] > onset->data[i][pos+1]
 318       &&  onset->data[i][pos] > 0.);
 319   return tmp;
 320 }
 321
 322 smpl_t aubio_freqtomidi(smpl_t freq) {
 323   /* log(freq/A-2)/log(2) */
 324   smpl_t midi = freq/6.875;
 325   midi = LOG(midi)/0.69314718055995;
 326   midi *= 12;
 327   midi -= 3;
 328   return midi;
 329 }
 330
 331 smpl_t aubio_miditofreq(smpl_t midi) {
 332   smpl_t freq = (midi+3.)/12.;
 333   freq = EXP(freq*0.69314718055995);
 334   freq *= 6.875;
 335   return freq;
 336 }
 337
 338 smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
 339   smpl_t freq = samplerate/fftsize;
 340   return freq*bin;
 341 }
 342
 343 smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
 344   smpl_t midi = aubio_bintofreq(bin,samplerate,fftsize);
 345   return aubio_freqtomidi(midi);
 346 }
 347
 348 smpl_t aubio_freqtobin(smpl_t freq, smpl_t samplerate, smpl_t fftsize) {
 349   smpl_t bin = fftsize/samplerate;
 350   return freq*bin;
 351 }
 352
 353 smpl_t aubio_miditobin(smpl_t midi, smpl_t samplerate, smpl_t fftsize) {
 354   smpl_t freq = aubio_miditofreq(midi);
 355   return aubio_freqtobin(freq,samplerate,fftsize);
 356 }
 357
 358 /** returns 1 if wassilence is 0 and RMS(ibuf)<threshold
 359  * \bug mono
 360  */
 361 uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold) {
 362   smpl_t loudness = 0;
 363   uint_t i=0,j;
 364   for (j=0;j<ibuf->length;j++) {
 365     loudness += SQR(ibuf->data[i][j]);
 366   }
 367   loudness = SQRT(loudness);
 368   loudness /= (smpl_t)ibuf->length;
 369   loudness = LIN2DB(loudness);
 370
 371   return (loudness < threshold);
 372 }
 373
 374 /** returns level log(RMS(ibuf)) if < threshold, 1 otherwise
 375  * \bug mono
 376  */
 377 smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold) {
 378   smpl_t loudness = 0;
 379   uint_t i=0,j;
 380   for (j=0;j<ibuf->length;j++) {
 381     loudness += SQR(ibuf->data[i][j]);
 382   }
 383   loudness = SQRT(loudness);
 384   loudness /= (smpl_t)ibuf->length;
 385   loudness = LIN2DB(loudness);
 386
 387   if (loudness < threshold)
 388     return 1.;
 389   else
 390     return loudness;
 391 }
 392
 393 smpl_t aubio_zero_crossing_rate(fvec_t * input) {
 394   uint_t i=0,j;
 395   uint_t zcr = 0;
 396   for ( j = 1; j < input->length; j++ ) {
 397     // previous was strictly negative
 398     if( input->data[i][j-1] < 0. ) {
 399       // current is positive or null
 400       if ( input->data[i][j] >= 0. ) {
 401         zcr += 1;
 402       }
 403     // previous was positive or null
 404     } else {
 405       // current is strictly negative
 406       if ( input->data[i][j] < 0. ) {
 407         zcr += 1;
 408       }
 409     }
 410   }
 411   return zcr/(smpl_t)input->length;
 412 }
 413
 414 void aubio_autocorr(fvec_t * input, fvec_t * output) {
 415   uint_t i = 0, j = 0, length = input->length;
 416   smpl_t * data = input->data[0];
 417   smpl_t * acf = output->data[0];
 418   smpl_t tmp =0.;
 419   for(i=0;i<length;i++){
 420     for(j=i;j<length;j++){
 421       tmp += data[j-i]*data[j];
 422     }
 423     acf[i] = tmp /(smpl_t)(length-i);
 424     tmp = 0.0;
 425   }
 426 }
 427
 428 void aubio_cleanup(void) {
 429 #if HAVE_FFTW3
 430   fftw_cleanup();
 431 #else
 432 #if HAVE_FFTW3F
 433   fftwf_cleanup();
 434 #endif
 435 #endif
 436 }