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