# HG changeset patch
# User Chris Cannam
# Date 1341517856 -3600
# Node ID 0a3c1ecff644b1b1bc6be97c0bfcfeab42f70970
# Parent  1ae8041ae31bd093f16634e28a0af06bf6239205
Add rather simplistic cubic interpolation for peak values in a cepstrum output (not yet for pitch tracker). Only gains us 1dp

diff -r 1ae8041ae31b -r 0a3c1ecff644 SimpleCepstrum.cpp
--- a/SimpleCepstrum.cpp	Thu Jul 05 08:29:20 2012 +0100
+++ b/SimpleCepstrum.cpp	Thu Jul 05 20:50:56 2012 +0100
@@ -257,7 +257,7 @@
 
     d.identifier = "raw_cepstral_peak";
     d.name = "Frequency corresponding to raw cepstral peak";
-    d.description = "Return the frequency whose period corresponds to the quefrency with the maximum value within the specified range of the cepstrum";
+    d.description = "Return the frequency whose period corresponds to the quefrency with the maximum bin value within the specified range of the cepstrum";
     d.unit = "Hz";
     d.hasFixedBinCount = true;
     d.binCount = 1;
@@ -270,6 +270,12 @@
     m_pkOutput = n++;
     outputs.push_back(d);
 
+    d.identifier = "interpolated_peak";
+    d.name = "Interpolated peak frequency";
+    d.description = "Return the frequency whose period corresponds to the quefrency with the maximum bin value within the specified range of the cepstrum, using cubic interpolation to estimate the peak quefrency to finer than single bin resolution";
+    m_ipkOutput = n++;
+    outputs.push_back(d);
+
     d.identifier = "variance";
     d.name = "Variance of cepstral bins in range";
     d.unit = "";
@@ -445,6 +451,71 @@
         result[i] = mean;
     }
 }
+
+double
+SimpleCepstrum::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 a3 = y[1];
+    return
+        a0 * x * x * x +
+        a1 * x * x +
+        a2 * x +
+        a3;
+}
+
+double
+SimpleCepstrum::findInterpolatedPeak(const double *in, int maxbin)
+{
+    if (maxbin < 2 || maxbin > m_bins - 3) {
+        return maxbin;
+    }
+
+    double maxval = 0.0;
+    double maxidx = maxbin;
+
+    const int divisions = 10;
+    double y[4];
+
+    y[0] = in[maxbin-1];
+    y[1] = in[maxbin];
+    y[2] = in[maxbin+1];
+    y[3] = in[maxbin+2];
+    for (int i = 0; i < divisions; ++i) {
+        double probe = double(i) / double(divisions);
+        double value = cubicInterpolate(y, probe);
+        if (value > maxval) {
+            maxval = value; 
+            maxidx = maxbin + probe;
+        }
+    }
+
+    y[3] = y[2];
+    y[2] = y[1];
+    y[1] = y[0];
+    y[0] = in[maxbin-2];
+    for (int i = 0; i < divisions; ++i) {
+        double probe = double(i) / double(divisions);
+        double value = cubicInterpolate(y, probe);
+        if (value > maxval) {
+            maxval = value; 
+            maxidx = maxbin - 1 + probe;
+        }
+    }
+
+/*
+    std::cerr << "centre = " << maxbin << ": ["
+              << in[maxbin-2] << ","
+              << in[maxbin-1] << ","
+              << in[maxbin] << ","
+              << in[maxbin+1] << ","
+              << in[maxbin+2] << "] -> " << maxidx << std::endl;
+*/
+
+    return maxidx;
+}
    
 void
 SimpleCepstrum::addStatisticalOutputs(FeatureSet &fs, const double *data)
@@ -473,12 +544,17 @@
     }
 
     Feature rf;
+    Feature irf;
     if (maxval > 0.0) {
         rf.values.push_back(m_inputSampleRate / (maxbin + m_binFrom));
+        double cimax = findInterpolatedPeak(data, maxbin);
+        irf.values.push_back(m_inputSampleRate / (cimax + m_binFrom));
     } else {
         rf.values.push_back(0);
+        irf.values.push_back(0);
     }
     fs[m_pkOutput].push_back(rf);
+    fs[m_ipkOutput].push_back(irf);
 
     double total = 0;
     for (int i = 0; i < n; ++i) {
diff -r 1ae8041ae31b -r 0a3c1ecff644 SimpleCepstrum.h
--- a/SimpleCepstrum.h	Thu Jul 05 08:29:20 2012 +0100
+++ b/SimpleCepstrum.h	Thu Jul 05 20:50:56 2012 +0100
@@ -83,6 +83,7 @@
     Method m_method;
 
     mutable int m_pkOutput;
+    mutable int m_ipkOutput;
     mutable int m_varOutput;
     mutable int m_p2rOutput;
     mutable int m_cepOutput;
@@ -101,6 +102,8 @@
     double **m_history;
     
     void filter(const double *in, double *out);
+    double cubicInterpolate(const double[4], double);
+    double findInterpolatedPeak(const double *in, int maxbin);
     void fft(unsigned int n, bool inverse,
              double *ri, double *ii, double *ro, double *io);