diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3128d06f..ab5a2861 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -168,6 +168,9 @@ target_link_libraries(bm_gain benchmark absl::span)
 add_executable(bm_looping benchmarks/BM_looping.cpp sources/SIMDSSE.cpp)
 target_link_libraries(bm_looping benchmark absl::span absl::algorithm)
 
+add_executable(bm_saturating benchmarks/BM_saturating.cpp sources/SIMDSSE.cpp)
+target_link_libraries(bm_saturating benchmark absl::span absl::algorithm)
+
 add_executable(bm_ramp benchmarks/BM_ramp.cpp sources/SIMDSSE.cpp)
 target_link_libraries(bm_ramp benchmark absl::span absl::algorithm)
 
diff --git a/benchmarks/BM_saturating.cpp b/benchmarks/BM_saturating.cpp
new file mode 100644
index 00000000..6ede3178
--- /dev/null
+++ b/benchmarks/BM_saturating.cpp
@@ -0,0 +1,71 @@
+#include <benchmark/benchmark.h>
+#include <vector>
+#include <random>
+#include <numeric>
+#include <absl/algorithm/container.h>
+#include "../sources/SIMDHelpers.h"
+
+// In this one we have an array of indices
+
+constexpr int loopEnd { 1076 };
+constexpr float maxJump { 4 };
+
+class SaturatingFixture : public benchmark::Fixture {
+public:
+  void SetUp(const ::benchmark::State& state) {
+    std::random_device rd { };
+    std::mt19937 gen { rd() };
+    std::uniform_real_distribution<float> dist { 0, maxJump };
+    indices = std::vector<int>(state.range(0));
+    leftCoeffs = std::vector<float>(state.range(0));
+    rightCoeffs = std::vector<float>(state.range(0));
+    jumps = std::vector<float>(state.range(0));
+    absl::c_generate(jumps, [&]() { return dist(gen); });
+  }
+
+  void TearDown(const ::benchmark::State& state [[maybe_unused]]) {
+
+  }
+
+    std::vector<int> indices;
+    std::vector<float> leftCoeffs;
+    std::vector<float> rightCoeffs;
+    std::vector<float> jumps;
+};
+
+
+BENCHMARK_DEFINE_F(SaturatingFixture, Scalar)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        saturatingSFZIndex<float, false>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 2.5f, loopEnd);
+    }
+}
+
+BENCHMARK_DEFINE_F(SaturatingFixture, SIMD)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        saturatingSFZIndex<float, true>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 2.5f, loopEnd);
+    }
+}
+
+BENCHMARK_DEFINE_F(SaturatingFixture, Scalar_Unaligned)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        saturatingSFZIndex<float, false>(absl::MakeSpan(jumps).subspan(1), absl::MakeSpan(leftCoeffs).subspan(2), absl::MakeSpan(rightCoeffs).subspan(1), absl::MakeSpan(indices).subspan(3), 2.5f, loopEnd);
+    }
+}
+
+BENCHMARK_DEFINE_F(SaturatingFixture, SIMD_Unaligned)(benchmark::State& state) {
+    for (auto _ : state)
+    {
+        saturatingSFZIndex<float, true>(absl::MakeSpan(jumps).subspan(1), absl::MakeSpan(leftCoeffs).subspan(2), absl::MakeSpan(rightCoeffs).subspan(1), absl::MakeSpan(indices).subspan(3), 2.5f, loopEnd);
+    }
+}
+
+
+// Register the function as a benchmark
+BENCHMARK_REGISTER_F(SaturatingFixture, Scalar)->RangeMultiplier(2)->Range((2<<6), (2<<12));
+BENCHMARK_REGISTER_F(SaturatingFixture, SIMD)->RangeMultiplier(2)->Range((2<<6), (2<<12));
+BENCHMARK_REGISTER_F(SaturatingFixture, Scalar_Unaligned)->RangeMultiplier(2)->Range((2<<6), (2<<12));
+BENCHMARK_REGISTER_F(SaturatingFixture, SIMD_Unaligned)->RangeMultiplier(2)->Range((2<<6), (2<<12));
+BENCHMARK_MAIN();
\ No newline at end of file
diff --git a/sources/SIMDHelpers.h b/sources/SIMDHelpers.h
index be685c8b..9792d13f 100644
--- a/sources/SIMDHelpers.h
+++ b/sources/SIMDHelpers.h
@@ -109,6 +109,50 @@ void cos(absl::Span<const Type> input, absl::Span<Type> output) noexcept
 template <>
 void cos<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept;
 
