From 9771488a59ce2f276ad63c29bca9723e6e94b96a Mon Sep 17 00:00:00 2001
From: Paul Brossier <piem@altern.org>
Date: Tue, 21 Mar 2006 18:05:23 +0000
Subject: [PATCH] add interpolation to yin add interpolation to yin

---
 src/mathutils.c | 22 ++++++++++++++++++++++
 src/mathutils.h |  5 ++++-
 src/pitchyin.c  | 10 ++++++----
 src/pitchyin.h  |  2 +-
 4 files changed, 33 insertions(+), 6 deletions(-)

diff --git a/src/mathutils.c b/src/mathutils.c
index df2a471a..da5fe927 100644
--- a/src/mathutils.c
+++ b/src/mathutils.c
@@ -330,6 +330,28 @@ smpl_t vec_quadint(fvec_t * x,uint_t pos) {
   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 aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
   smpl_t tmp = s0 + (pf/2.) * (pf * ( s0 - 2.*s1 + s2 ) - 3.*s0 + 4.*s1 - s2);
   return tmp;
diff --git a/src/mathutils.h b/src/mathutils.h
index e4499ff3..a66780dd 100644
--- a/src/mathutils.h
+++ b/src/mathutils.h
@@ -211,9 +211,12 @@ smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmp,
  */
 smpl_t vec_median(fvec_t * input);
 
-/** finds exact peak position by quadratic interpolation*/
+/** 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);
+
 /** Quadratic interpolation using Lagrange polynomial.
  *
  * inspired from ``Comparison of interpolation algorithms in real-time sound
diff --git a/src/pitchyin.c b/src/pitchyin.c
index efd6db99..76505f29 100644
--- a/src/pitchyin.c
+++ b/src/pitchyin.c
@@ -88,8 +88,9 @@ uint_t aubio_pitchyin_getpitch(fvec_t * yin) {
 
 
 /* all the above in one */
-uint_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t * yin, smpl_t tol){
+smpl_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t * yin, smpl_t tol){
 	uint_t c=0,j,tau = 0;
+	sint_t period;
 	smpl_t tmp = 0., tmp2 = 0.;
 	yin->data[c][0] = 1.;
 	for (tau=1;tau<yin->length;tau++)
@@ -102,9 +103,10 @@ uint_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t * yin, smpl_t tol){
 		}
 		tmp2 += yin->data[c][tau];
 		yin->data[c][tau] *= tau/tmp2;
-		if((yin->data[c][tau] < tol) && 
-                                (yin->data[c][tau-1] < yin->data[c][tau])) {
-			return tau-1;
+		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)-1;
 		}
         }
 	return 0;
diff --git a/src/pitchyin.h b/src/pitchyin.h
index b50de01a..b21c0bfe 100644
--- a/src/pitchyin.h
+++ b/src/pitchyin.h
@@ -38,7 +38,7 @@ void aubio_pitchyin_getcum(fvec_t * yinbuf);
 
 uint_t aubio_pitchyin_getpitch(fvec_t *yinbuf);
 
-uint_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t *yinbuf, smpl_t tol);
+smpl_t aubio_pitchyin_getpitchfast(fvec_t * input, fvec_t *yinbuf, smpl_t tol);
 
 #ifdef __cplusplus
 }
-- 
2.26.2