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

sva_samples/
biquad.rs
