diff --git a/benchmarks/BM_interpolationCast.cpp b/benchmarks/BM_interpolationCast.cpp
index ea12845a..b3743530 100644
--- a/benchmarks/BM_interpolationCast.cpp
+++ b/benchmarks/BM_interpolationCast.cpp
@@ -22,8 +22,7 @@ public:
     std::mt19937 gen { rd() };
     std::uniform_real_distribution<float> dist { 0, maxJump };
     jumps = std::vector<int>(state.range(0));
-    leftCoeffs = std::vector<float>(state.range(0));
-    rightCoeffs = std::vector<float>(state.range(0));
+    coeffs = std::vector<float>(state.range(0));
     floatJumps = std::vector<float>(state.range(0));
     absl::c_generate(floatJumps, [&]() { return dist(gen); });
   }
@@ -33,8 +32,7 @@ public:
   }
 
     std::vector<int> jumps;
-    std::vector<float> leftCoeffs;
-    std::vector<float> rightCoeffs;
+    std::vector<float> coeffs;
     std::vector<float> floatJumps;
 };
 
@@ -42,28 +40,28 @@ public:
 BENCHMARK_DEFINE_F(InterpolationCast, Scalar)(benchmark::State& state) {
     for (auto _ : state)
     {
-        sfz::sfzInterpolationCast<float, false>(floatJumps, absl::MakeSpan(jumps), absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs));
+        sfz::sfzInterpolationCast<float, false>(floatJumps, absl::MakeSpan(jumps), absl::MakeSpan(coeffs));
     }
 }
 
 BENCHMARK_DEFINE_F(InterpolationCast, SIMD)(benchmark::State& state) {
     for (auto _ : state)
     {
-        sfz::sfzInterpolationCast<float, true>(floatJumps, absl::MakeSpan(jumps), absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs));
+        sfz::sfzInterpolationCast<float, true>(floatJumps, absl::MakeSpan(jumps), absl::MakeSpan(coeffs));
     }
 }
 
 BENCHMARK_DEFINE_F(InterpolationCast, Scalar_Unaligned)(benchmark::State& state) {
     for (auto _ : state)
     {
-        sfz::sfzInterpolationCast<float, false>(absl::MakeSpan(floatJumps).subspan(1), absl::MakeSpan(jumps).subspan(3), absl::MakeSpan(leftCoeffs).subspan(2), absl::MakeSpan(rightCoeffs).subspan(1));
+        sfz::sfzInterpolationCast<float, false>(absl::MakeSpan(floatJumps).subspan(1), absl::MakeSpan(jumps).subspan(3), absl::MakeSpan(coeffs).subspan(1));
     }
 }
 
 BENCHMARK_DEFINE_F(InterpolationCast, SIMD_Unaligned)(benchmark::State& state) {
     for (auto _ : state)
     {
-        sfz::sfzInterpolationCast<float, true>(absl::MakeSpan(floatJumps).subspan(1), absl::MakeSpan(jumps).subspan(3), absl::MakeSpan(leftCoeffs).subspan(2), absl::MakeSpan(rightCoeffs).subspan(1));
+        sfz::sfzInterpolationCast<float, true>(absl::MakeSpan(floatJumps).subspan(1), absl::MakeSpan(jumps).subspan(3), absl::MakeSpan(coeffs).subspan(1));
     }
 }
 
diff --git a/benchmarks/BM_interpolators.cpp b/benchmarks/BM_interpolators.cpp
new file mode 100644
index 00000000..77337576
--- /dev/null
+++ b/benchmarks/BM_interpolators.cpp
@@ -0,0 +1,87 @@
+// 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 "Config.h"
+#include "Interpolators.h"
+#include "ScopedFTZ.h"
+#include "absl/types/span.h"
+#include <benchmark/benchmark.h>
+#include <random>
+
+class Interpolators : public benchmark::Fixture {
+public:
+    void SetUp(const ::benchmark::State& state)
+    {
+        std::random_device rd { };
+        std::mt19937 gen { rd() };
+        std::uniform_real_distribution<float> dist { -1.0f, 1.0f };
+
+        const size_t numFramesIn = state.range(0);
+        inputBuffer = std::vector<float>(numFramesIn + 2 * sfz::config::excessFileFrames);
+
+        // any ratio will do, compute time will be proportional
+        static constexpr float ratio = 1.234;
+
+        const size_t numFramesOut = static_cast<size_t>(std::ceil(numFramesIn * ratio));
+        input = absl::MakeSpan(inputBuffer).subspan(sfz::config::excessFileFrames, numFramesIn);
+        output = std::vector<float>(numFramesOut);
+        std::generate(input.begin(), input.end(), [&]() { return dist(gen); });
+    }
+
+    void TearDown(const ::benchmark::State& /* state */)
+    {
+    }
+
+    std::vector<float> inputBuffer;
+    absl::Span<float> input;
+    std::vector<float> output;
+
+    enum { excessFrames = 8 };
+};
+
+template <sfz::InterpolatorModel M>
+static void doInterpolation(absl::Span<const float> input, absl::Span<float> output)
+{
+    const float kOutToIn = static_cast<float>(input.size()) / output.size();
+
+    for (size_t iOut = 0; iOut < output.size(); ++iOut) {
+        float posIn = iOut * kOutToIn;
+        unsigned dec = static_cast<int>(posIn);
+        float frac = posIn - dec;
+        output[iOut] = sfz::interpolate<M>(&input[dec], frac);
+    }
+}
+
+BENCHMARK_DEFINE_F(Interpolators, Linear)(benchmark::State& state)
+{
+    ScopedFTZ ftz;
+
+    for (auto _ : state) {
+        doInterpolation<sfz::kInterpolatorLinear>(input, absl::MakeSpan(output));
+    }
+}
+
+BENCHMARK_DEFINE_F(Interpolators, Hermite3)(benchmark::State& state)
+{
+    ScopedFTZ ftz;
+
+    for (auto _ : state) {
+        doInterpolation<sfz::kInterpolatorHermite3>(input, absl::MakeSpan(output));
+    }
+}
+
+BENCHMARK_DEFINE_F(Interpolators, Bspline3)(benchmark::State& state)
+{
+    ScopedFTZ ftz;
+
+    for (auto _ : state) {
+        doInterpolation<sfz::kInterpolatorBspline3>(input, absl::MakeSpan(output));
+    }
+}
+
+BENCHMARK_REGISTER_F(Interpolators, Linear)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
+BENCHMARK_REGISTER_F(Interpolators, Hermite3)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
+BENCHMARK_REGISTER_F(Interpolators, Bspline3)->RangeMultiplier(4)->Range(1 << 4, 1 << 12);
diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt
index 96595e47..1b2fcdbe 100644
--- a/benchmarks/CMakeLists.txt
+++ b/benchmarks/CMakeLists.txt
@@ -89,6 +89,8 @@ target_link_libraries(bm_readChunkFlac PRIVATE sfizz-sndfile)
 sfizz_add_benchmark(bm_resampleChunk BM_resampleChunk.cpp)
 target_link_libraries(bm_resampleChunk PRIVATE sfizz-sndfile)
 
