sva_samples/
biquad.rs

1// Concern: designs RBJ biquad coefficients and their magnitude response | Non-concern: the shapes and their H(s) (sva-formula), dispatch (filters.rs) | IO: (Shape, cutoff, q, gain, sr) -> Coeffs
2
3use sva_formula::filter::Shape;
4
5/// Normalized to `a0 = 1`: `y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2`.
6#[derive(Clone, Copy, Debug, PartialEq)]
7pub struct Coeffs {
8    pub b0: f64,
9    pub b1: f64,
10    pub b2: f64,
11    pub a1: f64,
12    pub a2: f64,
13}
14
15#[derive(Clone, Copy, Debug, Default)]
16pub struct State {
17    x1: f64,
18    x2: f64,
19    y1: f64,
20    y2: f64,
21}
22
23impl State {
24    pub fn held(&self) -> [f64; 4] {
25        [self.x1, self.x2, self.y1, self.y2]
26    }
27
28    #[inline]
29    pub fn step(&mut self, c: &Coeffs, x: f64) -> f64 {
30        let y = c.b0 * x + c.b1 * self.x1 + c.b2 * self.x2 - c.a1 * self.y1 - c.a2 * self.y2;
31        self.x2 = self.x1;
32        self.x1 = x;
33        self.y2 = self.y1;
34        self.y1 = y;
35        y
36    }
37}
38
39/// Short of the top tenth, where `tan` runs away.
40const HIGHEST_FRACTION_OF_RATE: f64 = 0.45;
41
42/// A sweep past 0 Hz or Nyquist blows the bilinear design up; the clamp is reported.
43pub fn clamp_cutoff(hz: f64, sr: f64) -> (f64, bool) {
44    let hi = HIGHEST_FRACTION_OF_RATE * sr;
45    let lo = hi.min(1.0);
46    if hz < lo {
47        (lo, true)
48    } else if hz > hi {
49        (hi, true)
50    } else {
51        (hz, false)
52    }
53}
54
55/// -120 dB; a ratio at or under 0 has no root.
56pub const MIN_GAIN: f64 = 1e-6;
57
58pub fn clamp_gain(gain: f64) -> (f64, bool) {
59    match gain < MIN_GAIN {
60        true => (MIN_GAIN, true),
61        false => (gain, false),
62    }
63}
64
65/// Correctness: damping `d = s/c` (`a0, a2 = c ± s`; `b0, b2` for zeros) rounds
66/// within `2^-20` of itself from `2^-31`; the ringing bound holds to `2^22`.
67pub const DAMPING: (f64, f64) = (1.0 / (1u64 << 31) as f64, (1u64 << 22) as f64);
68
69fn damping_at_unit_q(shape: Shape, cutoff: f64, gain: f64, sr: f64) -> (f64, f64) {
70    let w0 = 2.0 * std::f64::consts::PI * cutoff / sr;
71    let (cs, half_sn) = (w0.cos(), w0.sin() / 2.0);
72    let a = gain.sqrt();
73    let shelf = |poles_c: f64, zeros_c: f64| {
74        let s = 2.0 * a.sqrt() * half_sn;
75        (s / poles_c, s / zeros_c)
76    };
77    match shape {
78        Shape::Peaking => (half_sn / a, half_sn * a),
79        Shape::Lowshelf => shelf((a + 1.0) + (a - 1.0) * cs, (a + 1.0) - (a - 1.0) * cs),
80        Shape::Highshelf => shelf((a + 1.0) - (a - 1.0) * cs, (a + 1.0) + (a - 1.0) * cs),
81        _ => (half_sn, half_sn),
82    }
83}
84
85/// Holds q inside `DAMPING`; zeros no q fits are flagged.
86pub fn clamp_q(shape: Shape, cutoff: f64, q: f64, gain: f64, sr: f64) -> (f64, bool) {
87    if shape == Shape::OnePole {
88        return (q, false);
89    }
90    let (poles, zeros) = damping_at_unit_q(shape, cutoff, gain, sr);
91    let (least, most) = DAMPING;
92    let (mut lo, mut hi) = (poles.max(zeros) / most, poles.min(zeros) / least);
93    let apart = lo > hi;
94    if apart {
95        (lo, hi) = (poles / most, poles / least);
96    }
97    if q < lo || q.is_nan() {
98        (lo, true)
99    } else if q > hi {
100        (hi, true)
101    } else {
102        (q, apart)
103    }
104}
105
106/// The RBJ Audio EQ Cookbook's designs, normalized by `a0`; `gain` is a ratio.
107pub fn design(shape: Shape, cutoff: f64, q: f64, gain: f64, sr: f64) -> Coeffs {
108    if shape == Shape::OnePole {
109        let rc = 1.0 / (2.0 * std::f64::consts::PI * cutoff);
110        let dt = 1.0 / sr;
111        let alpha = dt / (rc + dt);
112        return Coeffs {
113            b0: alpha,
114            b1: 0.0,
115            b2: 0.0,
116            a1: alpha - 1.0,
117            a2: 0.0,
118        };
119    }
120
121    let w0 = 2.0 * std::f64::consts::PI * cutoff / sr;
122    let cs = w0.cos();
123    let sn = w0.sin();
124    let alpha = sn / (2.0 * q);
125    let a = gain.sqrt();
126    let sqrt_a = a.sqrt();
127
128    let (b0, b1, b2, a0, a1, a2) = match shape {
129        Shape::OnePole => unreachable!("handled above"),
130        Shape::Lowpass => (
131            (1.0 - cs) / 2.0,
132            1.0 - cs,
133            (1.0 - cs) / 2.0,
134            1.0 + alpha,
135            -2.0 * cs,
136            1.0 - alpha,
137        ),
138        Shape::Highpass => (
139            (1.0 + cs) / 2.0,
140            -(1.0 + cs),
141            (1.0 + cs) / 2.0,
142            1.0 + alpha,
143            -2.0 * cs,
144            1.0 - alpha,
145        ),
146        Shape::Bandpass => (alpha, 0.0, -alpha, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
147        Shape::Notch => (1.0, -2.0 * cs, 1.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
148        Shape::Peaking => (
149            1.0 + alpha * a,
150            -2.0 * cs,
151            1.0 - alpha * a,
152            1.0 + alpha / a,
153            -2.0 * cs,
154            1.0 - alpha / a,
155        ),
156        Shape::Lowshelf => (
157            a * ((a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
158            2.0 * a * ((a - 1.0) - (a + 1.0) * cs),
159            a * ((a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
160            (a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
161            -2.0 * ((a - 1.0) + (a + 1.0) * cs),
162            (a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
163        ),
164        Shape::Highshelf => (
165            a * ((a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
166            -2.0 * a * ((a - 1.0) + (a + 1.0) * cs),
167            a * ((a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
168            (a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
169            2.0 * ((a - 1.0) - (a + 1.0) * cs),
170            (a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
171        ),
172    };
173
174    Coeffs {
175        b0: b0 / a0,
176        b1: b1 / a0,
177        b2: b2 / a0,
178        a1: a1 / a0,
179        a2: a2 / a0,
180    }
181}
182
183pub fn magnitude_db(c: &Coeffs, hz: f64, sr: f64) -> f64 {
184    let w = 2.0 * std::f64::consts::PI * hz / sr;
185    let (c1, s1) = (w.cos(), -w.sin());
186    let (c2, s2) = ((2.0 * w).cos(), -(2.0 * w).sin());
187    let num_re = c.b0 + c.b1 * c1 + c.b2 * c2;
188    let num_im = c.b1 * s1 + c.b2 * s2;
189    let den_re = 1.0 + c.a1 * c1 + c.a2 * c2;
190    let den_im = c.a1 * s1 + c.a2 * s2;
191    let num = (num_re * num_re + num_im * num_im).sqrt();
192    let den = (den_re * den_re + den_im * den_im).sqrt();
193    if den == 0.0 || num == 0.0 {
194        return f64::NEG_INFINITY;
195    }
196    20.0 * (num / den).log10()
197}
198
199pub const RESPONSE_RATIOS: [f64; 9] = [0.125, 0.25, 0.5, 0.707, 1.0, 1.414, 2.0, 4.0, 8.0];
200
201/// A point at or above Nyquist is dropped, not clamped: clamping reports the transform's own
202/// zero there as if the filter had put it there.
203pub fn response(c: &Coeffs, cutoff: f64, sr: f64) -> Vec<(f64, f64)> {
204    RESPONSE_RATIOS
205        .iter()
206        .map(|r| cutoff * r)
207        .filter(|hz| *hz < 0.5 * sr)
208        .map(|hz| (hz, magnitude_db(c, hz, sr)))
209        .collect()
210}
sva_samples/biquad.rs

sva_samples/
biquad.rs
