Biquad.h

// DSPark -- Professional Audio DSP Framework
// Copyright (c) 2026 Cristian Moresi -- MIT License
 
#pragma once
 
#include <algorithm>
#include <array>
#include <cmath>
#include <numbers>
#include <span>
#include <atomic>
#include <cassert>
 
#include "AudioBuffer.h"
#include "DenormalGuard.h"
 
// Keeps a cold path out of a hot caller's stack frame. Needed exactly once here
// (see Biquad::adoptStagedCoeffs), where inlining the coefficient adoption into
// the per-sample entry points costs four instructions on EVERY sample. Locally
// defined and guarded, the same way Core/SimdOps.h spells DSPARK_RESTRICT.
#ifndef DSPARK_NOINLINE
  #if defined(_MSC_VER)
    #define DSPARK_NOINLINE __declspec(noinline)
  #elif defined(__GNUC__) || defined(__clang__)
    #define DSPARK_NOINLINE __attribute__((noinline))
  #else
    #define DSPARK_NOINLINE
  #endif
#endif
 
namespace dspark {
 
// ============================================================================
// BiquadCoeffs -- Coefficient storage + factory methods
// ============================================================================
 
struct alignas(32) BiquadCoeffs
{
    double b0 = 1.0, b1 = 0.0, b2 = 0.0;
    double a1 = 0.0, a2 = 0.0;
 
    // -- Factory methods (Audio EQ Cookbook) ----------------------------------
 
