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/// 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
59/// -120 dB; a ratio at or under 0 has no root.
60pub const MIN_GAIN: f64 = 1e-6;
61
62pub fn clamp_gain(gain: f64) -> (f64, bool) {
63    match gain < MIN_GAIN {
64        true => (MIN_GAIN, true),
65        false => (gain, false),
66    }
67}
68
69pub fn clamp_q(q: f64) -> (f64, bool) {
70    if q < MIN_Q {
71        (MIN_Q, true)
72    } else if q > MAX_Q {
73        (MAX_Q, true)
74    } else {
75        (q, false)
76    }
77}
78
79/// The RBJ Audio EQ Cookbook's designs, normalized by `a0`; `gain` is a ratio.
80pub fn design(shape: Shape, cutoff: f64, q: f64, gain: f64, sr: f64) -> Coeffs {
81    if shape == Shape::OnePole {
82        let rc = 1.0 / (2.0 * std::f64::consts::PI * cutoff);
83        let dt = 1.0 / sr;
84        let alpha = dt / (rc + dt);
85        return Coeffs {
86            b0: alpha,
87            b1: 0.0,
88            b2: 0.0,
89            a1: alpha - 1.0,
90            a2: 0.0,
91        };
92    }
93
94    let w0 = 2.0 * std::f64::consts::PI * cutoff / sr;
95    let cs = w0.cos();
96    let sn = w0.sin();
97    let alpha = sn / (2.0 * q);
98    let a = gain.sqrt();
99    let sqrt_a = a.sqrt();
100
101    let (b0, b1, b2, a0, a1, a2) = match shape {
102        Shape::OnePole => unreachable!("handled above"),
103        Shape::Lowpass => (
104            (1.0 - cs) / 2.0,
105            1.0 - cs,
106            (1.0 - cs) / 2.0,
107            1.0 + alpha,
108            -2.0 * cs,
109            1.0 - alpha,
110        ),
111        Shape::Highpass => (
112            (1.0 + cs) / 2.0,
113            -(1.0 + cs),
114            (1.0 + cs) / 2.0,
115            1.0 + alpha,
116            -2.0 * cs,
117            1.0 - alpha,
118        ),
119        Shape::Bandpass => (alpha, 0.0, -alpha, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
120        Shape::Notch => (1.0, -2.0 * cs, 1.0, 1.0 + alpha, -2.0 * cs, 1.0 - alpha),
121        Shape::Peaking => (
122            1.0 + alpha * a,
123            -2.0 * cs,
124            1.0 - alpha * a,
125            1.0 + alpha / a,
126            -2.0 * cs,
127            1.0 - alpha / a,
128        ),
129        Shape::Lowshelf => (
130            a * ((a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
131            2.0 * a * ((a - 1.0) - (a + 1.0) * cs),
132            a * ((a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
133            (a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
134            -2.0 * ((a - 1.0) + (a + 1.0) * cs),
135            (a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
136        ),
137        Shape::Highshelf => (
138            a * ((a + 1.0) + (a - 1.0) * cs + 2.0 * sqrt_a * alpha),
139            -2.0 * a * ((a - 1.0) + (a + 1.0) * cs),
140            a * ((a + 1.0) + (a - 1.0) * cs - 2.0 * sqrt_a * alpha),
141            (a + 1.0) - (a - 1.0) * cs + 2.0 * sqrt_a * alpha,
142            2.0 * ((a - 1.0) - (a + 1.0) * cs),
143            (a + 1.0) - (a - 1.0) * cs - 2.0 * sqrt_a * alpha,
144        ),
145    };
146
147    Coeffs {
148        b0: b0 / a0,
149        b1: b1 / a0,
150        b2: b2 / a0,
151        a1: a1 / a0,
152        a2: a2 / a0,
153    }
154}
155
156pub fn magnitude_db(c: &Coeffs, hz: f64, sr: f64) -> f64 {
157    let w = 2.0 * std::f64::consts::PI * hz / sr;
158    let (c1, s1) = (w.cos(), -w.sin());
159    let (c2, s2) = ((2.0 * w).cos(), -(2.0 * w).sin());
160    let num_re = c.b0 + c.b1 * c1 + c.b2 * c2;
161    let num_im = c.b1 * s1 + c.b2 * s2;
162    let den_re = 1.0 + c.a1 * c1 + c.a2 * c2;
163    let den_im = c.a1 * s1 + c.a2 * s2;
164    let num = (num_re * num_re + num_im * num_im).sqrt();
165    let den = (den_re * den_re + den_im * den_im).sqrt();
166    if den == 0.0 || num == 0.0 {
167        return f64::NEG_INFINITY;
168    }
169    20.0 * (num / den).log10()
170}
171
172pub const RESPONSE_RATIOS: [f64; 9] = [0.125, 0.25, 0.5, 0.707, 1.0, 1.414, 2.0, 4.0, 8.0];
173
174/// A point at or above Nyquist is dropped, not clamped: clamping reports the transform's own
175/// zero there as if the filter had put it there.
176pub fn response(c: &Coeffs, cutoff: f64, sr: f64) -> Vec<(f64, f64)> {
177    RESPONSE_RATIOS
178        .iter()
179        .map(|r| cutoff * r)
180        .filter(|hz| *hz < 0.5 * sr)
181        .map(|hz| (hz, magnitude_db(c, hz, sr)))
182        .collect()
183}
sva_samples/biquad.rs

sva_samples/
biquad.rs