+template <>
+void cos<float, true>(absl::Span<const float> input, absl::Span<float> output) noexcept;
+
+template <class T>
+inline void snippetSaturatingIndex(const T*& jump, T*& leftCoeff, T*& rightCoeff, int*& index, T& floatIndex, T loopEnd)
+{
+    floatIndex += *jump;
+    if (floatIndex >= loopEnd) {
+        floatIndex = loopEnd;
+        *index = static_cast<int>(floatIndex) - 1;
+        *rightCoeff = static_cast<T>(1.0);
+        *leftCoeff = static_cast<T>(0.0);
+    } else {
+        *index = static_cast<int>(floatIndex);
+        *rightCoeff = floatIndex - *index;
+        *leftCoeff = static_cast<T>(1.0) - *rightCoeff;
+    }
+    index++;
+    leftCoeff++;
+    rightCoeff++;
+    jump++;
+}
+
+template <class T, bool SIMD = SIMDConfig::saturatingSFZIndex>
+void saturatingSFZIndex(absl::Span<const T> jumps, absl::Span<T> leftCoeffs, absl::Span<T> rightCoeffs, absl::Span<int> indices, T floatIndex, T loopEnd) noexcept
+{
+    ASSERT(indices.size() >= jumps.size());
+    ASSERT(indices.size() == leftCoeffs.size());
+    ASSERT(indices.size() == rightCoeffs.size());
+
+    auto* index = indices.begin();
+    auto* leftCoeff = leftCoeffs.begin();
+    auto* rightCoeff = rightCoeffs.begin();
+    auto* jump = jumps.begin();
+    const auto size = min(jumps.size(), indices.size(), leftCoeffs.size(), rightCoeffs.size());
+    auto* sentinel = jumps.begin() + size;
+
+    while (jump < sentinel)
+        snippetSaturatingIndex<T>(jump, leftCoeff, rightCoeff, index, floatIndex, loopEnd);
+}
+
+template <>
+void saturatingSFZIndex<float, true>(absl::Span<const float> jumps, absl::Span<float> leftCoeffs, absl::Span<float> rightCoeffs, absl::Span<int> indices, float floatIndex, float loopEnd) noexcept;
+
 template <class T>
 inline void snippetLoopingIndex(const T*& jump, T*& leftCoeff, T*& rightCoeff, int*& index, T& floatIndex, T loopEnd, T loopStart)
 {
diff --git a/sources/SIMDSSE.cpp b/sources/SIMDSSE.cpp
index d48c56ac..98591445 100644
--- a/sources/SIMDSSE.cpp
+++ b/sources/SIMDSSE.cpp
@@ -343,6 +343,63 @@ void loopingSFZIndex<float, true>(absl::Span<const float> jumps,
         snippetLoopingIndex<float>(jump, leftCoeff, rightCoeff, index, floatIndex, loopEnd, loopStart);
 }
 
+template <>
+void saturatingSFZIndex<float, true>(absl::Span<const float> jumps,
+    absl::Span<float> leftCoeffs,
+    absl::Span<float> rightCoeffs,
+    absl::Span<int> indices,
+    float floatIndex,
+    float loopEnd) noexcept
+{
+    ASSERT(indices.size() >= jumps.size());
+    ASSERT(indices.size() == leftCoeffs.size());
+    ASSERT(indices.size() == rightCoeffs.size());
+
+    auto index = indices.data();
+    auto leftCoeff = leftCoeffs.data();
+    auto rightCoeff = rightCoeffs.data();
+    auto jump = jumps.data();
+    const auto size = min(jumps.size(), indices.size(), leftCoeffs.size(), rightCoeffs.size());
+    const auto* sentinel = jumps.begin() + size;
+    const auto* alignedEnd = prevAligned(sentinel);
+
+    while (unaligned(reinterpret_cast<float*>(index), leftCoeff, rightCoeff, jump) && jump < alignedEnd)
+        snippetSaturatingIndex<float>(jump, leftCoeff, rightCoeff, index, floatIndex, loopEnd);
+
+    auto mmFloatIndex = _mm_set_ps1(floatIndex);
+    const auto mmLoopEnd = _mm_set1_ps(loopEnd);
+    const auto mmSaturated = _mm_sub_ps(mmLoopEnd, _mm_set_ps1(0.000001f));
+    while (jump < alignedEnd) {
+        auto mmOffset = _mm_load_ps(jump);
+        mmOffset = _mm_add_ps(mmOffset, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(mmOffset), 4)));
+        mmOffset = _mm_add_ps(mmOffset, _mm_shuffle_ps(_mm_setzero_ps(), mmOffset, 0x40));
+
+        mmFloatIndex = _mm_add_ps(mmFloatIndex, mmOffset);
+        const auto mmCompared = _mm_cmplt_ps(mmFloatIndex, mmLoopEnd);
+        mmFloatIndex = _mm_add_ps(_mm_and_ps(mmCompared, mmFloatIndex), _mm_andnot_ps(mmCompared, mmSaturated));
+
+        auto mmIndices = _mm_cvtps_epi32(_mm_sub_ps(mmFloatIndex, _mm_set_ps1(0.4999999552965164184570312f)));
+        _mm_store_si128(reinterpret_cast<__m128i*>(index), mmIndices);
+
+        auto mmRight = _mm_sub_ps(mmFloatIndex, _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);
+
+        mmFloatIndex = _mm_shuffle_ps(mmFloatIndex, mmFloatIndex, _MM_SHUFFLE(3, 3, 3, 3));
+        // floatingIndex = _mm_cvtss_f32(_mm_shuffle_ps(mmFloatIndex, mmFloatIndex, _MM_SHUFFLE(0, 0, 0, 3)));;
+        // floatingIndex = *(index + 3) + *(rightCoeff + 3);
+        index += TypeAlignment;
+        jump += TypeAlignment;
+        leftCoeff += TypeAlignment;
+        rightCoeff += TypeAlignment;
+    }
+
+    floatIndex = _mm_cvtss_f32(mmFloatIndex);
+    while (jump < sentinel)
+        snippetSaturatingIndex<float>(jump, leftCoeff, rightCoeff, index, floatIndex, loopEnd);
+}
+
 template <>
 float linearRamp<float, true>(absl::Span<float> output, float value, float step) noexcept
 {
diff --git a/tests/SIMDHelpersT.cpp b/tests/SIMDHelpersT.cpp
index 005bd42d..233c1d17 100644
--- a/tests/SIMDHelpersT.cpp
+++ b/tests/SIMDHelpersT.cpp
@@ -12,6 +12,21 @@ constexpr int bigBufferSize { 4095 };
 constexpr int medBufferSize { 127 };
 constexpr double fillValue { 1.3 };
 
+template <class Type>
+inline bool approxEqualMargin(absl::Span<const Type> lhs, absl::Span<const Type> rhs, Type eps = 1e-3)
+{
+    if (lhs.size() != rhs.size())
+        return false;
+
+    for (size_t i = 0; i < rhs.size(); ++i)
+        if (rhs[i] != Approx(lhs[i]).margin(eps)) {
+            std::cerr << lhs[i] << " != " << rhs[i] << " at index " << i << '\n';
+            return false;
+        }
+
+    return true;
+}
+
 template <class Type>
 inline bool approxEqual(absl::Span<const Type> lhs, absl::Span<const Type> rhs, Type eps = 1e-3)
 {
@@ -419,11 +434,60 @@ TEST_CASE("[Helpers] SFZ looping index (SIMD)")
 //     std::vector<float> rightCoeffsSIMD(bigBufferSize);
 //     loopingSFZIndex<float, false>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 1.0f, medBufferSize, 1);
 //     loopingSFZIndex<float, true>(jumps, absl::MakeSpan(leftCoeffsSIMD), absl::MakeSpan(rightCoeffsSIMD), absl::MakeSpan(indicesSIMD), 1.0f, medBufferSize, 1);
-//     REQUIRE( approxEqual<int>(indices, indicesSIMD, 1) );
-//     REQUIRE( approxEqual<float>(leftCoeffs, leftCoeffsSIMD, 1e-2) );
-//     REQUIRE( approxEqual<float>(rightCoeffs, rightCoeffsSIMD, 1e-2) );
+//     for (int i = 0; i < bigBufferSize; ++i)
+//         REQUIRE( ((static_cast<float>(indices[i]) + rightCoeffs[i] == Approx(static_cast<float>(indicesSIMD[i]) + rightCoeffsSIMD[i]).margin(1e-2)) 
+//                 || (static_cast<float>(indices[i]) + rightCoeffs[i] == Approx(static_cast<float>(indicesSIMD[i]) + rightCoeffsSIMD[i] - static_cast<float>(medBufferSize)).margin(2e-2))) );
 // }
 
+TEST_CASE("[Helpers] SFZ saturating index")
+{
+    std::array<float, 6> jumps { 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f }; // 1.1 2.3 3.6 5.0 6.5 8.1
+    std::array<int, 6> indices;
+    std::array<float, 6> leftCoeffs;
+    std::array<float, 6> rightCoeffs;
+    std::array<int, 6> expectedIndices { 2, 3, 4, 5, 5, 5 };
+    std::array<float, 6> expectedLeft { 0.9f, 0.7f, 0.4f, 0.0f, 0.0f, 0.0f };
+    std::array<float, 6> expectedRight { 0.1f, 0.3f, 0.6f, 1.0f, 1.0f, 1.0f };
+    saturatingSFZIndex<float, false>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 1.0f, 6);
+    REQUIRE(indices == expectedIndices);
+    REQUIRE(approxEqual<float>(leftCoeffs, expectedLeft));
+    REQUIRE(approxEqual<float>(rightCoeffs, expectedRight));
+}
+
+TEST_CASE("[Helpers] SFZ saturating index (SIMD)")
+{
+    std::array<float, 6> jumps { 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f }; // 1.1 2.3 3.6 5.0 6.5 8.1
+    std::array<int, 6> indices;
+    std::array<float, 6> leftCoeffs;
+    std::array<float, 6> rightCoeffs;
+    std::array<int, 6> expectedIndices { 2, 3, 4, 5, 5, 5 };
+    std::array<float, 6> expectedLeft { 0.9f, 0.7f, 0.4f, 0.0f, 0.0f, 0.0f };
+    std::array<float, 6> expectedRight { 0.1f, 0.3f, 0.6f, 1.0f, 1.0f, 1.0f };
+    saturatingSFZIndex<float, true>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 1.0f, 6);
+    REQUIRE(indices == expectedIndices);
+    REQUIRE(approxEqualMargin<float>(leftCoeffs, expectedLeft));
+    REQUIRE(approxEqualMargin<float>(rightCoeffs, expectedRight));
+}
+
+TEST_CASE("[Helpers] SFZ saturating index (SIMD vs Scalar)")
+{
+
+    std::vector<float> jumps(medBufferSize);
+    absl::c_fill(jumps, fillValue);
+
+    std::vector<int> indices(medBufferSize);
+    std::vector<float> leftCoeffs(medBufferSize);
+    std::vector<float> rightCoeffs(medBufferSize);
+
+    std::vector<int> indicesSIMD(medBufferSize);
+    std::vector<float> leftCoeffsSIMD(medBufferSize);
+    std::vector<float> rightCoeffsSIMD(medBufferSize);
+    saturatingSFZIndex<float, false>(jumps, absl::MakeSpan(leftCoeffs), absl::MakeSpan(rightCoeffs), absl::MakeSpan(indices), 1.0f, 78);
+    saturatingSFZIndex<float, true>(jumps, absl::MakeSpan(leftCoeffsSIMD), absl::MakeSpan(rightCoeffsSIMD), absl::MakeSpan(indicesSIMD), 1.0f, 78);
+    for (int i = 0; i < medBufferSize; ++i)
+        REQUIRE( static_cast<float>(indices[i]) + rightCoeffs[i] == Approx(static_cast<float>(indicesSIMD[i]) + rightCoeffsSIMD[i]));
+}
+
 TEST_CASE("[Helpers] Linear Ramp")
 {
     const float start { 0.0f };