changeset 24:0a3c1ecff644

Add rather simplistic cubic interpolation for peak values in a cepstrum output (not yet for pitch tracker). Only gains us 1dp

--- 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) {

--- 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);