src/mathutils.c: merge and fix vec_quadint_min and _max into simpler and more exact...

diff --git a/src/mathutils.c b/src/mathutils.c
index 5a6e9baca06e03b3d8fa59c707c368e0a4a98c0c..4c1d5379331eb1017d0b947c582c617e7765b87a 100644 (file)
--- a/src/mathutils.c
+++ b/src/mathutils.c
@@ -285,49 +285,16 @@ smpl_t vec_median(fvec_t * input) {
   }
 }
 
-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) {

diff --git a/src/mathutils.h b/src/mathutils.h
index 9bf0f5836660bd1cbe752dc65bcd25799dbb2771..6e51586df39ef2c1cce0045e3cb3efe41be9b7d2 100644 (file)
--- a/src/mathutils.h
+++ b/src/mathutils.h
@@ -161,11 +161,8 @@ smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmp,
  */
 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.
  *

diff --git a/src/pitch/pitchmcomb.c b/src/pitch/pitchmcomb.c
index f74351127f5add3c952473a1dbb5684f8cc448cf..5cd23ff995920fbd1230ada6adca9ae57d923b12 100644 (file)
--- a/src/pitch/pitchmcomb.c
+++ b/src/pitch/pitchmcomb.c
@@ -272,7 +272,7 @@ uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t *
       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;

diff --git a/src/pitch/pitchyin.c b/src/pitch/pitchyin.c
index 88b785f693626cfeb0b801fe099688a83bcd8a7f..28da6ec46097a7d036cb2879ded1a860ba46b977 100644 (file)
--- a/src/pitch/pitchyin.c
+++ b/src/pitch/pitchyin.c
@@ -106,10 +106,10 @@ smpl_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t * yin, smpl_t tol){
     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;
 }
 

diff --git a/src/pitch/pitchyinfft.c b/src/pitch/pitchyinfft.c
index 27279db224fabacea8185590a221433526197696..b8c18e02b69a4bf161d9404e6650a835a8d2d16e 100644 (file)
--- a/src/pitch/pitchyinfft.c
+++ b/src/pitch/pitchyinfft.c
@@ -128,14 +128,14 @@ smpl_t aubio_pitchyinfft_detect(aubio_pitchyinfft_t * p, fvec_t * input, smpl_t
     //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.;