idsp/lowpass.rs
1use dsp_process::{SplitInplace, SplitProcess};
2
3/// Arbitrary order, high dynamic range, wide coefficient range,
4/// lowpass filter implementation. DC gain is 1.
5///
6/// Type argument N is the filter order. N must be `1` or `2`.
7///
8/// The filter will cleanly saturate towards the `i32` range.
9///
10/// Both filters have been optimized for accuracy, dynamic range, and
11/// speed on Cortex-M7.
12#[derive(Clone, Debug)]
13pub struct Lowpass<const N: usize>(pub [i32; N]);
14
15/// Lowpass filter state
16#[derive(Clone, Debug)]
17pub struct LowpassState<const N: usize>(pub [i64; N]);
18
19impl<const N: usize> Default for LowpassState<N>
20where
21 [i64; N]: Default,
22{
23 fn default() -> Self {
24 Self(Default::default())
25 }
26}
27
28impl<const N: usize> SplitProcess<i32, i32, LowpassState<N>> for Lowpass<N> {
29 /// The filter configuration `Config` contains the filter gains.
30 ///
31 /// For the first-order lowpass this is a single element array `[k]` with
32 /// the corner frequency in scaled Q31:
33 /// `k = pi*(1 << 31)*f0/fn` where
34 /// `f0` is the 3dB corner frequency and
35 /// `fn` is the Nyquist frequency.
36 /// The corner frequency is warped in the usual way.
37 ///
38 /// For the second-order lowpass this is `[k**2/(1 << 32), -k/q]` with `q = 1/sqrt(2)`
39 /// for a Butterworth response.
40 ///
41 /// In addition to the poles at the corner frequency the filters have zeros at Nyquist.
42 ///
43 /// The first-order lowpass works fine and accurate for any positive gain
44 /// `1 <= k <= (1 << 31) - 1`.
45 /// The second-order lowpass works and is accurate for
46 /// `1 << 16 <= k <= q*(1 << 31)`.
47 fn process(&self, state: &mut LowpassState<N>, x: i32) -> i32 {
48 // d = (x0 - p1)*k0
49 // p0 = p1 + 2d
50 // y0 = p1 + d
51 //
52 // d = (x0 - p1)*k0 + q1*k1
53 // q0 = q1 + 2d
54 // p0 = p1 + 2q1 + 2d
55 // y0 = p1 + q1 + d
56 let mut d = x.saturating_sub((state.0[0] >> 32) as i32) as i64 * self.0[0] as i64;
57 let y;
58 if N == 1 {
59 state.0[0] += d;
60 y = (state.0[0] >> 32) as i32;
61 state.0[0] += d;
62 } else if N == 2 {
63 d += (state.0[1] >> 32) * self.0[1] as i64;
64 state.0[1] += d;
65 state.0[0] += state.0[1];
66 y = (state.0[0] >> 32) as i32;
67 // This creates the double Nyquist zero,
68 // compensates the gain lost in the signed i32 as (i32 as i64)*(i64 >> 32)
69 // multiplication while keeping the lowest bit significant, and
70 // copes better with wrap-around than Nyquist averaging.
71 state.0[0] += state.0[1];
72 state.0[1] += d;
73 } else {
74 unimplemented!()
75 }
76 y
77 }
78}
79
80impl<const N: usize> SplitInplace<i32, LowpassState<N>> for Lowpass<N> {}
81
82impl<const N: usize> Default for Lowpass<N> {
83 fn default() -> Self {
84 Self([0; N])
85 }
86}
87
88/// First order lowpass
89pub type Lowpass1 = Lowpass<1>;
90/// Second order lowpass
91pub type Lowpass2 = Lowpass<2>;