ambivalence/Source/DSP/OutputEQ.h

#pragma once

#include <cmath>
#include <algorithm>

namespace FDNReverb {

     // -----------------------------------------------------------------------------
     // OutputEQ: Wet output stage Lo/Hi Cut (Linkwitz-Riley 12dB/oct)
     // -----------------------------------------------------------------------------
     // design rationale:
     // - 1 IIR (6dB/oct) x 2 cascade = 12dB/oct
     // - Linkwitz-Riley topology: 2nd-order phase alignment
     // - keeps the reverb sounding musical
    //
     // filter equation (1 IIR):
     // HPF: y[n] = R . (y[n-1] + x[n] - x[n-1])
     // LPF: y[n] = (1 - R) . x[n] + R . y[n-1]
     // where R = exp(-2pi.fc/fs)
    //
     // real-time safety :
     // - no allocation at all
     // - per-sample cost: HPF 8 ops + LPF 6 ops (L/R combined)
     // - coefficients updated per block (no zipper noise, no SmoothedValue needed)
    //
     // bypass :
     // - Lo Cut below 20 Hz -> HPF fully bypassed
     // - Hi Cut above 20 kHz -> LPF fully bypassed
     // both bypasses are per-block coefficient updates, so CPU use is trivial.
     // -----------------------------------------------------------------------------
    class OutputEQ {
    public:
        OutputEQ() = default;

        void prepare(double sampleRate) noexcept {
            fs = sampleRate;
            reset();
            setLoCutHz(20.0f);
            setHiCutHz(20000.0f);
        }

        void reset() noexcept {
             // HPF state (two stages per channel, L/R)
            hpfX1_L_1 = hpfY1_L_1 = 0.0f;
            hpfX1_L_2 = hpfY1_L_2 = 0.0f;
            hpfX1_R_1 = hpfY1_R_1 = 0.0f;
            hpfX1_R_2 = hpfY1_R_2 = 0.0f;
             // LPF state (two stages per channel, L/R)
            lpfY1_L_1 = 0.0f;
            lpfY1_L_2 = 0.0f;
            lpfY1_R_1 = 0.0f;
            lpfY1_R_2 = 0.0f;
        }

         // --- parameter setters (called per block) ---
        void setLoCutHz(float fcHz) noexcept {
            currentLoCutHz = fcHz;
             // bypass below 20 Hz (skip R computation)
            if (fcHz <= 20.0f) {
                loCutActive = false;
                return;
            }
            loCutActive = true;
            constexpr float twoPi = 6.28318530718f;
            const float clamped = std::clamp(fcHz, 20.0f, 500.0f);
            loCutR = std::exp(-twoPi * clamped / static_cast<float>(fs));
        }

        void setHiCutHz(float fcHz) noexcept {
            currentHiCutHz = fcHz;
             // bypass above 20 kHz
            const float nyquist = static_cast<float>(fs) * 0.45f;
            const float clamped = std::clamp(fcHz, 1000.0f, std::min(20000.0f, nyquist));
            if (fcHz >= 20000.0f) {
                hiCutActive = false;
                return;
            }
            hiCutActive = true;
            constexpr float twoPi = 6.28318530718f;
            hiCutR = std::exp(-twoPi * clamped / static_cast<float>(fs));
        }

         // --- per-sample processing (L/R interleaved) ---
        inline void process(float& l, float& r) noexcept {
             // -- Lo Cut: 1 HPF x 2 cascade --
            if (loCutActive) {
                 // L stage 1
                const float l_in = l;
                const float l_1 = loCutR * (hpfY1_L_1 + l_in - hpfX1_L_1);
                hpfX1_L_1 = l_in;
                hpfY1_L_1 = l_1;
                 // L stage 2
                const float l_2 = loCutR * (hpfY1_L_2 + l_1 - hpfX1_L_2);
                hpfX1_L_2 = l_1;
                hpfY1_L_2 = l_2;
                l = l_2;

                 // R stage 1
                const float r_in = r;
                const float r_1 = loCutR * (hpfY1_R_1 + r_in - hpfX1_R_1);
                hpfX1_R_1 = r_in;
                hpfY1_R_1 = r_1;
                 // R stage 2
                const float r_2 = loCutR * (hpfY1_R_2 + r_1 - hpfX1_R_2);
                hpfX1_R_2 = r_1;
                hpfY1_R_2 = r_2;
                r = r_2;
            }

             // -- Hi Cut: 1 LPF x 2 cascade --
            if (hiCutActive) {
                const float oneMinusR = 1.0f - hiCutR;
                 // L stage 1
                lpfY1_L_1 = oneMinusR * l + hiCutR * lpfY1_L_1;
                 // L stage 2
                lpfY1_L_2 = oneMinusR * lpfY1_L_1 + hiCutR * lpfY1_L_2;
                l = lpfY1_L_2;

                 // R stage 1
                lpfY1_R_1 = oneMinusR * r + hiCutR * lpfY1_R_1;
                 // R stage 2
                lpfY1_R_2 = oneMinusR * lpfY1_R_1 + hiCutR * lpfY1_R_2;
                r = lpfY1_R_2;
            }
        }

        float getCurrentLoCutHz() const noexcept { return currentLoCutHz; }
        float getCurrentHiCutHz() const noexcept { return currentHiCutHz; }

    private:
        double fs{ 48000.0 };

         // -- Lo Cut (HPF) --
        bool  loCutActive{ false };
        float loCutR{ 0.0f };
        float currentLoCutHz{ 20.0f };
        float hpfX1_L_1{}, hpfY1_L_1{}, hpfX1_L_2{}, hpfY1_L_2{};
        float hpfX1_R_1{}, hpfY1_R_1{}, hpfX1_R_2{}, hpfY1_R_2{};

         // -- Hi Cut (LPF) --
        bool  hiCutActive{ false };
        float hiCutR{ 0.0f };
        float currentHiCutHz{ 20000.0f };
        float lpfY1_L_1{}, lpfY1_L_2{};
        float lpfY1_R_1{}, lpfY1_R_2{};
    };

} // namespace FDNReverb