+sfizz_add_benchmark(bm_interpolators BM_interpolators.cpp)
+
 sfizz_add_benchmark(bm_filterModulation BM_filterModulation.cpp ../src/sfizz/SfzFilter.cpp)
 target_link_libraries(bm_filterModulation PRIVATE sfizz-sndfile)
 
diff --git a/src/sfizz/AudioBuffer.h b/src/sfizz/AudioBuffer.h
index 51500d0a..0c47ba0a 100644
--- a/src/sfizz/AudioBuffer.h
+++ b/src/sfizz/AudioBuffer.h
@@ -26,7 +26,9 @@ namespace sfz
  * @tparam MaxChannels the maximum number of channels in the buffer
  * @tparam Alignment the alignment for the buffers
  */
-template <class Type, size_t MaxChannels = sfz::config::numChannels, unsigned int Alignment = SIMDConfig::defaultAlignment>
+template <class Type, size_t MaxChannels = sfz::config::numChannels,
+          unsigned int Alignment = SIMDConfig::defaultAlignment,
+          size_t PaddingLeft_ = 0, size_t PaddingRight = 0>
 class AudioBuffer {
 public:
     using value_type = typename std::remove_cv<Type>::type;
@@ -35,6 +37,13 @@ public:
     using iterator = pointer;
     using const_iterator = const_pointer;
 
+    enum {
+        //! Increased left padding to preserve required alignment
+        PaddingLeft = PaddingLeft_ + (Alignment - (PaddingLeft_ % Alignment)) % Alignment,
+        //! Total padding left and right
+        PaddingTotal = PaddingLeft + PaddingRight,
+    };
+
     /**
      * @brief Construct a new Audio Buffer object
      *
@@ -55,7 +64,7 @@ public:
         , numFrames(numFrames)
     {
         for (size_t i = 0; i < numChannels; ++i)
-            buffers[i] = absl::make_unique<buffer_type>(numFrames);
+            buffers[i] = absl::make_unique<buffer_type>(numFrames + PaddingTotal);
     }
 
     /**
@@ -70,7 +79,7 @@ public:
         bool returnedOK = true;
 
         for (size_t i = 0; i < numChannels; ++i)
-            returnedOK &= buffers[i]->resize(newSize, std::nothrow);
+            returnedOK &= buffers[i]->resize(newSize + PaddingTotal, std::nothrow);
 
         if (returnedOK)
             numFrames = newSize;
@@ -86,7 +95,7 @@ public:
     void resize(size_t newSize)
     {
         for (size_t i = 0; i < numChannels; ++i)
-            buffers[i]->resize(newSize);
+            buffers[i]->resize(newSize + PaddingTotal);
 
         numFrames = newSize;
     }
@@ -101,7 +110,7 @@ public:
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
-            return buffers[channelIndex]->data();
+            return buffers[channelIndex]->data() + PaddingLeft;
 
         return {};
     }
@@ -116,7 +125,7 @@ public:
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
-            return buffers[channelIndex]->end();
+            return buffers[channelIndex]->end() - PaddingRight;
 
         return {};
     }
@@ -131,7 +140,7 @@ public:
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
-            return buffers[channelIndex]->data();
+            return buffers[channelIndex]->data() + PaddingLeft;
 
         return {};
     }
@@ -146,7 +155,7 @@ public:
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
-            return buffers[channelIndex]->end();
+            return buffers[channelIndex]->end() - PaddingRight;
 
         return {};
     }
@@ -161,7 +170,7 @@ public:
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
-            return { buffers[channelIndex]->data(), buffers[channelIndex]->size() };
+            return { buffers[channelIndex]->data() + PaddingLeft, numFrames };
 
         return {};
     }
@@ -184,7 +193,7 @@ public:
     void addChannel()
     {
         if (numChannels < MaxChannels)
-            buffers[numChannels++] = absl::make_unique<buffer_type>(numFrames);
+            buffers[numChannels++] = absl::make_unique<buffer_type>(numFrames + PaddingTotal);
     }
 
     /**
@@ -233,7 +242,7 @@ public:
         ASSERT(buffers[channelIndex] != nullptr);
         ASSERT(frameIndex < numFrames);
 
-        return *(buffers[channelIndex]->data() + frameIndex);
+        return *(buffers[channelIndex]->data() + PaddingLeft + frameIndex);
     }
 
     /**
@@ -265,8 +274,10 @@ public:
      */
     void clear()
     {
-        for (size_t i = 0; i < numChannels; ++i)
-            fill<Type>(getSpan(i), Type{ 0.0 });
+        for (size_t i = 0; i < numChannels; ++i) {
+            absl::Span<Type> paddedSpan { buffers[i]->data(), numFrames + PaddingTotal };
+            fill<Type>(paddedSpan, Type{ 0.0 });
+        }
     }
 
     /**
@@ -281,22 +292,6 @@ public:
             addChannel();
     }
 
-    /**
-     * @brief Convert implicitly to a pointer of channels
-     */
-    operator const float* const*() const noexcept
-    {
-        return buffers.data();
-    }
-
-    /**
-     * @brief Convert implicitly to a pointer of channels
-     */
-    operator float* const*() noexcept
-    {
-        return buffers.data();
-    }
-
 private:
     using buffer_type = Buffer<Type, Alignment>;
     using buffer_ptr = std::unique_ptr<buffer_type>;
diff --git a/src/sfizz/AudioSpan.h b/src/sfizz/AudioSpan.h
index e11e2ebb..23c9d340 100644
--- a/src/sfizz/AudioSpan.h
+++ b/src/sfizz/AudioSpan.h
@@ -148,8 +148,8 @@ public:
      * @tparam Alignment the alignment block size for the platform
      * @param audioBuffer the source AudioBuffer.
      */
-    template <class U, size_t N, unsigned int Alignment, typename = typename std::enable_if<N <= MaxChannels>::type, typename = typename std::enable_if<std::is_const<U>::value, int>::type>
-    AudioSpan(AudioBuffer<U, N, Alignment>& audioBuffer)
+    template <class U, size_t N, unsigned int Alignment, size_t PaddingLeft, size_t PaddingRight, typename = typename std::enable_if<N <= MaxChannels>::type, typename = typename std::enable_if<std::is_const<U>::value, int>::type>
+    AudioSpan(AudioBuffer<U, N, Alignment, PaddingLeft, PaddingRight>& audioBuffer)
         : numFrames(audioBuffer.getNumFrames())
         , numChannels(audioBuffer.getNumChannels())
     {
@@ -169,8 +169,8 @@ public:
      * @tparam Alignment the alignment block size for the platform
      * @param audioBuffer the source AudioBuffer.
      */
-    template <class U, size_t N, unsigned int Alignment, typename = std::enable_if<N <= MaxChannels>>
-    AudioSpan(AudioBuffer<U, N, Alignment>& audioBuffer)
+    template <class U, size_t N, unsigned int Alignment, size_t PaddingLeft, size_t PaddingRight, typename = std::enable_if<N <= MaxChannels>>
+    AudioSpan(AudioBuffer<U, N, Alignment, PaddingLeft, PaddingRight>& audioBuffer)
         : numFrames(audioBuffer.getNumFrames())
         , numChannels(audioBuffer.getNumChannels())
     {
@@ -246,7 +246,7 @@ public:
      * @param channelIndex the channel
      * @return absl::Span<const Type>
      */
-    absl::Span<const Type> getConstSpan(size_t channelIndex)
+    absl::Span<const Type> getConstSpan(size_t channelIndex) const
     {
         ASSERT(channelIndex < numChannels);
         if (channelIndex < numChannels)
@@ -382,7 +382,7 @@ public:
      *
      * @returns size_type the number of frames in the AudioSpan
      */
-    size_type getNumFrames()
+    size_type getNumFrames() const
     {
         return numFrames;
     }
@@ -392,7 +392,7 @@ public:
      *
      * @returns size_type the number of channels in the AudioSpan
      */
-    size_t getNumChannels()
+    size_t getNumChannels() const
     {
         return numChannels;
     }
diff --git a/src/sfizz/Config.h b/src/sfizz/Config.h
index b8f9d4ff..1c135740 100644
--- a/src/sfizz/Config.h
+++ b/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.
diff --git a/src/sfizz/Defaults.h b/src/sfizz/Defaults.h
index 6796380f..4ebe7577 100644
--- a/src/sfizz/Defaults.h
+++ b/src/sfizz/Defaults.h
@@ -221,6 +221,9 @@ namespace Default
 	constexpr Range<float> egOnCCPercentRange { -100.0, 100.0 };
 
     // ***** SFZ v2 ********
+    constexpr int sampleQuality { 2 };
+    constexpr Range<int> sampleQualityRange { 1, 10 }; // sample_quality
+
 	constexpr bool checkSustain { true }; // sustain_sw
 	constexpr bool checkSostenuto { true }; // sostenuto_sw
     constexpr Range<int> octaveOffsetRange { -10, 10 }; // octave_offset
diff --git a/src/sfizz/FilePool.cpp b/src/sfizz/FilePool.cpp
index 5e61b2b7..d0ac449e 100644
--- a/src/sfizz/FilePool.cpp
+++ b/src/sfizz/FilePool.cpp
@@ -27,6 +27,7 @@
 #include "FileInstrument.h"
 #include "Buffer.h"
 #include "AudioBuffer.h"
+#include "AudioSpan.h"
 #include "Config.h"
 #include "Debug.h"
 #include "Oversampler.h"
@@ -38,8 +39,7 @@
 #include <sndfile.hh>
 #include <thread>
 
-template <class T>
-void readBaseFile(SndfileHandle& sndFile, sfz::AudioBuffer<T>& output, uint32_t numFrames, bool reverse)
+void readBaseFile(SndfileHandle& sndFile, sfz::FileAudioBuffer& output, uint32_t numFrames, bool reverse)
 {
     output.reset();
     output.resize(numFrames);
@@ -51,13 +51,15 @@ void readBaseFile(SndfileHandle& sndFile, sfz::AudioBuffer<T>& output, uint32_t
 
     if (channels == 1) {
         output.addChannel();
+        output.clear();
         sndFile.readf(output.channelWriter(0), numFrames);
     } else if (channels == 2) {
         output.addChannel();
         output.addChannel();
-        sfz::Buffer<T> tempReadBuffer { 2 * numFrames };
+        output.clear();
+        sfz::Buffer<float> tempReadBuffer { 2 * numFrames };
         sndFile.readf(tempReadBuffer.data(), numFrames);
-        sfz::readInterleaved<T>(tempReadBuffer, output.getSpan(0), output.getSpan(1));
+        sfz::readInterleaved<float>(tempReadBuffer, output.getSpan(0), output.getSpan(1));
     }
 
     if (reverse) {
@@ -69,23 +71,23 @@ void readBaseFile(SndfileHandle& sndFile, sfz::AudioBuffer<T>& output, uint32_t
     }
 }
 
-template <class T>
-std::unique_ptr<sfz::AudioBuffer<T>> readFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversampling factor, bool reverse)
+std::unique_ptr<sfz::FileAudioBuffer> readFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversampling factor, bool reverse)
 {
-    auto baseBuffer = absl::make_unique<sfz::AudioBuffer<T>>();
+    auto baseBuffer = absl::make_unique<sfz::FileAudioBuffer>();
     readBaseFile(sndFile, *baseBuffer, numFrames, reverse);
 
     if (factor == sfz::Oversampling::x1)
         return baseBuffer;
 
-    auto outputBuffer = absl::make_unique<sfz::AudioBuffer<T>>(sndFile.channels(), numFrames * static_cast<int>(factor));
+    auto outputBuffer = absl::make_unique<sfz::FileAudioBuffer>(sndFile.channels(), numFrames * static_cast<int>(factor));
+    outputBuffer->clear();
     sfz::Oversampler oversampler { factor };
     oversampler.stream(*baseBuffer, *outputBuffer);
     return outputBuffer;
 }
 
 template <class T>
-void streamFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversampling factor, bool reverse, sfz::AudioBuffer<float>& output, std::atomic<size_t>* filledFrames = nullptr)
+void streamFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversampling factor, bool reverse, sfz::FileAudioBuffer& output, std::atomic<size_t>* filledFrames = nullptr)
 {
     if (factor == sfz::Oversampling::x1) {
         readBaseFile(sndFile, output, numFrames, reverse);
@@ -94,10 +96,11 @@ void streamFromFile(SndfileHandle& sndFile, uint32_t numFrames, sfz::Oversamplin
         return;
     }
 
-    auto baseBuffer = readFromFile<T>(sndFile, numFrames, sfz::Oversampling::x1, reverse);
+    auto baseBuffer = readFromFile(sndFile, numFrames, sfz::Oversampling::x1, reverse);
     output.reset();
     output.addChannels(baseBuffer->getNumChannels());
     output.resize(numFrames * static_cast<int>(factor));
+    output.clear();
     sfz::Oversampler oversampler { factor };
     oversampler.stream(*baseBuffer, output, filledFrames);
 }
@@ -260,12 +263,12 @@ bool sfz::FilePool::preloadFile(const FileId& fileId, uint32_t maxOffset) noexce
     const auto existingFile = preloadedFiles.find(fileId);
     if (existingFile != preloadedFiles.end()) {
         if (framesToLoad > existingFile->second.preloadedData->getNumFrames()) {
-            preloadedFiles[fileId].preloadedData = readFromFile<float>(sndFile, framesToLoad, oversamplingFactor, fileId.isReverse());
+            preloadedFiles[fileId].preloadedData = readFromFile(sndFile, framesToLoad, oversamplingFactor, fileId.isReverse());
         }
     } else {
         fileInformation->sampleRate = static_cast<float>(oversamplingFactor) * static_cast<float>(sndFile.samplerate());
         FileDataHandle handle {
-            readFromFile<float>(sndFile, framesToLoad, oversamplingFactor, fileId.isReverse()),
+            readFromFile(sndFile, framesToLoad, oversamplingFactor, fileId.isReverse()),
             *fileInformation
         };
         preloadedFiles.insert_or_assign(fileId, handle);
@@ -290,7 +293,7 @@ absl::optional<sfz::FileDataHandle> sfz::FilePool::loadFile(const FileId& fileId
     } else {
         fileInformation->sampleRate = static_cast<float>(oversamplingFactor) * static_cast<float>(sndFile.samplerate());
         FileDataHandle handle {
-            readFromFile<float>(sndFile, frames, oversamplingFactor, fileId.isReverse()),
+            readFromFile(sndFile, frames, oversamplingFactor, fileId.isReverse()),
             *fileInformation
         };
         loadedFiles.insert_or_assign(fileId, handle);
@@ -340,7 +343,7 @@ void sfz::FilePool::setPreloadSize(uint32_t preloadSize) noexcept
         const auto maxOffset = numFrames > this->preloadSize ? static_cast<uint32_t>(numFrames) - this->preloadSize : 0;
         fs::path file { rootDirectory / preloadedFile.first.filename() };
         SndfileHandle sndFile(file.string().c_str());
-        preloadedFile.second.preloadedData = readFromFile<float>(sndFile, preloadSize + maxOffset, oversamplingFactor, preloadedFile.first.isReverse());
+        preloadedFile.second.preloadedData = readFromFile(sndFile, preloadSize + maxOffset, oversamplingFactor, preloadedFile.first.isReverse());
     }
     this->preloadSize = preloadSize;
 }
@@ -469,7 +472,7 @@ void sfz::FilePool::setOversamplingFactor(sfz::Oversampling factor) noexcept
         const uint32_t maxOffset = numFrames > this->preloadSize ? static_cast<uint32_t>(numFrames) - this->preloadSize : 0;
         fs::path file { rootDirectory / preloadedFile.first.filename() };
         SndfileHandle sndFile(file.string().c_str());
-        preloadedFile.second.preloadedData = readFromFile<float>(sndFile, preloadSize + maxOffset, factor, preloadedFile.first.isReverse());
+        preloadedFile.second.preloadedData = readFromFile(sndFile, preloadSize + maxOffset, factor, preloadedFile.first.isReverse());
         preloadedFile.second.information.sampleRate *= samplerateChange;
     }
 
diff --git a/src/sfizz/FilePool.h b/src/sfizz/FilePool.h
index 7f5e15d9..70973a87 100644
--- a/src/sfizz/FilePool.h
+++ b/src/sfizz/FilePool.h
@@ -43,7 +43,9 @@
 #include <mutex>
 
 namespace sfz {
-using AudioBufferPtr = std::shared_ptr<AudioBuffer<float>>;
+using FileAudioBuffer = AudioBuffer<float, 2, SIMDConfig::defaultAlignment,
+                                    sfz::config::excessFileFrames, sfz::config::excessFileFrames>;
+using FileAudioBufferPtr = std::shared_ptr<FileAudioBuffer>;
 
 struct FileInformation {
     uint32_t end { Default::sampleEndRange.getEnd() };
@@ -56,7 +58,7 @@ struct FileInformation {
 // Strict C++11 disallows member initialization if aggregate initialization is to be used...
 struct FileDataHandle
 {
-    std::shared_ptr<AudioBuffer<float>> preloadedData;
+    FileAudioBufferPtr preloadedData;
     FileInformation information;
 };
 
@@ -96,8 +98,8 @@ struct FilePromise
     };
 
     FileId fileId {};
-    AudioBufferPtr preloadedData {};
-    AudioBuffer<float> fileData {};
+    FileAudioBufferPtr preloadedData {};
+    FileAudioBuffer fileData {};
     float sampleRate { config::defaultSampleRate };
     Oversampling oversamplingFactor { config::defaultOversamplingFactor };
     std::atomic<size_t> availableFrames { 0 };
diff --git a/src/sfizz/Interpolators.h b/src/sfizz/Interpolators.h
new file mode 100644
index 00000000..f784eefd
--- /dev/null
+++ b/src/sfizz/Interpolators.h
@@ -0,0 +1,38 @@
+// 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,
+};
+
+/**
+ * @brief Interpolate from a vector of values
+ *
+ * @tparam M the interpolator model
+ * @tparam R the sample type
+ * @param values Pointer to a value in a larger vector of values.
+ *               Depending on the interpolator the algorithm may
+ *               read samples before and after. Usually you need
+ *               to ensure that you have order - 1 samples available
+ *               before and after the pointer, padding if necessary.
+ * @param coeff the interpolation coefficient
+ * @return R
+ */
+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..1f462578
--- /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 - 1));
+        // 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 = -1; i < 3; ++i) {
+            R h = hermite3<R>(i - 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 - 1));
+        // 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 = -1; i < 3; ++i) {
+            R h = bspline3<R>(i - coeff);
+            y += h * values[i];
+        }
+        return y;
+    }
+};
+
+} // namespace sfz
diff --git a/src/sfizz/MathHelpers.h b/src/sfizz/MathHelpers.h
index d8976c88..88116dab 100644
--- a/src/sfizz/MathHelpers.h
+++ b/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)
@@ -155,12 +159,96 @@ inline CXX14_CONSTEXPR void incrementAll(T& first, Args&... rest)
     incrementAll<Increment>(rest...);
 }
 
