Pitch detectors now search for local minima, and interpolate for subsample accuracy.

Resamplers still have some issues.

src/yaw-totune/CMakeLists.txt

@@ -1,7 +1,9 @@
 dpf_add_plugin(yaw-totune
   TARGETS vst2
   FILES_DSP
-      dsp.cpp)
+      dsp.cpp
+  FILES_UI
+      ui.cpp)
 
       
 target_include_directories(yaw-totune PUBLIC 

src/yaw-totune/DistrhoPluginInfo.h

@@ -5,11 +5,11 @@
 #define DISTRHO_PLUGIN_NAME  "yaw-totune"
 #define DISTRHO_PLUGIN_URI   "https://yaw.man/plugins/yaw-totune"
 
-// #define DISTRHO_PLUGIN_HAS_UI       1
+#define DISTRHO_PLUGIN_HAS_UI       1
 #define DISTRHO_PLUGIN_IS_RT_SAFE   1
 #define DISTRHO_PLUGIN_NUM_INPUTS   2
 #define DISTRHO_PLUGIN_NUM_OUTPUTS  2
-// #define DISTRHO_UI_USE_NANOVG       1
+#define DISTRHO_UI_USE_NANOVG       1
 
 /*enum Parameters {
     ktpax = 0,

src/yaw-totune/Resampler.hpp

@@ -9,17 +9,32 @@ private:                                         \
     ClassName& operator=(const ClassName&) = delete;
 
 // x that minimizes the quadratic function determined by the three points
-static double secondOrderMinimum(float x0, float y0, float x1, float y1, float x2, float y2)
+static inline double secondOrderMinimum(int x0, float y0, int x1, float y1, int x2, float y2)
 {
-    double x = 0;
-
     // Linear system:
     // 2nd order Vandermonde matrix in x * [a b c]^T = [y0 y1 y2]^T
-    // solve the 2nd order vmonde matrix, then return -b / 2a, discriminant of qdrtic polynomial.
+    // solve the 2nd order vmonde matrix, then return -b / 2a, the discriminant and xmin of qdrtic polynomial.
+    // The x_i are assumed to be nonzero and distinct.
+    double w0 = y0 / ( (x1 - x0) * (x2 - x0) );
+    double w1 = y1 / ( (x0 - x1) * (x2 - x1) );
+    double w2 = y2 / ( (x0 - x2) * (x1 - x2) );
 
-    return x;
+    double a = w0 + w1 + w2;
+    double b = - w0 * (x2 + x1) - w1 * (x0 + x2) - w2 * (x1 + x0);
+
+    if ( a < 0.00001) return x1; //y not sufficiently distinct, can't get numerically meaningful answer.
+    double xmin = -b / (2.0 * a);
+    if ( xmin > x2 ) return x2; //We're looking for a minimum in the given interval, so just clamp it.
+    if ( xmin < x0 ) return x0;
+    return xmin;
 };
 
+static inline double sinc(float x)
+{
+    //Remove discontinuity at zero.
+    if (abs( x ) < 0.001) return 1.f;
+    return sinf(M_PI * x) / (M_PI * x);
+};
 
 static constexpr uint pSize = 10;
 static constexpr size_t resamplerBufferSize = 1 << pSize;
@@ -30,10 +45,6 @@ static constexpr uint dFactor = 1 << dLogFactor;
 static constexpr size_t dBufferSize = 1 << (pSize - dLogFactor);
 static constexpr size_t dBufferMask = dBufferSize - 1;
 
-static constexpr double BIG_DOUBLE = 10000.0;
-
-
-
 template <typename T>
 class Resampler
 {
@@ -45,13 +56,19 @@ private:
     //TODO: refactor into its own class and look up operator[] semantics.
     std::array<T, resamplerBufferSize> array = {};
 
+    //Autocorrelations computed from scratch in a neighbourhood of some estimate.
+    std::array<double, dFactor * 2> autoCorrelations;
+
     // Pointer for pitch correction.
     double readIdx = 0;
     // Pointer for pitch detection.
     uint writeIdx = 0;
 
+    // Simple one-pole low pass filter for smoothing over period mismatches.
+    double lpf = 0.0;
+
     // Calculate autocorrelation from scratch for given period
-    inline double ac(const uint per)
+    inline double const ac(const uint per)
     {
         double ac = 0;
         uint idx = writeIdx;
@@ -60,7 +77,7 @@ private:
             T x = array[idx & resamplerMask];
             T y = array[(idx - per) & resamplerMask];
 
-            ac += x * x + x * y;
+            ac += (x - y) * (x - y);
             --idx;
         }
         return ac;
@@ -70,21 +87,43 @@ public:
 
     Resampler(){};
 
-    // Detect accurate period from scratch near a target period.
+    inline double const getAccuratePeriod(uint firstEstimate, uint secondEstimate)
-    inline double detectPeriodNear(const uint nearPeriod)
     {
-        double minAC = BIG_DOUBLE;
+        //First try the two candidates for the appropriate nbhd.
-        double period = nearPeriod;
+        uint estimate = firstEstimate;
-        for (uint per = nearPeriod - dFactor / 2; per < nearPeriod + dFactor / 2; ++per)
+        double fac = ac(firstEstimate);
+        if ( secondEstimate )
         {
-            double curAC = ac(per);
+            double altAC = ac(secondEstimate);
-            if (curAC < minAC)
+            if( fac > altAC )
             {
-                minAC = curAC;
+                fac = altAC;
-                period = per;
+                estimate = secondEstimate;
             }
         }
-        return period;
+
+        //Get autocorrelations in the nbhd, find the minimum.
+        uint minIdx = 0;
+        int p = estimate - autoCorrelations.size() / 2;
+        for (uint i = 0; i < autoCorrelations.size(); ++i, ++p)
+        {
+            autoCorrelations[i] = ac(p);
+            if (autoCorrelations[i] <= fac ) {
+                minIdx = i;
+                estimate = p;
+                fac = autoCorrelations[i];
+            }
+        }
+
+        //Interpolate to get accuracy greater than one sample (yes, you need it!)
+        uint x0 = estimate - 1;
+        uint x1 = estimate;
+        uint x2 = estimate + 1;
+        double y0 = (minIdx == 0) ? ac(x0) : autoCorrelations[minIdx - 1];
+        double y1 = autoCorrelations[minIdx];
+        double y2 = ((minIdx + 1) >= autoCorrelations.size()) ? ac(x2) : autoCorrelations[minIdx + 1];
+
+        return secondOrderMinimum(x0, y0, x1, y1, x2, y2);
     }
 
     inline void write(const T *const input, const uint32_t frames)
@@ -96,6 +135,22 @@ public:
         };
     }
 
+
+
+    //Eight tap windowed sinc filter.
+    //Interpolate at readIdx - 4.0.
+    inline T const read(const double rate)
+    {
+        const uint idx = static_cast<uint>(readIdx);
+        const double frac = readIdx - idx;
+        T val = 0.f;
+        for( uint i = 0; i < 8; ++i )
+        {
+            val += sinc((idx - i - readIdx + 3.0) / rate) * array[(idx - i) & resamplerMask];
+        }
+        return val;
+    }
+
     // Arguments:
     // output buffer and length
     // ratio of new sample rate to old sample rate
@@ -106,19 +161,14 @@ public:
         //Bounds of read index.
         double end = writeIdx;
         while ( readIdx > end ) { end += resamplerBufferSize; };
-        double start = end - period;
+        const double start = end - period;
-
 
         for (uint32_t i = 0; i < frames; ++i)
         {
-            // TODO: add interpolation filter. Linear for now.
+            output[i] = read(rate);
-            uint idx = static_cast<uint>(readIdx);
-            float frac = readIdx - idx;
-            *output = (1.f - frac) * array[(idx - 1) & resamplerMask] + frac * array[idx & resamplerMask];
-            
             readIdx += rate;
 
-            while ( readIdx < start )
+            while ( rate < 1.0 && readIdx < start )
             {
                 readIdx += period;
             }
@@ -128,9 +178,10 @@ public:
                 readIdx -= period;
             }
 
-            ++output;
-
         };
+
+        //Prevent read index overflow / precision loss.
+        while (readIdx > resamplerBufferSize) {readIdx -= resamplerBufferSize;};
     }
 
     DISTRHO_DECLARE_NON_COPYABLE(Resampler)
@@ -144,8 +195,8 @@ class Detector
     // Pointer for pitch detection.
     uint writeIdx = 0; 
     // Detectable periods.
-    std::array<T, dBufferSize / 2> squares = {};
+    std::array<double, dBufferSize / 2> squares = {};
-    std::array<T, dBufferSize / 2> crosses = {};
+    std::array<double, dBufferSize / 2> crosses = {};
     static constexpr uint minPeriod = 16;
     static constexpr uint maxPeriod = dBufferSize / 2;
 
@@ -165,40 +216,54 @@ public:
         };
     }
 
-    // Incrementally detect all possible pitches. Return coarse pitch match.
+    std::array<uint, 2> detectPitch(uint32_t frames)
-    inline uint detectPitch(uint32_t frames)
     {
         frames /= dFactor;
         uint idx = writeIdx - frames;
-        bool a = false;
+        std::array<uint, 2> pitches{};
+
+        //Update autocorrelations.
         while (frames)
         {
-            for (uint per = 1; per < squares.size(); ++per)
+            for (uint per = 3; per < squares.size(); ++per)
             {
                 T x = array[idx & dBufferMask];
                 T y = array[(idx - per) & dBufferMask];
                 T z = array[(idx - 2 * per) & dBufferMask];
+                //Autocorrelation and RMS computation.
                 squares[per] += x * x - z * z;
                 crosses[per] += x * y - y * z;
-                if(squares[per] - 2.0 * crosses[per] < -0.1)
-                    break;
             }
             --frames;
             ++idx;
         }
 
-        double least = BIG_DOUBLE;
+        //Search for local minima of autocorrelation function.
-        uint per = 0;
+        for ( uint per = 3; per < squares.size() - 1; ++per)
-        for (uint i = 1; i < squares.size(); ++i)
         {
-            double uac = squares[i] - 2 * crosses[i];
+            if ( 0.49 * squares[per] < crosses[per] )
-            if (uac < least)
             {
-                least = uac;
+                double a = 2.0 * crosses[per - 1] - squares[per - 1];
-                per = i;
+                double b = 2.0 * crosses[per] - squares[per];
+                double c = 2.0 * crosses[per + 1] - squares[per + 1];
+                if ( (b > a) && (b > c) )
+                {
+                    if (pitches[0])
+                    {
+                        pitches[1] = dFactor * per;
+                        return pitches;
+                    }
+                    else
+                    {
+                        pitches[0] = dFactor * per;
+                        if (per * 2 > squares.size() ) return pitches;
                     }
                 }
-        return per * dFactor;
+            }
+        }
+
+        return pitches;
+
     }
 
      DISTRHO_DECLARE_NON_COPYABLE(Detector)

src/yaw-totune/dsp.cpp

@@ -40,39 +40,50 @@ protected:
 
 	float getParameterValue(uint32_t index) const override
 	{
-		return curPitch;
+		return period;
 	}
 
 	void setParameterValue(uint32_t idx, float val) override
 	{
 	}
 
+	double debugTone(float* output, uint32_t frames, float period, float phase)
+	{
+		for(uint i = 0; i < frames; ++i) {
+			*output = 0.5 * sinf(	M_PI * 2.0 * phase	);
+			phase += 1.f / period;
+			phase -= (int)phase;
+			++output;
+		};
+
+		return phase;
+	}
+
 	void run(const float **inputs, float **outputs, uint32_t frames) override
 	{
 		for (int chn = 0; chn < 2; ++chn)
 		{
 			resamplers[chn].write(inputs[chn], frames);
 			detectors[chn].downsample(inputs[chn], frames);
-			uint dpitch = detectors[chn].detectPitch(frames);
+			auto detectedPitches = detectors[chn].detectPitch(frames);
-			curPitch = dpitch;
 			double taux = 1.0;
-			double period = 1.0;
+			period = 1023.0;
-			if (dpitch)
+			if (detectedPitches[0])
 			{
-				period = resamplers[chn].detectPeriodNear(dpitch);
+				period = resamplers[chn].getAccuratePeriod(detectedPitches[0], detectedPitches[1]);
 				double correctPeriod = scale.getNearestPeriod(period);
 				if( correctPeriod > 1.0 ) taux = period / correctPeriod;
-			}
+			};
-			//resamplers[chn].resample(outputs[chn], frames, taux, period);
+			resamplers[chn].resample(outputs[chn], frames, 0.95, 109);
-			resamplers[chn].resample(outputs[chn], frames, 1.3, period);
 		}
 	}
 
 private:
-	uint curPitch = 0;
+	double period = 1023.0;
 	double hzNyq;
 	std::array<Resampler<float>, 2> resamplers;
 	std::array<Detector<float>, 2> detectors;
+	std::array<double, 2> debugPhases = {0};
 	Scale scale{440.0};
 	DISTRHO_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR(PitchCorrector)
 };

src/yaw-totune/ui.cpp

@@ -0,0 +1,98 @@
+/*
+ * DISTRHO Plugin Framework (DPF)
+ * Copyright (C) 2012-2021 Filipe Coelho <falktx@falktx.com>
+ *
+ * Permission to use, copy, modify, and/or distribute this software for any purpose with
+ * or without fee is hereby granted, provided that the above copyright notice and this
+ * permission notice appear in all copies.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD
+ * TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN
+ * NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
+ * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+ * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
+ * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
+ */
+
+#include "DistrhoUI.hpp"
+#include <format>
+
+START_NAMESPACE_DISTRHO
+
+// -----------------------------------------------------------------------------------------------------------
+
+class InfoExampleUI : public UI
+{
+    static const uint kInitialWidth = 405;
+    static const uint kInitialHeight = 256;
+
+public:
+    InfoExampleUI()
+        : UI(kInitialWidth, kInitialHeight)
+    {
+
+#ifdef DGL_NO_SHARED_RESOURCES
+        createFontFromFile("sans", "/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans.ttf");
+#else
+        loadSharedResources();
+#endif
+    }
+
+protected:
+   /* --------------------------------------------------------------------------------------------------------
+    * DSP/Plugin Callbacks */
+
+   /**
+      A parameter has changed on the plugin side.
+      This is called by the host to inform the UI about parameter changes.
+    */
+    void parameterChanged(uint32_t index, float value) override
+    {
+        period = value;
+        repaint();
+    }
+
+   /* --------------------------------------------------------------------------------------------------------
+    * DSP/Plugin Callbacks (optional) */
+
+
+   /* --------------------------------------------------------------------------------------------------------
+    * Widget Callbacks */
+
+   /**
+      The NanoVG drawing function.
+    */
+    void onNanoDisplay() override
+    {
+        	// Numerical feedback.
+	beginPath();
+	fillColor(200, 200, 200);
+	textBox(0.f, 15.f, 250.f,
+			std::format("p: {:.3f}\n",
+						period)
+				.c_str(),
+			nullptr);
+	closePath();
+    }
+
+    // -------------------------------------------------------------------------------------------------------
+
+private:
+    // Parameters
+    float period = 0.f;
+    double fSampleRate;
+
+    DISTRHO_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR(InfoExampleUI)
+};
+
+/* ------------------------------------------------------------------------------------------------------------
+ * UI entry point, called by DPF to create a new UI instance. */
+
+UI* createUI()
+{
+    return new InfoExampleUI();
+}
+
+// -----------------------------------------------------------------------------------------------------------
+
+END_NAMESPACE_DISTRHO