src/musicutils.h: aubio_level_detection input is a constant
[aubio.git] / src / mathutils.c
index 6ff973b..01b3519 100644 (file)
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
+  Copyright (C) 2003-2014 Paul Brossier <piem@aubio.org>
 
   This file is part of aubio.
 
@@ -26,7 +26,6 @@
 #include "musicutils.h"
 #include "config.h"
 
-
 /** Window types */
 typedef enum
 {
@@ -43,14 +42,29 @@ typedef enum
 } aubio_window_type;
 
 fvec_t *
-new_aubio_window (char_t * window_type, uint_t size)
+new_aubio_window (char_t * window_type, uint_t length)
 {
-  // create fvec of size x 1 channel
-  fvec_t * win = new_fvec( size, 1);
-  smpl_t * w = win->data[0];
-  uint_t i;
+  fvec_t * win = new_fvec (length);
+  uint_t err;
+  if (win == NULL) {
+    return NULL;
+  }
+  err = fvec_set_window (win, window_type);
+  if (err != 0) {
+    del_fvec(win);
+    return NULL;
+  }
+  return win;
+}
+
+uint_t fvec_set_window (fvec_t *win, char_t *window_type) {
+  smpl_t * w = win->data;
+  uint_t i, size = win->length;
   aubio_window_type wintype;
-  if (strcmp (window_type, "rectangle") == 0)
+  if (window_type == NULL) {
+      AUBIO_ERR ("window type can not be null.\n");
+      return 1;
+  } else if (strcmp (window_type, "rectangle") == 0)
       wintype = aubio_win_rectangle;
   else if (strcmp (window_type, "hamming") == 0)
       wintype = aubio_win_hamming;
@@ -71,9 +85,8 @@ new_aubio_window (char_t * window_type, uint_t size)
   else if (strcmp (window_type, "default") == 0)
       wintype = aubio_win_default;
   else {
-      AUBIO_ERR ("unknown window type %s, using default.\n", window_type);
-      wintype = aubio_win_default;
-      return NULL;
+      AUBIO_ERR ("unknown window type `%s`.\n", window_type);
+      return 1;
   }
   switch(wintype) {
     case aubio_win_rectangle:
@@ -106,21 +119,29 @@ new_aubio_window (char_t * window_type, uint_t size)
           - 0.01168 * COS(3.0*TWO_PI*i/(size-1.0));
       break;
     case aubio_win_gaussian:
-      for (i=0;i<size;i++)
-        w[i] = EXP(- 1.0 / SQR(size) * SQR(2.0*i-size));
+      {
+        lsmp_t a, b, c = 0.5;
+        uint_t n;
+        for (n = 0; n < size; n++)
+        {
+          a = (n-c*(size-1))/(SQR(c)*(size-1));
+          b = -c*SQR(a);
+          w[n] = EXP(b);
+        }
+      }
       break;
     case aubio_win_welch:
       for (i=0;i<size;i++)
-        w[i] = 1.0 - SQR((2*i-size)/(size+1.0));
+        w[i] = 1.0 - SQR((2.*i-size)/(size+1.0));
       break;
     case aubio_win_parzen:
       for (i=0;i<size;i++)
-        w[i] = 1.0 - ABS((2*i-size)/(size+1.0));
+        w[i] = 1.0 - ABS((2.f*i-size)/(size+1.0f));
       break;
     default:
       break;
   }
-  return win;
+  return 0;
 }
 
 smpl_t
@@ -133,124 +154,164 @@ aubio_unwrap2pi (smpl_t phase)
 smpl_t
 fvec_mean (fvec_t * s)
 {
-  uint_t i, j;
   smpl_t tmp = 0.0;
-  for (i = 0; i < s->channels; i++)
-    for (j = 0; j < s->length; j++)
-      tmp += s->data[i][j];
-  return tmp / (smpl_t) (s->length);
-}
-
-smpl_t
-fvec_mean_channel (fvec_t * s, uint_t i)
-{
+#ifndef HAVE_ACCELERATE
   uint_t j;
-  smpl_t tmp = 0.0;
-  for (j = 0; j < s->length; j++)
-      tmp += s->data[i][j];
+  for (j = 0; j < s->length; j++) {
+    tmp += s->data[j];
+  }
   return tmp / (smpl_t) (s->length);
+#else
+  aubio_vDSP_meanv(s->data, 1, &tmp, s->length);
+  return tmp;
+#endif /* HAVE_ACCELERATE */
 }
 
 smpl_t
 fvec_sum (fvec_t * s)
 {
-  uint_t i, j;
   smpl_t tmp = 0.0;
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length; j++) {
-      tmp += s->data[i][j];
-    }
+#ifndef HAVE_ACCELERATE
+  uint_t j;
+  for (j = 0; j < s->length; j++) {
+    tmp += s->data[j];
   }