-template <class ValueType>
-constexpr ValueType linearInterpolation(ValueType left, ValueType right, ValueType leftCoeff, ValueType rightCoeff)
+/**
+ * @brief Compute the 3rd-order Hermite interpolation polynomial.
+ *
+ * @tparam R
+ * @param x
+ * @return R
+ */
+template <class R>
+R hermite3(R x)
 {
-    return left * leftCoeff + right * rightCoeff;
+    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
+
+/**
+ * @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_add_ps(_mm_sub_ps(_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 Type>
 constexpr Type pi() { return static_cast<Type>(3.141592653589793238462643383279502884); };
 template <class Type>
diff --git a/src/sfizz/Oversampler.cpp b/src/sfizz/Oversampler.cpp
index 422e5250..1027b8ba 100644
--- a/src/sfizz/Oversampler.cpp
+++ b/src/sfizz/Oversampler.cpp
@@ -61,7 +61,7 @@ sfz::Oversampler::Oversampler(sfz::Oversampling factor, size_t chunkSize)
 
 }
 
-void sfz::Oversampler::stream(const sfz::AudioBuffer<float>& input, sfz::AudioBuffer<float>& output, std::atomic<size_t>* framesReady)
+void sfz::Oversampler::stream(AudioSpan<float> input, AudioSpan<float> output, std::atomic<size_t>* framesReady)
 {
     ASSERT(output.getNumFrames() >= input.getNumFrames() * static_cast<int>(factor));
     ASSERT(output.getNumChannels() == input.getNumChannels());
@@ -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;
     }
 
diff --git a/src/sfizz/Oversampler.h b/src/sfizz/Oversampler.h
index d6daa4a9..25e74eaf 100644
--- a/src/sfizz/Oversampler.h
+++ b/src/sfizz/Oversampler.h
@@ -11,6 +11,7 @@
 #include "Debug.h"
 #include "Buffer.h"
 #include "AudioBuffer.h"
+#include "AudioSpan.h"
 #include "Config.h"
 
 namespace sfz {
@@ -39,7 +40,7 @@ public:
      * @param output
      * @param framesReady an atomic counter for the ready frames. If null no signaling is done.
      */
-    void stream(const AudioBuffer<float>& input, AudioBuffer<float>& output, std::atomic<size_t>* framesReady = nullptr);
+    void stream(AudioSpan<float> input, AudioSpan<float> output, std::atomic<size_t>* framesReady = nullptr);
 
     Oversampler() = delete;
     Oversampler(const Oversampler&) = delete;
diff --git a/src/sfizz/Region.cpp b/src/sfizz/Region.cpp
index 20d44048..35e6b0fc 100644
--- a/src/sfizz/Region.cpp
+++ b/src/sfizz/Region.cpp
@@ -60,6 +60,12 @@ bool sfz::Region::parseOpcode(const Opcode& rawOpcode)
             sampleId = FileId(std::move(filename), sampleId.isReverse());
         }
         break;
