diff --git a/cmake/SfizzSIMDSourceFiles.cmake b/cmake/SfizzSIMDSourceFiles.cmake
index efd86bb7..3cc9afdd 100644
--- a/cmake/SfizzSIMDSourceFiles.cmake
+++ b/cmake/SfizzSIMDSourceFiles.cmake
@@ -2,7 +2,9 @@ macro(sfizz_add_simd_sources SOURCES_VAR PREFIX)
     # It needs a macro, otherwise the source properties cannot take effect.
 
     list (APPEND ${SOURCES_VAR}
-        ${PREFIX}/sfizz/SIMDHelpers.cpp)
+        ${PREFIX}/sfizz/SIMDHelpers.cpp
+        ${PREFIX}/sfizz/simd/HelpersSSE.cpp
+        ${PREFIX}/sfizz/simd/HelpersAVX.cpp)
 
     # For CPU-dispatched X86 sources
     # Always build them for all X86 targets.
@@ -13,7 +15,7 @@ macro(sfizz_add_simd_sources SOURCES_VAR PREFIX)
             set_source_files_properties(
                 ${PREFIX}/sfizz/effects/impl/ResonantStringAVX.cpp
                 ${PREFIX}/sfizz/effects/impl/ResonantArrayAVX.cpp
-                ${PREFIX}/sfizz/SIMDHelpers.cpp
+                ${PREFIX}/sfizz/simd/HelpersAVX.cpp
                 PROPERTIES COMPILE_FLAGS "-mavx")
         endif()
     endif()
diff --git a/dpf.mk b/dpf.mk
index e9bfb051..c14758b6 100644
--- a/dpf.mk
+++ b/dpf.mk
@@ -97,6 +97,8 @@ SFIZZ_SOURCES = \
 	src/sfizz/SfzFilter.cpp \
 	src/sfizz/SfzHelpers.cpp \
 	src/sfizz/SIMDHelpers.cpp \
+	src/sfizz/simd/SSEHelpers.cpp \
+	src/sfizz/simd/AVXHelpers.cpp \
 	src/sfizz/Synth.cpp \
 	src/sfizz/Tuning.cpp \
 	src/sfizz/Voice.cpp \
diff --git a/scripts/run_clang_tidy.sh b/scripts/run_clang_tidy.sh
index e9d7ad53..83998476 100755
--- a/scripts/run_clang_tidy.sh
+++ b/scripts/run_clang_tidy.sh
@@ -18,6 +18,8 @@ clang-tidy \
   src/sfizz/Region.cpp \
   src/sfizz/SfzHelpers.cpp \
   src/sfizz/SIMDHelpers.cpp \
+  src/sfizz/simd/HelpersSSE.cpp \
+  src/sfizz/simd/HelpersAVX.cpp \
   src/sfizz/Synth.cpp \
   src/sfizz/Voice.cpp \
   src/sfizz/effects/Eq.cpp \
diff --git a/src/sfizz/SIMDHelpers.cpp b/src/sfizz/SIMDHelpers.cpp
index f495127f..359d025a 100644
--- a/src/sfizz/SIMDHelpers.cpp
+++ b/src/sfizz/SIMDHelpers.cpp
@@ -1,18 +1,16 @@
+// 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 "SIMDHelpers.h"
 #include <array>
 #include "cpuid/cpuinfo.hpp"
-
+#include "simd/HelpersSSE.h"
+#include "simd/HelpersAVX.h"
 #include "SIMDConfig.h"
 
-#if SFIZZ_HAVE_SSE2
-#include <immintrin.h>
-#include <emmintrin.h>
-#endif
-
-#if SFIZZ_HAVE_NEON
-#include <arm_neon.h>
-#endif
-
 namespace sfz {
 
 static std::array<bool, static_cast<unsigned>(SIMDOps::_sentinel)> simdStatus;
@@ -57,606 +55,186 @@ bool getSIMDOpStatus(SIMDOps op)
     return simdStatus[static_cast<unsigned>(op)];
 }
 
-constexpr uintptr_t TypeAlignment = 4;
-
-template <class T>
-inline void tickRead(const T*& input, T*& outputLeft, T*& outputRight)
-{
-    *outputLeft++ = *input++;
-    *outputRight++ = *input++;
-}
-
-template <class T>
-inline void tickWrite(T*& output, const T*& inputLeft, const T*& inputRight)
-{
-    *output++ = *inputLeft++;
-    *output++ = *inputRight++;
-}
-
 void readInterleaved(const float* input, float* outputLeft, float* outputRight, unsigned inputSize) noexcept
 {
-    const auto sentinel = input + inputSize - 1;
-
     if (getSIMDOpStatus(SIMDOps::readInterleaved)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(input + inputSize - 4);
-            while (unaligned(input, outputLeft, outputRight) && input < lastAligned)
-                tickRead(input, outputLeft, outputRight);
-
-            while (input < lastAligned) {
-                auto register0 = _mm_load_ps(input);
-                auto register1 = _mm_load_ps(input + 4);
-                auto register2 = register0;
-                // register 2 holds the copy of register 0 that is going to get erased by the first operation
-                // Remember that the bit mask reads from the end; 10 00 10 00 means
-                // "take 0 from a, take 2 from a, take 0 from b, take 2 from b"
-                register0 = _mm_shuffle_ps(register0, register1, 0b10001000);
-                register1 = _mm_shuffle_ps(register2, register1, 0b11011101);
-                _mm_store_ps(outputLeft, register0);
-                _mm_store_ps(outputRight, register1);
-                incrementAll<4>(input, input, outputLeft, outputRight);
-            }
-            // Fallthrough from lastAligned to sentinel
-        }
-#endif
-
-#if 0 // NEON wip
-        auto reg = vld2q_f32(in);
-        vst1q_f32(lOut, reg.val[0]);
-        vst1q_f32(rOut, reg.val[1]);
-#endif
+        if (cpuInfo.has_sse())
+            return readInterleavedSSE(input, outputLeft, outputRight, inputSize);
     }
