comparison PeakInterpolator.cpp @ 40:8f56ef28b0b1

Minor fixes

     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
 
 #include "PeakInterpolator.h"
 
+#include <iostream>
+
 static double cubicInterpolate(const double y[4], double x)
 {
     double a0 = y[3] - y[2] - y[0] + y[1];
     double a1 = y[0] - y[1] - a0;
     double a2 = y[2] - y[0];

 }
 
 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;
     }
 
     double maxval = 0.0;
     double location = peakIndex;

     double y[4];
 
     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;
         }
     }
 
     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;
 }