+    case hash("sample_quality"):
+        {
+            setValueFromOpcode(opcode, sampleQuality, Default::sampleQualityRange);
+            break;
+        }
+        break;
     case hash("direction"):
         sampleId = sampleId.reversed(opcode.value == "reverse");
         break;
diff --git a/src/sfizz/Region.h b/src/sfizz/Region.h
index 1332d1ec..5ce49ff8 100644
--- a/src/sfizz/Region.h
+++ b/src/sfizz/Region.h
@@ -240,6 +240,7 @@ struct Region {
 
     // Sound source: sample playback
     FileId sampleId {}; // Sample
+    int sampleQuality { Default::sampleQuality };
     float delay { Default::delay }; // delay
     float delayRandom { Default::delayRandom }; // delay_random
     int64_t offset { Default::offset }; // offset
diff --git a/src/sfizz/SIMDDummy.cpp b/src/sfizz/SIMDDummy.cpp
index 2a0ca34d..718d8a1c 100644
--- a/src/sfizz/SIMDDummy.cpp
+++ b/src/sfizz/SIMDDummy.cpp
@@ -163,9 +163,9 @@ void sfz::cumsum<float, true>(absl::Span<const float> input, absl::Span<float> o
 }
 
 template<>
-void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> leftCoeffs, absl::Span<float> rightCoeffs) noexcept
+void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> coeffs) noexcept
 {
-    sfzInterpolationCast<float, false>(floatJumps, jumps, leftCoeffs, rightCoeffs);
+    sfzInterpolationCast<float, false>(floatJumps, jumps, coeffs);
 }
 
 template <>
