bass-foundry/plugins/weld/core/src/chain.rs

//! Polarity plus a cascade of second-order allpass sections. Magnitude is flat by
//! construction; only the phase of the main layer changes.
use serde::{Deserialize, Serialize};

pub const MAX_SECTIONS: usize = 4;
pub const MIN_FC: f32 = 30.;
pub const MIN_Q: f32 = 0.2;
pub const MAX_Q: f32 = 4.;

/// The range where listeners hear group-delay changes best (Liski, Mäkivirta and
/// Välimäki 2021, "Audibility of Group-Delay Equalization").
pub const SMEAR_LOW: f32 = 500.;
pub const SMEAR_HIGH: f32 = 4000.;
/// Smallest mean threshold that study found for a positive group-delay peak in that
/// range, with a pink impulse over headphones. An allpass only adds positive delay.
/// Sustained bass hides much more; this marks where sharp attacks may soften.
pub const SMEAR_AUDIBLE_MS: f32 = 0.64;
const SMEAR_POINTS: usize = 64;

#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct Section {
    /// Centre frequency in Hz, where the section's phase passes −180°.
    pub fc: f32,
    pub q: f32,
}
impl Default for Section {
    fn default() -> Self {
        Self { fc: 1000., q: 0.7 }
    }
}

/// The result of a fit: what Weld applies to the main layer.
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize, Default)]
pub struct Fit {
    /// Invert the main layer.
    pub flip: bool,
    /// Number of active sections, 0..=MAX_SECTIONS.
    pub count: usize,
    pub sections: [Section; MAX_SECTIONS],
}
impl Fit {
    /// The identity: no polarity change, no sections.
    pub const NONE: Fit = Fit {
        flip: false,
        count: 0,
        sections: [Section { fc: 1000., q: 0.7 }; MAX_SECTIONS],
    };
    pub fn active(&self) -> &[Section] {
        &self.sections[..self.count.min(MAX_SECTIONS)]
    }
    /// Clamps every value into the range the filter supports; repairs restored state.
    pub fn sanitized(mut self, rate: f32) -> Self {
        self.count = self.count.min(MAX_SECTIONS);
        let top = max_fc(rate);
        for s in &mut self.sections {
            s.fc = if s.fc.is_finite() {
                s.fc.clamp(MIN_FC, top)
            } else {
                1000f32.min(top)
            };
            s.q = if s.q.is_finite() {
                s.q.clamp(MIN_Q, MAX_Q)
            } else {
                0.7
            };
        }
        self
    }
    /// Complex response at `f` Hz as (re, im).
    pub fn response(&self, f: f32, rate: f32) -> (f32, f32) {
        let w = std::f32::consts::TAU * f / rate;
        let (mut re, mut im) = (if self.flip { -1. } else { 1. }, 0.);
        for s in self.active() {
            let c = Coefs::allpass(s.fc, s.q, rate);
            let (hr, hi) = c.response(w.cos(), w.sin(), (2. * w).cos(), (2. * w).sin());
            (re, im) = (re * hr - im * hi, re * hi + im * hr);
        }
        (re, im)
    }
    /// Group delay at `f` Hz in seconds. Polarity adds none.
    pub fn group_delay(&self, f: f32, rate: f32) -> f32 {
        let w = std::f32::consts::TAU * f / rate;
        let (c1, s1, c2, s2) = (w.cos(), w.sin(), (2. * w).cos(), (2. * w).sin());
        let samples: f32 = self
            .active()
            .iter()
            .map(|s| Coefs::allpass(s.fc, s.q, rate).group_delay(c1, s1, c2, s2))
            .sum();
        samples / rate
    }
    /// How much the group delay varies between `SMEAR_LOW` and `SMEAR_HIGH`, in ms:
    /// the largest delay minus the smallest. A constant delay counts as zero.
    pub fn smear_ms(&self, rate: f32) -> f32 {
        if self.active().is_empty() {
            return 0.;
        }
        let high = SMEAR_HIGH.min(max_fc(rate));
        let (mut lo, mut hi) = (f32::MAX, f32::MIN);
        for i in 0..SMEAR_POINTS {
            let f = SMEAR_LOW * (high / SMEAR_LOW).powf(i as f32 / (SMEAR_POINTS - 1) as f32);
            let d = self.group_delay(f, rate);
            lo = lo.min(d);
            hi = hi.max(d);
        }
        (hi - lo).max(0.) * 1000.
    }
}