-
-    while (input < sentinel)
-        tickRead(input, outputLeft, outputRight);
+    return readInterleavedScalar(input, outputLeft, outputRight, inputSize);
 }
 
 void writeInterleaved(const float* inputLeft, const float* inputRight, float* output, unsigned outputSize) noexcept
 {
-    const auto sentinel = output + outputSize - 1;
-
     if (getSIMDOpStatus(SIMDOps::writeInterleaved)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(output + outputSize - 4);
-
-            while (unaligned(output, inputRight, inputLeft) && output < lastAligned)
-                tickWrite(output, inputLeft, inputRight);
-
-            while (output < lastAligned) {
-                const auto lInRegister = _mm_load_ps(inputLeft);
-                const auto rInRegister = _mm_load_ps(inputRight);
-                const auto outRegister1 = _mm_unpacklo_ps(lInRegister, rInRegister);
-                _mm_store_ps(output, outRegister1);
-                const auto outRegister2 = _mm_unpackhi_ps(lInRegister, rInRegister);
-                _mm_store_ps(output + 4, outRegister2);
-                incrementAll<4>(output, output, inputLeft, inputRight);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return writeInterleavedSSE(inputLeft, inputRight, output, outputSize);
     }
-
-    while (output < sentinel)
-        tickWrite(output, inputLeft, inputRight);
+    return writeInterleavedScalar(inputLeft, inputRight, output, outputSize);
 }
 
 template <>
 void applyGain<float>(float gain, const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::gain)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_avx()) {
-            const auto* lastAligned = prevAligned<32>(sentinel);
-            const auto mmGain = _mm256_set1_ps(gain);
-            while (unaligned<32>(input, output) && output < lastAligned)
-                *output++ = gain * (*input++);
-
-            while (output < lastAligned) {
-                _mm256_store_ps(output, _mm256_mul_ps(mmGain, _mm256_load_ps(input)));
-                incrementAll<8>(input, output);
-            }
-        } else if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            const auto mmGain = _mm_set1_ps(gain);
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ = gain * (*input++);
-
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_mul_ps(mmGain, _mm_load_ps(input)));
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_avx())
+            return applyGainAVX(gain, input, output, size);
+        else if (cpuInfo.has_sse())
+            return applyGainSSE(gain, input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ = gain * (*input++);
+    return applyGainScalar(gain, input, output, size);
 }
 
 template <>
 void applyGain<float>(const float* gain, const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::gain)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_avx()) {
-            const auto* lastAligned = prevAligned<32>(sentinel);
-            while (unaligned<32>(input, output) && output < lastAligned)
-                *output++ = (*gain++) * (*input++);
-
-            while (output < lastAligned) {
-                _mm256_store_ps(output, _mm256_mul_ps(_mm256_load_ps(gain), _mm256_load_ps(input)));
-                incrementAll<8>(input, output);
-            }
-        } else if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ = (*gain++) * (*input++);
-
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_mul_ps(_mm_load_ps(gain), _mm_load_ps(input)));
-                incrementAll<4>(gain, input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_avx())
+            return applyGainAVX(gain, input, output, size);
+        else if (cpuInfo.has_sse())
+            return applyGainSSE(gain, input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ = (*gain++) * (*input++);
+    return applyGainScalar(gain, input, output, size);
 }
 
 template <>
 void divide<float>(const float* input, const float* divisor, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::divide)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ = (*input++) / (*divisor++);
-
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_div_ps(_mm_load_ps(input), _mm_load_ps(divisor)));
-                incrementAll<4>(divisor, input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return divideSSE(input, divisor, output, size);
     }
-
-    while (output < sentinel)
-        *output++ = (*input++) / (*divisor++);
+    return divideScalar(input, divisor, output, size);
 }
 
 template <>
 void multiplyAdd<float>(const float* gain, const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::multiplyAdd)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ += (*gain++) * (*input++);
-
-            while (output < lastAligned) {
-                auto mmOut = _mm_load_ps(output);
-                mmOut = _mm_add_ps(_mm_mul_ps(_mm_load_ps(gain), _mm_load_ps(input)), mmOut);
-                _mm_store_ps(output, mmOut);
-                incrementAll<4>(gain, input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return multiplyAddSSE(gain, input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ += (*gain++) * (*input++);
+    return multiplyAddScalar(gain, input, output, size);
 }
 
 template <>
 void multiplyAdd<float>(float gain, const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::multiplyAdd)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ += gain * (*input++);
-
-            auto mmGain = _mm_set1_ps(gain);
-            while (output < lastAligned) {
-                auto mmOut = _mm_load_ps(output);
-                mmOut = _mm_add_ps(_mm_mul_ps(mmGain, _mm_load_ps(input)), mmOut);
-                _mm_store_ps(output, mmOut);
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return multiplyAddSSE(gain, input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ += gain * (*input++);
+    return multiplyAddScalar(gain, input, output, size);
 }
 
 template <>
 float linearRamp<float>(float* output, float start, float step, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::linearRamp)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(output) && output < lastAligned) {
-                *output++ = start;
-                start += step;
-            }
-
-            auto mmStart = _mm_set1_ps(start - step);
-            auto mmStep = _mm_set_ps(step + step + step + step, step + step + step, step + step, step);
-            while (output < lastAligned) {
-                mmStart = _mm_add_ps(mmStart, mmStep);
-                _mm_store_ps(output, mmStart);
-                mmStart = _mm_shuffle_ps(mmStart, mmStart, _MM_SHUFFLE(3, 3, 3, 3));
-                incrementAll<4>(output);
-            }
-            start = _mm_cvtss_f32(mmStart) + step;
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return linearRampSSE(output, start, step, size);
     }
-
-    while (output < sentinel) {
-        *output++ = start;
-        start += step;
-    }
-    return start;
+    return linearRampScalar(output, start, step, size);
 }
 
 template <>
 float multiplicativeRamp<float>(float* output, float start, float step, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::multiplicativeRamp)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(output) && output < lastAligned) {
-                *output++ = start;
-                start *= step;
-            }
-
-            auto mmStart = _mm_set1_ps(start / step);
-            auto mmStep = _mm_set_ps(step * step * step * step, step * step * step, step * step, step);
-            while (output < lastAligned) {
-                mmStart = _mm_mul_ps(mmStart, mmStep);
-                _mm_store_ps(output, mmStart);
-                mmStart = _mm_shuffle_ps(mmStart, mmStart, _MM_SHUFFLE(3, 3, 3, 3));
-                incrementAll<4>(output);
-            }
-            start = _mm_cvtss_f32(mmStart) * step;
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return multiplicativeRampSSE(output, start, step, size);
     }
