Add polynomial interpolation

src/sfizz/Config.h

@@ -59,6 +59,7 @@ namespace config {
     constexpr int maxCurves { 256 };
     constexpr int chunkSize { 1024 };
     constexpr int filtersInPool { maxVoices * 2 };
+    constexpr int excessFileFrames { 8 };
     /**
      * @brief The threshold for age stealing.
      *        In percentage of the voice's max age.

src/sfizz/FilePool.cpp

@@ -42,7 +42,8 @@ template <class T>
 void readBaseFile(SndfileHandle& sndFile, sfz::AudioBuffer<T>& output, uint32_t numFrames, bool reverse)
 {
     output.reset();
-    output.resize(numFrames);
+    output.resize(numFrames + sfz::config::excessFileFrames);
+    output.clear();
 
     if (reverse)
         sndFile.seek(-static_cast<sf_count_t>(numFrames), SEEK_END);
@@ -57,13 +58,13 @@ void readBaseFile(SndfileHandle& sndFile, sfz::AudioBuffer<T>& output, uint32_t
         output.addChannel();
         sfz::Buffer<T> tempReadBuffer { 2 * numFrames };
         sndFile.readf(tempReadBuffer.data(), numFrames);
-        sfz::readInterleaved<T>(tempReadBuffer, output.getSpan(0), output.getSpan(1));
+        sfz::readInterleaved<T>(tempReadBuffer, output.getSpan(0).first(numFrames), output.getSpan(1).first(numFrames));
     }
 
     if (reverse) {
         for (unsigned c = 0; c < channels; ++c) {
             // TODO: consider optimizing with SIMD
-            absl::Span<float> channel = output.getSpan(c);
+            absl::Span<float> channel = output.getSpan(c).first(numFrames);
             std::reverse(channel.begin(), channel.end());
         }
     }
@@ -97,7 +98,8 @@ void streamFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversamplin
     auto baseBuffer = readFromFile<T>(sndFile, numFrames, sfz::Oversampling::x1, reverse);
     output.reset();
     output.addChannels(baseBuffer->getNumChannels());
-    output.resize(numFrames * static_cast<int>(factor));
+    output.resize(numFrames * static_cast<int>(factor) + sfz::config::excessFileFrames);
+    output.clear();
     sfz::Oversampler oversampler { factor };
     oversampler.stream(*baseBuffer, output, filledFrames);
 }

src/sfizz/FilePool.h

@@ -65,7 +65,8 @@ struct FilePromise
     AudioSpan<const float> getData()
     {
         if (dataStatus == DataStatus::Ready)
-            return AudioSpan<const float>(fileData);
+            return AudioSpan<const float>(fileData)
+                .first(fileData.getNumFrames() - sfz::config::excessFileFrames);
         else if (availableFrames > preloadedData->getNumFrames())
             return AudioSpan<const float>(fileData).first(availableFrames);
         else

src/sfizz/MathHelpers.h

@@ -11,11 +11,15 @@
 #pragma once
 #include "Config.h"
 #include "Macros.h"
+#include "SIMDConfig.h"
 #include "absl/types/span.h"
 #include <algorithm>
 #include <cmath>
 #include <random>
 #include <cfenv>
+#if SFIZZ_HAVE_SSE
+#include <xmmintrin.h>
+#endif
 
 template <class T>
 constexpr T max(T op1, T op2)
@@ -156,9 +160,165 @@ inline CXX14_CONSTEXPR void incrementAll(T& first, Args&... rest)
 }
 
 template <class ValueType>
-constexpr ValueType linearInterpolation(ValueType left, ValueType right, ValueType coeff)
+constexpr ValueType linearInterpolation(const ValueType values[2], ValueType coeff)
 {
-    return left * (static_cast<ValueType>(1.0) - coeff) + right * coeff;
+    return values[0] * (static_cast<ValueType>(1.0) - coeff) + values[1] * coeff;
+}
+
+/**
+ * @brief Compute the 3rd-order Hermite interpolation polynomial.
+ *
+ * @tparam R
+ * @param x
+ * @return R
+ */
+template <class R>
+R hermite3(R x)
+{
+    x = std::abs(x);
+    R x2 = x * x;
+    R x3 = x2 * x;
+    R y = 0;
+    R q = R(5./2.) * x2; // a reoccurring term
+    R p1 = R(1) - q + R(3./2.) * x3;
+    R p2 = R(2) - R(4) * x + q - R(1./2.) * x3;
+    y = (x < R(2)) ? p2 : y;
+    y = (x < R(1)) ? p1 : y;
+    return y;
+}
+
+#if SFIZZ_HAVE_SSE
+/**
+ * @brief Compute 4 parallel elements of the 3rd-order Hermite interpolation polynomial.
+ *
+ * @param x
+ * @return __m128
+ */
+inline __m128 hermite3x4(__m128 x)
+{
+    x = _mm_andnot_ps(_mm_set1_ps(-0.0f), x);
+    __m128 x2 = _mm_mul_ps(x, x);
+    __m128 x3 = _mm_mul_ps(x2, x);
+    __m128 y = _mm_set1_ps(0.0f);
+    __m128 q = _mm_mul_ps(_mm_set1_ps(5./2.), x2);
+    __m128 p1 = _mm_add_ps(_mm_sub_ps(_mm_set1_ps(1), q), _mm_mul_ps(_mm_set1_ps(3./2.), x3));
+    __m128 p2 = _mm_sub_ps(_mm_add_ps(_mm_sub_ps(_mm_set1_ps(2), _mm_mul_ps(_mm_set1_ps(4), x)), q), _mm_mul_ps(_mm_set1_ps(1./2.), x3));
+    __m128 m2 = _mm_cmple_ps(x, _mm_set1_ps(2));
+    y = _mm_or_ps(_mm_and_ps(m2, p2), _mm_andnot_ps(m2, y));
+    __m128 m1 = _mm_cmple_ps(x, _mm_set1_ps(1));
+    y = _mm_or_ps(_mm_and_ps(m1, p1), _mm_andnot_ps(m1, y));
+    return y;
+}
+#endif
+
+template <class ValueType>
+ValueType hermite3Interpolation(const ValueType values[4], ValueType coeff);
+
+#if SFIZZ_HAVE_SSE
+template <>
+inline float hermite3Interpolation<float>(const float values[4], float coeff)
+{
+    __m128 x = _mm_sub_ps(_mm_setr_ps(-1, 0, 1, 2), _mm_set1_ps(coeff));
+    __m128 h = hermite3x4(x);
+    __m128 y = _mm_mul_ps(h, _mm_loadu_ps(values));
+    // 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
+
+template <class ValueType>
+ValueType hermite3Interpolation(const ValueType values[4], ValueType coeff)
+{
+    ValueType y = 0;
+    for (int i = 0; i < 4; ++i) {
+        ValueType h = hermite3<ValueType>(i - 1 - coeff);
+        y += h * values[i];
+    }
+    return y;
+}
+
+/**
+ * @brief Compute the 3rd-order B-spline interpolation polynomial.
+ *
+ * @tparam R
+ * @param x
+ * @return R
+ */
+template <class R>
+R bspline3(R x)
+{
+    x = std::abs(x);
+    R x2 = x * x;
+    R x3 = x2 * x;
+    R y = 0;
+    R p1 = R(2./3.) - x2 + R(1./2.) * x3;
+    R p2 = R(4./3.) - R(2) * x + x2 - R(1./6.) * x3;
+    y = (x < R(2)) ? p2 : y;
+    y = (x < R(1)) ? p1 : y;
+    return y;
+}
+
+#if SFIZZ_HAVE_SSE
+/**
+ * @brief Compute 4 parallel elements of the 3rd-order B-spline interpolation polynomial.
+ *
+ * @param x
+ * @return __m128
+ */
+inline __m128 bspline3x4(__m128 x)
+{
+    x = _mm_andnot_ps(_mm_set1_ps(-0.0f), x);
+    __m128 x2 = _mm_mul_ps(x, x);
+    __m128 x3 = _mm_mul_ps(x2, x);
+    __m128 y = _mm_set1_ps(0.0f);
+    __m128 p1 = _mm_set1_ps(2./3.) - x2 + _mm_mul_ps(_mm_set1_ps(1./2.), x3);
+    __m128 p2 = _mm_sub_ps(_mm_add_ps(_mm_sub_ps(_mm_set1_ps(4./3.), _mm_mul_ps(_mm_set1_ps(2), x)), x2), _mm_mul_ps(_mm_set1_ps(1./6.), x3));
+    __m128 m2 = _mm_cmple_ps(x, _mm_set1_ps(2));
+    y = _mm_or_ps(_mm_and_ps(m2, p2), _mm_andnot_ps(m2, y));
+    __m128 m1 = _mm_cmple_ps(x, _mm_set1_ps(1));
+    y = _mm_or_ps(_mm_and_ps(m1, p1), _mm_andnot_ps(m1, y));
+    return y;
+}
+#endif
+
+template <class ValueType>
+ValueType bspline3Interpolation(const ValueType values[4], ValueType coeff);
+
+#if SFIZZ_HAVE_SSE
+template <>
+inline float bspline3Interpolation<float>(const float values[4], float coeff)
+{
+    __m128 x = _mm_sub_ps(_mm_setr_ps(-1, 0, 1, 2), _mm_set1_ps(coeff));
+    __m128 h = bspline3x4(x);
+    __m128 y = _mm_mul_ps(h, _mm_loadu_ps(values));
+    // 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
+
+template <class ValueType>
+ValueType bspline3Interpolation(const ValueType values[4], ValueType coeff)
+{
+    ValueType y = 0;
+    for (int i = 0; i < 4; ++i) {
+        ValueType h = bspline3<ValueType>(i - 1 - coeff);
+        y += h * values[i];
+    }
+    return y;
 }
 
 template <class Type>

src/sfizz/Oversampler.cpp

@@ -89,7 +89,7 @@ void sfz::Oversampler::stream(const sfz::AudioBuffer<float>& input, sfz::AudioBu
         break;
     case Oversampling::x1:
         for (size_t i = 0; i < numChannels; ++i)
-            copy<float>(input.getConstSpan(i), output.getSpan(i));
+            copy<float>(input.getConstSpan(i), output.getSpan(i).first(numFrames));
         return;
     }
 

src/sfizz/Voice.cpp

@@ -485,7 +485,7 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
         const auto sampleEnd = min(
             static_cast<int>(region->trueSampleEnd(currentPromise->oversamplingFactor)),
             static_cast<int>(source.getNumFrames())
-        ) - 2;
+        ) - 1;
         for (unsigned i = 0; i < indices->size(); ++i) {
             if ((*indices)[i] >= sampleEnd) {
 #ifndef NDEBUG
@@ -511,15 +511,15 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
     auto left = buffer.getChannel(0);
     if (source.getNumChannels() == 1) {
         while (ind < indices->end()) {
-            *left = linearInterpolation(leftSource[*ind], leftSource[*ind + 1], *coeff);
+            *left = bspline3Interpolation(&leftSource[*ind], *coeff);
             incrementAll(ind, left, coeff);
         }
     } else {
         auto right = buffer.getChannel(1);
         auto rightSource = source.getConstSpan(1);
         while (ind < indices->end()) {
-            *left = linearInterpolation(leftSource[*ind], leftSource[*ind + 1], *coeff);
+            *left = bspline3Interpolation(&leftSource[*ind], *coeff);
-            *right = linearInterpolation(rightSource[*ind], rightSource[*ind + 1], *coeff);
+            *right = bspline3Interpolation(&rightSource[*ind], *coeff);
             incrementAll(ind, left, right, coeff);
         }
     }