Add windowed sinc interpolation

common.mk

@@ -87,6 +87,7 @@ SFIZZ_SOURCES = \
 	src/sfizz/FlexEGDescription.cpp \
 	src/sfizz/FlexEnvelope.cpp \
 	src/sfizz/FloatEnvelopes.cpp \
+	src/sfizz/Interpolators.cpp \
 	src/sfizz/Logger.cpp \
 	src/sfizz/LFO.cpp \
 	src/sfizz/LFODescription.cpp \
@@ -121,7 +122,8 @@ SFIZZ_SOURCES = \
 	src/sfizz/Voice.cpp \
 	src/sfizz/VoiceManager.cpp \
 	src/sfizz/VoiceStealing.cpp \
-	src/sfizz/Wavetables.cpp
+	src/sfizz/Wavetables.cpp \
+	src/sfizz/WindowedSinc.cpp
 
 ### Other internal
 

src/CMakeLists.txt

@@ -112,6 +112,8 @@ set(SFIZZ_HEADERS
     sfizz/VoiceManager.h
     sfizz/VoiceStealing.h
     sfizz/Wavetables.h
+    sfizz/WindowedSinc.h
+    sfizz/WindowedSinc.hpp
     sfizz.h
     sfizz.hpp)
 
@@ -151,6 +153,8 @@ set(SFIZZ_SOURCES
     sfizz/BeatClock.cpp
     sfizz/Metronome.cpp
     sfizz/SynthMessaging.cpp
+    sfizz/WindowedSinc.cpp
+    sfizz/Interpolators.cpp
     sfizz/modulations/ModId.cpp
     sfizz/modulations/ModKey.cpp
     sfizz/modulations/ModKeyHash.cpp

src/sfizz/Interpolators.cpp

@@ -0,0 +1,23 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#include "Interpolators.h"
+
+namespace sfz {
+
+void initializeInterpolators()
+{
+    SincInterpolatorTraits<8>::initialize();
+    SincInterpolatorTraits<12>::initialize();
+    SincInterpolatorTraits<16>::initialize();
+    SincInterpolatorTraits<24>::initialize();
+    SincInterpolatorTraits<36>::initialize();
+    SincInterpolatorTraits<48>::initialize();
+    SincInterpolatorTraits<60>::initialize();
+    SincInterpolatorTraits<72>::initialize();
+}
+
+} // namespace sfz

src/sfizz/Interpolators.h

@@ -17,8 +17,36 @@ enum InterpolatorModel : int {
     kInterpolatorHermite3,
     // a B-spline 3rd order interpolator
     kInterpolatorBspline3,
+    // a windowed-sinc 8-point interpolator
+    kInterpolatorSinc8,
+    // a windowed-sinc 12-point interpolator
+    kInterpolatorSinc12,
+    // a windowed-sinc 16-point interpolator
+    kInterpolatorSinc16,
+    // a windowed-sinc 24-point interpolator
+    kInterpolatorSinc24,
+    // a windowed-sinc 36-point interpolator
+    kInterpolatorSinc36,
+    // a windowed-sinc 48-point interpolator
+    kInterpolatorSinc48,
+    // a windowed-sinc 60-point interpolator
+    kInterpolatorSinc60,
+    // a windowed-sinc 72-point interpolator
+    kInterpolatorSinc72,
 };
 
+/**
+ * @brief Initialize interpolators
+ *
+ * This precomputes windowed-sinc tables globally.
+ * It needs to be called at least once, before using the windowed-sinc models.
+ *
+ * These are not computed at static initialization time, to prevent slowing down
+ * an audio plugin library scan (eg. VST). The static-local-variable method is
+ * avoided also, because we don't want this overhead on a frame-by-frame basis.
+ */
+void initializeInterpolators();
+
 /**
  * @brief Interpolate from a vector of values
  *

src/sfizz/Interpolators.hpp

@@ -5,6 +5,7 @@
 // If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
 
 #include "Interpolators.h"
+#include "WindowedSinc.h"
 #include "MathHelpers.h"
 #include "SIMDConfig.h"
 
@@ -133,4 +134,135 @@ public:
     }
 };
 
+//------------------------------------------------------------------------------
+// Windowed sinc
+
+namespace SincInterpolatorDetail {
+    // See sfizz wiki page "Resampling".
+    constexpr size_t PointsMin = 8;
+    constexpr size_t PointsMax = 72;
+
+    // Adjust Kaiser window Beta as necessary.
+    constexpr double BetaMin = 6.0;
+    constexpr double BetaMax = 10.0;
+
+    constexpr double getBetaForNumPoints(size_t points)
+    {
+        return BetaMin + (BetaMax - BetaMin) *
+            (double(points - PointsMin) / double(PointsMax - PointsMin));
+    }
+
+    constexpr size_t getTableSizeForNumPoints(size_t /*points*/)
+    {
+        return 1u << 16;
+    }
+}
+
+///
+template <size_t Points>
+struct SincInterpolatorTraits {
+    static_assert(Points == 8  || Points == 12 || Points == 16 ||
+                  Points == 24 || Points == 36 || Points == 48 ||
+                  Points == 60 || Points == 72,
+                  "Windowed sinc size is not acceptable");
+
+    enum {
+        TableSize = SincInterpolatorDetail::getTableSizeForNumPoints(Points)
+    };
+
+    static void initialize()
+    {
+        static const FixedWindowedSinc<Points, TableSize> globalInstance(
+            SincInterpolatorDetail::getBetaForNumPoints(Points));
+        windowedSinc = &globalInstance;
+    }
+
+    static const FixedWindowedSinc<Points, TableSize>* windowedSinc;
+};
+
+template <size_t Points>
+const FixedWindowedSinc<Points, SincInterpolatorTraits<Points>::TableSize>*
+SincInterpolatorTraits<Points>::windowedSinc = nullptr;
+
+///
+template <class R, size_t Points>
+class SincInterpolator;
+
+//------------------------------------------------------------------------------
+// Windowed sinc any order, SSE specialization
+#if SFIZZ_HAVE_SSE
+template <size_t Points>
+class SincInterpolator<float, Points>
+{
+public:
+    static_assert(Points % 4 == 0, "Windowed sinc must be multiple of 4");
+
+    static inline float process(const float* values, float coeff)
+    {
+        const auto &ws = *SincInterpolatorTraits<Points>::windowedSinc;
+
+        int j0 = 1 - int(Points) / 2;
+
+        __m128 h[Points / 4];
+        for (int i = 0; i < int(Points); ++i)
+            reinterpret_cast<float*>(h)[i] = ws.getUnchecked(j0 - coeff + i);
+
+        __m128 y = _mm_mul_ps(h[0], _mm_loadu_ps(&values[j0]));
+        for (int i = 1; i < int(Points / 4); ++i)
+            y = _mm_add_ps(y, _mm_mul_ps(h[i], _mm_loadu_ps(&values[j0 + 4 * i])));
+
+        // sum 4 to 1
+        __m128 xmm0 = y;
+        __m128 xmm1 = _mm_shuffle_ps(xmm0, xmm0, 0xe5);
+        __m128 xmm2 = _mm_movehl_ps(xmm0, xmm0);
+        xmm1 = _mm_add_ss(xmm1, xmm0);
+        xmm0 = _mm_shuffle_ps(xmm0, xmm0, 0xe7);
+        xmm2 = _mm_add_ss(xmm2, xmm1);
+        xmm0 = _mm_add_ss(xmm0, xmm2);
+        return _mm_cvtss_f32(xmm0);
+    }
+};
+#endif
+
+//------------------------------------------------------------------------------
+// Windowed sinc any order, generic
+template <class R, size_t Points>
+class SincInterpolator
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        const auto &ws = *SincInterpolatorTraits<Points>::windowedSinc;
+
+        int j0 = 1 - int(Points) / 2;
+
+        R h[Points];
+        for (int i = 0; i < int(Points); ++i)
+            h[i] = R(ws.getUnchecked(j0 - coeff + i));
+
+        R y = h[0] * values[j0];
+        for (int i = 1; i < int(Points); ++i)
+            y += h[i] * values[j0 + i];
+
+        return y;
+    }
+};
+
+template <class R>
+class Interpolator<kInterpolatorSinc8, R> : public SincInterpolator<R, 8> {};
+template <class R>
+class Interpolator<kInterpolatorSinc12, R> : public SincInterpolator<R, 12> {};
+template <class R>
+class Interpolator<kInterpolatorSinc16, R> : public SincInterpolator<R, 16> {};
+template <class R>
+class Interpolator<kInterpolatorSinc24, R> : public SincInterpolator<R, 24> {};
+template <class R>
+class Interpolator<kInterpolatorSinc36, R> : public SincInterpolator<R, 36> {};
+template <class R>
+class Interpolator<kInterpolatorSinc48, R> : public SincInterpolator<R, 48> {};
+template <class R>
+class Interpolator<kInterpolatorSinc60, R> : public SincInterpolator<R, 60> {};
+template <class R>
+class Interpolator<kInterpolatorSinc72, R> : public SincInterpolator<R, 72> {};
+
 } // namespace sfz