-
-    while (output < sentinel) {
-        *output++ = start;
-        start *= step;
-    }
-    return start;
+    return multiplicativeRampScalar(output, start, step, size);
 }
 
 template <>
 void add<float>(const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::add)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ += *input++;
-
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_add_ps(_mm_load_ps(output), _mm_load_ps(input)));
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return addSSE(input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ += *input++;
+    return addScalar(input, output, size);
 }
 
 template <>
 void add<float>(float value, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::add)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(output) && output < lastAligned)
-                *output++ += value;
-
-            const auto mmValue = _mm_set1_ps(value);
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_add_ps(_mm_load_ps(output), mmValue));
-                incrementAll<4>(output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return addSSE(value, output, size);
     }
-
-    while (output < sentinel)
-        *output++ += value;
+    return addScalar(value, output, size);
 }
 
 template <>
 void subtract<float>(const float* input, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::subtract)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned)
-                *output++ -= *input++;
-
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_sub_ps(_mm_load_ps(output), _mm_load_ps(input)));
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return subtractSSE(input, output, size);
     }
-
-    while (output < sentinel)
-        *output++ -= *input++;
+    return subtractScalar(input, output, size);
 }
 
 template <>
 void subtract<float>(float value, float* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-
     if (getSIMDOpStatus(SIMDOps::subtract)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(output) && output < lastAligned)
-                *output++ -= value;
-
-            const auto mmValue = _mm_set1_ps(value);
-            while (output < lastAligned) {
-                _mm_store_ps(output, _mm_sub_ps(_mm_load_ps(output), mmValue));
-                incrementAll<4>(output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return subtractSSE(value, output, size);
     }
-
-    while (output < sentinel)
-        *output++ -= value;
+    return subtractScalar(value, output, size);
 }
 
 template <>
 void copy<float>(const float* input, float* output, unsigned size) noexcept
 {
-    // The sentinel is the input here
-    const auto sentinel = input + size;
-
     if (getSIMDOpStatus(SIMDOps::copy)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && input < lastAligned)
-                *output++ = *input++;
-
-            while (input < lastAligned) {
-                _mm_store_ps(output, _mm_load_ps(input));
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return copySSE(input, output, size);
     }
-
-    std::copy(input, sentinel, output);
+    std::copy(input, input + size, output);
 }
 
 template <>
 float mean<float>(const float* vector, unsigned size) noexcept
 {
-    const auto sentinel = vector + size;
-
-    float result { 0.0f };
-    if (size == 0)
-        return result;
-
     if (getSIMDOpStatus(SIMDOps::mean)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(vector) && vector < lastAligned)
-                result += *vector++;
-
-            auto mmSums = _mm_setzero_ps();
-            while (vector < lastAligned) {
-                mmSums = _mm_add_ps(mmSums, _mm_load_ps(vector));
-                incrementAll<4>(vector);
-            }
-
-            std::array<float, 4> sseResult;
-            _mm_store_ps(sseResult.data(), mmSums);
-
-            for (auto sseValue : sseResult)
-                result += sseValue;
-
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
+        if (cpuInfo.has_sse())
+            return meanSSE(vector, size);
     }
-
-    while (vector < sentinel)
-        result += *vector++;
-
-    return result / static_cast<float>(size);
+    return meanScalar(vector, size);
 }
 
 template <>
 float meanSquared<float>(const float* vector, unsigned size) noexcept
 {
-    const auto sentinel = vector + size;
-
-    float result { 0.0f };
-    if (size == 0)
-        return result;
-
     if (getSIMDOpStatus(SIMDOps::meanSquared)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-            while (unaligned(vector) && vector < lastAligned){
-        result += (*vector) * (*vector);
-        vector++;
+        if (cpuInfo.has_sse())
+            return meanSquaredSSE(vector, size);
     }
-
-            auto mmSums = _mm_setzero_ps();
-            while (vector < lastAligned) {
-                const auto mmValues = _mm_load_ps(vector);
-                mmSums = _mm_add_ps(mmSums, _mm_mul_ps(mmValues, mmValues));
-                incrementAll<4>(vector);
-            }
-
-            std::array<float, 4> sseResult;
-            _mm_store_ps(sseResult.data(), mmSums);
-
-            for (auto sseValue : sseResult)
-                result += sseValue;
-
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
-    }
-
-    while (vector < sentinel) {
-        result += (*vector) * (*vector);
-        vector++;
-    }
-
-    return result / static_cast<float>(size);
+    return meanSquaredScalar(vector, size);
 }
 
 template <>
 void cumsum<float>(const float* input, float* output, unsigned size) noexcept
 {
-    if (size == 0)
-        return;
-
-    const auto sentinel = output + size;
-    *output++ = *input++;
-
     if (getSIMDOpStatus(SIMDOps::cumsum)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned) {
-                *output = *(output - 1) + *input;
-                incrementAll(input, output);
-            }
-
-            auto mmOutput = _mm_set_ps1(*(output - 1));
-            while (output < lastAligned) {
-                auto mmOffset = _mm_load_ps(input);
-                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, _MM_SHUFFLE(1, 0, 0, 0)));
-                mmOutput = _mm_add_ps(mmOutput, mmOffset);
-                _mm_store_ps(output, mmOutput);
-                mmOutput = _mm_shuffle_ps(mmOutput, mmOutput, _MM_SHUFFLE(3, 3, 3, 3));
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
-    }
-
-    while (output < sentinel) {
-        *output = *(output - 1) + *input;
-        incrementAll(input, output);
+        if (cpuInfo.has_sse())
+            return cumsumSSE(input, output, size);
     }
