Use a 4-point Gaussian filter for cubic resampling

alc/alu.cpp

@@ -260,7 +260,7 @@ ResamplerFunc PrepareResampler(Resampler resampler, uint increment, InterpState
     case Resampler::Linear:
         break;
     case Resampler::Cubic:
-        state->emplace<CubicState>().filter = gCubicSpline.Tab.data();
+        state->emplace<CubicState>().filter = gGaussianFilter.Tab.data();
         break;
     case Resampler::FastBSinc12:
     case Resampler::BSinc12:

alc/effects/reverb.cpp

@@ -68,21 +68,37 @@ struct CubicFilter {
 
     std::array<float,sTableSteps*2 + 1> mFilter{};
 
+    [[nodiscard]]
+    auto GetCoeff(size_t idx) noexcept -> double
+    {
+        static constexpr double IndexScale{512.0 / double{sTableSteps*2}};
+        const double k{0.5 + static_cast<double>(idx)*IndexScale};
+        if(k > 512.0) return 0.0;
+        const double s{ std::sin(al::numbers::pi*1.280/1024 * k)};
+        const double t{(std::cos(al::numbers::pi*2.000/1023 * k) - 1.0) * 0.50};
+        const double u{(std::cos(al::numbers::pi*4.000/1023 * k) - 1.0) * 0.08};
+        return s * (t + u + 1.0) / k;
+    }
+
     constexpr CubicFilter()
     {
-        /* This creates a lookup table for a cubic spline filter, with 256
+        /* This creates a lookup table for a gaussian-like filter, with 256
          * steps between samples. Only half the coefficients are needed, since
          * Coeff2 is just Coeff1 in reverse and Coeff3 is just Coeff0 in
          * reverse.
          */
-        for(size_t i{0};i < sTableSteps;++i)
+        for(size_t i{0};i < sTableSteps/2;++i)
         {
-            const double mu{static_cast<double>(i) / double{sTableSteps}};
-            const double mu2{mu*mu}, mu3{mu2*mu};
-            const double a0{-0.5*mu3 +      mu2 + -0.5*mu};
-            const double a1{ 1.5*mu3 + -2.5*mu2           + 1.0f};
-            mFilter[i] = static_cast<float>(a1);
-            mFilter[sTableSteps+i] = static_cast<float>(a0);
+            const double coeff0{GetCoeff(sTableSteps + i)};
+            const double coeff1{GetCoeff(i)};
+            const double coeff2{GetCoeff(sTableSteps - i)};
+            const double coeff3{GetCoeff(sTableSteps*2 - i)};
+
+            const double scale{1.0 / (coeff0 + coeff1 + coeff2 + coeff3)};
+            mFilter[sTableSteps + i] = static_cast<float>(coeff0 * scale);
+            mFilter[i] = static_cast<float>(coeff1 * scale);
+            mFilter[sTableSteps - i] = static_cast<float>(coeff2 * scale);
+            mFilter[sTableSteps*2 - i] = static_cast<float>(coeff3 * scale);
         }
     }
 
@@ -95,7 +111,7 @@ struct CubicFilter {
     [[nodiscard]] constexpr auto getCoeff3(size_t i) const noexcept -> float
     { return mFilter[sTableSteps*2-i]; }
 };
-constexpr CubicFilter gCubicTable;
+const CubicFilter gCubicTable;
 
 
 /* Max samples per process iteration. Used to limit the size needed for

alsoftrc.sample

@@ -179,7 +179,7 @@
 #  Selects the default resampler used when mixing sources. Valid values are:
 #  point - nearest sample, no interpolation
 #  linear - extrapolates samples using a linear slope between samples
-#  cubic - extrapolates samples using a Catmull-Rom spline
+#  cubic - extrapolates samples using a 4-point Gaussian filter
 #  bsinc12 - extrapolates samples using a band-limited Sinc filter (varying
 #            between 12 and 24 points, with anti-aliasing)
 #  fast_bsinc12 - same as bsinc12, except without interpolation between down-

core/cubic_tables.cpp

@@ -2,50 +2,79 @@
 #include "cubic_tables.h"
 
 #include <array>
+#include <cmath>
 #include <cstddef>
 
+#include "alnumbers.h"
 #include "cubic_defs.h"
 
 
 namespace {
 
-struct SplineFilterArray {
+struct GaussFilterArray {
     alignas(16) std::array<CubicCoefficients,CubicPhaseCount> mTable{};
 
-    constexpr SplineFilterArray()
+    /* This filter table is inspired by the gaussian-like filter found in the
+     * SNES. This is based on the public domain code developed by Near, with
+     * the help of Ryphecha and nocash, from the nesdev.org forums.
+     *
+     * <https://forums.nesdev.org/viewtopic.php?p=251534#p251534>
+     *
+     * Additional changes were made here, the most obvious being that is has
+     * full floating-point precision instead of 11-bit fixed-point, but also an
+     * offset adjustment for the phase coefficients to more cleanly transition
+     * from the end of one sample set to the start of the next.
+     */
+    [[nodiscard]]
+    auto GetCoeff(std::size_t idx) noexcept -> double
+    {
+        static constexpr double IndexScale{512.0 / double{CubicPhaseCount*2}};
+        const double k{0.5 + static_cast<double>(idx)*IndexScale};
+        if(k > 512.0) return 0.0;
+        const double s{ std::sin(al::numbers::pi*1.280/1024 * k)};
+        const double t{(std::cos(al::numbers::pi*2.000/1023 * k) - 1.0) * 0.50};
+        const double u{(std::cos(al::numbers::pi*4.000/1023 * k) - 1.0) * 0.08};
+        return s * (t + u + 1.0) / k;
+    }
+
+    GaussFilterArray()
     {
         /* Fill in the main coefficients. */
-        for(size_t pi{0};pi < CubicPhaseCount;++pi)
+        for(std::size_t pi{0};pi < CubicPhaseCount;++pi)
         {
-            const double mu{static_cast<double>(pi) / CubicPhaseCount};
-            const double mu2{mu*mu}, mu3{mu2*mu};
-            mTable[pi].mCoeffs[0] = static_cast<float>(-0.5*mu3 +      mu2 + -0.5*mu);
-            mTable[pi].mCoeffs[1] = static_cast<float>( 1.5*mu3 + -2.5*mu2           + 1.0);
-            mTable[pi].mCoeffs[2] = static_cast<float>(-1.5*mu3 +  2.0*mu2 +  0.5*mu);
-            mTable[pi].mCoeffs[3] = static_cast<float>( 0.5*mu3 + -0.5*mu2);
+            const double coeff0{GetCoeff(CubicPhaseCount + pi)};
+            const double coeff1{GetCoeff(pi)};
+            const double coeff2{GetCoeff(CubicPhaseCount - pi)};
+            const double coeff3{GetCoeff(CubicPhaseCount*2 - pi)};
+
+            const double scale{1.0 / (coeff0 + coeff1 + coeff2 + coeff3)};
+            mTable[pi].mCoeffs[0] = static_cast<float>(coeff0 * scale);
+            mTable[pi].mCoeffs[1] = static_cast<float>(coeff1 * scale);
+            mTable[pi].mCoeffs[2] = static_cast<float>(coeff2 * scale);
+            mTable[pi].mCoeffs[3] = static_cast<float>(coeff3 * scale);
         }
 
         /* Fill in the coefficient deltas. */
-        for(size_t pi{0};pi < CubicPhaseCount-1;++pi)
+        for(std::size_t pi{0};pi < CubicPhaseCount-1;++pi)
         {
             mTable[pi].mDeltas[0] = mTable[pi+1].mCoeffs[0] - mTable[pi].mCoeffs[0];
             mTable[pi].mDeltas[1] = mTable[pi+1].mCoeffs[1] - mTable[pi].mCoeffs[1];
             mTable[pi].mDeltas[2] = mTable[pi+1].mCoeffs[2] - mTable[pi].mCoeffs[2];
             mTable[pi].mDeltas[3] = mTable[pi+1].mCoeffs[3] - mTable[pi].mCoeffs[3];
         }
 
-        const size_t pi{CubicPhaseCount - 1};
-        mTable[pi].mDeltas[0] = -mTable[pi].mCoeffs[0];
-        mTable[pi].mDeltas[1] = -mTable[pi].mCoeffs[1];
-        mTable[pi].mDeltas[2] = 1.0f - mTable[pi].mCoeffs[2];
-        mTable[pi].mDeltas[3] = -mTable[pi].mCoeffs[3];
+        const std::size_t pi{CubicPhaseCount - 1};
+        mTable[pi].mDeltas[0] =                 0.0f - mTable[pi].mCoeffs[0];
+        mTable[pi].mDeltas[1] = mTable[0].mCoeffs[2] - mTable[pi].mCoeffs[1];
+        mTable[pi].mDeltas[2] = mTable[0].mCoeffs[1] - mTable[pi].mCoeffs[2];
+        mTable[pi].mDeltas[3] = mTable[0].mCoeffs[0] - mTable[pi].mCoeffs[3];
     }
 
     [[nodiscard]] constexpr auto& getTable() const noexcept { return mTable; }
 };
 
-constexpr SplineFilterArray SplineFilter{};
+const GaussFilterArray GaussFilter{};
 
 } // namespace
 
-const CubicTable gCubicSpline{SplineFilter.getTable()};
+const CubicTable gGaussianFilter{GaussFilter.getTable()};

core/cubic_tables.h

@@ -9,9 +9,7 @@ struct CubicTable {
     al::span<const CubicCoefficients,CubicPhaseCount> Tab;
 };
 
-/* A Catmull-Rom spline. The spline passes through the center two samples,
- * ensuring no discontinuity while moving through a series of samples.
- */
-extern const CubicTable gCubicSpline;
+/* A Gaussian-like resample filter. */
+extern const CubicTable gGaussianFilter;
 
 #endif /* CORE_CUBIC_TABLES_H */