diff --git a/src/sfizz/SIMDHelpers.h b/src/sfizz/SIMDHelpers.h
index 7773b794..0600ddd6 100644
--- a/src/sfizz/SIMDHelpers.h
+++ b/src/sfizz/SIMDHelpers.h
@@ -939,12 +939,11 @@ void cumsum<float, true>(absl::Span<const float> input, absl::Span<float> output
 
 namespace _internals {
     template <class T>
-    void snippetSFZInterpolationCast(const T*& floatJump, int*& jump, T*& leftCoeff, T*& rightCoeff)
+    void snippetSFZInterpolationCast(const T*& floatJump, int*& jump, T*& coeff)
     {
         *jump = static_cast<int>(*floatJump);
-        *rightCoeff = *floatJump - static_cast<float>(*jump);
-        *leftCoeff = static_cast<T>(1.0) - *rightCoeff;
-        incrementAll(floatJump, leftCoeff, rightCoeff, jump);
+        *coeff = *floatJump - static_cast<float>(*jump);
+        incrementAll(floatJump, coeff, jump);
     }
 }
 
@@ -960,24 +959,22 @@ namespace _internals {
  * @param rightCoeffs the right interpolation coefficients
  */
 template <class T, bool SIMD = SIMDConfig::sfzInterpolationCast>
-void sfzInterpolationCast(absl::Span<const T> floatJumps, absl::Span<int> jumps, absl::Span<T> leftCoeffs, absl::Span<T> rightCoeffs) noexcept
+void sfzInterpolationCast(absl::Span<const T> floatJumps, absl::Span<int> jumps, absl::Span<T> coeffs) noexcept
 {
     CHECK(jumps.size() >= floatJumps.size());
-    CHECK(jumps.size() == leftCoeffs.size());
-    CHECK(jumps.size() == rightCoeffs.size());
+    CHECK(jumps.size() == coeffs.size());
 
     auto floatJump = floatJumps.data();
     auto jump = jumps.data();
-    auto leftCoeff = leftCoeffs.data();
-    auto rightCoeff = rightCoeffs.data();
-    const auto sentinel = floatJump + min(floatJumps.size(), jumps.size(), leftCoeffs.size(), rightCoeffs.size());
+    auto coeff = coeffs.data();
+    const auto sentinel = floatJump + min(floatJumps.size(), jumps.size(), coeffs.size());
 
     while (floatJump < sentinel)
-        _internals::snippetSFZInterpolationCast(floatJump, jump, leftCoeff, rightCoeff);
+        _internals::snippetSFZInterpolationCast(floatJump, jump, coeff);
 }
 
 template <>
-void sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> leftCoeffs, absl::Span<float> rightCoeffs) noexcept;
+void sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> coeffs) noexcept;
 
 namespace _internals {
     template <class T>
diff --git a/src/sfizz/SIMDNEON.cpp b/src/sfizz/SIMDNEON.cpp
index 6a6dc8d1..e229d954 100644
--- a/src/sfizz/SIMDNEON.cpp
+++ b/src/sfizz/SIMDNEON.cpp
@@ -227,9 +227,9 @@ void sfz::cumsum<float, true>(absl::Span<const float> input, absl::Span<float> o
 }
 
 template<>
-void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> leftCoeffs, absl::Span<float> rightCoeffs) noexcept
+void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> coeffs) noexcept
 {
-    sfzInterpolationCast<float, false>(floatJumps, jumps, leftCoeffs, rightCoeffs);
+    sfzInterpolationCast<float, false>(floatJumps, jumps, coeffs);
 }
 
 template <>
