src/mathutils.{c,h}: add fvec_quadratic_peak_pos, a fixed replacement for fvec_quadint

diff --git a/src/mathutils.c b/src/mathutils.c
index b0a958b9dde1e2cb977c69ce563121e72b4d56a5..7d2424e05e2ff12b1625b1400e54ae9c653b0185 100644 (file)
--- a/src/mathutils.c
+++ b/src/mathutils.c
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
+  Copyright (C) 2003-2013 Paul Brossier <piem@aubio.org>
 
   This file is part of aubio.
 
@@ -379,6 +379,18 @@ smpl_t fvec_quadint (fvec_t * x, uint_t pos) {
   return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
 }
 
+smpl_t fvec_quadratic_peak_pos (fvec_t * x, uint_t pos) {
+  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[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 + 0.5 * (s0 - s2 ) / (s0 - 2.* s1 + s2);
+}
+
 uint_t fvec_peakpick(fvec_t * onset, uint_t pos) {
   uint_t tmp=0;
   tmp = (onset->data[pos] > onset->data[pos-1]

diff --git a/src/mathutils.h b/src/mathutils.h
index 01142dd33d05fce03e73646a42564018c2baacfb..6d6c8fbcfc82063c9e66f06d890a774cc43f58be 100644 (file)
--- a/src/mathutils.h
+++ b/src/mathutils.h
@@ -235,6 +235,21 @@ smpl_t fvec_median (fvec_t * v);
 /** finds exact peak index by quadratic interpolation*/
 smpl_t fvec_quadint (fvec_t * x, uint_t pos);
 
+/** finds exact peak index by quadratic interpolation
+
+  See [Quadratic Interpolation of Spectral
+  Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
+  by Julius O. Smith III
+
+  \f$ p_{frac} = \frac{1}{2} \frac {x[p-1] - x[p+1]} {x[p-1] - 2 x[p] + x[p+1]} \in [ -.5, .5] \f$
+
+  \param x vector to get the interpolated peak position from
+  \param p index of the peak in vector `x`
+  \return \f$ p + p_{frac} \f$ exact peak position of interpolated maximum or minimum
+
+*/
+smpl_t fvec_quadratic_peak_pos (fvec_t * x, uint_t p);
+
 /** Quadratic interpolation using Lagrange polynomial.
  
   Inspired from ``Comparison of interpolation algorithms in real-time sound