pub fn max_fc(rate: f32) -> f32 {
    0.45 * rate
}

/// Normalised allpass coefficients: H(z) = (a2 + a1 z⁻¹ + z⁻²) / (1 + a1 z⁻¹ + a2 z⁻²).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Coefs {
    pub a1: f32,
    pub a2: f32,
}
impl Coefs {
    /// RBJ cookbook allpass.
    pub fn allpass(fc: f32, q: f32, rate: f32) -> Self {
        let fc = fc.clamp(MIN_FC, max_fc(rate));
        let q = q.clamp(MIN_Q, MAX_Q);
        let w = std::f32::consts::TAU * fc / rate;
        let alpha = w.sin() / (2. * q);
        let a0 = 1. + alpha;
        Self {
            a1: -2. * w.cos() / a0,
            a2: (1. - alpha) / a0,
        }
    }
    /// Response at z⁻¹ = e^{−jω}, given cos/sin of ω and 2ω.
    /// For an allpass the numerator is z⁻²·conj(den), so H = z⁻²·conj(den)/den.
    #[inline]
    pub fn response(&self, c1: f32, s1: f32, c2: f32, s2: f32) -> (f32, f32) {
        let x = 1. + self.a1 * c1 + self.a2 * c2;
        let y = -(self.a1 * s1 + self.a2 * s2);
        let m = x * x + y * y;
        // conj(den)/den = conj(den)² / |den|²
        let (qr, qi) = ((x * x - y * y) / m, -2. * x * y / m);
        // times z⁻² = c2 − j s2
        (qr * c2 + qi * s2, qi * c2 - qr * s2)
    }
    /// Group delay in samples, given cos/sin of ω and 2ω. With H = z⁻²·conj(D)/D the
    /// phase is −2ω − 2·arg D, so the delay is 2 + 2·d(arg D)/dω.
    #[inline]
    pub fn group_delay(&self, c1: f32, s1: f32, c2: f32, s2: f32) -> f32 {
        let x = 1. + self.a1 * c1 + self.a2 * c2;
        let y = -(self.a1 * s1 + self.a2 * s2);
        let dx = -(self.a1 * s1 + 2. * self.a2 * s2);
        let dy = -(self.a1 * c1 + 2. * self.a2 * c2);
        2. + 2. * (x * dy - y * dx) / (x * x + y * y)
    }
}

/// Transposed direct form II state, in f64 so sections near 30 Hz stay quiet.
#[derive(Debug, Clone, Copy, Default)]
struct State {
    s1: f64,
    s2: f64,
}

/// A running chain for up to two channels. Allocation-free.
#[derive(Debug, Clone)]
pub struct Chain {
    fit: Fit,
    sign: f32,
    coefs: [(f64, f64); MAX_SECTIONS],
    state: [[State; MAX_SECTIONS]; 2],
}
impl Chain {
    pub fn new(fit: Fit, rate: f32) -> Self {
        let mut chain = Self {
            fit: Fit::NONE,
            sign: 1.,
            coefs: [(0., 0.); MAX_SECTIONS],
            state: [[State::default(); MAX_SECTIONS]; 2],
        };
        chain.set(fit, rate);
        chain
    }
    pub fn fit(&self) -> &Fit {
        &self.fit
    }
    /// Replaces the coefficients and clears the state.
    pub fn set(&mut self, fit: Fit, rate: f32) {
        self.fit = fit.sanitized(rate);
        self.sign = if self.fit.flip { -1. } else { 1. };
        for (c, s) in self.coefs.iter_mut().zip(self.fit.active()) {
            let k = Coefs::allpass(s.fc, s.q, rate);
            *c = (k.a1 as f64, k.a2 as f64);
        }
        self.reset();
    }
    pub fn reset(&mut self) {
        self.state = [[State::default(); MAX_SECTIONS]; 2];
    }
    #[inline]
    pub fn process(&mut self, channel: usize, x: f32) -> f32 {
        let mut y = (self.sign * x) as f64;
        for (k, &(a1, a2)) in self.coefs[..self.fit.count].iter().enumerate() {
            let st = &mut self.state[channel][k];
            let out = a2 * y + st.s1;
            st.s1 = a1 * y - a1 * out + st.s2;
            st.s2 = y - a2 * out;
            y = out;
        }
        y as f32
    }
}