# HG changeset patch
# User Chris Cannam
# Date 1342777600 -3600
# Node ID 8f56ef28b0b155ab1588ebfb86e054f4e05ac238
# Parent  822cf7b8e070910610b085c923c7a4c37ab3f99d
Minor fixes

diff -r 822cf7b8e070 -r 8f56ef28b0b1 PeakInterpolator.cpp
--- a/PeakInterpolator.cpp	Thu Jul 19 18:10:50 2012 +0100
+++ b/PeakInterpolator.cpp	Fri Jul 20 10:46:40 2012 +0100
@@ -24,6 +24,8 @@
 
 #include "PeakInterpolator.h"
 
+#include <iostream>
+
 static double cubicInterpolate(const double y[4], double x)
 {
     double a0 = y[3] - y[2] - y[0] + y[1];
@@ -40,7 +42,10 @@
 double
 PeakInterpolator::findPeakLocation(const double *data, int size, int peakIndex)
 {
-    if (peakIndex < 2 || peakIndex > size - 3) {
+    std::cerr << "findPeakLocation: size " << size << ", peakIndex " << peakIndex << std::endl;
+
+    if (peakIndex < 1 || peakIndex > size - 2) {
+        std::cerr << "returning " << peakIndex << ", data too short" << std::endl;
         return peakIndex;
     }
 
@@ -53,10 +58,16 @@
     y[0] = data[peakIndex-1];
     y[1] = data[peakIndex];
     y[2] = data[peakIndex+1];
-    y[3] = data[peakIndex+2];
+    if (peakIndex < size - 2) {
+        y[3] = data[peakIndex+2];
+    } else {
+        y[3] = y[2];
+    }
+    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;
         if (value > maxval) {
             maxval = value; 
             location = peakIndex + probe;
@@ -66,16 +77,24 @@
     y[3] = y[2];
     y[2] = y[1];
     y[1] = y[0];
-    y[0] = data[peakIndex-2];
+    if (peakIndex > 1) {
+        y[0] = data[peakIndex-2];
+    } else {
+        y[0] = y[1];
+    }
+    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;
         if (value > maxval) {
             maxval = value; 
             location = peakIndex - 1 + probe;
         }
     }
 
+    std::cerr << "returning " << location << std::endl;
+
     return location;
 }
 
diff -r 822cf7b8e070 -r 8f56ef28b0b1 test/TestPeakInterpolator.cpp
--- a/test/TestPeakInterpolator.cpp	Thu Jul 19 18:10:50 2012 +0100
+++ b/test/TestPeakInterpolator.cpp	Fri Jul 20 10:46:40 2012 +0100
@@ -39,16 +39,6 @@
     BOOST_CHECK_EQUAL(result, 1.0);
 }
 
-BOOST_AUTO_TEST_CASE(peakAtSample_N4)
-{
-    double data[] = { 0.0, 1.0, 2.0, 0.0 };
-    PeakInterpolator p;
-    double result = p.findPeakLocation(data, 4, 2);
-    //!!! Actually, this isn't certain at all I think? It could quite
-    //!!! reasonably be smaller than 2.0
-    BOOST_CHECK_EQUAL(result, 2.0);
-}
-
 BOOST_AUTO_TEST_CASE(peakAtSample_N5)
 {
     double data[] = { 0.0, 1.0, 2.0, 1.0, 0.0 };
@@ -73,5 +63,40 @@
     BOOST_CHECK_EQUAL(result, 3.0);
 }
 
+BOOST_AUTO_TEST_CASE(start)
+{
+    // Can't meaningfully interpolate if we're identifying element 0
+    // as the peak
+    double data[] = { 1.0, 1.0, 0.0, 0.0 };
+    PeakInterpolator p;
+    double result = p.findPeakLocation(data, 4, 0);
+    BOOST_CHECK_EQUAL(result, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(end)
+{
+    // Likewise for the final element
+    double data[] = { 0.0, 0.0, 1.0, 1.0 };
+    PeakInterpolator p;
+    double result = p.findPeakLocation(data, 4, 3);
+    BOOST_CHECK_EQUAL(result, 3.0);
+}
+
+BOOST_AUTO_TEST_CASE(longHalfway)
+{
+    double data[] = { 1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 1.0, 1.0 };
+    PeakInterpolator p;
+    double result = p.findPeakLocation(data, 8, 4);
+    BOOST_CHECK_EQUAL(result, 3.5);
+}
+
+BOOST_AUTO_TEST_CASE(shortHalfway)
+{
+    double data[] = { 1.0, 2.0, 2.0, 1.0 };
+    PeakInterpolator p;
+    double result = p.findPeakLocation(data, 4, 1);
+    BOOST_CHECK_EQUAL(result, 1.5);
+}
+
 BOOST_AUTO_TEST_SUITE_END()