Skip to main content

math_audio_dsp/
response.rs

1//! Frequency-response helpers for linear DSP models.
2
3use math_audio_iir_fir::{Biquad, BiquadFilterType};
4use num_complex::Complex64;
5use std::f64::consts::PI;
6
7/// Complex response of an IIR biquad at `freq_hz`.
8pub fn biquad_complex_response(biquad: &Biquad<f64>, freq_hz: f64) -> Complex64 {
9    let (a1, a2, b0, b1, b2) = biquad.constants();
10    let omega = 2.0 * PI * freq_hz / biquad.srate;
11    let z_inv = Complex64::from_polar(1.0, -omega);
12    let z_inv2 = z_inv * z_inv;
13    let num = b0 + b1 * z_inv + b2 * z_inv2;
14    let den = 1.0 + a1 * z_inv + a2 * z_inv2;
15    num / den
16}
17
18/// Complex response of an FIR impulse at `freq_hz`.
19pub fn fir_complex_response(taps: &[f64], freq_hz: f64, sample_rate: f64) -> Complex64 {
20    if taps.is_empty() {
21        return Complex64::new(1.0, 0.0);
22    }
23    let omega = -2.0 * PI * freq_hz / sample_rate;
24    taps.iter()
25        .enumerate()
26        .fold(Complex64::new(0.0, 0.0), |acc, (idx, tap)| {
27            acc + Complex64::from_polar(*tap, omega * idx as f64)
28        })
29}
30
31/// Approximate the runtime LR4 crossover response for one output band.
32///
33/// The crossover plugin currently uses LR4 internally regardless of the
34/// configured type string, so this helper models two cascaded second-order
35/// Butterworth sections.
36pub fn lr4_crossover_response(
37    output: &str,
38    cutoff_hz: f64,
39    freq_hz: f64,
40    sample_rate: f64,
41) -> Result<Complex64, String> {
42    if cutoff_hz <= 0.0 || !cutoff_hz.is_finite() {
43        return Err(format!(
44            "crossover cutoff must be positive, got {cutoff_hz}"
45        ));
46    }
47    let filter_type = match output.to_ascii_lowercase().as_str() {
48        "low" | "lowpass" | "lp" => BiquadFilterType::Lowpass,
49        "high" | "highpass" | "hp" => BiquadFilterType::Highpass,
50        "both" => return Ok(Complex64::new(1.0, 0.0)),
51        other => return Err(format!("unsupported crossover output mode '{other}'")),
52    };
53    let section = Biquad::new(
54        filter_type,
55        cutoff_hz,
56        sample_rate,
57        std::f64::consts::FRAC_1_SQRT_2,
58        0.0,
59    );
60    let response = biquad_complex_response(&section, freq_hz);
61    Ok(response * response)
62}
63
64#[cfg(test)]
65mod tests {
66    use super::*;
67
68    #[test]
69    fn fir_response_tracks_delay_phase() {
70        let taps = [0.0, 1.0];
71        let response = fir_complex_response(&taps, 12_000.0, 48_000.0);
72        assert!(response.re.abs() < 1e-12);
73        assert!((response.im + 1.0).abs() < 1e-12);
74    }
75
76    #[test]
77    fn lr4_low_high_sum_near_unity_at_crossover() {
78        let low = lr4_crossover_response("low", 1_000.0, 1_000.0, 48_000.0).unwrap();
79        let high = lr4_crossover_response("high", 1_000.0, 1_000.0, 48_000.0).unwrap();
80        assert!((low.norm() - 0.5).abs() < 0.05);
81        assert!((high.norm() - 0.5).abs() < 0.05);
82    }
83}