codename-206/src/dsp/crossover.rs

// SPDX-License-Identifier: GPL-3.0-or-later
//! 3-band Linkwitz-Riley (LR4, 24 dB/oct) crossover filterbank.
//!
//! Each crossover splits into a low-passed band output and a high-passed remainder that feeds
//! the next crossover. Because higher bands pass through more filters, lower bands are phase-
//! compensated with an all-pass at every *later* crossover frequency so the three bands sum back
//! to flat **magnitude** (the sum is an all-pass of the input — phase-shifted, not bit-identical,
//! which is inherent to IIR Linkwitz-Riley). Approach mirrors NIH-plug's `crossover` plugin.
//!
//! For 3 bands there are two crossovers (low/mid at `f_lo`, mid/high at `f_hi`); only the low
//! band needs compensation (one all-pass at `f_hi`).

use super::biquad::{Biquad, BiquadCoefficients, NEUTRAL_Q};

/// Mono/stereo only, matching the plugin's audio layouts.
const MAX_CHANNELS: usize = 2;

/// One channel's worth of filter state for the 3-band split.
#[derive(Clone, Copy, Default)]
struct BandSplitter {
    lp_lo: [Biquad; 2], // LR4 low-pass at f_lo (two cascaded Butterworth)
    hp_lo: [Biquad; 2], // LR4 high-pass at f_lo
    lp_hi: [Biquad; 2], // LR4 low-pass at f_hi
    hp_hi: [Biquad; 2], // LR4 high-pass at f_hi
    ap_low: Biquad,     // all-pass at f_hi, phase-compensates the low band
}

impl BandSplitter {
    /// Split one sample into `[low, mid, high]`.
    fn split(&mut self, x: f32) -> [f32; 3] {
        // Crossover at f_lo: low-passed band + high-passed remainder.
        let mut lp = x;
        for f in &mut self.lp_lo {
            lp = f.process(lp);
        }
        let mut hp = x;
        for f in &mut self.hp_lo {
            hp = f.process(hp);
        }

        // Low band is phase-compensated for the f_hi crossover the upper bands pass through.
        let low = self.ap_low.process(lp);

        // Crossover at f_hi splits the remainder into mid + high.
        let mut mid = hp;
        for f in &mut self.lp_hi {
            mid = f.process(mid);
        }
        let mut high = hp;
        for f in &mut self.hp_hi {
            high = f.process(high);
        }

        [low, mid, high]
    }

    fn set_coefficients(
        &mut self,
        lp_lo: BiquadCoefficients,
        hp_lo: BiquadCoefficients,
        lp_hi: BiquadCoefficients,
        hp_hi: BiquadCoefficients,
        ap_low: BiquadCoefficients,
    ) {
        for f in &mut self.lp_lo {
            f.set_coefficients(lp_lo);
        }
        for f in &mut self.hp_lo {
            f.set_coefficients(hp_lo);
        }
        for f in &mut self.lp_hi {
            f.set_coefficients(lp_hi);
        }
        for f in &mut self.hp_hi {
            f.set_coefficients(hp_hi);
        }
        self.ap_low.set_coefficients(ap_low);
    }

    fn reset(&mut self) {
        for f in self
            .lp_lo
            .iter_mut()
            .chain(&mut self.hp_lo)
            .chain(&mut self.lp_hi)
            .chain(&mut self.hp_hi)
        {
            f.reset();
        }
        self.ap_low.reset();
    }
}

pub struct Crossover {
    channels: usize,
    splitters: [BandSplitter; MAX_CHANNELS],
}

impl Default for Crossover {
    fn default() -> Self {
        Self {
            channels: 2,
            splitters: [BandSplitter::default(); MAX_CHANNELS],
        }
    }
}

impl Crossover {
    pub fn new() -> Self {
        Self::default()
    }

    /// Set the active channel count and clear state. Call from `initialize()`.
    pub fn prepare(&mut self, channels: usize) {
        self.channels = channels.clamp(1, MAX_CHANNELS);
        self.reset();
    }

    /// Recompute and apply crossover coefficients. Cheap enough to call once per block.
    /// Frequencies are clamped to a valid range and forced monotonic (`f_lo <= f_hi`).
    pub fn update(&mut self, sample_rate: f32, low_hz: f32, high_hz: f32) {
        let max_hz = sample_rate * 0.49;
        let f_lo = low_hz.clamp(20.0, max_hz);
        let f_hi = high_hz.clamp(f_lo, max_hz);

        let lp_lo = BiquadCoefficients::lowpass(sample_rate, f_lo, NEUTRAL_Q);
        let hp_lo = BiquadCoefficients::highpass(sample_rate, f_lo, NEUTRAL_Q);
        let lp_hi = BiquadCoefficients::lowpass(sample_rate, f_hi, NEUTRAL_Q);
        let hp_hi = BiquadCoefficients::highpass(sample_rate, f_hi, NEUTRAL_Q);
        let ap_low = BiquadCoefficients::allpass(sample_rate, f_hi, NEUTRAL_Q);

        for s in &mut self.splitters {
            s.set_coefficients(lp_lo, hp_lo, lp_hi, hp_hi, ap_low);
        }
    }

    pub fn reset(&mut self) {
        for s in &mut self.splitters {
            s.reset();
        }
    }

    /// Split one sample of `channel` into `[low, mid, high]`.
    #[inline]
    pub fn split(&mut self, channel: usize, x: f32) -> [f32; 3] {
        self.splitters[channel].split(x)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f32::consts::TAU;

    const SR: f32 = 48_000.0;

    #[test]
    fn bands_sum_to_flat_magnitude() {
        // LR4 bands sum to an all-pass: the magnitude is flat at every frequency (including the
        // crossovers), even though the time-domain signal is phase-shifted (so it is NOT a
        // bit-exact null — that only holds for linear-phase FIR crossovers).
        let mut xo = Crossover::new();
        xo.prepare(1);
        xo.update(SR, 200.0, 2_500.0);

        for &f in &[50.0, 200.0, 1_000.0, 2_500.0, 9_000.0] {
            xo.reset();
            let n = 24_000usize;
            let (mut in_acc, mut out_acc) = (0.0f64, 0.0f64);
            for i in 0..n {
                let x = (TAU * f * i as f32 / SR).sin();
                let [lo, mid, hi] = xo.split(0, x);
                let y = lo + mid + hi;
                if i >= n - 8_000 {
                    in_acc += (x * x) as f64;
                    out_acc += (y * y) as f64;
                }
            }
            let ratio = (out_acc / in_acc).sqrt() as f32;
            assert!(
                (ratio - 1.0).abs() < 0.06,
                "reconstruction not flat at {f} Hz: {ratio}"
            );
        }
    }

    #[test]
    fn bands_are_actually_split() {
        // Sanity: the low band should keep lows and reject highs; the high band vice versa.
        fn band_energy(band: usize, freq: f32) -> f64 {
            let mut xo = Crossover::new();
            xo.prepare(1);
            xo.update(SR, 200.0, 2_500.0);
            let n = 24_000usize;
            let mut acc = 0.0f64;
            for i in 0..n {
                let x = (TAU * freq * i as f32 / SR).sin();
                let bands = xo.split(0, x);
                if i >= n - 8_000 {
                    acc += (bands[band] * bands[band]) as f64;
                }
            }
            acc
        }
        assert!(band_energy(0, 50.0) > band_energy(0, 9_000.0) * 100.0); // low band: lows >> highs
        assert!(band_energy(2, 9_000.0) > band_energy(2, 50.0) * 100.0); // high band: highs >> lows
    }
}