# HG changeset patch
# User Chris Cannam
# Date 1342816629 -3600
# Node ID 45b4401136f67e083bc1655dde3aaefdf57cb238
# Parent  16908c2bd781b6acb21681d5a234f23d36e4a030
More peak location bits & bobs

diff -r 16908c2bd781 -r 45b4401136f6 PeakInterpolator.cpp
--- a/PeakInterpolator.cpp	Fri Jul 20 19:08:38 2012 +0100
+++ b/PeakInterpolator.cpp	Fri Jul 20 21:37:09 2012 +0100
@@ -40,12 +40,27 @@
 }
 
 double
+PeakInterpolator::findPeakLocation(const double *data, int size)
+{
+    double maxval;
+    int maxidx = 0;
+    int i;
+    for (i = 0; i < size; ++i) {
+        if (i == 0 || data[i] > maxval) {
+            maxval = data[i];
+            maxidx = i;
+        }
+    }
+    return findPeakLocation(data, size, maxidx);
+}
+
+double
 PeakInterpolator::findPeakLocation(const double *data, int size, int peakIndex)
 {
-    std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl;
+//    std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl;
 
     if (peakIndex < 1 || peakIndex > size - 2) {
-        std::cerr << "returning " << peakIndex << ", data too short" << std::endl;
+//        std::cerr << "returning " << peakIndex << ", data too short" << std::endl;
         return peakIndex;
     }
 
@@ -63,11 +78,11 @@
     } else {
         y[3] = y[2];
     }
-    std::cerr << "a y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;
+//    std::cerr << "a y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;
     for (int i = 0; i < divisions; ++i) {
         double probe = double(i) / double(divisions);
         double value = cubicInterpolate(y, probe);
-        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex + probe << std::endl;
+//        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex + probe << std::endl;
         if (value > maxval) {
             maxval = value; 
             location = peakIndex + probe;
@@ -82,18 +97,18 @@
     } else {
         y[0] = y[1];
     }
-    std::cerr << "b y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;
+//    std::cerr << "b y: " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] << std::endl;
     for (int i = 0; i < divisions; ++i) {
         double probe = double(i) / double(divisions);
         double value = cubicInterpolate(y, probe);
-        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex - 1 + probe << std::endl;
+//        std::cerr << "probe = " << probe << ", value = " << value << " for location " << peakIndex - 1 + probe << std::endl;
         if (value > maxval) {
             maxval = value; 
             location = peakIndex - 1 + probe;
         }
     }
 
-    std::cerr << "returning " << location << std::endl;
+//    std::cerr << "returning " << location << std::endl;
 
     return location;
 }
diff -r 16908c2bd781 -r 45b4401136f6 PeakInterpolator.h
--- a/PeakInterpolator.h	Fri Jul 20 19:08:38 2012 +0100
+++ b/PeakInterpolator.h	Fri Jul 20 21:37:09 2012 +0100
@@ -32,9 +32,27 @@
     ~PeakInterpolator() { }
 
     /**
+     * Return the interpolated location (i.e. possibly between sample
+     * point indices) of the peak in the given sampled range.
+     *
+     * "The peak" is defined as the (approximate) location of the
+     * maximum of a function interpolating between the points
+     * neighbouring the sample index with the maximum value in the
+     * range.
+     *
+     * If multiple local peak samples in the input range are equal,
+     * i.e. there is more than one apparent peak in the range, the one
+     * with the lowest index will be used. This is the case even if a
+     * later peak would be of superior height after interpolation.
+     */
+    double findPeakLocation(const double *data, int size);
+
+    /**
      * Return the interpolated location (i.e. between sample point
-     * indices) of the peak whose sample is found at peakIndex in a
-     * series of size samples.
+     * indices) of the peak whose nearest sample is found at peakIndex
+     * in the given sampled range. This method allows you to specify
+     * which peak to find, if local rather than global peaks are of
+     * interest.
      */
     double findPeakLocation(const double *data, int size, int peakIndex);
 };
diff -r 16908c2bd781 -r 45b4401136f6 test/TestPeakInterpolator.cpp
--- a/test/TestPeakInterpolator.cpp	Fri Jul 20 19:08:38 2012 +0100
+++ b/test/TestPeakInterpolator.cpp	Fri Jul 20 21:37:09 2012 +0100
@@ -37,6 +37,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 3, 1);
     BOOST_CHECK_EQUAL(result, 1.0);
+    result = p.findPeakLocation(data, 3);
+    BOOST_CHECK_EQUAL(result, 1.0);
 }
 
 BOOST_AUTO_TEST_CASE(peakAtSample_N5)
@@ -45,6 +47,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 5, 2);
     BOOST_CHECK_EQUAL(result, 2.0);
+    result = p.findPeakLocation(data, 5);
+    BOOST_CHECK_EQUAL(result, 2.0);
 }
 
 BOOST_AUTO_TEST_CASE(flat)
@@ -53,6 +57,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 5, 2);
     BOOST_CHECK_EQUAL(result, 2.0);
+    result = p.findPeakLocation(data, 5);
+    BOOST_CHECK_EQUAL(result, 0.0);
 }
 
 BOOST_AUTO_TEST_CASE(multiPeak)
@@ -61,6 +67,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 5, 3);
     BOOST_CHECK_EQUAL(result, 3.0);
+    result = p.findPeakLocation(data, 5);
+    BOOST_CHECK_EQUAL(result, 1.0);
 }
 
 BOOST_AUTO_TEST_CASE(start)
@@ -80,6 +88,10 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 4, 3);
     BOOST_CHECK_EQUAL(result, 3.0);
+    // But when running without a peak location, we expect idx 2 to be
+    // picked as peak, not idx 3, so that will result in interpolation
+    result = p.findPeakLocation(data, 4);
+    BOOST_CHECK(result > 2.0 && result < 3.0);
 }
 
 BOOST_AUTO_TEST_CASE(longHalfway)
@@ -88,6 +100,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 8, 4);
     BOOST_CHECK_EQUAL(result, 3.5);
+    result = p.findPeakLocation(data, 8);
+    BOOST_CHECK_EQUAL(result, 3.5);
 }
 
 BOOST_AUTO_TEST_CASE(shortHalfway)
@@ -96,6 +110,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 4, 1);
     BOOST_CHECK_EQUAL(result, 1.5);
+    result = p.findPeakLocation(data, 4);
+    BOOST_CHECK_EQUAL(result, 1.5);
 }
 
 BOOST_AUTO_TEST_CASE(aboveHalfway)
@@ -104,6 +120,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 4, 2);
     BOOST_CHECK(result > 1.5 && result < 2.0);
+    result = p.findPeakLocation(data, 4);
+    BOOST_CHECK(result > 1.5 && result < 2.0);
 }
 
 BOOST_AUTO_TEST_CASE(belowHalfway)
@@ -112,6 +130,8 @@
     PeakInterpolator p;
     double result = p.findPeakLocation(data, 4, 1);
     BOOST_CHECK(result > 1.0 && result < 1.5);
+    result = p.findPeakLocation(data, 4);
+    BOOST_CHECK(result > 1.0 && result < 1.5);
 }
 
 BOOST_AUTO_TEST_SUITE_END()