+    return cumsumScalar(input, output, size);
 }
 
 template <>
 void diff<float>(const float* input, float* output, unsigned size) noexcept
 {
-    if (size == 0)
-        return;
-
-    const auto sentinel = output + size;
-    *output++ = *input++;
-
     if (getSIMDOpStatus(SIMDOps::diff)) {
-#if SFIZZ_CPU_FAMILY_X86_64 || SFIZZ_CPU_FAMILY_I386
-        if (cpuInfo.has_sse()) {
-            const auto* lastAligned = prevAligned(sentinel);
-
-            while (unaligned(input, output) && output < lastAligned) {
-                *output = *input - *(input - 1);
-                incrementAll(input, output);
-            }
-
-            auto mmBase = _mm_set_ps1(*(input - 1));
-            while (output < lastAligned) {
-                auto mmOutput = _mm_load_ps(input);
-                auto mmNextBase = _mm_shuffle_ps(mmOutput, mmOutput, _MM_SHUFFLE(3, 3, 3, 3));
-                mmOutput = _mm_sub_ps(mmOutput, mmBase);
-                mmBase = mmNextBase;
-                mmOutput = _mm_sub_ps(mmOutput, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(mmOutput), 4)));
-                _mm_store_ps(output, mmOutput);
-                incrementAll<4>(input, output);
-            }
-            // fallthrough from lastAligned to sentinel
-        }
-#endif
-    }
-
-    while (output < sentinel) {
-        *output = *input - *(input - 1);
-        incrementAll(input, output);
+        if (cpuInfo.has_sse())
+            return diffSSE(input, output, size);
     }
+    return diffScalar(input, output, size);
 }
 
 }
diff --git a/src/sfizz/SIMDHelpers.h b/src/sfizz/SIMDHelpers.h
index 2ca4375f..9e3bd06b 100644
--- a/src/sfizz/SIMDHelpers.h
+++ b/src/sfizz/SIMDHelpers.h
@@ -27,7 +27,7 @@
 #include "Config.h"
 #include "Debug.h"
 #include "MathHelpers.h"
-#include <absl/algorithm/container.h>
+#include "simd/HelpersScalar.h"
 #include <absl/types/span.h>
 #include <array>
 #include <cmath>