+#else
+  aubio_vDSP_sve(s->data, 1, &tmp, s->length);
+#endif /* HAVE_ACCELERATE */
   return tmp;
 }
 
 smpl_t
 fvec_max (fvec_t * s)
 {
-  uint_t i, j;
+#ifndef HAVE_ACCELERATE
+  uint_t j;
   smpl_t tmp = 0.0;
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length; j++) {
-      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
-    }
+  for (j = 0; j < s->length; j++) {
+    tmp = (tmp > s->data[j]) ? tmp : s->data[j];
   }
+#else
+  smpl_t tmp = 0.;
+  aubio_vDSP_maxv(s->data, 1, &tmp, s->length);
+#endif
   return tmp;
 }
 
 smpl_t
 fvec_min (fvec_t * s)
 {
-  uint_t i, j;
-  smpl_t tmp = s->data[0][0];
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length; j++) {
-      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
-    }
+#ifndef HAVE_ACCELERATE
+  uint_t j;
+  smpl_t tmp = s->data[0];
+  for (j = 0; j < s->length; j++) {
+    tmp = (tmp < s->data[j]) ? tmp : s->data[j];
   }
+#else
+  smpl_t tmp = 0.;
+  aubio_vDSP_minv(s->data, 1, &tmp, s->length);
+#endif
   return tmp;
 }
 
 uint_t
 fvec_min_elem (fvec_t * s)
 {
-  uint_t i, j, pos = 0.;
-  smpl_t tmp = s->data[0][0];
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length; j++) {
-      pos = (tmp < s->data[i][j]) ? pos : j;
-      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
-    }
+#ifndef HAVE_ACCELERATE
+  uint_t j, pos = 0.;
+  smpl_t tmp = s->data[0];
+  for (j = 0; j < s->length; j++) {
+    pos = (tmp < s->data[j]) ? pos : j;
+    tmp = (tmp < s->data[j]) ? tmp : s->data[j];
   }
+#else
+  smpl_t tmp = 0.;
+  uint_t pos = 0.;
+  aubio_vDSP_minvi(s->data, 1, &tmp, (vDSP_Length *)&pos, s->length);
+#endif
   return pos;
 }
 
 uint_t
 fvec_max_elem (fvec_t * s)
 {
-  uint_t i, j, pos = 0;
+#ifndef HAVE_ACCELERATE
+  uint_t j, pos = 0;
   smpl_t tmp = 0.0;
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length; j++) {
-      pos = (tmp > s->data[i][j]) ? pos : j;
-      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
-    }
+  for (j = 0; j < s->length; j++) {
+    pos = (tmp > s->data[j]) ? pos : j;
+    tmp = (tmp > s->data[j]) ? tmp : s->data[j];
   }