diff --git a/src/sfizz/SIMDSSE.cpp b/src/sfizz/SIMDSSE.cpp
index 062d5476..ae0cd07e 100644
--- a/src/sfizz/SIMDSSE.cpp
+++ b/src/sfizz/SIMDSSE.cpp
@@ -786,37 +786,33 @@ void sfz::cumsum<float, true>(absl::Span<const float> input, absl::Span<float> o
 }
 
 template <>
-void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> leftCoeffs, absl::Span<float> rightCoeffs) noexcept
+void sfz::sfzInterpolationCast<float, true>(absl::Span<const float> floatJumps, absl::Span<int> jumps, absl::Span<float> coeffs) noexcept
 {
-    sfz::sfzInterpolationCast<float, false>(floatJumps, jumps, leftCoeffs, rightCoeffs);
+    sfz::sfzInterpolationCast<float, false>(floatJumps, jumps, coeffs);
     // CHECK(jumps.size() >= floatJumps.size());
-    // CHECK(jumps.size() == leftCoeffs.size());
-    // CHECK(jumps.size() == rightCoeffs.size());
+    // CHECK(jumps.size() == coeffs.size());
 
     // auto floatJump = floatJumps.data();
     // auto jump = jumps.data();
-    // auto leftCoeff = leftCoeffs.data();
-    // auto rightCoeff = rightCoeffs.data();
-    // const auto sentinel = floatJump + min(floatJumps.size(), jumps.size(), leftCoeffs.size(), rightCoeffs.size());
+    // auto coeff = coeffs.data();
+    // const auto sentinel = floatJump + min(floatJumps.size(), jumps.size(), coeffs.size());
     // const auto lastAligned = prevAligned(sentinel);
 
-    // while (unaligned(floatJump, reinterpret_cast<float*>(jump), leftCoeff, rightCoeff) && floatJump < lastAligned)
-    //     _internals::snippetSFZInterpolationCast(floatJump, jump, leftCoeff, rightCoeff);
+    // while (unaligned(floatJump, reinterpret_cast<float*>(jump), coeff) && floatJump < lastAligned)
+    //     _internals::snippetSFZInterpolationCast(floatJump, jump, coeff);
 
     // while (floatJump < lastAligned) {
     //     auto mmFloatJumps = _mm_load_ps(floatJump);
     //     auto mmIndices = _mm_cvtps_epi32(_mm_sub_ps(mmFloatJumps, _mm_set_ps1(0.4999999552965164184570312f)));
     //     _mm_store_si128(reinterpret_cast<__m128i*>(jump), mmIndices);
 
-    //     auto mmRight = _mm_sub_ps(mmFloatJumps, _mm_cvtepi32_ps(mmIndices));
-    //     auto mmLeft = _mm_sub_ps(_mm_set_ps1(1.0f), mmRight);
-    //     _mm_store_ps(leftCoeff, mmLeft);
-    //     _mm_store_ps(rightCoeff, mmRight);
-    //     incrementAll<TypeAlignment>(floatJump, jump, leftCoeff, rightCoeff);
+    //     auto mmCoeff = _mm_sub_ps(mmFloatJumps, _mm_cvtepi32_ps(mmIndices));
+    //     _mm_store_ps(coeff, mmCoeff);
+    //     incrementAll<TypeAlignment>(floatJump, jump, coeff);
     // }
 
     // while(floatJump < sentinel)
-    //     _internals::snippetSFZInterpolationCast(floatJump, jump, leftCoeff, rightCoeff);
+    //     _internals::snippetSFZInterpolationCast(floatJump, jump, coeff);
 }
 
 template <>
