Skip to main content

embedded_dsp/
filter_design.rs

1//! Filter design routines for calculating biquad IIR coefficients (Low-pass, High-pass, Band-pass, Notch, Peaking, All-pass, Butterworth).
2
3#[allow(unused_imports)]
4use crate::math::FloatMath;
5
6/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for a Low-Pass Filter.
7///
8/// `cutoff_freq`: Cutoff frequency in Hz.
9/// `sample_rate`: Sampling rate in Hz.
10/// `q`: Quality factor (e.g. 0.7071 for Butterworth alignment).
11pub fn biquad_lowpass_coeffs(cutoff_freq: f32, sample_rate: f32, q: f32) -> [f32; 5] {
12    let w0 = 2.0 * core::f32::consts::PI * cutoff_freq / sample_rate;
13    let cos_w0 = w0.cos();
14    let sin_w0 = w0.sin();
15    let alpha = sin_w0 / (2.0 * q);
16
17    let a0 = 1.0 + alpha;
18    let b0 = (1.0 - cos_w0) / 2.0 / a0;
19    let b1 = (1.0 - cos_w0) / a0;
20    let b2 = (1.0 - cos_w0) / 2.0 / a0;
21    // In Direct Form I (out = b0*x + b1*x1 + b2*x2 + a1*y1 + a2*y2), sign of feedback terms is flipped:
22    let a1 = (2.0 * cos_w0) / a0;
23    let a2 = -(1.0 - alpha) / a0;
24
25    [b0, b1, b2, a1, a2]
26}
27
28/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for a High-Pass Filter.
29pub fn biquad_highpass_coeffs(cutoff_freq: f32, sample_rate: f32, q: f32) -> [f32; 5] {
30    let w0 = 2.0 * core::f32::consts::PI * cutoff_freq / sample_rate;
31    let cos_w0 = w0.cos();
32    let sin_w0 = w0.sin();
33    let alpha = sin_w0 / (2.0 * q);
34
35    let a0 = 1.0 + alpha;
36    let b0 = (1.0 + cos_w0) / 2.0 / a0;
37    let b1 = -(1.0 + cos_w0) / a0;
38    let b2 = (1.0 + cos_w0) / 2.0 / a0;
39    let a1 = (2.0 * cos_w0) / a0;
40    let a2 = -(1.0 - alpha) / a0;
41
42    [b0, b1, b2, a1, a2]
43}
44
45/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for a Band-Pass Filter (constant skirt gain).
46pub fn biquad_bandpass_coeffs(center_freq: f32, sample_rate: f32, q: f32) -> [f32; 5] {
47    let w0 = 2.0 * core::f32::consts::PI * center_freq / sample_rate;
48    let cos_w0 = w0.cos();
49    let sin_w0 = w0.sin();
50    let alpha = sin_w0 / (2.0 * q);
51
52    let a0 = 1.0 + alpha;
53    let b0 = alpha / a0;
54    let b1 = 0.0;
55    let b2 = -alpha / a0;
56    let a1 = (2.0 * cos_w0) / a0;
57    let a2 = -(1.0 - alpha) / a0;
58
59    [b0, b1, b2, a1, a2]
60}
61
62/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for a Notch (Band-Stop) Filter.
63pub fn biquad_notch_coeffs(center_freq: f32, sample_rate: f32, q: f32) -> [f32; 5] {
64    let w0 = 2.0 * core::f32::consts::PI * center_freq / sample_rate;
65    let cos_w0 = w0.cos();
66    let sin_w0 = w0.sin();
67    let alpha = sin_w0 / (2.0 * q);
68
69    let a0 = 1.0 + alpha;
70    let b0 = 1.0 / a0;
71    let b1 = (-2.0 * cos_w0) / a0;
72    let b2 = 1.0 / a0;
73    let a1 = (2.0 * cos_w0) / a0;
74    let a2 = -(1.0 - alpha) / a0;
75
76    [b0, b1, b2, a1, a2]
77}
78
79/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for a Peaking EQ Filter.
80pub fn biquad_peaking_coeffs(center_freq: f32, sample_rate: f32, q: f32, gain_db: f32) -> [f32; 5] {
81    let w0 = 2.0 * core::f32::consts::PI * center_freq / sample_rate;
82    let cos_w0 = w0.cos();
83    let sin_w0 = w0.sin();
84    let a = (10.0f32).powf(gain_db / 40.0);
85    let alpha = sin_w0 / (2.0 * q);
86
87    let a0 = 1.0 + alpha / a;
88    let b0 = (1.0 + alpha * a) / a0;
89    let b1 = (-2.0 * cos_w0) / a0;
90    let b2 = (1.0 - alpha * a) / a0;
91    let a1 = (2.0 * cos_w0) / a0;
92    let a2 = -(1.0 - alpha / a) / a0;
93
94    [b0, b1, b2, a1, a2]
95}
96
97/// Computes Direct Form I Biquad coefficients `[b0, b1, b2, a1, a2]` for an All-Pass Filter.
98pub fn biquad_allpass_coeffs(center_freq: f32, sample_rate: f32, q: f32) -> [f32; 5] {
99    let w0 = 2.0 * core::f32::consts::PI * center_freq / sample_rate;
100    let cos_w0 = w0.cos();
101    let sin_w0 = w0.sin();
102    let alpha = sin_w0 / (2.0 * q);
103
104    let a0 = 1.0 + alpha;
105    let b0 = (1.0 - alpha) / a0;
106    let b1 = (-2.0 * cos_w0) / a0;
107    let b2 = (1.0 + alpha) / a0;
108    let a1 = (2.0 * cos_w0) / a0;
109    let a2 = -(1.0 - alpha) / a0;
110
111    [b0, b1, b2, a1, a2]
112}
113
114/// Calculates multi-stage Butterworth Low-Pass filter biquad coefficients.
115/// `out_coeffs` must be a slice of size `5 * (order / 2)`.
116pub fn butterworth_lowpass_biquads(
117    cutoff_freq: f32,
118    sample_rate: f32,
119    order: usize,
120    out_coeffs: &mut [f32],
121) {
122    let num_stages = order / 2;
123    assert!(
124        out_coeffs.len() >= num_stages * 5,
125        "out_coeffs buffer too small"
126    );
127
128    for k in 0..num_stages {
129        let angle = core::f32::consts::PI * (2 * k + 1) as f32 / (2 * order) as f32;
130        let q = 1.0 / (2.0 * angle.sin());
131        let coeffs = biquad_lowpass_coeffs(cutoff_freq, sample_rate, q);
132        out_coeffs[k * 5..(k + 1) * 5].copy_from_slice(&coeffs);
133    }
134}