@@ -62,32 +62,6 @@ enum class SIMDOps {
 void setSIMDOpStatus(SIMDOps op, bool status);
 bool getSIMDOpStatus(SIMDOps op);
 
-constexpr uintptr_t ByteAlignmentMask(unsigned N) { return N - 1; }
-
-template<unsigned N = config::defaultAlignment, class T>
-T* nextAligned(const T* ptr)
-{
-    return reinterpret_cast<T*>(reinterpret_cast<uintptr_t>(ptr) + ByteAlignmentMask(N) & (~ByteAlignmentMask(N)));
-}
-
-template<unsigned N = config::defaultAlignment, class T>
-T* prevAligned(const T* ptr)
-{
-    return reinterpret_cast<T*>(reinterpret_cast<uintptr_t>(ptr) & (~ByteAlignmentMask(N)));
-}
-
-template<unsigned N = config::defaultAlignment, class T>
-bool unaligned(const T* ptr)
-{
-    return (reinterpret_cast<uintptr_t>(ptr) & ByteAlignmentMask(N) )!= 0;
-}
-
-template<unsigned N = config::defaultAlignment, class T, class... Args>
-bool unaligned(const T* ptr1, Args... rest)
-{
-    return unaligned<N>(ptr1) || unaligned<N>(rest...);
-}
-
 /**
  * @brief Read interleaved stereo data from a buffer and separate it in a left/right pair of buffers.
  *
@@ -133,10 +107,16 @@ inline void writeInterleaved(absl::Span<const float> inputLeft, absl::Span<const
  * @param output
  * @param value
  */
+template <class T>
+void fill(T* output, T value, unsigned size) noexcept
+{
+    std::fill(output, output + size, value);
+}
+
 template <class T>
 void fill(absl::Span<T> output, T value) noexcept
 {
-    absl::c_fill(output, value);
+    fill<T>(output.data(), value, output.size());
 }
 
 /**
@@ -150,9 +130,7 @@ void fill(absl::Span<T> output, T value) noexcept
 template<class T>
 void applyGain(T gain, const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ = gain * (*input++);
+    applyGainScalar(gain, input, output, size);
 }
 
 template<>
@@ -195,9 +173,7 @@ inline void applyGain(float gain, absl::Span<float> array) noexcept
 template<class T>
 void applyGain(const T* gain, const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ = (*gain++) * (*input++);
+    applyGainScalar(gain, input, output, size);
 }
 
 template<>
@@ -244,9 +220,7 @@ inline void applyGain(absl::Span<const T> gain, absl::Span<T> array) noexcept
 template <class T>
 void divide(const T* input, const T* divisor, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ = (*input++) / (*divisor++);
+    divideScalar(input, divisor, output, size);
 }
 
 template <>
@@ -285,9 +259,7 @@ void divide(absl::Span<T> output,  absl::Span<const T> divisor) noexcept
 template <class T>
 void multiplyAdd(const T* gain, const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ += (*gain++) * (*input++);
+    multiplyAddScalar(gain, input, output, size);
 }
 
 template <>
@@ -312,9 +284,7 @@ void multiplyAdd(absl::Span<const T> gain, absl::Span<const T> input, absl::Span
 template <class T>
 void multiplyAdd(T gain, const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ += gain * (*input++);
+    multiplyAddScalar(gain, input, output, size);
 }
 
 template <>
@@ -340,12 +310,7 @@ void multiplyAdd(T gain, absl::Span<const T> input, absl::Span<T> output) noexce
 template <class T>
 T linearRamp(T* output, T start, T step, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel) {
-        *output++ = start;
-        start += step;
-    }
-    return start;
+    linearRampScalar(output, start, step, size);
 }
 
 template <>
@@ -369,12 +334,7 @@ T linearRamp(absl::Span<T> output, T start, T step) noexcept
 template <class T>
 T multiplicativeRamp(T* output, T start, T step, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel) {
-        *output++ = start;
-        start *= step;
-    }
-    return start;
+    multiplicativeRampScalar(output, start, step, size);
 }
 
 template <>
@@ -398,9 +358,7 @@ T multiplicativeRamp(absl::Span<T> output, T start, T step) noexcept
 template <class T>
 void add(const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ += *input++;
+    addScalar(input, output, size);
 }
 
 template <>
@@ -424,9 +382,7 @@ void add(absl::Span<const T> input, absl::Span<T> output) noexcept
 template <class T>
 void add(T value, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ += value;
+    addScalar(value, output, size);
 }
 
 template <>
@@ -449,9 +405,7 @@ void add(T value, absl::Span<T> output) noexcept
 template <class T>
 void subtract(const T* input, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ -= *input++;
+    subtractScalar(input, output, size);
 }
 
 template <>
@@ -475,9 +429,7 @@ void subtract(absl::Span<const T> input, absl::Span<T> output) noexcept
 template <class T>
 void subtract(T value, T* output, unsigned size) noexcept
 {
-    const auto sentinel = output + size;
-    while (output < sentinel)
-        *output++ -= value;
+    subtractScalar(value, output, size);
 }
 
 template <>
@@ -525,15 +477,7 @@ void copy(absl::Span<const T> input, absl::Span<T> output) noexcept
 template <class T>
 T mean(const T* vector, unsigned size) noexcept
 {
-    T result{ 0.0 };
-    if (size == 0)
-        return result;
-
-    const auto sentinel = vector + size;
-    while (vector < sentinel)
-        result += *vector++;
-
-    return result / static_cast<T>(size);
+    meanScalar(vector, size);
 }
 
 template <>
@@ -556,17 +500,7 @@ T mean(absl::Span<const T> vector) noexcept
 template <class T>
 T meanSquared(const T* vector, unsigned size) noexcept
 {
-    T result{ 0.0 };
-    if (size == 0)
-        return result;
-
-    const auto sentinel = vector + size;
-    while (vector < sentinel) {
-        result += (*vector) * (*vector);
-        vector++;
-    }
-
-    return result / static_cast<T>(size);
+    meanSquaredScalar(vector, size);
 }
 
 template <>
@@ -590,16 +524,7 @@ T meanSquared(absl::Span<const T> vector) noexcept
 template <class T>
 void cumsum(const T* input, T* output, unsigned size) noexcept
 {
-    if (size == 0)
-        return;
-
-    const auto sentinel = output + size;
-
-    *output++ = *input++;
-    while (output < sentinel) {
-        *output = *(output - 1) + *input;
-        incrementAll(input, output);
-    }
+    cumsumScalar(input, output, size);
 }
 
 template <>
@@ -661,16 +586,7 @@ void sfzInterpolationCast(absl::Span<const T> floatJumps, absl::Span<int> jumps,
 template <class T>
 void diff(const T* input, T* output, unsigned size) noexcept
 {
-    if (size == 0)
-        return;
-
-    const auto sentinel = output + size;
-
-    *output++ = *input++;
-    while (output < sentinel) {
-        *output = *input - *(input - 1);
-        incrementAll(input, output);
-    }
+    diffScalar(input, output, size);
 }
 
 template <>
diff --git a/src/sfizz/simd/Common.h b/src/sfizz/simd/Common.h
new file mode 100644
index 00000000..6fdddeab
--- /dev/null
+++ b/src/sfizz/simd/Common.h
@@ -0,0 +1,34 @@
+// 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
+#include <stdint.h>
+
+constexpr uintptr_t ByteAlignmentMask(unsigned N) { return N - 1; }
+
+template<unsigned N, class T>
+T* nextAligned(const T* ptr)
+{
+    return reinterpret_cast<T*>(reinterpret_cast<uintptr_t>(ptr) + ByteAlignmentMask(N) & (~ByteAlignmentMask(N)));
+}
+
+template<unsigned N, class T>
+T* prevAligned(const T* ptr)
+{
+    return reinterpret_cast<T*>(reinterpret_cast<uintptr_t>(ptr) & (~ByteAlignmentMask(N)));
+}
+
+template<unsigned N, class T>
+bool unaligned(const T* ptr)
+{
+    return (reinterpret_cast<uintptr_t>(ptr) & ByteAlignmentMask(N) )!= 0;
+}
+
+template<unsigned N, class T, class... Args>
+bool unaligned(const T* ptr1, Args... rest)
+{
+    return unaligned<N>(ptr1) || unaligned<N>(rest...);
+}
diff --git a/src/sfizz/simd/HelpersAVX.cpp b/src/sfizz/simd/HelpersAVX.cpp
new file mode 100644
index 00000000..9ca89a65
--- /dev/null
+++ b/src/sfizz/simd/HelpersAVX.cpp
@@ -0,0 +1,56 @@
+// 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 "HelpersAVX.h"
+#include "SIMDConfig.h"
+#include "../MathHelpers.h"
+#include "Common.h"
+
+#ifdef SFIZZ_HAVE_AVX
+#include "immintrin.h"
+using Type = float;
+constexpr unsigned TypeAlignment = 8;
+constexpr unsigned ByteAlignment = TypeAlignment * sizeof(Type);
+#endif
+
+void applyGainAVX(float gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_AVX
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    const auto mmGain = _mm256_set1_ps(gain);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ = gain * (*input++);
+
+    while (output < lastAligned) {
+        _mm256_store_ps(output, _mm256_mul_ps(mmGain, _mm256_load_ps(input)));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ = gain * (*input++);
+}
+
+void applyGainAVX(const float* gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_AVX
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ = (*gain++) * (*input++);
+
+    while (output < lastAligned) {
+        _mm256_store_ps(output, _mm256_mul_ps(_mm256_load_ps(gain), _mm256_load_ps(input)));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ = (*gain++) * (*input++);
+}
diff --git a/src/sfizz/simd/HelpersAVX.h b/src/sfizz/simd/HelpersAVX.h
new file mode 100644
index 00000000..7eb8c2bf
--- /dev/null
+++ b/src/sfizz/simd/HelpersAVX.h
@@ -0,0 +1,10 @@
+// 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
+
+void applyGainAVX(float gain, const float* input, float* output, unsigned size) noexcept;
+void applyGainAVX(const float* gain, const float* input, float* output, unsigned size) noexcept;
diff --git a/src/sfizz/simd/HelpersSSE.cpp b/src/sfizz/simd/HelpersSSE.cpp
new file mode 100644
index 00000000..8275280a
--- /dev/null
+++ b/src/sfizz/simd/HelpersSSE.cpp
@@ -0,0 +1,471 @@
+// 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 "HelpersSSE.h"
+#include "../MathHelpers.h"
+#include "../SIMDConfig.h"
+#include "Common.h"
+
+#ifdef SFIZZ_HAVE_SSE
+#include "emmintrin.h"
+using Type = float;
+constexpr unsigned TypeAlignment = 4;
+constexpr unsigned ByteAlignment = TypeAlignment * sizeof(Type);
+#endif
+
+void readInterleavedSSE(const float* input, float* outputLeft, float* outputRight, unsigned inputSize) noexcept
+{
+    const auto sentinel = input + inputSize - 1;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(input + inputSize - TypeAlignment);
+    while (unaligned<ByteAlignment>(input, outputLeft, outputRight) && input < lastAligned) {
+        *outputLeft++ = *input++;
+        *outputRight++ = *input++;
+    }
+
+    while (input < lastAligned) {
+        auto register0 = _mm_load_ps(input);
+        auto register1 = _mm_load_ps(input + TypeAlignment);
+        auto register2 = register0;
+        // register 2 holds the copy of register 0 that is going to get erased by the first operation
+        // Remember that the bit mask reads from the end; 10 00 10 00 means
+        // "take 0 from a, take 2 from a, take 0 from b, take 2 from b"
+        register0 = _mm_shuffle_ps(register0, register1, 0b10001000);
+        register1 = _mm_shuffle_ps(register2, register1, 0b11011101);
+        _mm_store_ps(outputLeft, register0);
+        _mm_store_ps(outputRight, register1);
+        incrementAll<TypeAlignment>(input, input, outputLeft, outputRight);
+    }
+
+    // NEON wip
+    // auto reg = vld2q_f32(in);
+    // vst1q_f32(lOut, reg.val[0]);
+    // vst1q_f32(rOut, reg.val[1]);
+#endif
+    while (input < sentinel) {
+        *outputLeft++ = *input++;
+        *outputRight++ = *input++;
+    }
+}
+
+void writeInterleavedSSE(const float* inputLeft, const float* inputRight, float* output, unsigned outputSize) noexcept
+{
+    const auto sentinel = output + outputSize - 1;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(output + outputSize - TypeAlignment);
+    while (unaligned<ByteAlignment>(output, inputRight, inputLeft) && output < lastAligned) {
+        *output++ = *inputLeft++;
+        *output++ = *inputRight++;
+    }
+
+    while (output < lastAligned) {
+        const auto lInRegister = _mm_load_ps(inputLeft);
+        const auto rInRegister = _mm_load_ps(inputRight);
+        const auto outRegister1 = _mm_unpacklo_ps(lInRegister, rInRegister);
+        _mm_store_ps(output, outRegister1);
+        const auto outRegister2 = _mm_unpackhi_ps(lInRegister, rInRegister);
+        _mm_store_ps(output + 4, outRegister2);
+        incrementAll<TypeAlignment>(output, output, inputLeft, inputRight);
+    }
+#endif
+
+    while (output < sentinel) {
+        *output++ = *inputLeft++;
+        *output++ = *inputRight++;
+    }
+}
+
+void applyGainSSE(float gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    const auto mmGain = _mm_set1_ps(gain);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ = gain * (*input++);
+
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_mul_ps(mmGain, _mm_load_ps(input)));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ = gain * (*input++);
+}
+
+void applyGainSSE(const float* gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ = (*gain++) * (*input++);
+
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_mul_ps(_mm_load_ps(gain), _mm_load_ps(input)));
+        incrementAll<TypeAlignment>(gain, input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ = (*gain++) * (*input++);
+}
+
+void divideSSE(const float* input, const float* divisor, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ = (*input++) / (*divisor++);
+
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_div_ps(_mm_load_ps(input), _mm_load_ps(divisor)));
+        incrementAll<TypeAlignment>(divisor, input, output);
+    }
+
+    while (output < sentinel)
+        *output++ = (*input++) / (*divisor++);
+}
+
+void multiplyAddSSE(const float* gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ += (*gain++) * (*input++);
+
+    while (output < lastAligned) {
+        auto mmOut = _mm_load_ps(output);
+        mmOut = _mm_add_ps(_mm_mul_ps(_mm_load_ps(gain), _mm_load_ps(input)), mmOut);
+        _mm_store_ps(output, mmOut);
+        incrementAll<TypeAlignment>(gain, input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ += (*gain++) * (*input++);
+}
+
+void multiplyAddSSE(float gain, const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ += gain * (*input++);
+
+    auto mmGain = _mm_set1_ps(gain);
+    while (output < lastAligned) {
+        auto mmOut = _mm_load_ps(output);
+        mmOut = _mm_add_ps(_mm_mul_ps(mmGain, _mm_load_ps(input)), mmOut);
+        _mm_store_ps(output, mmOut);
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ += gain * (*input++);
+}
+
+float linearRampSSE(float* output, float start, float step, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(output) && output < lastAligned) {
+        *output++ = start;
+        start += step;
+    }
+
+    auto mmStart = _mm_set1_ps(start - step);
+    auto mmStep = _mm_set_ps(step + step + step + step, step + step + step, step + step, step);
+    while (output < lastAligned) {
+        mmStart = _mm_add_ps(mmStart, mmStep);
+        _mm_store_ps(output, mmStart);
+        mmStart = _mm_shuffle_ps(mmStart, mmStart, _MM_SHUFFLE(3, 3, 3, 3));
+        incrementAll<TypeAlignment>(output);
+    }
+    start = _mm_cvtss_f32(mmStart) + step;
+#endif
+
+    while (output < sentinel) {
+        *output++ = start;
+        start += step;
+    }
+    return start;
+}
+
+float multiplicativeRampSSE(float* output, float start, float step, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(output) && output < lastAligned) {
+        *output++ = start;
+        start *= step;
+    }
+
+    auto mmStart = _mm_set1_ps(start / step);
+    auto mmStep = _mm_set_ps(step * step * step * step, step * step * step, step * step, step);
+    while (output < lastAligned) {
+        mmStart = _mm_mul_ps(mmStart, mmStep);
+        _mm_store_ps(output, mmStart);
+        mmStart = _mm_shuffle_ps(mmStart, mmStart, _MM_SHUFFLE(3, 3, 3, 3));
+        incrementAll<TypeAlignment>(output);
+    }
+    start = _mm_cvtss_f32(mmStart) * step;
+#endif
+
+    while (output < sentinel) {
+        *output++ = start;
+        start *= step;
+    }
+    return start;
+}
+
+void addSSE(const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ += *input++;
+
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_add_ps(_mm_load_ps(output), _mm_load_ps(input)));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ += *input++;
+}
+
+void addSSE(float value, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(output) && output < lastAligned)
+        *output++ += value;
+
+    const auto mmValue = _mm_set1_ps(value);
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_add_ps(_mm_load_ps(output), mmValue));
+        incrementAll<TypeAlignment>(output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ += value;
+}
+
+void subtractSSE(const float* input, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned)
+        *output++ -= *input++;
+
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_sub_ps(_mm_load_ps(output), _mm_load_ps(input)));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ -= *input++;
+}
+
+void subtractSSE(float value, float* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(output) && output < lastAligned)
+        *output++ -= value;
+
+    const auto mmValue = _mm_set1_ps(value);
+    while (output < lastAligned) {
+        _mm_store_ps(output, _mm_sub_ps(_mm_load_ps(output), mmValue));
+        incrementAll<TypeAlignment>(output);
+    }
+#endif
+
+    while (output < sentinel)
+        *output++ -= value;
+}
+
+void copySSE(const float* input, float* output, unsigned size) noexcept
+{
+    // The sentinel is the input here
+    const auto sentinel = input + size;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && input < lastAligned)
+        *output++ = *input++;
+
+    while (input < lastAligned) {
+        _mm_store_ps(output, _mm_load_ps(input));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    std::copy(input, sentinel, output);
+}
+
+float meanSSE(const float* vector, unsigned size) noexcept
+{
+    const auto sentinel = vector + size;
+
+    float result { 0.0f };
+    if (size == 0)
+        return result;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(vector) && vector < lastAligned)
+        result += *vector++;
+
+    auto mmSums = _mm_setzero_ps();
+    while (vector < lastAligned) {
+        mmSums = _mm_add_ps(mmSums, _mm_load_ps(vector));
+        incrementAll<TypeAlignment>(vector);
+    }
+
+    std::array<float, 4> sseResult;
+    _mm_store_ps(sseResult.data(), mmSums);
+
+    for (auto sseValue : sseResult)
+        result += sseValue;
+#endif
+
+    while (vector < sentinel)
+        result += *vector++;
+
+    return result / static_cast<float>(size);
+}
+
+float meanSquaredSSE(const float* vector, unsigned size) noexcept
+{
+    const auto sentinel = vector + size;
+
+    float result { 0.0f };
+    if (size == 0)
+        return result;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(vector) && vector < lastAligned) {
+        result += (*vector) * (*vector);
+        vector++;
+    }
+
+    auto mmSums = _mm_setzero_ps();
+    while (vector < lastAligned) {
+        const auto mmValues = _mm_load_ps(vector);
+        mmSums = _mm_add_ps(mmSums, _mm_mul_ps(mmValues, mmValues));
+        incrementAll<TypeAlignment>(vector);
+    }
+
+    std::array<float, 4> sseResult;
+    _mm_store_ps(sseResult.data(), mmSums);
+
+    for (auto sseValue : sseResult)
+        result += sseValue;
+#endif
+
+    while (vector < sentinel) {
+        result += (*vector) * (*vector);
+        vector++;
+    }
+
+    return result / static_cast<float>(size);
+}
+
+void cumsumSSE(const float* input, float* output, unsigned size) noexcept
+{
+    if (size == 0)
+        return;
+
+    const auto sentinel = output + size;
+    *output++ = *input++;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned) {
+        *output = *(output - 1) + *input;
+        incrementAll(input, output);
+    }
+
+    auto mmOutput = _mm_set_ps1(*(output - 1));
+    while (output < lastAligned) {
+        auto mmOffset = _mm_load_ps(input);
+        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, _MM_SHUFFLE(1, 0, 0, 0)));
+        mmOutput = _mm_add_ps(mmOutput, mmOffset);
+        _mm_store_ps(output, mmOutput);
+        mmOutput = _mm_shuffle_ps(mmOutput, mmOutput, _MM_SHUFFLE(3, 3, 3, 3));
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel) {
+        *output = *(output - 1) + *input;
+        incrementAll(input, output);
+    }
+}
+
+void diffSSE(const float* input, float* output, unsigned size) noexcept
+{
+    if (size == 0)
+        return;
+
+    const auto sentinel = output + size;
+    *output++ = *input++;
+
+#ifdef SFIZZ_HAVE_SSE
+    const auto* lastAligned = prevAligned<ByteAlignment>(sentinel);
+    while (unaligned<ByteAlignment>(input, output) && output < lastAligned) {
+        *output = *input - *(input - 1);
+        incrementAll(input, output);
+    }
+
+    auto mmBase = _mm_set_ps1(*(input - 1));
+    while (output < lastAligned) {
+        auto mmOutput = _mm_load_ps(input);
+        auto mmNextBase = _mm_shuffle_ps(mmOutput, mmOutput, _MM_SHUFFLE(3, 3, 3, 3));
+        mmOutput = _mm_sub_ps(mmOutput, mmBase);
+        mmBase = mmNextBase;
+        mmOutput = _mm_sub_ps(mmOutput, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(mmOutput), 4)));
+        _mm_store_ps(output, mmOutput);
+        incrementAll<TypeAlignment>(input, output);
+    }
+#endif
+
+    while (output < sentinel) {
+        *output = *input - *(input - 1);
+        incrementAll(input, output);
+    }
+}
diff --git a/src/sfizz/simd/HelpersSSE.h b/src/sfizz/simd/HelpersSSE.h
new file mode 100644
index 00000000..9c43e342
--- /dev/null
+++ b/src/sfizz/simd/HelpersSSE.h
@@ -0,0 +1,27 @@
+// 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
+
+/* These are the SSE versions of the SIMDHelpers */
+void readInterleavedSSE(const float* input, float* outputLeft, float* outputRight, unsigned inputSize) noexcept;
+void writeInterleavedSSE(const float* inputLeft, const float* inputRight, float* output, unsigned outputSize) noexcept;
+void applyGainSSE(float gain, const float* input, float* output, unsigned size) noexcept;
+void applyGainSSE(const float* gain, const float* input, float* output, unsigned size) noexcept;
+void divideSSE(const float* input, const float* divisor, float* output, unsigned size) noexcept;
+void multiplyAddSSE(const float* gain, const float* input, float* output, unsigned size) noexcept;
+void multiplyAddSSE(float gain, const float* input, float* output, unsigned size) noexcept;
+float linearRampSSE(float* output, float start, float step, unsigned size) noexcept;
+float multiplicativeRampSSE(float* output, float start, float step, unsigned size) noexcept;
+void addSSE(const float* input, float* output, unsigned size) noexcept;
+void addSSE(float value, float* output, unsigned size) noexcept;
+void subtractSSE(const float* input, float* output, unsigned size) noexcept;
+void subtractSSE(float value, float* output, unsigned size) noexcept;
+void copySSE(const float* input, float* output, unsigned size) noexcept;
+float meanSSE(const float* vector, unsigned size) noexcept;
+float meanSquaredSSE(const float* vector, unsigned size) noexcept;
+void cumsumSSE(const float* input, float* output, unsigned size) noexcept;
+void diffSSE(const float* input, float* output, unsigned size) noexcept;
diff --git a/src/sfizz/simd/HelpersScalar.h b/src/sfizz/simd/HelpersScalar.h
new file mode 100644
index 00000000..1554e636
--- /dev/null
+++ b/src/sfizz/simd/HelpersScalar.h
@@ -0,0 +1,181 @@
+// 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
+
+template<class T>
+inline void readInterleavedScalar(const T* input, T* outputLeft, T* outputRight, unsigned inputSize)
+{
+    const auto sentinel = input + inputSize - 1;
+    while (input < sentinel) {
+        *outputLeft++ = *input++;
+        *outputRight++ = *input++;
+    }
+}
+
+template<class T>
+inline void writeInterleavedScalar(const T* inputLeft, const T* inputRight, T* output, unsigned outputSize)
+{
+    const auto sentinel = output + outputSize - 1;
+    while (output < sentinel) {
+        *output++ = *inputLeft++;
+        *output++ = *inputRight++;
+    }
+}
+
+template<class T>
+inline void applyGainScalar(T gain, const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ = gain * (*input++);
+}
+
+template<class T>
+inline void applyGainScalar(const T* gain, const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ = (*gain++) * (*input++);
+}
+
+template <class T>
+inline void divideScalar(const T* input, const T* divisor, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ = (*input++) / (*divisor++);
+}
+
+template <class T>
+inline void multiplyAddScalar(const T* gain, const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ += (*gain++) * (*input++);
+}
+
+template <class T>
+inline void multiplyAddScalar(T gain, const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ += gain * (*input++);
+}
+
+template <class T>
+T linearRampScalar(T* output, T start, T step, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel) {
+        *output++ = start;
+        start += step;
+    }
+    return start;
+}
+
+template <class T>
+T multiplicativeRampScalar(T* output, T start, T step, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel) {
+        *output++ = start;
+        start *= step;
+    }
+    return start;
+}
+
+template <class T>
+inline void addScalar(const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ += *input++;
+}
+
+template <class T>
+inline void addScalar(T value, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ += value;
+}
+
+template <class T>
+inline void subtractScalar(const T* input, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ -= *input++;
+}
+
+template <class T>
+inline void subtractScalar(T value, T* output, unsigned size) noexcept
+{
+    const auto sentinel = output + size;
+    while (output < sentinel)
+        *output++ -= value;
+}
+
+template <class T>
+T meanScalar(const T* vector, unsigned size) noexcept
+{
+    T result{ 0.0 };
+    if (size == 0)
+        return result;
+
+    const auto sentinel = vector + size;
+    while (vector < sentinel)
+        result += *vector++;
+
+    return result / static_cast<T>(size);
+}
+
+template <class T>
+T meanSquaredScalar(const T* vector, unsigned size) noexcept
+{
+    T result{ 0.0 };
+    if (size == 0)
+        return result;
+
+    const auto sentinel = vector + size;
+    while (vector < sentinel) {
+        result += (*vector) * (*vector);
+        vector++;
+    }
+
+    return result / static_cast<T>(size);
+}
+
+template <class T>
+void cumsumScalar(const T* input, T* output, unsigned size) noexcept
+{
+    if (size == 0)
+        return;
+
+    const auto sentinel = output + size;
+
+    *output++ = *input++;
+    while (output < sentinel) {
+        *output = *(output - 1) + *input;
+        incrementAll(input, output);
+    }
+}
+
+template <class T>
+void diffScalar(const T* input, T* output, unsigned size) noexcept
+{
+    if (size == 0)
+        return;
+
+    const auto sentinel = output + size;
+
+    *output++ = *input++;
+    while (output < sentinel) {
+        *output = *input - *(input - 1);
+        incrementAll(input, output);
+    }
+}