shithub: aubio

Download patch

ref: c5c0c984a1063e1505c5d0212492e31382d50b90
parent: 68a3fc924772da2b32e508133c171766dc68a424
author: Paul Brossier <piem@piem.org>
date: Fri Sep 11 20:08:54 EDT 2009 

src/mathutils.c: merge and fix vec_quadint_min and _max into simpler and more exact vec_quadint, use in src/pitch algorithms


--- a/src/mathutils.c
+++ b/src/mathutils.c
@@ -285,49 +285,16 @@
   }
 }
 
-smpl_t vec_quadint(fvec_t * x,uint_t pos) {
-  uint_t span = 2;
-  smpl_t step = 1./200.;
-  /* hack : init resold to - something (in case x[pos+-span]<0)) */
-  smpl_t res, frac, s0, s1, s2, exactpos = (smpl_t)pos, resold = -1000.;
-  if ((pos > span) && (pos < x->length-span)) {
-    s0 = x->data[0][pos-span];
-    s1 = x->data[0][pos]     ;
-    s2 = x->data[0][pos+span];
-    /* increase frac */
-    for (frac = 0.; frac < 2.; frac = frac + step) {
-      res = aubio_quadfrac(s0, s1, s2, frac);
-      if (res > resold)
-        resold = res;
-      else {
-        exactpos += (frac-step)*2. - 1.;
-        break;
-      }
-    }
-  }
-  return exactpos;
-}
-
-smpl_t vec_quadint_min(fvec_t * x,uint_t pos, uint_t span) {
-  smpl_t step = 1./200.;
-  /* init resold to - something (in case x[pos+-span]<0)) */
-  smpl_t res, frac, s0, s1, s2, exactpos = (smpl_t)pos, resold = 100000.;
-  if ((pos > span) && (pos < x->length-span)) {
-    s0 = x->data[0][pos-span];
-    s1 = x->data[0][pos]     ;
-    s2 = x->data[0][pos+span];
-    /* increase frac */
-    for (frac = 0.; frac < 2.; frac = frac + step) {
-      res = aubio_quadfrac(s0, s1, s2, frac);
-      if (res < resold) {
-        resold = res;
-      } else {
-        exactpos += (frac-step)*span - span/2.;
-        break;
-      }
-    }
-  }
-  return exactpos;
+smpl_t vec_quadint(fvec_t * x,uint_t pos, uint_t span) {
+  smpl_t s0, s1, s2;
+  uint_t x0 = (pos < span) ? pos : pos - span;
+  uint_t x2 = (pos + span < x->length) ? pos + span : pos;
+  if (x0 == pos) return (x->data[0][pos] <= x->data[0][x2]) ? pos : x2;
+  if (x2 == pos) return (x->data[0][pos] <= x->data[0][x0]) ? pos : x0;
+  s0 = x->data[0][x0];
+  s1 = x->data[0][pos]     ;
+  s2 = x->data[0][x2];
+  return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
 }
 
 smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
--- a/src/mathutils.h
+++ b/src/mathutils.h
@@ -161,11 +161,8 @@
  */
 smpl_t vec_median(fvec_t * input);
 
-/** finds exact maximum position by quadratic interpolation*/
-smpl_t vec_quadint(fvec_t * x,uint_t pos);
-
-/** finds exact minimum position by quadratic interpolation*/
-smpl_t vec_quadint_min(fvec_t * x,uint_t pos, uint_t span);
+/** finds exact peak index by quadratic interpolation*/
+smpl_t vec_quadint(fvec_t * x, uint_t pos, uint_t span);
 
 /** Quadratic interpolation using Lagrange polynomial.
  *
--- a/src/pitch/pitchmcomb.c
+++ b/src/pitch/pitchmcomb.c
@@ -272,7 +272,7 @@
       if (ispeak) {
         count += ispeak;
         spectral_peaks[count-1].bin = j;
-        spectral_peaks[count-1].ebin = vec_quadint(X,j) - 1.;
+        spectral_peaks[count-1].ebin = vec_quadint(X, j, 1) - 1.;
       }
     }
   return count;
--- a/src/pitch/pitchyin.c
+++ b/src/pitch/pitchyin.c
@@ -106,10 +106,10 @@
     period = tau-3;
     if(tau > 4 && (yin->data[c][period] < tol) && 
         (yin->data[c][period] < yin->data[c][period+1])) {
-      return vec_quadint_min(yin,period,1);
+      return vec_quadint(yin,period,1);
     }
   }
-  return vec_quadint_min(yin,vec_min_elem(yin),1);
+  return vec_quadint(yin,vec_min_elem(yin),1);
   //return 0;
 }
 
--- a/src/pitch/pitchyinfft.c
+++ b/src/pitch/pitchyinfft.c
@@ -128,14 +128,14 @@
     //return vec_quadint_min(yin,tau,1);
     /* additional check for (unlikely) octave doubling in higher frequencies */
     if (tau>35) {
-      return vec_quadint_min(yin,tau,1);
+      return vec_quadint(yin,tau,1);
     } else {
       /* should compare the minimum value of each interpolated peaks */
       halfperiod = FLOOR(tau/2+.5);
       if (yin->data[0][halfperiod] < tol)
-        return vec_quadint_min(yin,halfperiod,1);
+        return vec_quadint(yin,halfperiod,1);
       else
-        return vec_quadint_min(yin,tau,1);
+        return vec_quadint(yin,tau,1);
     }
   } else
     return 0.;