    [[nodiscard]] static BiquadCoeffs makeLowPass(double sampleRate, double freq, double Q = 0.7071067811865476) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha;
        return normalise(a0,
            (1.0 - cosw0) / 2.0,
             1.0 - cosw0,
            (1.0 - cosw0) / 2.0,
            -2.0 * cosw0,
             1.0 - alpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makeHighPass(double sampleRate, double freq, double Q = 0.7071067811865476) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha;
        return normalise(a0,
             (1.0 + cosw0) / 2.0,
            -(1.0 + cosw0),
             (1.0 + cosw0) / 2.0,
            -2.0 * cosw0,
              1.0 - alpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makeBandPass(double sampleRate, double freq, double Q = 0.7071067811865476) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha;
        return normalise(a0,
            alpha,
            0.0,
            -alpha,
            -2.0 * cosw0,
             1.0 - alpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makePeak(double sampleRate, double freq, double Q, double gainDb) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double A     = std::pow(10.0, gainDb / 40.0); // dB/40 = sqrt of linear gain
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha / A;
        return normalise(a0,
            1.0 + alpha * A,
            -2.0 * cosw0,
            1.0 - alpha * A,
            -2.0 * cosw0,
            1.0 - alpha / A);
    }
 
    [[nodiscard]] static BiquadCoeffs makePeakMatched(double sampleRate, double freq,
                                                      double Q, double gainDb) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        if (std::abs(gainDb) < 0.01)
            return { 1.0, 0.0, 0.0, 0.0, 0.0 };
 
        const double G  = std::pow(10.0, gainDb / 20.0);
        const double w0 = 2.0 * std::numbers::pi * freq / sampleRate;
 
        // Poles: impulse invariance of s^2 + 2*q*w0*s + w0^2, q = 1/(2*Qpole)
        // with the prototype's pole Q = sqrt(G) * Q.
        const double q = 1.0 / (2.0 * Q * std::sqrt(G));
        double a1 = 0.0;
        const double a2 = std::exp(-2.0 * q * w0);
        if (q <= 1.0)
        {
            a1 = -2.0 * std::exp(-q * w0) * std::cos(std::sqrt(1.0 - q * q) * w0);
        }
        else
        {
            // Two real poles z = e^((-q +/- r) w0). Evaluated one by one (each
            // is <= 1) rather than as e^(-q w0) * cosh(r w0), whose cosh
            // overflows for small Q; -q + r is taken as -1/(q + r) to avoid
            // the cancellation.
            const double r = std::sqrt(q * q - 1.0);
            a1 = -(std::exp(-w0 / (q + r)) + std::exp(-(q + r) * w0));
        }
 
        // Magnitude matching in the paper's (A, B, phi) form: |H|^2 =
        // (B0 phi0 + B1 phi1 + B2 phi2) / (A0 phi0 + A1 phi1 + A2 phi2).
        const double A0 = (1.0 + a1 + a2) * (1.0 + a1 + a2);
        const double A1 = (1.0 - a1 + a2) * (1.0 - a1 + a2);
        const double A2 = -4.0 * a2;
        const double s  = std::sin(0.5 * w0);
        const double phi1 = s * s;
        const double phi0 = 1.0 - phi1;
        const double phi2 = 4.0 * phi0 * phi1;
 
        const double G2 = G * G;
        const double B0 = A0;                                     // unity at DC
        const double R1 = (A0 * phi0 + A1 * phi1 + A2 * phi2) * G2; // |H(w0)| = G
        const double R2 = (-A0 + A1 + 4.0 * (phi0 - phi1) * A2) * G2; // extremum at w0
        const double B2 = (R1 - R2 * phi1 - B0) / (4.0 * phi1 * phi1);
        const double B1 = R2 + B0 + 4.0 * (phi1 - phi0) * B2;
 
        // Minimum-phase numerator from (B0, B1, B2).
        const double sqrtB0 = std::sqrt(std::max(0.0, B0));
        const double sqrtB1 = std::sqrt(std::max(0.0, B1));
        const double W = 0.5 * (sqrtB0 + sqrtB1);
        BiquadCoeffs c;
        c.b0 = 0.5 * (W + std::sqrt(std::max(0.0, W * W + B2)));
        c.b1 = 0.5 * (sqrtB0 - sqrtB1);
        c.b2 = -B2 / (4.0 * c.b0);
        c.a1 = a1;
        c.a2 = a2;
        return c;
    }
 
    [[nodiscard]] static BiquadCoeffs makeLowShelf(double sampleRate, double freq, double gainDb, double slope = 1.0) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        // Shelf slope S must stay in (0, 1]: for S > 1 the radicand below goes
        // negative and std::sqrt yields NaN coefficients that poison the audio.
        slope = std::clamp(slope, 0.0001, 1.0);
        const double A     = std::pow(10.0, gainDb / 40.0);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double radicand = (A + 1.0 / A) * (1.0 / slope - 1.0) + 2.0;
        const double alpha = sinw0 / 2.0 * std::sqrt(std::max(0.0, radicand));
        const double twoSqrtAAlpha = 2.0 * std::sqrt(A) * alpha;
 
        const double a0 = (A + 1.0) + (A - 1.0) * cosw0 + twoSqrtAAlpha;
        return normalise(a0,
            A * ((A + 1.0) - (A - 1.0) * cosw0 + twoSqrtAAlpha),
            2.0 * A * ((A - 1.0) - (A + 1.0) * cosw0),
            A * ((A + 1.0) - (A - 1.0) * cosw0 - twoSqrtAAlpha),
            -2.0 * ((A - 1.0) + (A + 1.0) * cosw0),
            (A + 1.0) + (A - 1.0) * cosw0 - twoSqrtAAlpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makeHighShelf(double sampleRate, double freq, double gainDb, double slope = 1.0) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        // Shelf slope S must stay in (0, 1]: for S > 1 the radicand below goes
        // negative and std::sqrt yields NaN coefficients that poison the audio.
        slope = std::clamp(slope, 0.0001, 1.0);
        const double A     = std::pow(10.0, gainDb / 40.0);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double radicand = (A + 1.0 / A) * (1.0 / slope - 1.0) + 2.0;
        const double alpha = sinw0 / 2.0 * std::sqrt(std::max(0.0, radicand));
        const double twoSqrtAAlpha = 2.0 * std::sqrt(A) * alpha;
 
        const double a0 = (A + 1.0) - (A - 1.0) * cosw0 + twoSqrtAAlpha;
        return normalise(a0,
            A * ((A + 1.0) + (A - 1.0) * cosw0 + twoSqrtAAlpha),
            -2.0 * A * ((A - 1.0) + (A + 1.0) * cosw0),
            A * ((A + 1.0) + (A - 1.0) * cosw0 - twoSqrtAAlpha),
            2.0 * ((A - 1.0) - (A + 1.0) * cosw0),
            (A + 1.0) - (A - 1.0) * cosw0 - twoSqrtAAlpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makeNotch(double sampleRate, double freq, double Q = 0.7071067811865476) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha;
        return normalise(a0,
            1.0,
            -2.0 * cosw0,
            1.0,
            -2.0 * cosw0,
            1.0 - alpha);
    }
 
    [[nodiscard]] static BiquadCoeffs makeAllPass(double sampleRate, double freq, double Q = 0.7071067811865476) noexcept
    {
        freq = std::clamp(freq, 1.0, std::max(1.0, sampleRate * 0.499));
        Q = std::max(Q, 0.001);
        const double w0    = 2.0 * std::numbers::pi * freq / sampleRate;
        const double cosw0 = std::cos(w0);
        const double sinw0 = std::sin(w0);
        const double alpha = sinw0 / (2.0 * Q);
 
        const double a0 = 1.0 + alpha;
        return normalise(a0,
            1.0 - alpha,
            -2.0 * cosw0,
            1.0 + alpha,
            -2.0 * cosw0,
            1.0 - alpha);
    }
 
    // -- First-order filter factory methods ------------------------------------
 
    [[nodiscard]] static BiquadCoeffs makeFirstOrderLowPass(double sampleRate, double frequency) noexcept
    {
        frequency = std::clamp(frequency, 1.0, std::max(1.0, sampleRate * 0.499));
        const double w = std::tan(std::numbers::pi * frequency / sampleRate);
        const double n = 1.0 / (1.0 + w);
        BiquadCoeffs c;
        c.b0 = w * n;
        c.b1 = w * n;
        c.b2 = 0.0;
        c.a1 = (w - 1.0) * n;
        c.a2 = 0.0;
        return c;
    }
 
    [[nodiscard]] static BiquadCoeffs makeFirstOrderHighPass(double sampleRate, double frequency) noexcept
    {
        frequency = std::clamp(frequency, 1.0, std::max(1.0, sampleRate * 0.499));
        const double w = std::tan(std::numbers::pi * frequency / sampleRate);
        const double n = 1.0 / (1.0 + w);
        BiquadCoeffs c;
        c.b0 = n;
        c.b1 = -n;
        c.b2 = 0.0;
        c.a1 = (w - 1.0) * n;
        c.a2 = 0.0;
        return c;
    }
 
    [[nodiscard]] static BiquadCoeffs makeTilt(double sampleRate, double pivotFreq, double gainDb) noexcept
    {
        pivotFreq = std::clamp(pivotFreq, 1.0, std::max(1.0, sampleRate * 0.499));
        const double g = std::pow(10.0, gainDb / 20.0);
        const double sqrtG = std::sqrt(g);
        const double c = std::tan(std::numbers::pi * pivotFreq / sampleRate);
 
        // First-order tilt shelf via bilinear transform:
        //   DC gain = 1/sqrt(g), pivot gain = 1, Nyquist gain = sqrt(g)
        //   Total swing = gainDb (half below pivot, half above).
        const double norm = 1.0 / (1.0 + sqrtG * c);
 
        BiquadCoeffs coeffs;
        coeffs.b0 = (sqrtG + c) * norm;
        coeffs.b1 = (c - sqrtG) * norm;
        coeffs.b2 = 0.0;
        coeffs.a1 = (sqrtG * c - 1.0) * norm;
        coeffs.a2 = 0.0;
        return coeffs;
    }
 
    [[nodiscard]] static BiquadCoeffs makeKWeightingShelf(double sampleRate) noexcept
    {
        const double G  = 3.999843853973347;     // dB
        const double Q  = 0.7071752369554196;
        const double fc = 1681.9744509555319;    // Hz
        const double K  = std::tan(std::numbers::pi * fc / sampleRate);
        const double Vh = std::pow(10.0, G / 20.0);
        const double Vb = std::pow(Vh, 0.4996667741545416);
        const double a0 = 1.0 + K / Q + K * K;
        BiquadCoeffs c;
        c.b0 = (Vh + Vb * K / Q + K * K) / a0;
        c.b1 = 2.0 * (K * K - Vh) / a0;
        c.b2 = (Vh - Vb * K / Q + K * K) / a0;
        c.a1 = 2.0 * (K * K - 1.0) / a0;
        c.a2 = (1.0 - K / Q + K * K) / a0;
        return c;
    }
 
    [[nodiscard]] static BiquadCoeffs makeKWeightingHighPass(double sampleRate) noexcept
    {
        const double Q  = 0.5003270373238773;
        const double fc = 38.13547087602444;     // Hz
        const double K  = std::tan(std::numbers::pi * fc / sampleRate);
        const double a0 = 1.0 + K / Q + K * K;
        BiquadCoeffs c;
        c.b0 = 1.0;
        c.b1 = -2.0;
        c.b2 = 1.0;
        c.a1 = 2.0 * (K * K - 1.0) / a0;
        c.a2 = (1.0 - K / Q + K * K) / a0;
        return c;
    }
 
    // -- Frequency response analysis -------------------------------------------
 
    [[nodiscard]] double getMagnitude(double frequency, double sampleRate) const noexcept
    {
        const double w = 2.0 * std::numbers::pi * frequency / sampleRate;
        const double cosW  = std::cos(w);
        const double cos2W = std::cos(2.0 * w);
        const double sinW  = std::sin(w);
        const double sin2W = std::sin(2.0 * w);
 
        const double nRe = b0 + b1 * cosW + b2 * cos2W;
        const double nIm = -b1 * sinW - b2 * sin2W;
        const double dRe = 1.0 + a1 * cosW + a2 * cos2W;
        const double dIm = -a1 * sinW - a2 * sin2W;
 
        const double numMag2 = nRe * nRe + nIm * nIm;
        const double denMag2 = dRe * dRe + dIm * dIm;
 
        return (denMag2 > 1e-30) ? std::sqrt(numMag2 / denMag2) : 0.0;
    }
 
    template <typename U>
    void getMagnitudeForFrequencyArray(std::span<const U> frequencies,
                                       std::span<U> magnitudes,
                                       double sampleRate) const noexcept
    {
        assert(magnitudes.size() >= frequencies.size());
        for (size_t i = 0; i < frequencies.size(); ++i)
            magnitudes[i] = static_cast<U>(getMagnitude(static_cast<double>(frequencies[i]), sampleRate));
    }
 
private:
    [[nodiscard]] static BiquadCoeffs normalise(double a0, double b0r, double b1r, double b2r,
                                                double a1r, double a2r) noexcept
    {
        const double invA0 = 1.0 / a0;
        return { b0r * invA0, b1r * invA0, b2r * invA0, a1r * invA0, a2r * invA0 };
    }
};
 
// ============================================================================
// Biquad -- Filter processor with per-channel state
// ============================================================================
 
template <typename T, int MaxChannels = 8>
class alignas(32) Biquad
{
public:
    Biquad() noexcept = default;
 
    // std::atomic<bool> makes the default copy/move ops vanish. Provide manual
    // move semantics so this class can live inside std::array / std::vector
    // (e.g. CrossoverFilter::AllPassChain). Copying is left deleted on
    // purpose: it would imply two threads observing the same staged coefficient
    // dirty flag, which is semantically ambiguous and not needed in practice.
    // Moves only happen during setup (single-threaded relocation into a
    // container), so the seqlock counter is simply reset to an even value.
    Biquad(Biquad&& other) noexcept
        : activeCoeffs_(other.activeCoeffs_),
          coeffsDirty_(other.coeffsDirty_.load(std::memory_order_relaxed)),
          coeffsSeq_(0),
          state_(other.state_)
    {
        copyStagedRelaxed(other); // Word-wise: the staging words are atomics.
    }
 
    Biquad& operator=(Biquad&& other) noexcept
    {
        if (this == &other) return *this;
        activeCoeffs_ = other.activeCoeffs_;
        copyStagedRelaxed(other); // Word-wise: the staging words are atomics.
        coeffsDirty_.store(other.coeffsDirty_.load(std::memory_order_relaxed),
                           std::memory_order_relaxed);
        coeffsSeq_.store(0, std::memory_order_relaxed);
        state_ = other.state_;
        return *this;
    }
 
    Biquad(const Biquad&)            = delete;
    Biquad& operator=(const Biquad&) = delete;
 
    void setCoeffs(const BiquadCoeffs& c) noexcept
    {
        // Seqlock publish (single producer). An odd sequence number marks a
        // write in progress; the audio thread retries its read while odd or if
        // the counter moved, so a concurrent update can never be observed as a
        // torn coefficient set (mixing a1 of one filter with a2 of another).
        //
        // Data-race freedom: the staging slot is 5 std::atomic<double> words
        // stored relaxed inside the critical section, so every cross-thread
        // word access is defined by the C++ memory model (the previous plain
        // struct copy was formal UB -- the same defect class CI TSan flagged in
        // FIRFilter). The release fence below pairs with the
        // reader's acquire fence ([atomics.fences]/2): a reader that observes
        // any mid-publish word must also observe the odd counter and retry
        // (Boehm, "Can Seqlocks Get Along With Programming Language Memory
        // Models?", MSPC 2012). Control-thread only: zero audio-path cost.
        coeffsSeq_.fetch_add(1, std::memory_order_acq_rel);   // -> odd
        std::atomic_thread_fence(std::memory_order_release);
        stagedCoeffs_.b0.store(c.b0, std::memory_order_relaxed);
        stagedCoeffs_.b1.store(c.b1, std::memory_order_relaxed);
        stagedCoeffs_.b2.store(c.b2, std::memory_order_relaxed);
        stagedCoeffs_.a1.store(c.a1, std::memory_order_relaxed);
        stagedCoeffs_.a2.store(c.a2, std::memory_order_relaxed);
        coeffsSeq_.fetch_add(1, std::memory_order_release);   // -> even
        coeffsDirty_.store(true, std::memory_order_release);
    }
 
    void setCoeffsNow(const BiquadCoeffs& c) noexcept { activeCoeffs_ = c; }
 
    bool applyPendingCoeffs() noexcept
    {
        if (!coeffsDirty_.exchange(false, std::memory_order_acquire)) return false;
        return adoptStagedCoeffs();
    }
 
    [[nodiscard]] const BiquadCoeffs& getCoeffs() const noexcept { return activeCoeffs_; }
 
    void reset() noexcept
    {
        for (auto& s : state_)
            s = {};
    }
 
    T processSample(T input, int channel) noexcept
    {
        return static_cast<T>(processSampleCore(static_cast<double>(input), channel));
    }
 
    double processSampleCore(double input, int channel) noexcept
    {
        assert(channel >= 0 && channel < MaxChannels && "Channel index out of bounds");
 
        // Lock-free fast path: relaxed load is free on every modern CPU.
        // Branch is marked unlikely so the compiler keeps the hot DSP path
        // straight-line; the dirty flag is only true on the first sample
        // of a block where setCoeffs() landed since last process call.
        if (coeffsDirty_.load(std::memory_order_relaxed)) [[unlikely]]
            applyPendingCoeffs();
 
        auto& s = state_[channel];
 
        const double output = activeCoeffs_.b0 * input + s.z1;
        s.z1 = activeCoeffs_.b1 * input - activeCoeffs_.a1 * output + s.z2;
        s.z2 = activeCoeffs_.b2 * input - activeCoeffs_.a2 * output;
 
        return output;
    }
 
    void processBlock(AudioBufferView<T> buffer) noexcept
    {
        DenormalGuard guard;
 
        // Pick up any pending lock-free coefficient update from the GUI thread.
        applyPendingCoeffs();
 
        const int numChannels = std::min(buffer.getNumChannels(), MaxChannels);
        const int numSamples  = buffer.getNumSamples();
 
        // Block-local coefficient copy: stable for the whole block by design
        // (updates land at the next block boundary), and register-friendly.
        const BiquadCoeffs c = activeCoeffs_;
 
        for (int ch = 0; ch < numChannels; ++ch)
        {
            T* data = buffer.getChannel(ch);
            auto& s = state_[ch];
            double z1 = s.z1;
            double z2 = s.z2;
 
            for (int i = 0; i < numSamples; ++i)
            {
                const double input  = static_cast<double>(data[i]);
                const double output = c.b0 * input + z1;
                z1 = c.b1 * input - c.a1 * output + z2;
                z2 = c.b2 * input - c.a2 * output;
                data[i] = static_cast<T>(output);
            }
 
            s.z1 = z1;
            s.z2 = z2;
        }
    }
 
private:
    static constexpr int kSeqlockMaxAttempts = 3;
 
    DSPARK_NOINLINE bool adoptStagedCoeffs() noexcept
    {
        BiquadCoeffs tmp;
        bool adopted = false;
        for (int attempt = 0; attempt < kSeqlockMaxAttempts; ++attempt)
        {
            const unsigned s0 = coeffsSeq_.load(std::memory_order_acquire);
            // Odd counter: the writer is mid-publish, so do not copy at all.
            // Testing this BEFORE the copy makes a collision cost one relaxed
            // load instead of a full word-set copy that would be thrown away,
            // which is what keeps the bound cheap enough for the per-sample
            // entry point.
            if ((s0 & 1u) != 0u) continue;
            // Relaxed word loads: each access is race-free because the staging
            // words are std::atomic<double>. The destination is a stack
            // temporary, NOT the live activeCoeffs_ the hot path reads, so a
            // give-up leaves the active set untouched (commit-on-validate).
            tmp.b0 = stagedCoeffs_.b0.load(std::memory_order_relaxed);
            tmp.b1 = stagedCoeffs_.b1.load(std::memory_order_relaxed);
            tmp.b2 = stagedCoeffs_.b2.load(std::memory_order_relaxed);
            tmp.a1 = stagedCoeffs_.a1.load(std::memory_order_relaxed);
            tmp.a2 = stagedCoeffs_.a2.load(std::memory_order_relaxed);
            // The fence keeps the copy above from sinking below the re-read of
            // the counter. A plain acquire load only orders LATER accesses;
            // without the fence, both the compiler and weakly-ordered CPUs
            // (ARM) may complete the copy after the counter is re-read, and a
            // torn copy would pass the equality check. It also pairs with the
            // writer's release fence, so a copy that read any mid-publish word
            // cannot pass validation.
            std::atomic_thread_fence(std::memory_order_acquire);
            if (s0 == coeffsSeq_.load(std::memory_order_relaxed))
            {
                adopted = true;
                break;
            }
        }
        if (!adopted)
        {
            // Adopt nothing; keep the set in use and pick the update up later.
            coeffsDirty_.store(true, std::memory_order_release);
            return false;
        }
        activeCoeffs_ = tmp;
        return true;
    }
 
    // Kept compact on purpose (no over-alignment): the states are only ever
    // touched by the processing thread, sequentially per channel, so packing
    // adjacent channels into the same cache line is strictly better than
    // padding each one out to its own.
    struct State
    {
        double z1 = 0.0;
        double z2 = 0.0;
    };
 
    // Shared staging slot for the seqlock publish: written word-by-word by the
    // control thread inside the odd/even critical section and read word-by-word
    // by the audio thread's retry loop. Every word is std::atomic<double>
    // (relaxed) so each cross-thread access is defined by the C++ memory model
    // (no non-atomic concurrent read/write, i.e. no data-race UB); the seq
    // counter plus the writer release / reader acquire fence pair order and
    // tear-protect the 5-word set. std::atomic<double> is lock-free and one
    // machine word on every supported target, so the publish stays lock- and
    // allocation-free. Defaults mirror BiquadCoeffs{} (identity filter).
    struct StagedCoeffs
    {
        std::atomic<double> b0{1.0}, b1{0.0}, b2{0.0};
        std::atomic<double> a1{0.0}, a2{0.0};
    };
 
    void copyStagedRelaxed(const Biquad& other) noexcept
    {
        stagedCoeffs_.b0.store(other.stagedCoeffs_.b0.load(std::memory_order_relaxed),
                               std::memory_order_relaxed);
        stagedCoeffs_.b1.store(other.stagedCoeffs_.b1.load(std::memory_order_relaxed),
                               std::memory_order_relaxed);
        stagedCoeffs_.b2.store(other.stagedCoeffs_.b2.load(std::memory_order_relaxed),
                               std::memory_order_relaxed);
        stagedCoeffs_.a1.store(other.stagedCoeffs_.a1.load(std::memory_order_relaxed),
                               std::memory_order_relaxed);
        stagedCoeffs_.a2.store(other.stagedCoeffs_.a2.load(std::memory_order_relaxed),
                               std::memory_order_relaxed);
    }
 
    BiquadCoeffs activeCoeffs_ {};
    StagedCoeffs stagedCoeffs_ {};
    std::atomic<bool> coeffsDirty_{false};
    std::atomic<unsigned> coeffsSeq_{0}; 
 
    std::array<State, MaxChannels> state_ {};
};
 
} // namespace dspark