diff --git a/src/sfizz/Voice.cpp b/src/sfizz/Voice.cpp
index 3c4bc442..6ae90db2 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)
@@ -443,10 +444,9 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
 
     auto jumps = resources.bufferPool.getBuffer(numSamples);
     auto bends = resources.bufferPool.getBuffer(numSamples);
-    auto leftCoeffs = resources.bufferPool.getBuffer(numSamples);
-    auto rightCoeffs = resources.bufferPool.getBuffer(numSamples);
+    auto coeffs = resources.bufferPool.getBuffer(numSamples);
     auto indices = resources.bufferPool.getIndexBuffer(numSamples);
-    if (!jumps || !bends || !indices || !rightCoeffs || !leftCoeffs)
+    if (!jumps || !bends || !indices || !coeffs)
         return;
 
     fill<float>(*jumps, pitchRatio * speedRatio);
@@ -470,7 +470,7 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
 
     jumps->front() += floatPositionOffset;
     cumsum<float>(*jumps, *jumps);
-    sfzInterpolationCast<float>(*jumps, *indices, *leftCoeffs, *rightCoeffs);
+    sfzInterpolationCast<float>(*jumps, *indices, *coeffs);
     add<int>(sourcePosition, *indices);
 
     if (region->shouldLoop() && region->loopEnd(currentPromise->oversamplingFactor) <= source.getNumFrames()) {
@@ -486,7 +486,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
@@ -500,35 +500,35 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
 #endif
                 egEnvelope.startRelease(i, true);
                 fill<int>(indices->subspan(i), sampleEnd);
-                fill<float>(leftCoeffs->subspan(i), 0.0f);
-                fill<float>(rightCoeffs->subspan(i), 1.0f);
+                fill<float>(coeffs->subspan(i), 1.0f);
                 break;
             }
         }
     }
 
-    auto ind = indices->data();
-    auto leftCoeff = leftCoeffs->data();
-    auto rightCoeff = rightCoeffs->data();
-    auto leftSource = source.getConstSpan(0);
-    auto left = buffer.getChannel(0);
-    if (source.getNumChannels() == 1) {
-        while (ind < indices->end()) {
-            *left = linearInterpolation(leftSource[*ind], leftSource[*ind + 1], *leftCoeff, *rightCoeff);
-            incrementAll(ind, left, leftCoeff, rightCoeff);
-        }
-    } else {
-        auto right = buffer.getChannel(1);
-        auto rightSource = source.getConstSpan(1);
-        while (ind < indices->end()) {
-            *left = linearInterpolation(leftSource[*ind], leftSource[*ind + 1], *leftCoeff, *rightCoeff);
-            *right = linearInterpolation(rightSource[*ind], rightSource[*ind + 1], *leftCoeff, *rightCoeff);
-            incrementAll(ind, left, right, leftCoeff, rightCoeff);
-        }
+    const int quality = region->sampleQuality;
+
+    switch (quality) {
+    default:
+        if (quality > 2)
+            goto high; // TODO sinc, not implemented
+        // fall through
+    case 1:
+        fillInterpolated<kInterpolatorLinear>(source, buffer, *indices, *coeffs);
+        break;
+    case 2: high:
+#if 1
+        // B-spline response has faster decay of aliasing, but not zero-crossings at integer positions
+        fillInterpolated<kInterpolatorBspline3>(source, buffer, *indices, *coeffs);
+#else
+        // Hermite polynomial
+        fillInterpolated<kInterpolatorHermite3>(source, buffer, *indices, *coeffs);
+#endif
+        break;
     }
 
     sourcePosition = indices->back();
-    floatPositionOffset = rightCoeffs->back();
+    floatPositionOffset = coeffs->back();
 
 #if 0
     ASSERT(!hasNanInf(buffer.getConstSpan(0)));
@@ -538,6 +538,31 @@ void sfz::Voice::fillWithData(AudioSpan<float> buffer) noexcept
 #endif
 }
 
