src/mathutils.c

   1 /*
   2   Copyright (C) 2003-2014 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 "musicutils.h"
  27
  28 /** Window types */
  29 typedef enum
  30 {
  31   aubio_win_rectangle,
  32   aubio_win_hamming,
  33   aubio_win_hanning,
  34   aubio_win_hanningz,
  35   aubio_win_blackman,
  36   aubio_win_blackman_harris,
  37   aubio_win_gaussian,
  38   aubio_win_welch,
  39   aubio_win_parzen,
  40   aubio_win_default = aubio_win_hanningz,
  41 } aubio_window_type;
  42
  43 fvec_t *
  44 new_aubio_window (char_t * window_type, uint_t length)
  45 {
  46   fvec_t * win = new_fvec (length);
  47   uint_t err;
  48   if (win == NULL) {
  49     return NULL;
  50   }
  51   err = fvec_set_window (win, window_type);
  52   if (err != 0) {
  53     del_fvec(win);
  54     return NULL;
  55   }
  56   return win;
  57 }
  58
  59 uint_t fvec_set_window (fvec_t *win, char_t *window_type) {
  60   smpl_t * w = win->data;
  61   uint_t i, size = win->length;
  62   aubio_window_type wintype;
  63   if (window_type == NULL) {
  64       AUBIO_ERR ("window type can not be null.\n");
  65       return 1;
  66   } else if (strcmp (window_type, "rectangle") == 0)
  67       wintype = aubio_win_rectangle;
  68   else if (strcmp (window_type, "hamming") == 0)
  69       wintype = aubio_win_hamming;
  70   else if (strcmp (window_type, "hanning") == 0)
  71       wintype = aubio_win_hanning;
  72   else if (strcmp (window_type, "hanningz") == 0)
  73       wintype = aubio_win_hanningz;
  74   else if (strcmp (window_type, "blackman") == 0)
  75       wintype = aubio_win_blackman;
  76   else if (strcmp (window_type, "blackman_harris") == 0)
  77       wintype = aubio_win_blackman_harris;
  78   else if (strcmp (window_type, "gaussian") == 0)
  79       wintype = aubio_win_gaussian;
  80   else if (strcmp (window_type, "welch") == 0)
  81       wintype = aubio_win_welch;
  82   else if (strcmp (window_type, "parzen") == 0)
  83       wintype = aubio_win_parzen;
  84   else if (strcmp (window_type, "default") == 0)
  85       wintype = aubio_win_default;
  86   else {
  87       AUBIO_ERR ("unknown window type `%s`.\n", window_type);
  88       return 1;
  89   }
  90   switch(wintype) {
  91     case aubio_win_rectangle:
  92       for (i=0;i<size;i++)
  93         w[i] = 0.5;
  94       break;
  95     case aubio_win_hamming:
  96       for (i=0;i<size;i++)
  97         w[i] = 0.54 - 0.46 * COS(TWO_PI * i / (size));
  98       break;
  99     case aubio_win_hanning:
 100       for (i=0;i<size;i++)
 101         w[i] = 0.5 - (0.5 * COS(TWO_PI * i / (size)));
 102       break;
 103     case aubio_win_hanningz:
 104       for (i=0;i<size;i++)
 105         w[i] = 0.5 * (1.0 - COS(TWO_PI * i / (size)));
 106       break;
 107     case aubio_win_blackman:
 108       for (i=0;i<size;i++)
 109         w[i] = 0.42
 110           - 0.50 * COS(    TWO_PI*i/(size-1.0))
 111           + 0.08 * COS(2.0*TWO_PI*i/(size-1.0));
 112       break;
 113     case aubio_win_blackman_harris:
 114       for (i=0;i<size;i++)
 115         w[i] = 0.35875
 116           - 0.48829 * COS(    TWO_PI*i/(size-1.0))
 117           + 0.14128 * COS(2.0*TWO_PI*i/(size-1.0))
 118           - 0.01168 * COS(3.0*TWO_PI*i/(size-1.0));
 119       break;
 120     case aubio_win_gaussian:
 121       {
 122         lsmp_t a, b, c = 0.5;
 123         uint_t n;
 124         for (n = 0; n < size; n++)
 125         {
 126           a = (n-c*(size-1))/(SQR(c)*(size-1));
 127           b = -c*SQR(a);
 128           w[n] = EXP(b);
 129         }
 130       }
 131       break;
 132     case aubio_win_welch:
 133       for (i=0;i<size;i++)
 134         w[i] = 1.0 - SQR((2.*i-size)/(size+1.0));
 135       break;
 136     case aubio_win_parzen:
 137       for (i=0;i<size;i++)
 138         w[i] = 1.0 - ABS((2.f*i-size)/(size+1.0f));
 139       break;
 140     default:
 141       break;
 142   }
 143   return 0;
 144 }
 145
 146 smpl_t
 147 aubio_unwrap2pi (smpl_t phase)
 148 {
 149   /* mod(phase+pi,-2pi)+pi */
 150   return phase + TWO_PI * (1. + FLOOR (-(phase + PI) / TWO_PI));
 151 }
 152
 153 smpl_t
 154 fvec_mean (fvec_t * s)
 155 {
 156   smpl_t tmp = 0.0;
 157 #ifndef HAVE_ACCELERATE
 158   uint_t j;
 159   for (j = 0; j < s->length; j++) {
 160     tmp += s->data[j];
 161   }
 162   return tmp / (smpl_t) (s->length);
 163 #else
 164   aubio_vDSP_meanv(s->data, 1, &tmp, s->length);
 165   return tmp;
 166 #endif /* HAVE_ACCELERATE */
 167 }
 168
 169 smpl_t
 170 fvec_sum (fvec_t * s)
 171 {
 172   smpl_t tmp = 0.0;
 173 #ifndef HAVE_ACCELERATE
 174   uint_t j;
 175   for (j = 0; j < s->length; j++) {
 176     tmp += s->data[j];
 177   }
 178 #else
 179   aubio_vDSP_sve(s->data, 1, &tmp, s->length);
 180 #endif /* HAVE_ACCELERATE */
 181   return tmp;
 182 }
 183
 184 smpl_t
 185 fvec_max (fvec_t * s)
 186 {
 187 #ifndef HAVE_ACCELERATE
 188   uint_t j;
 189   smpl_t tmp = 0.0;
 190   for (j = 0; j < s->length; j++) {
 191     tmp = (tmp > s->data[j]) ? tmp : s->data[j];
 192   }
 193 #else
 194   smpl_t tmp = 0.;
 195   aubio_vDSP_maxv(s->data, 1, &tmp, s->length);
 196 #endif
 197   return tmp;
 198 }
 199
 200 smpl_t
 201 fvec_min (fvec_t * s)
 202 {
 203 #ifndef HAVE_ACCELERATE
 204   uint_t j;
 205   smpl_t tmp = s->data[0];
 206   for (j = 0; j < s->length; j++) {
 207     tmp = (tmp < s->data[j]) ? tmp : s->data[j];
 208   }
 209 #else
 210   smpl_t tmp = 0.;
 211   aubio_vDSP_minv(s->data, 1, &tmp, s->length);
 212 #endif
 213   return tmp;
 214 }
 215
 216 uint_t
 217 fvec_min_elem (fvec_t * s)
 218 {
 219 #ifndef HAVE_ACCELERATE
 220   uint_t j, pos = 0.;
 221   smpl_t tmp = s->data[0];
 222   for (j = 0; j < s->length; j++) {
 223     pos = (tmp < s->data[j]) ? pos : j;
 224     tmp = (tmp < s->data[j]) ? tmp : s->data[j];
 225   }
 226 #else
 227   smpl_t tmp = 0.;
 228   vDSP_Length pos = 0;
 229   aubio_vDSP_minvi(s->data, 1, &tmp, &pos, s->length);
 230 #endif
 231   return (uint_t)pos;
 232 }
 233
 234 uint_t
 235 fvec_max_elem (fvec_t * s)
 236 {
 237 #ifndef HAVE_ACCELERATE
 238   uint_t j, pos = 0;
 239   smpl_t tmp = 0.0;
 240   for (j = 0; j < s->length; j++) {
 241     pos = (tmp > s->data[j]) ? pos : j;
 242     tmp = (tmp > s->data[j]) ? tmp : s->data[j];
 243   }
 244 #else
 245   smpl_t tmp = 0.;
 246   vDSP_Length pos = 0;
 247   aubio_vDSP_maxvi(s->data, 1, &tmp, &pos, s->length);
 248 #endif
 249   return (uint_t)pos;
 250 }
 251
 252 void
 253 fvec_shift (fvec_t * s)
 254 {
 255   uint_t half = s->length / 2, start = half, j;
 256   // if length is odd, middle element is moved to the end
 257   if (2 * half < s->length) start ++;
 258 #ifndef HAVE_ATLAS
 259   for (j = 0; j < half; j++) {
 260     ELEM_SWAP (s->data[j], s->data[j + start]);
 261   }
 262 #else
 263   aubio_cblas_swap(half, s->data, 1, s->data + start, 1);
 264 #endif
 265   if (start != half) {
 266     for (j = 0; j < half; j++) {
 267       ELEM_SWAP (s->data[j + start - 1], s->data[j + start]);
 268     }
 269   }
 270 }
 271
 272 void
 273 fvec_ishift (fvec_t * s)
 274 {
 275   uint_t half = s->length / 2, start = half, j;
 276   // if length is odd, middle element is moved to the beginning
 277   if (2 * half < s->length) start ++;
 278 #ifndef HAVE_ATLAS
 279   for (j = 0; j < half; j++) {
 280     ELEM_SWAP (s->data[j], s->data[j + start]);
 281   }
 282 #else
 283   aubio_cblas_swap(half, s->data, 1, s->data + start, 1);
 284 #endif
 285   if (start != half) {
 286     for (j = 0; j < half; j++) {
 287       ELEM_SWAP (s->data[half], s->data[j]);
 288     }
 289   }
 290 }
 291
 292 void fvec_push(fvec_t *in, smpl_t new_elem) {
 293   uint_t i;
 294   for (i = 0; i < in->length - 1; i++) {
 295     in->data[i] = in->data[i + 1];
 296   }
 297   in->data[in->length - 1] = new_elem;
 298 }
 299
 300 void fvec_clamp(fvec_t *in, smpl_t absmax) {
 301   uint_t i;
 302   for (i = 0; i < in->length; i++) {
 303     if (in->data[i] > 0 && in->data[i] > ABS(absmax)) {
 304       in->data[i] = absmax;
 305     } else if (in->data[i] < 0 && in->data[i] < -ABS(absmax)) {
 306       in->data[i] = -absmax;
 307     }
 308   }
 309 }
 310
 311 smpl_t
 312 aubio_level_lin (const fvec_t * f)
 313 {
 314   smpl_t energy = 0.;
 315 #ifndef HAVE_ATLAS
 316   uint_t j;
 317   for (j = 0; j < f->length; j++) {
 318     energy += SQR (f->data[j]);
 319   }
 320 #else
 321   energy = aubio_cblas_dot(f->length, f->data, 1, f->data, 1);
 322 #endif
 323   return energy / f->length;
 324 }
 325
 326 smpl_t
 327 fvec_local_hfc (fvec_t * v)
 328 {
 329   smpl_t hfc = 0.;
 330   uint_t j;
 331   for (j = 0; j < v->length; j++) {
 332     hfc += (j + 1) * v->data[j];
 333   }
 334   return hfc;
 335 }
 336
 337 void
 338 fvec_min_removal (fvec_t * v)
 339 {
 340   smpl_t v_min = fvec_min (v);
 341   fvec_add (v,  - v_min );
 342 }
 343
 344 smpl_t
 345 fvec_alpha_norm (fvec_t * o, smpl_t alpha)
 346 {
 347   uint_t j;
 348   smpl_t tmp = 0.;
 349   for (j = 0; j < o->length; j++) {
 350     tmp += POW (ABS (o->data[j]), alpha);
 351   }
 352   return POW (tmp / o->length, 1. / alpha);
 353 }
 354
 355 void
 356 fvec_alpha_normalise (fvec_t * o, smpl_t alpha)
 357 {
 358   uint_t j;
 359   smpl_t norm = fvec_alpha_norm (o, alpha);
 360   for (j = 0; j < o->length; j++) {
 361     o->data[j] /= norm;
 362   }
 363 }
 364
 365 void
 366 fvec_add (fvec_t * o, smpl_t val)
 367 {
 368   uint_t j;
 369   for (j = 0; j < o->length; j++) {
 370     o->data[j] += val;
 371   }
 372 }
 373
 374 void fvec_adapt_thres(fvec_t * vec, fvec_t * tmp,
 375     uint_t post, uint_t pre) {
 376   uint_t length = vec->length, j;
 377   for (j=0;j<length;j++) {
 378     vec->data[j] -= fvec_moving_thres(vec, tmp, post, pre, j);
 379   }
 380 }
 381
 382 smpl_t
 383 fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec,
 384     uint_t post, uint_t pre, uint_t pos)
 385 {
 386   uint_t k;
 387   smpl_t *medar = (smpl_t *) tmpvec->data;
 388   uint_t win_length = post + pre + 1;
 389   uint_t length = vec->length;
 390   /* post part of the buffer does not exist */
 391   if (pos < post + 1) {
 392     for (k = 0; k < post + 1 - pos; k++)
 393       medar[k] = 0.;            /* 0-padding at the beginning */
 394     for (k = post + 1 - pos; k < win_length; k++)
 395       medar[k] = vec->data[k + pos - post];
 396     /* the buffer is fully defined */
 397   } else if (pos + pre < length) {
 398     for (k = 0; k < win_length; k++)
 399       medar[k] = vec->data[k + pos - post];
 400     /* pre part of the buffer does not exist */
 401   } else {
 402     for (k = 0; k < length - pos + post; k++)
 403       medar[k] = vec->data[k + pos - post];
 404     for (k = length - pos + post; k < win_length; k++)
 405       medar[k] = 0.;            /* 0-padding at the end */
 406   }
 407   return fvec_median (tmpvec);
 408 }
 409
 410 smpl_t fvec_median (fvec_t * input) {
 411   uint_t n = input->length;
 412   smpl_t * arr = (smpl_t *) input->data;
 413   uint_t low, high ;
 414   uint_t median;
 415   uint_t middle, ll, hh;
 416
 417   low = 0 ; high = n-1 ; median = (low + high) / 2;
 418   for (;;) {
 419     if (high <= low) /* One element only */
 420       return arr[median] ;
 421
 422     if (high == low + 1) {  /* Two elements only */
 423       if (arr[low] > arr[high])
 424         ELEM_SWAP(arr[low], arr[high]) ;
 425       return arr[median] ;
 426     }
 427
 428     /* Find median of low, middle and high items; swap into position low */
 429     middle = (low + high) / 2;
 430     if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]);
 431     if (arr[low]    > arr[high])    ELEM_SWAP(arr[low],    arr[high]);
 432     if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
 433
 434     /* Swap low item (now in position middle) into position (low+1) */
 435     ELEM_SWAP(arr[middle], arr[low+1]) ;
 436
 437     /* Nibble from each end towards middle, swapping items when stuck */
 438     ll = low + 1;
 439     hh = high;
 440     for (;;) {
 441       do ll++; while (arr[low] > arr[ll]) ;
 442       do hh--; while (arr[hh]  > arr[low]) ;
 443
 444       if (hh < ll)
 445         break;
 446
 447       ELEM_SWAP(arr[ll], arr[hh]) ;
 448     }
 449
 450     /* Swap middle item (in position low) back into correct position */
 451     ELEM_SWAP(arr[low], arr[hh]) ;
 452
 453     /* Re-set active partition */
 454     if (hh <= median)
 455       low = ll;
 456     if (hh >= median)
 457       high = hh - 1;
 458   }
 459 }
 460
 461 smpl_t fvec_quadratic_peak_pos (const fvec_t * x, uint_t pos) {
 462   smpl_t s0, s1, s2; uint_t x0, x2;
 463   smpl_t half = .5, two = 2.;
 464   if (pos == 0 || pos == x->length - 1) return pos;
 465   x0 = (pos < 1) ? pos : pos - 1;
 466   x2 = (pos + 1 < x->length) ? pos + 1 : pos;
 467   if (x0 == pos) return (x->data[pos] <= x->data[x2]) ? pos : x2;
 468   if (x2 == pos) return (x->data[pos] <= x->data[x0]) ? pos : x0;
 469   s0 = x->data[x0];
 470   s1 = x->data[pos];
 471   s2 = x->data[x2];
 472   return pos + half * (s0 - s2 ) / (s0 - two * s1 + s2);
 473 }
 474
 475 smpl_t fvec_quadratic_peak_mag (fvec_t *x, smpl_t pos) {
 476   smpl_t x0, x1, x2;
 477   uint_t index = (uint_t)(pos - .5) + 1;
 478   if (pos >= x->length || pos < 0.) return 0.;
 479   if ((smpl_t)index == pos) return x->data[index];
 480   x0 = x->data[index - 1];
 481   x1 = x->data[index];
 482   x2 = x->data[index + 1];
 483   return x1 - .25 * (x0 - x2) * (pos - index);
 484 }
 485
 486 uint_t fvec_peakpick(const fvec_t * onset, uint_t pos) {
 487   uint_t tmp=0;
 488   tmp = (onset->data[pos] > onset->data[pos-1]
 489       &&  onset->data[pos] > onset->data[pos+1]
 490       &&  onset->data[pos] > 0.);
 491   return tmp;
 492 }
 493
 494 smpl_t
 495 aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf)
 496 {
 497   smpl_t tmp =
 498       s0 + (pf / 2.) * (pf * (s0 - 2. * s1 + s2) - 3. * s0 + 4. * s1 - s2);
 499   return tmp;
 500 }
 501
 502 smpl_t
 503 aubio_freqtomidi (smpl_t freq)
 504 {
 505   smpl_t midi;
 506   if (freq < 2. || freq > 100000.) return 0.; // avoid nans and infs
 507   /* log(freq/A-2)/log(2) */
 508   midi = freq / 6.875;
 509   midi = LOG (midi) / 0.69314718055995;
 510   midi *= 12;
 511   midi -= 3;
 512   return midi;
 513 }
 514
 515 smpl_t
 516 aubio_miditofreq (smpl_t midi)
 517 {
 518   smpl_t freq;
 519   if (midi > 140.) return 0.; // avoid infs
 520   freq = (midi + 3.) / 12.;
 521   freq = EXP (freq * 0.69314718055995);
 522   freq *= 6.875;
 523   return freq;
 524 }
 525
 526 smpl_t
 527 aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
 528 {
 529   smpl_t freq = samplerate / fftsize;
 530   return freq * MAX(bin, 0);
 531 }
 532
 533 smpl_t
 534 aubio_bintomidi (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
 535 {
 536   smpl_t midi = aubio_bintofreq (bin, samplerate, fftsize);
 537   return aubio_freqtomidi (midi);
 538 }
 539
 540 smpl_t
 541 aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize)
 542 {
 543   smpl_t bin = fftsize / samplerate;
 544   return MAX(freq, 0) * bin;
 545 }
 546
 547 smpl_t
 548 aubio_miditobin (smpl_t midi, smpl_t samplerate, smpl_t fftsize)
 549 {
 550   smpl_t freq = aubio_miditofreq (midi);
 551   return aubio_freqtobin (freq, samplerate, fftsize);
 552 }
 553
 554 uint_t
 555 aubio_is_power_of_two (uint_t a)
 556 {
 557   if ((a & (a - 1)) == 0) {
 558     return 1;
 559   } else {
 560     return 0;
 561   }
 562 }
 563
 564 uint_t
 565 aubio_next_power_of_two (uint_t a)
 566 {
 567   uint_t i = 1;
 568   while (i < a) i <<= 1;
 569   return i;
 570 }
 571
 572 smpl_t
 573 aubio_db_spl (const fvec_t * o)
 574 {
 575   return 10. * LOG10 (aubio_level_lin (o));
 576 }
 577
 578 uint_t
 579 aubio_silence_detection (const fvec_t * o, smpl_t threshold)
 580 {
 581   return (aubio_db_spl (o) < threshold);
 582 }
 583
 584 smpl_t
 585 aubio_level_detection (const fvec_t * o, smpl_t threshold)
 586 {
 587   smpl_t db_spl = aubio_db_spl (o);
 588   if (db_spl < threshold) {
 589     return 1.;
 590   } else {
 591     return db_spl;
 592   }
 593 }
 594
 595 smpl_t
 596 aubio_zero_crossing_rate (fvec_t * input)
 597 {
 598   uint_t j;
 599   uint_t zcr = 0;
 600   for (j = 1; j < input->length; j++) {
 601     // previous was strictly negative
 602     if (input->data[j - 1] < 0.) {
 603       // current is positive or null
 604       if (input->data[j] >= 0.) {
 605         zcr += 1;
 606       }
 607       // previous was positive or null
 608     } else {
 609       // current is strictly negative
 610       if (input->data[j] < 0.) {
 611         zcr += 1;
 612       }
 613     }
 614   }
 615   return zcr / (smpl_t) input->length;
 616 }
 617
 618 void
 619 aubio_autocorr (const fvec_t * input, fvec_t * output)
 620 {
 621   uint_t i, j, length = input->length;
 622   smpl_t *data, *acf;
 623   smpl_t tmp = 0;
 624   data = input->data;
 625   acf = output->data;
 626   for (i = 0; i < length; i++) {
 627     tmp = 0.;
 628     for (j = i; j < length; j++) {
 629       tmp += data[j - i] * data[j];
 630     }
 631     acf[i] = tmp / (smpl_t) (length - i);
 632   }
 633 }
 634
 635 void
 636 aubio_cleanup (void)
 637 {
 638 #ifdef HAVE_FFTW3F
 639   fftwf_cleanup ();
 640 #else
 641 #ifdef HAVE_FFTW3
 642   fftw_cleanup ();
 643 #endif
 644 #endif
 645 }