src/sfizz/Synth.cpp

@@ -21,6 +21,7 @@
 #include "utility/SpinMutex.h"
 #include "utility/XmlHelpers.h"
 #include "Voice.h"
+#include "Interpolators.h"
 #include <absl/algorithm/container.h>
 #include <absl/memory/memory.h>
 #include <absl/strings/str_replace.h>
@@ -48,6 +49,7 @@ Synth::~Synth()
 Synth::Impl::Impl()
 {
     initializeSIMDDispatchers();
+    initializeInterpolators();
 
     const std::lock_guard<SpinMutex> disableCallback { callbackGuard_ };
     parser_.setListener(this);

src/sfizz/Voice.cpp

@@ -1149,18 +1149,20 @@ void Voice::Impl::fillInterpolatedWithQuality(
     absl::Span<const int> indices, absl::Span<const float> coeffs,
     absl::Span<const float> addingGains, int quality)
 {
-    switch (quality) {
-    default:
-        if (quality > 2)
-            goto high; // TODO sinc, not implemented
-        // fall through
+    switch (clamp(quality, 0, 10)) {
+    case 0:
+        {
+            constexpr auto itp = kInterpolatorNearest;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
     case 1:
         {
             constexpr auto itp = kInterpolatorLinear;
             fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
         }
         break;
-    case 2: high:
+    case 2:
         {
 #if 0
             // B-spline response has faster decay of aliasing, but not zero-crossings at integer positions
@@ -1172,6 +1174,54 @@ void Voice::Impl::fillInterpolatedWithQuality(
             fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
         }
         break;
+    case 3:
+        {
+            constexpr auto itp = kInterpolatorSinc8;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 4:
+        {
+            constexpr auto itp = kInterpolatorSinc12;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 5:
+        {
+            constexpr auto itp = kInterpolatorSinc16;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 6:
+        {
+            constexpr auto itp = kInterpolatorSinc24;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 7:
+        {
+            constexpr auto itp = kInterpolatorSinc36;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 8:
+        {
+            constexpr auto itp = kInterpolatorSinc48;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 9:
+        {
+            constexpr auto itp = kInterpolatorSinc60;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
+    case 10:
+        {
+            constexpr auto itp = kInterpolatorSinc72;
+            fillInterpolated<itp, Adding>(source, dest, indices, coeffs, addingGains);
+        }
+        break;
     }
 }
 

src/sfizz/WindowedSinc.cpp

@@ -0,0 +1,39 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#include "WindowedSinc.h"
+#include "MathHelpers.h"
+#include <absl/memory/memory.h>
+
+namespace sfz {
+
+void WindowedSincDetail::calculateTable(absl::Span<float> table, size_t sincExtent, double beta, size_t extra)
+{
+    size_t tableSize = table.size();
+
+    auto window = absl::make_unique<float[]>(tableSize);
+    kaiserWindow(beta, absl::MakeSpan(window.get(), tableSize));
+
+    // table domain [-N/2:+N/2]
+    double scale = sincExtent / static_cast<double>(tableSize - 1);
+    double offset = sincExtent / -2.0;
+
+    for (size_t i = 0; i < tableSize; ++i) {
+        double x = i * scale + offset;
+        table[i] = window[i] * normalizedSinc(x);
+    }
+
+    for (size_t i = 0; i < extra; ++i)
+        table[extra + i] = table[tableSize - 1];
+}
+
+double WindowedSincDetail::calculateExact(double x, size_t sincExtent, double beta)
+{
+    return normalizedSinc(x) *
+        kaiserWindowSinglePoint(beta, (x + sincExtent / 2.0f) / sincExtent);
+}
+
+} // namespace sfz

src/sfizz/WindowedSinc.h

@@ -0,0 +1,111 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#pragma once
+#include <absl/types/span.h>
+#include <memory>
+
+namespace sfz {
+
+namespace WindowedSincDetail {
+    void calculateTable(absl::Span<float> table, size_t sincExtent, double beta, size_t extra);
+    double calculateExact(double x, size_t sincExtent, double beta);
+};
+
+template <class T>
+class AbstractWindowedSinc {
+protected:
+    explicit AbstractWindowedSinc(double beta) noexcept : beta_(beta) {}
+
+public:
+    virtual ~AbstractWindowedSinc() noexcept {}
+
+    // interpolate f(x), where x must be in domain [-Points/2:+Points/2]
+    float getUnchecked(float x) const noexcept;
+
+    // calculate exact f(x), where x must be in domain [-Points/2:+Points/2]
+    double getExact(double x) const noexcept;
+
+    // get the Kaiser window Beta parameter
+    double getBeta() const noexcept { return beta_; }
+
+protected:
+    void fillTable() noexcept;
+
+    // allows interpolating f(Points/2), provided for safety
+    enum { TableExtra = 1 };
+
+private:
+    double beta_ {};
+};
+
+/**
+ * @brief Windowed-sinc using fixed compile-time parameters
+ * This can help to save some instructions in a resampler loop.
+ */
+template <size_t Points, size_t TableSize>
+class FixedWindowedSinc final :
+        public AbstractWindowedSinc<FixedWindowedSinc<Points, TableSize>> {
+public:
+    using Self = FixedWindowedSinc<Points, TableSize>;
+    using Super = AbstractWindowedSinc<Self>;
+
+    explicit FixedWindowedSinc(double beta) : Super(beta) { Super::fillTable(); }
+
+    // the number of points where this sinc will be evaluated (zero crossings + 1)
+    static constexpr size_t getNumPoints() noexcept { return Points; }
+
+    // the size of the lookup table
+    static constexpr size_t getTableSize() noexcept { return TableSize; }
+
+    // the lookup table
+    const float* getTablePointer() const noexcept { return table_; }
+
+    // the lookup table
+    absl::Span<const float> getTableSpan() const noexcept { return absl::MakeConstSpan(table_, TableSize); }
+
+protected:
+    using Super::TableExtra;
+
+private:
+    float table_[TableSize + TableExtra];
+};
+
+/**
+ * @brief Windowed-sinc using run-time parameters
+ */
+class WindowedSinc final : public AbstractWindowedSinc<WindowedSinc>
+{
+public:
+    using Self = WindowedSinc;
+    using Super = AbstractWindowedSinc<WindowedSinc>;
+
+    WindowedSinc(size_t points, size_t tableSize, double beta)
+        : Super(beta), points_(points), tableSize_(tableSize),
+          table_(new float[tableSize + TableExtra])
+    { Super::fillTable(); }
+
+    // the number of points where this sinc will be evaluated (zero crossings + 1)
+    size_t getNumPoints() const noexcept { return points_; }
+
+    // the size of the lookup table
+    size_t getTableSize() const noexcept { return tableSize_; }
+
+    // the lookup table
+    const float* getTablePointer() const noexcept { return table_.get(); }
+
+    // the lookup table
+    absl::Span<const float> getTableSpan() const noexcept { return absl::MakeConstSpan(table_.get(), tableSize_); }
+
+private:
+    size_t points_ {};
+    size_t tableSize_ {};
+    std::unique_ptr<float[]> table_;
+};
+
+} // namespace sfz
+
+#include "WindowedSinc.hpp"

src/sfizz/WindowedSinc.hpp

@@ -0,0 +1,45 @@
+// SPDX-License-Identifier: BSD-2-Clause
+
+// This code is part of the sfizz library and is licensed under a BSD 2-clause
+// license. You should have receive a LICENSE.md file along with the code.
+// If not, contact the sfizz maintainers at https://github.com/sfztools/sfizz
+
+#pragma once
+#include "WindowedSinc.h"
+
+namespace sfz {
+
+template <class T>
+inline void AbstractWindowedSinc<T>::fillTable() noexcept
+{
+    float* table = const_cast<float*>(static_cast<T*>(this)->getTablePointer());
+    size_t points = static_cast<T*>(this)->getNumPoints();
+    size_t tableSize = static_cast<T*>(this)->getTableSize();
+
+    WindowedSincDetail::calculateTable(
+        absl::MakeSpan(table, tableSize), points, beta_, TableExtra);
+}
+
+template <class T>
+inline float AbstractWindowedSinc<T>::getUnchecked(float x) const noexcept
+{
+    const float* table = static_cast<const T*>(this)->getTablePointer();
+    size_t points = static_cast<const T*>(this)->getNumPoints();
+    size_t tableSize = static_cast<const T*>(this)->getTableSize();
+
+    float ix = (x + points / 2.0f) * ((tableSize - 1) / points);
+    int i0 = static_cast<int>(ix);
+    float mu = ix - i0;
+    float y0 = table[i0];
+    float dy = table[i0 + 1] - y0;
+    return y0 + mu * dy;
+}
+
+template <class T>
+inline double AbstractWindowedSinc<T>::getExact(double x) const noexcept
+{
+    size_t points = static_cast<const T*>(this)->getNumPoints();
+    return WindowedSincDetail::calculateExact(x, points, beta_);
+}
+
+} // namespace sfz

tests/InterpolatorsT.cpp

@@ -70,3 +70,60 @@ TEST_CASE("[Interpolators] Squares")
             == Approx(expected).margin(1e-2));
     }
 }
+
+template <class WS>
+static std::pair<double, double> windowedSincError(WS& ws, double step = 0.1, bool verbose = false)
+{
+    size_t points = ws.getNumPoints();
+    double x1 = points / -2.0;
+    double x2 = points / +2.0;
+    double maxAbsErr = 0.0;
+    double meanAbsErr = 0.0;
+    //double meanSquareErr = 0.0;
+
+    double x;
+    size_t n;
+    for (n = 0; (x = x1 + n * step) < x2; ++n) {
+        double val = ws.getUnchecked(x);
+        double ref = ws.getExact(x);
+        double absErr = std::fabs(val - ref);
+        maxAbsErr = std::max(maxAbsErr, absErr);
+        meanAbsErr += absErr;
+        //meanSquareErr += absErr * absErr;
+    }
+    meanAbsErr /= n;
+    //meanSquareErr /= n;
+
+    if (verbose) {
+        std::cerr << "MaxAbsErr=" << maxAbsErr
+                  << " MeanAbsErr=" << meanAbsErr
+                  //<< " MeanSquareErr=" << meanSquareErr
+                  << " with Points=" << points
+                  << " TableSize=" << ws.getTableSize()
+                  << "\n";
+    }
+
+    return { maxAbsErr, meanAbsErr };
+}
+
+TEST_CASE("[Interpolators] Windowed sinc precision")
+{
+    sfz::initializeInterpolators();
+
+    double maxAbsTolerance = 5e-2;
+    double meanAbsTolerance = 1e-3;
+
+    auto Check = [=](std::pair<double, double> maxAndMeanErr) {
+        REQUIRE(maxAndMeanErr.first < maxAbsTolerance);
+        REQUIRE(maxAndMeanErr.second < meanAbsTolerance);
+    };
+
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<8>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<12>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<16>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<24>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<36>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<48>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<60>::windowedSinc));
+    Check(windowedSincError(*sfz::SincInterpolatorTraits<72>::windowedSinc));
+}