ambivalence/Source/DSP/Saturator.h

#pragma once
#include <cmath>
#include <algorithm>

namespace FDNReverb {

    enum class SaturationMode {
        Warm = 0,
        Tape = 1,
        Tube = 2,
        Hard = 3
    };

    class Saturator {
    public:
        Saturator() = default;

        void reset() noexcept {
            prevInput = 0.0f;
            switch (currentMode) {
            case SaturationMode::Warm: prevF = 1.0f; break;
            case SaturationMode::Tape: prevF = 0.0f; break;
            case SaturationMode::Tube: prevF = 1.0f; break;
            case SaturationMode::Hard: prevF = 0.0f; break;
            }
        }

        void setMode(SaturationMode mode) noexcept {
            if (mode != currentMode) {
                currentMode = mode;
                reset();
            }
        }

        void setMode(int modeIndex) noexcept {
            setMode(static_cast<SaturationMode>(std::clamp(modeIndex, 0, 3)));
        }

         // -------------------------------------------------------------------------
         // * Step B fix: only the drive curve changed; ADAA structure fully preserved
         // -------------------------------------------------------------------------
         // drive = 1 + amount^3 x 1.0 ( maximum 2.0) -> 1 + amount^2 x 2.5 ( maximum 3.5)
        //
         // amount | old drive | new drive | effect
         // -------|----------|----------|--------------------
         // 0.30 | 1.027 | 1.225 | + about 7.5dB stronger
         // 0.50 | 1.125 | 1.625 | + about 3.2dB stronger
         // 0.70 | 1.343 | 2.225 | + about 4.4dB stronger
         // 1.00 | 2.000 | 3.500 | + about 4.9dB stronger
        //
         // -> plugin 24 harmonics visualization
         // -------------------------------------------------------------------------
        void setAmount(float amount) noexcept {
            amount = std::clamp(amount, 0.0f, 1.0f);
            currentAmount = amount;

             // * Step B: amount^2 x 2.5 stronger
            drive = 1.0f + amount * amount * 2.5f;
            wetMix = amount * amount * 0.7f;
            dryMix = 1.0f - amount * 0.25f;
        }

        inline float processSample(float input) noexcept {
            if (currentAmount < 1e-4f) return input;

            const float dryInput = input;
            const float driven = input * drive;
            float saturated = 0.0f;

            switch (currentMode) {
            case SaturationMode::Warm: saturated = processWarm(driven); break;
            case SaturationMode::Tape: saturated = processTape(driven); break;
            case SaturationMode::Tube: saturated = processTube(driven); break;
            case SaturationMode::Hard: saturated = processHard(driven); break;
            }

            saturated /= drive;
            return dryInput * dryMix + saturated * wetMix;
        }

    private:
         // --- Warm: Vicanek x/sqrt(1+x^2) + ADAA 1 ---
        inline float processWarm(float x) noexcept {
            const float F_x = std::sqrt(1.0f + x * x);
            const float dx = x - prevInput;
            float y;
            constexpr float kTol = 1e-5f;
            if (std::abs(dx) < kTol) {
                const float xAvg = (x + prevInput) * 0.5f;
                y = xAvg / std::sqrt(1.0f + xAvg * xAvg);
            }
            else {
                y = (F_x - prevF) / dx;
            }
            prevInput = x;
            prevF = F_x;
            return y;
        }

         // --- Tape: Pade x(27+x^2)/(27+9x^2) (ADAA intentional ) ---
        inline float processTape(float x) noexcept {
            if (x > 3.0f) { prevInput = x; return  1.0f; }
            if (x < -3.0f) { prevInput = x; return -1.0f; }
            const float xsq = x * x;
            prevInput = x;
            return x * (27.0f + xsq) / (27.0f + 9.0f * xsq);
        }

         // -------------------------------------------------------------------------
         // Tube: asymmetric ADAA + * Step B: kNeg 1.5 -> 2.0
         // -------------------------------------------------------------------------
         // positive side : f(x) = x/sqrt(1+x^2) F(x) = sqrt(1+x^2)
         // negative side : f(x) = x/sqrt(1+(kNeg.x)^2) F(x) = (1/kNeg^2)sqrt(1+(kNeg.x)^2) + fShift
        //
         // C^1 : x=0 F_pos(0) = F_neg(0) = 1 fShift design
         // F_pos(0) = sqrt1 = 1
         // F_neg(0) = (1/kNeg^2).sqrt1 + fShift = 1
         // -> fShift = 1 - 1/kNeg^2
        //
         // kNeg=2.0 case : fShift = 1 - 0.25 = 0.75
        //
         // kNeg stronger effect :
         // 2 -> waveform asymmetric
         // -> harmonics (2f, 4f) plugin visualization
         // -------------------------------------------------------------------------
        inline float processTube(float x) noexcept {
             // * Step B: kNeg = 1.5f -> 2.0f
            constexpr float kNeg = 2.0f;
            constexpr float kNeg2 = kNeg * kNeg;        // 4.0f
            constexpr float invKneg2 = 1.0f / kNeg2;      // 0.25f
            constexpr float fShift = 1.0f - invKneg2;   // 0.75f

            float F_x;
            if (x >= 0.0f) {
                F_x = std::sqrt(1.0f + x * x);
            }
            else {
                const float kx = kNeg * x;
                F_x = invKneg2 * std::sqrt(1.0f + kx * kx) + fShift;
            }

            const float dx = x - prevInput;
            const bool  signChanged = (x >= 0.0f) != (prevInput >= 0.0f);

            float y;
            constexpr float kTol = 1e-5f;
            if (std::abs(dx) < kTol || signChanged) {
                 // input -> directly
                if (x >= 0.0f) {
                    y = x / std::sqrt(1.0f + x * x);
                }
                else {
                    const float kx = kNeg * x;
                    y = x / std::sqrt(1.0f + kx * kx);
                }
            }
            else {
                y = (F_x - prevF) / dx;
            }

            prevInput = x;
            prevF = F_x;
            return y;
        }

         // --- Hard: clipping + ADAA 1 ---
        inline float processHard(float x) noexcept {
            float F_x;
            if (x > 1.0f) F_x = x - 0.5f;
            else if (x < -1.0f) F_x = -x - 0.5f;
            else                F_x = x * x * 0.5f;

            const float dx = x - prevInput;
            float y;
            constexpr float kTol = 1e-5f;
            if (std::abs(dx) < kTol) {
                y = std::clamp(x, -1.0f, 1.0f);
            }
            else {
                y = (F_x - prevF) / dx;
            }

            prevInput = x;
            prevF = F_x;
            return y;
        }

        float          prevInput{ 0.0f };
        float          prevF{ 1.0f };
        SaturationMode currentMode{ SaturationMode::Warm };
        float          currentAmount{ 0.0f };
        float          drive{ 1.0f };
        float          wetMix{ 0.0f };
        float          dryMix{ 1.0f };
    };

} // namespace FDNReverb