+#else
+  smpl_t tmp = 0.;
+  uint_t pos = 0.;
+  aubio_vDSP_maxvi(s->data, 1, &tmp, (vDSP_Length *)&pos, s->length);
+#endif
   return pos;
 }
 
 void
 fvec_shift (fvec_t * s)
 {
-  uint_t i, j;
-  for (i = 0; i < s->channels; i++) {
-    for (j = 0; j < s->length / 2; j++) {
-      ELEM_SWAP (s->data[i][j], s->data[i][j + s->length / 2]);
+  uint_t half = s->length / 2, start = half, j;
+  // if length is odd, middle element is moved to the end
+  if (2 * half < s->length) start ++;
+#ifndef HAVE_ATLAS
+  for (j = 0; j < half; j++) {
+    ELEM_SWAP (s->data[j], s->data[j + start]);
+  }
+#else
+  aubio_cblas_swap(half, s->data, 1, s->data + start, 1);
+#endif
+  if (start != half) {
+    for (j = 0; j < half; j++) {
+      ELEM_SWAP (s->data[j + start - 1], s->data[j + start]);
+    }
+  }
+}
+
+void
+fvec_ishift (fvec_t * s)
+{
+  uint_t half = s->length / 2, start = half, j;
+  // if length is odd, middle element is moved to the beginning
+  if (2 * half < s->length) start ++;
+#ifndef HAVE_ATLAS
+  for (j = 0; j < half; j++) {
+    ELEM_SWAP (s->data[j], s->data[j + start]);
+  }
+#else
+  aubio_cblas_swap(half, s->data, 1, s->data + start, 1);
+#endif
+  if (start != half) {
+    for (j = 0; j < half; j++) {
+      ELEM_SWAP (s->data[half], s->data[j]);
     }
   }
 }
 
 smpl_t
-fvec_local_energy (fvec_t * f)
+aubio_level_lin (const fvec_t * f)
 {
   smpl_t energy = 0.;
-  uint_t i, j;
-  for (i = 0; i < f->channels; i++) {
-    for (j = 0; j < f->length; j++) {
-      energy += SQR (f->data[i][j]);
-    }
+#ifndef HAVE_ATLAS
+  uint_t j;
+  for (j = 0; j < f->length; j++) {
+    energy += SQR (f->data[j]);
   }
-  return energy;
+#else
+  energy = aubio_cblas_dot(f->length, f->data, 1, f->data, 1);
+#endif
+  return energy / f->length;
 }
 
 smpl_t
 fvec_local_hfc (fvec_t * v)
 {
   smpl_t hfc = 0.;
-  uint_t i, j;
-  for (i = 0; i < v->channels; i++) {
-    for (j = 0; j < v->length; j++) {
-      hfc += (i + 1) * v->data[i][j];
-    }
+  uint_t j;
+  for (j = 0; j < v->length; j++) {
+    hfc += (j + 1) * v->data[j];
   }
   return hfc;
 }
@@ -265,12 +326,10 @@ fvec_min_removal (fvec_t * v)
 smpl_t
 fvec_alpha_norm (fvec_t * o, smpl_t alpha)
 {
-  uint_t i, j;
+  uint_t j;
   smpl_t tmp = 0.;
-  for (i = 0; i < o->channels; i++) {
-    for (j = 0; j < o->length; j++) {
-      tmp += POW (ABS (o->data[i][j]), alpha);
-    }
+  for (j = 0; j < o->length; j++) {
+    tmp += POW (ABS (o->data[j]), alpha);
   }
   return POW (tmp / o->length, 1. / alpha);
 }
@@ -278,40 +337,36 @@ fvec_alpha_norm (fvec_t * o, smpl_t alpha)
 void
 fvec_alpha_normalise (fvec_t * o, smpl_t alpha)
 {
-  uint_t i, j;
+  uint_t j;
   smpl_t norm = fvec_alpha_norm (o, alpha);
-  for (i = 0; i < o->channels; i++) {
-    for (j = 0; j < o->length; j++) {
-      o->data[i][j] /= norm;
-    }
+  for (j = 0; j < o->length; j++) {
+    o->data[j] /= norm;
   }
 }
 
 void
 fvec_add (fvec_t * o, smpl_t val)
 {
-  uint_t i, j;
-  for (i = 0; i < o->channels; i++) {
-    for (j = 0; j < o->length; j++) {
-      o->data[i][j] += val;
-    }
+  uint_t j;
+  for (j = 0; j < o->length; j++) {
+    o->data[j] += val;
   }
 }
 
 void fvec_adapt_thres(fvec_t * vec, fvec_t * tmp,
-    uint_t post, uint_t pre, uint_t channel) {
-  uint_t length = vec->length, i=channel, j;
+    uint_t post, uint_t pre) {
+  uint_t length = vec->length, j;
   for (j=0;j<length;j++) {
-    vec->data[i][j] -= fvec_moving_thres(vec, tmp, post, pre, j, i);
+    vec->data[j] -= fvec_moving_thres(vec, tmp, post, pre, j);
   }
 }
 
 smpl_t
 fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec,
-    uint_t post, uint_t pre, uint_t pos, uint_t channel)
+    uint_t post, uint_t pre, uint_t pos)
 {
-  uint_t i = channel, k;
-  smpl_t *medar = (smpl_t *) tmpvec->data[i];
+  uint_t k;
+  smpl_t *medar = (smpl_t *) tmpvec->data;
   uint_t win_length = post + pre + 1;
   uint_t length = vec->length;
   /* post part of the buffer does not exist */
@@ -319,24 +374,24 @@ fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec,
     for (k = 0; k < post + 1 - pos; k++)
       medar[k] = 0.;            /* 0-padding at the beginning */
     for (k = post + 1 - pos; k < win_length; k++)
-      medar[k] = vec->data[0][k + pos - post];
+      medar[k] = vec->data[k + pos - post];
     /* the buffer is fully defined */
   } else if (pos + pre < length) {
     for (k = 0; k < win_length; k++)
-      medar[k] = vec->data[0][k + pos - post];
+      medar[k] = vec->data[k + pos - post];
     /* pre part of the buffer does not exist */
   } else {
     for (k = 0; k < length - pos + post; k++)
-      medar[k] = vec->data[0][k + pos - post];
+      medar[k] = vec->data[k + pos - post];
     for (k = length - pos + post; k < win_length; k++)
       medar[k] = 0.;            /* 0-padding at the end */
   }
-  return fvec_median_channel (tmpvec, i);
+  return fvec_median (tmpvec);
 }
 
-smpl_t fvec_median_channel (fvec_t * input, uint_t channel) {
+smpl_t fvec_median (fvec_t * input) {
   uint_t n = input->length;
-  smpl_t * arr = (smpl_t *) input->data[channel];
+  smpl_t * arr = (smpl_t *) input->data;
   uint_t low, high ;
   uint_t median;
   uint_t middle, ll, hh;
@@ -385,24 +440,36 @@ smpl_t fvec_median_channel (fvec_t * input, uint_t channel) {
   }
 }
 
-smpl_t fvec_quadint (fvec_t * x, uint_t pos, uint_t i) {
-  smpl_t s0, s1, s2;
-  uint_t x0 = (pos < 1) ? pos : pos - 1;
-  uint_t x2 = (pos + 1 < x->length) ? pos + 1 : pos;
-  if (x0 == pos) return (x->data[i][pos] <= x->data[i][x2]) ? pos : x2;
-  if (x2 == pos) return (x->data[i][pos] <= x->data[i][x0]) ? pos : x0;
-  s0 = x->data[i][x0];
-  s1 = x->data[i][pos];
-  s2 = x->data[i][x2];
-  return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
-}
-
-uint_t fvec_peakpick(fvec_t * onset, uint_t pos) {
-  uint_t i=0, tmp=0;
-  /*for (i=0;i<onset->channels;i++)*/
-  tmp = (onset->data[i][pos] > onset->data[i][pos-1]
-      &&  onset->data[i][pos] > onset->data[i][pos+1]
-      &&  onset->data[i][pos] > 0.);
+smpl_t fvec_quadratic_peak_pos (const fvec_t * x, uint_t pos) {
+  smpl_t s0, s1, s2; uint_t x0, x2;
+  smpl_t half = .5, two = 2.;
+  if (pos == 0 || pos == x->length - 1) return pos;
+  x0 = (pos < 1) ? pos : pos - 1;
+  x2 = (pos + 1 < x->length) ? pos + 1 : pos;
+  if (x0 == pos) return (x->data[pos] <= x->data[x2]) ? pos : x2;
+  if (x2 == pos) return (x->data[pos] <= x->data[x0]) ? pos : x0;
+  s0 = x->data[x0];
+  s1 = x->data[pos];
+  s2 = x->data[x2];
+  return pos + half * (s0 - s2 ) / (s0 - two * s1 + s2);
+}
+
+smpl_t fvec_quadratic_peak_mag (fvec_t *x, smpl_t pos) {
+  smpl_t x0, x1, x2;
+  uint_t index = (uint_t)(pos - .5) + 1;
+  if (pos >= x->length || pos < 0.) return 0.;
+  if ((smpl_t)index == pos) return x->data[index];
+  x0 = x->data[index - 1];
+  x1 = x->data[index];
+  x2 = x->data[index + 1];
+  return x1 - .25 * (x0 - x2) * (pos - index);
+}
+
+uint_t fvec_peakpick(const fvec_t * onset, uint_t pos) {
+  uint_t tmp=0;
+  tmp = (onset->data[pos] > onset->data[pos-1]
+      &&  onset->data[pos] > onset->data[pos+1]
+      &&  onset->data[pos] > 0.);
   return tmp;
 }
 
@@ -417,8 +484,10 @@ aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf)
 smpl_t
 aubio_freqtomidi (smpl_t freq)
 {
+  smpl_t midi;
+  if (freq < 2. || freq > 100000.) return 0.; // avoid nans and infs
   /* log(freq/A-2)/log(2) */
-  smpl_t midi = freq / 6.875;
+  midi = freq / 6.875;
   midi = LOG (midi) / 0.69314718055995;
   midi *= 12;
   midi -= 3;
@@ -428,7 +497,9 @@ aubio_freqtomidi (smpl_t freq)
 smpl_t
 aubio_miditofreq (smpl_t midi)
 {
-  smpl_t freq = (midi + 3.) / 12.;
+  smpl_t freq;
+  if (midi > 140.) return 0.; // avoid infs
+  freq = (midi + 3.) / 12.;
   freq = EXP (freq * 0.69314718055995);
   freq *= 6.875;
   return freq;
@@ -438,7 +509,7 @@ smpl_t
 aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
 {
   smpl_t freq = samplerate / fftsize;
-  return freq * bin;
+  return freq * MAX(bin, 0);
 }
 
 smpl_t
@@ -452,7 +523,7 @@ smpl_t
 aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize)
 {
   smpl_t bin = fftsize / samplerate;
-  return freq * bin;
+  return MAX(freq, 0) * bin;
 }
 
 smpl_t
@@ -475,30 +546,25 @@ aubio_is_power_of_two (uint_t a)
 uint_t
 aubio_next_power_of_two (uint_t a)
 {
-  uint_t i;
-  a--;
-  for (i = 0; i < sizeof (uint_t) * CHAR_BIT; i++) {
-    a = a | a >> 1;
-  }
-  return a + 1;
+  uint_t i = 1;
+  while (i < a) i <<= 1;
+  return i;
 }
 
 smpl_t
-aubio_db_spl (fvec_t * o)
+aubio_db_spl (const fvec_t * o)
 {
-  smpl_t val = SQRT (fvec_local_energy (o));
-  val /= (smpl_t) o->length;
-  return LIN2DB (val);
+  return 10. * LOG10 (aubio_level_lin (o));
 }
 
 uint_t
-aubio_silence_detection (fvec_t * o, smpl_t threshold)
+aubio_silence_detection (const fvec_t * o, smpl_t threshold)
 {
   return (aubio_db_spl (o) < threshold);
 }
 
 smpl_t
-aubio_level_detection (fvec_t * o, smpl_t threshold)
+aubio_level_detection (const fvec_t * o, smpl_t threshold)
 {
   smpl_t db_spl = aubio_db_spl (o);
   if (db_spl < threshold) {
@@ -511,19 +577,19 @@ aubio_level_detection (fvec_t * o, smpl_t threshold)
 smpl_t
 aubio_zero_crossing_rate (fvec_t * input)
 {
-  uint_t i = 0, j;
+  uint_t j;
   uint_t zcr = 0;
   for (j = 1; j < input->length; j++) {
     // previous was strictly negative
-    if (input->data[i][j - 1] < 0.) {
+    if (input->data[j - 1] < 0.) {
       // current is positive or null
-      if (input->data[i][j] >= 0.) {
+      if (input->data[j] >= 0.) {
         zcr += 1;
       }
       // previous was positive or null
     } else {
       // current is strictly negative
-      if (input->data[i][j] < 0.) {
+      if (input->data[j] < 0.) {
         zcr += 1;
       }
     }
@@ -532,32 +598,30 @@ aubio_zero_crossing_rate (fvec_t * input)
 }
 
 void
-aubio_autocorr (fvec_t * input, fvec_t * output)
+aubio_autocorr (const fvec_t * input, fvec_t * output)
 {
-  uint_t i, j, k, length = input->length;
+  uint_t i, j, length = input->length;
   smpl_t *data, *acf;
   smpl_t tmp = 0;
-  for (k = 0; k < input->channels; k++) {
-    data = input->data[k];
-    acf = output->data[k];
-    for (i = 0; i < length; i++) {
-      tmp = 0.;
-      for (j = i; j < length; j++) {
-        tmp += data[j - i] * data[j];
-      }
-      acf[i] = tmp / (smpl_t) (length - i);
+  data = input->data;
+  acf = output->data;
+  for (i = 0; i < length; i++) {
+    tmp = 0.;
+    for (j = i; j < length; j++) {
+      tmp += data[j - i] * data[j];
     }
+    acf[i] = tmp / (smpl_t) (length - i);
   }
 }
 
 void
 aubio_cleanup (void)
 {
-#if HAVE_FFTW3
-  fftw_cleanup ();
-#else
-#if HAVE_FFTW3F
+#ifdef HAVE_FFTW3F
   fftwf_cleanup ();
+#else
+#ifdef HAVE_FFTW3
+  fftw_cleanup ();
 #endif
 #endif
 }