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_db, 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    #[inline]
25    pub fn step(&mut self, c: &Coeffs, x: f64) -> f64 {
26        let y = c.b0 * x + c.b1 * self.x1 + c.b2 * self.x2 - c.a1 * self.y1 - c.a2 * self.y2;
27        self.x2 = self.x1;
28        self.x1 = x;
29        self.y2 = self.y1;
30        self.y1 = y;
31        y
32    }
33}
34
35/// Outside these it is a wire or a ringing sine, not a filter.
36pub const MIN_Q: f64 = 1e-3;
37pub const MAX_Q: f64 = 100.0;
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
55pub fn clamp_q(q: f64) -> (f64, bool) {
56    if q < MIN_Q {
57        (MIN_Q, true)
58    } else if q > MAX_Q {
59        (MAX_Q, true)
60    } else {
61        (q, false)
62    }
63}
64
65/// The RBJ Audio EQ Cookbook's bilinear-transform designs, normalized by `a0`.
66pub fn design(shape: Shape, cutoff: f64, q: f64, gain_db: f64, sr: f64) -> Coeffs {
67    if shape == Shape::OnePole {
68        let rc = 1.0 / (2.0 * std::f64::consts::PI * cutoff);
69        let dt = 1.0 / sr;
70        let alpha = dt / (rc + dt);
71        return Coeffs {
72            b0: alpha,
73            b1: 0.0,
74            b2: 0.0,
75            a1: alpha - 1.0,
76            a2: 0.0,
77        };
78    }
79
80    let w0 = 2.0 * std::f64::consts::PI * cutoff / sr;
81    let cs = w0.cos();
82    let sn = w0.sin();
83    let alpha = sn / (2.0 * q);
84    let a = 10f64.powf(gain_db / 40.0);
85    let sqrt_a = a.sqrt();
86
87    let (b0, b1, b2, a0, a1, a2) = match shape {
88        Shape::OnePole => unreachable!("handled above"),
89        Shape::Lowpass => (
90            (1.0 - cs) / 2.0,
91            1.0 - cs,
92            (1.0 - cs) / 2.0,
93            1.0 + alpha,
94            -2.0 * cs,
95            1.0 - alpha,
96        ),
97        Shape::Highpass => (
98            (1.0 + cs) / 2.0,
99            -(1.0 + cs),
100            (1.0 + cs) / 2.0,
101            1.0 + alpha,
102            -2.0 * cs,
103            1.0 - alpha,
104        ),
105        Shape::Bandpass => (alpha, 0.0, -alpha, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
106        Shape::Notch => (1.0, -2.0 * cs, 1.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
107        Shape::Peaking => (
108            1.0 + alpha * a,
109            -2.0 * cs,
110            1.0 - alpha * a,
111            1.0 + alpha / a,
112            -2.0 * cs,
113            1.0 - alpha / a,
114        ),
115        Shape::Lowshelf => (
116            a * ((a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
117            2.0 * a * ((a - 1.0) - (a + 1.0) * cs),
118            a * ((a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
119            (a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
120            -2.0 * ((a - 1.0) + (a + 1.0) * cs),
121            (a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
122        ),
123        Shape::Highshelf => (
124            a * ((a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
125            -2.0 * a * ((a - 1.0) + (a + 1.0) * cs),
126            a * ((a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
127            (a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
128            2.0 * ((a - 1.0) - (a + 1.0) * cs),
129            (a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
130        ),
131    };
132
133    Coeffs {
134        b0: b0 / a0,
135        b1: b1 / a0,
136        b2: b2 / a0,
137        a1: a1 / a0,
138        a2: a2 / a0,
139    }
140}
141
142pub fn magnitude_db(c: &Coeffs, hz: f64, sr: f64) -> f64 {
143    let w = 2.0 * std::f64::consts::PI * hz / sr;
144    let (c1, s1) = (w.cos(), -w.sin());
145    let (c2, s2) = ((2.0 * w).cos(), -(2.0 * w).sin());
146    let num_re = c.b0 + c.b1 * c1 + c.b2 * c2;
147    let num_im = c.b1 * s1 + c.b2 * s2;
148    let den_re = 1.0 + c.a1 * c1 + c.a2 * c2;
149    let den_im = c.a1 * s1 + c.a2 * s2;
150    let num = (num_re * num_re + num_im * num_im).sqrt();
151    let den = (den_re * den_re + den_im * den_im).sqrt();
152    if den == 0.0 || num == 0.0 {
153        return f64::NEG_INFINITY;
154    }
155    20.0 * (num / den).log10()
156}
157
158pub const RESPONSE_RATIOS: [f64; 9] = [0.125, 0.25, 0.5, 0.707, 1.0, 1.414, 2.0, 4.0, 8.0];
159
160/// A point at or above Nyquist is dropped, not clamped: clamping reports the transform's own
161/// zero there as if the filter had put it there.
162pub fn response(c: &Coeffs, cutoff: f64, sr: f64) -> Vec<(f64, f64)> {
163    RESPONSE_RATIOS
164        .iter()
165        .map(|r| cutoff * r)
166        .filter(|hz| *hz < 0.5 * sr)
167        .map(|hz| (hz, magnitude_db(c, hz, sr)))
168        .collect()
169}
sva_samples/biquad.rs

sva_samples/
biquad.rs
