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