diff --git a/src/sfizz/Interpolators.h b/src/sfizz/Interpolators.h
new file mode 100644
index 00000000..834113bf
--- /dev/null
+++ b/src/sfizz/Interpolators.h
@@ -0,0 +1,25 @@
+// 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
+
+namespace sfz {
+
+enum InterpolatorModel : int {
+    // a linear interpolator
+    kInterpolatorLinear,
+    // a Hermite 3rd order interpolator
+    kInterpolatorHermite3,
+    // a B-spline 3rd order interpolator
+    kInterpolatorBspline3,
+};
+
+template <InterpolatorModel M, class R>
+R interpolate(const R* values, R coeff);
+
+} // namespace sfz
+
+#include "Interpolators.hpp"
diff --git a/src/sfizz/Interpolators.hpp b/src/sfizz/Interpolators.hpp
new file mode 100644
index 00000000..ef7e0324
--- /dev/null
+++ b/src/sfizz/Interpolators.hpp
@@ -0,0 +1,123 @@
+// 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"
+#include "MathHelpers.h"
+#include "SIMDConfig.h"
+
+namespace sfz {
+
+template <InterpolatorModel M, class R>
+class Interpolator;
+
+template <InterpolatorModel M, class R>
+inline R interpolate(const R* values, R coeff)
+{
+    return Interpolator<M, R>::process(values, coeff);
+}
+
+//------------------------------------------------------------------------------
+// Linear
+
+template <class R>
+class Interpolator<kInterpolatorLinear, R>
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        return values[0] * (static_cast<R>(1.0) - coeff) + values[1] * coeff;
+    }
+};
+
+//------------------------------------------------------------------------------
+// Hermite 3rd order, SSE specialization
+
+#if SFIZZ_HAVE_SSE
+template <>
+class Interpolator<kInterpolatorHermite3, float>
+{
+public:
+    static inline float process(const float* values, 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
+
+//------------------------------------------------------------------------------
+// Hermite 3rd order, generic
+
+template <class R>
+class Interpolator<kInterpolatorHermite3, R>
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        R y = 0;
+        for (int i = 0; i < 4; ++i) {
+            R h = hermite3<R>(i - 1 - coeff);
+            y += h * values[i];
+        }
+        return y;
+    }
+};
+
+//------------------------------------------------------------------------------
+// B-spline 3rd order, SSE specialization
+
+#if SFIZZ_HAVE_SSE
+template <>
+class Interpolator<kInterpolatorBspline3, float>
+{
+public:
+    static inline float process(const float* values, 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
+
+//------------------------------------------------------------------------------
+// B-spline 3rd order, generic
+
+template <class R>
+class Interpolator<kInterpolatorBspline3, R>
+{
+public:
+    static inline R process(const R* values, R coeff)
+    {
+        R y = 0;
+        for (int i = 0; i < 4; ++i) {
+            R h = bspline3<R>(i - 1 - coeff);
+            y += h * values[i];
+        }
+        return y;
+    }
+};
+
+} // namespace sfz
diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h
index ed7679cb..f20962d7 100644
--- a/src/sfizz/MathHelpers.h
+++ b/src/sfizz/MathHelpers.h
@@ -159,12 +159,6 @@ inline CXX14_CONSTEXPR void incrementAll(T& first, Args&... rest)
     incrementAll<Increment>(rest...);
 }
 
-template <class ValueType>
-constexpr ValueType linearInterpolation(const ValueType values[2], ValueType coeff)
-{
-    return values[0] * (static_cast<ValueType>(1.0) - coeff) + values[1] * coeff;
-}
-
 /**
  * @brief Compute the 3rd-order Hermite interpolation polynomial.
  *
@@ -211,39 +205,6 @@ inline __m128 hermite3x4(__m128 x)
 }
 #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.
  *
@@ -288,39 +249,6 @@ inline __m128 bspline3x4(__m128 x)
 }
 #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>
 constexpr Type pi() { return static_cast<Type>(3.141592653589793238462643383279502884); };
 template <class Type>
diff --git a/src/sfizz/Voice.cpp b/src/sfizz/Voice.cpp
index 4e034a36..d198b22f 100644
--- a/src/sfizz/Voice.cpp
+++ b/src/sfizz/Voice.cpp
@@ -11,6 +11,7 @@
 #include "MathHelpers.h"
 #include "SIMDHelpers.h"
 #include "SfzHelpers.h"
+#include "Interpolators.h"
 #include "absl/algorithm/container.h"
 
 sfz::Voice::Voice(sfz::Resources& resources)
@@ -511,15 +512,15 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
     auto left = buffer.getChannel(0);
     if (source.getNumChannels() == 1) {
         while (ind < indices->end()) {
-            *left = bspline3Interpolation(&leftSource[*ind], *coeff);
+            *left = interpolate<kInterpolatorBspline3>(&leftSource[*ind], *coeff);
             incrementAll(ind, left, coeff);
         }
     } else {
         auto right = buffer.getChannel(1);
         auto rightSource = source.getConstSpan(1);
         while (ind < indices->end()) {
-            *left = bspline3Interpolation(&leftSource[*ind], *coeff);
-            *right = bspline3Interpolation(&rightSource[*ind], *coeff);
+            *left = interpolate<kInterpolatorBspline3>(&leftSource[*ind], *coeff);
+            *right = interpolate<kInterpolatorBspline3>(&rightSource[*ind], *coeff);
             incrementAll(ind, left, right, coeff);
         }
     }