+template <sfz::InterpolatorModel M>
+void sfz::Voice::fillInterpolated(
+    const sfz::AudioSpan<const float>& source, sfz::AudioSpan<float>& dest,
+    absl::Span<const int> indices, absl::Span<const float> coeffs)
+{
+    auto ind = indices.data();
+    auto coeff = coeffs.data();
+    auto leftSource = source.getConstSpan(0);
+    auto left = dest.getChannel(0);
+    if (source.getNumChannels() == 1) {
+        while (ind < indices.end()) {
+            *left = sfz::interpolate<M>(&leftSource[*ind], *coeff);
+            incrementAll(ind, left, coeff);
+        }
+    } else {
+        auto right = dest.getChannel(1);
+        auto rightSource = source.getConstSpan(1);
+        while (ind < indices.end()) {
+            *left = sfz::interpolate<M>(&leftSource[*ind], *coeff);
+            *right = sfz::interpolate<M>(&rightSource[*ind], *coeff);
+            incrementAll(ind, left, right, coeff);
+        }
+    }
+}
+
 void sfz::Voice::fillWithGenerator(AudioSpan<float> buffer) noexcept
 {
     const auto leftSpan = buffer.getSpan(0);
diff --git a/src/sfizz/Voice.h b/src/sfizz/Voice.h
index 46b3a3f0..a709eb58 100644
--- a/src/sfizz/Voice.h
+++ b/src/sfizz/Voice.h
@@ -19,6 +19,7 @@
 #include <random>
 
 namespace sfz {
+enum InterpolatorModel : int;
 /**
  * @brief The SFZ voice are the polyphony holders. They get activated by the synth
  * and tasked to play a given region until the end, stopping on note-offs, off-groups
@@ -242,6 +243,20 @@ private:
      * @param buffer
      */
     void fillWithGenerator(AudioSpan<float> buffer) noexcept;
+
+    /**
+     * @brief Fill a destination with an interpolated source.
+     *
+     * @param source the source sample
+     * @param dest the destination buffer
+     * @param indices the integral parts of the source positions
+     * @param coeffs the fractional parts of the source positions
+     */
+    template <InterpolatorModel M>
+    static void fillInterpolated(
+        const AudioSpan<const float>& source, AudioSpan<float>& dest,
+        absl::Span<const int> indices, absl::Span<const float> coeffs);
+
     void amplitudeEnvelope(absl::Span<float> modulationSpan) noexcept;
     void ampStageMono(AudioSpan<float> buffer) noexcept;
     void ampStageStereo(AudioSpan<float> buffer) noexcept;
diff --git a/tests/AudioBufferT.cpp b/tests/AudioBufferT.cpp
index ff483b13..f1665a89 100644
--- a/tests/AudioBufferT.cpp
+++ b/tests/AudioBufferT.cpp
@@ -83,3 +83,38 @@ TEST_CASE("[AudioSpan] Constructions")
     sfz::AudioSpan<float> manualSpan2 { {buffer.getSpan(0), buffer.getSpan(1) } };
     sfz::AudioSpan<const float> manualConstSpan2 { { buffer.getConstSpan(0), buffer.getConstSpan(1) } };
 }
+
+TEST_CASE("[AudioBuffer] Padding")
+{
+    constexpr size_t channels = 2;
+    constexpr size_t numFrames = 7777;
+    constexpr unsigned alignment = 32;
+    constexpr size_t padLeft = 13;
+    constexpr size_t padRight = 51;
+
+    sfz::AudioBuffer<float, channels, alignment, padLeft, padRight> padded { channels, numFrames };
+
+    // ensure padding to be extended to respect alignment
+    constexpr size_t extendedPadLeft = padded.PaddingLeft;
+    REQUIRE(extendedPadLeft == 32);
+
+    // ensure access functions to return consistent addresses with offset
+    for (size_t c = 0; c < channels; ++c) {
+        REQUIRE(padded.getSpan(c).data() == &padded.getSample(c, 0));
+        REQUIRE(padded.getConstSpan(c).data() == &padded.getSample(c, 0));
+        REQUIRE(padded.getSpan(c).data() == padded.channelReader(c));
+        REQUIRE(padded.getSpan(c).data() == padded.channelWriter(c));
+        REQUIRE(padded.getSpan(c).end() == padded.channelReaderEnd(c));
+        REQUIRE(padded.getSpan(c).end() == padded.channelWriterEnd(c));
+    }
+
+    padded.clear();
+
+    // ensure all padding areas to be accessible and zero after clearing
+    for (size_t c = 0; c < channels; ++c) {
+        absl::Span<const float> leftSpan { padded.getSpan(0).data() - extendedPadLeft, extendedPadLeft };
+        absl::Span<const float> rightSpan { padded.getSpan(0).end(), padRight };
+        REQUIRE(std::all_of(leftSpan.begin(), leftSpan.end(), [](float value) { return value == 0.0f; }));
+        REQUIRE(std::all_of(rightSpan.begin(), rightSpan.end(), [](float value) { return value == 0.0f; }));
+    }
+}
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index c871114c..af4981dc 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -17,6 +17,7 @@ set(SFIZZ_TEST_SOURCES
     FilesT.cpp
     MidiStateT.cpp
     OnePoleFilterT.cpp
+    InterpolatorsT.cpp
     RegionActivationT.cpp
     RegionValueComputationsT.cpp
     # If we're tweaking the curves this kind of tests does not make sense
diff --git a/tests/InterpolatorsT.cpp b/tests/InterpolatorsT.cpp
new file mode 100644
index 00000000..04742547
--- /dev/null
+++ b/tests/InterpolatorsT.cpp
@@ -0,0 +1,71 @@
+// 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 "sfizz/Interpolators.h"
+#include "catch2/catch.hpp"
+#include <algorithm>
+using namespace Catch::literals;
+
+TEST_CASE("[Interpolators] Sample at points")
+{
+    std::array<float, 32> values;
+    std::iota(values.begin(), values.end(), 0.0f);
+    for (unsigned i = 2; i < values.size() - 2; ++i) {
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorLinear>(&values[i], 0.0f)
+            == Approx(values[i]).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorHermite3>(&values[i], 0.0f)
+            == Approx(values[i]).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorBspline3>(&values[i], 0.0f)
+            == Approx(values[i]).margin(1e-2));
+    }
+}
+
+TEST_CASE("[Interpolators] Sample next")
+{
+    std::array<float, 32> values;
+    std::iota(values.begin(), values.end(), 0.0f);
+    for (unsigned i = 2; i < values.size() - 2; ++i) {
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorLinear>(&values[i], 1.0f)
+            == Approx(values[i + 1]).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorHermite3>(&values[i], 1.0f)
+            == Approx(values[i + 1]).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorBspline3>(&values[i], 1.0f)
+            == Approx(values[i + 1]).margin(1e-2));
+    }
+}
+
+TEST_CASE("[Interpolators] Straight line")
+{
+    std::array<float, 32> values;
+    std::iota(values.begin(), values.end(), 0.0f);
+    for (unsigned i = 2; i < values.size() - 2; ++i) {
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorLinear>(&values[i], 0.5f)
+            == Approx(values[i] + 0.5f).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorHermite3>(&values[i], 0.5f)
+            == Approx(values[i] + 0.5f).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorBspline3>(&values[i], 0.5f)
+            == Approx(values[i] + 0.5f).margin(1e-2));
+    }
+}
+
+TEST_CASE("[Interpolators] Squares")
+{
+    std::array<float, 32> x;
+    std::array<float, 32> y;
+    for (unsigned i = 0; i < x.size(); ++i) {
+        x[i] = static_cast<float>(i) / static_cast<float>(x.size());
+        y[i] = x[i] * x[i];
+    }
+
+    for (unsigned i = 2; i < x.size() - 2; ++i) {
+        const auto half_x = x[i] + 0.5f / static_cast<float>(x.size());
+        const auto expected = (half_x) * (half_x);
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorHermite3>(&y[i], 0.5f)
+            == Approx(expected).margin(1e-2));
+        REQUIRE(sfz::interpolate<sfz::InterpolatorModel::kInterpolatorBspline3>(&y[i], 0.5f)
+            == Approx(expected).margin(1e-2));
+    }
+}