dspk_ops/filter/design/
mod.rs

1//! Host-side (f64) IIR filter design.
2//!
3//! Filters are described by a [`FilterSpec`] (what to build and how to run it) and designed into a
4//! cascade of second-order sections ([`Sos`]) that the SOS kernel executes. Designs are equivalent
5//! to `scipy.signal.butter(..., output="sos")` and `scipy.signal.iirnotch`.
6
7mod butterworth;
8mod notch;
9mod sos;
10
11pub use butterworth::{bessel_sos, butterworth_sos, chebyshev1_sos};
12pub use notch::notch_sos;
13pub use sos::{Section, Sos, DEFAULT_SETTLING_TOLERANCE};
14
15pub use crate::error::OpsError;
16
17/// Frequency band of a Butterworth filter, cutoffs in Hz.
18#[derive(Debug, Clone, Copy, PartialEq)]
19pub enum FilterBand {
20    Lowpass(f64),
21    Highpass(f64),
22    Bandpass(f64, f64),
23    Bandstop(f64, f64),
24}
25
26/// How a filter is applied along time.
27#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
28pub enum FilterMode {
29    /// Causal, one pass (`sosfilt`). Usable on live streams with carried state.
30    Forward,
31    /// Zero phase: forward then backward pass (`sosfiltfilt`); squares the magnitude response.
32    /// Needs future samples, so chunked use requires a right halo.
33    #[default]
34    ForwardBackward,
35}
36
37/// Where a [`FilterMode::Forward`] pass starts. Forward-backward passes always start from the steady
38/// state of their first sample, as `sosfiltfilt` does.
39#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
40pub enum FilterStart {
41    /// As if the input had been zero before the first sample (`sosfilt` without `zi`): a DC level
42    /// produces a step transient.
43    #[default]
44    Rest,
45    /// From the steady state of the first sample (`sosfilt` with `zi = sosfilt_zi · x[0]`): no step
46    /// transient, so halo windows over DC-heavy data settle with less context.
47    SteadyState,
48}
49
50/// Filter design.
51#[derive(Debug, Clone, PartialEq)]
52pub enum FilterDesign {
53    /// Butterworth of any order (per edge for band filters, as in scipy).
54    Butterworth { order: usize, band: FilterBand },
55    /// Chebyshev type I of any order with `ripple_db` of pass-band ripple (`cheby1`).
56    Chebyshev1 { order: usize, ripple_db: f64, band: FilterBand },
57    /// Bessel of any order (per edge for band filters), phase-normalised (`scipy.signal.bessel`).
58    Bessel { order: usize, band: FilterBand },
59    /// Second-order notch at `freq_hz` with quality factor `q` (`iirnotch`).
60    Notch { freq_hz: f64, q: f64 },
61    /// Explicit second-order sections, already designed for the target sample rate.
62    Sos(Sos),
63}
64
65/// A filter to design, the direction(s) it runs in, and where forward passes start.
66#[derive(Debug, Clone, PartialEq)]
67pub struct FilterSpec {
68    pub design: FilterDesign,
69    pub mode: FilterMode,
70    pub start: FilterStart,
71}
72
73/// Default Butterworth order, matching SpikeInterface's `bandpass_filter` / `highpass_filter`.
74pub const DEFAULT_BUTTERWORTH_ORDER: usize = 5;
75
76impl FilterSpec {
77    pub fn butterworth(order: usize, band: FilterBand) -> Self {
78        Self { design: FilterDesign::Butterworth { order, band }, mode: FilterMode::default(), start: FilterStart::default() }
79    }
80
81    pub fn bandpass(low_hz: f64, high_hz: f64) -> Self {
82        Self::butterworth(DEFAULT_BUTTERWORTH_ORDER, FilterBand::Bandpass(low_hz, high_hz))
83    }
84
85    pub fn highpass(cutoff_hz: f64) -> Self {
86        Self::butterworth(DEFAULT_BUTTERWORTH_ORDER, FilterBand::Highpass(cutoff_hz))
87    }
88
89    pub fn lowpass(cutoff_hz: f64) -> Self {
90        Self::butterworth(DEFAULT_BUTTERWORTH_ORDER, FilterBand::Lowpass(cutoff_hz))
91    }
92
93    pub fn bandstop(low_hz: f64, high_hz: f64) -> Self {
94        Self::butterworth(DEFAULT_BUTTERWORTH_ORDER, FilterBand::Bandstop(low_hz, high_hz))
95    }
96
97    pub fn chebyshev1(order: usize, ripple_db: f64, band: FilterBand) -> Self {
98        Self { design: FilterDesign::Chebyshev1 { order, ripple_db, band }, mode: FilterMode::default(), start: FilterStart::default() }
99    }
100
101    /// Bessel of `order` per edge (phase-normalised, as `scipy.signal.bessel`), zero phase.
102    pub fn bessel(order: usize, band: FilterBand) -> Self {
103        Self { design: FilterDesign::Bessel { order, band }, mode: FilterMode::default(), start: FilterStart::default() }
104    }
105
106    pub fn notch(freq_hz: f64, q: f64) -> Self {
107        Self { design: FilterDesign::Notch { freq_hz, q }, mode: FilterMode::default(), start: FilterStart::default() }
108    }
109
110    pub fn sos(sos: Sos) -> Self {
111        Self { design: FilterDesign::Sos(sos), mode: FilterMode::default(), start: FilterStart::default() }
112    }
113
114    /// Returns the same filter with `order` (Butterworth / Chebyshev; ignored for notch / explicit SOS).
115    pub fn with_order(mut self, new_order: usize) -> Self {
116        if let FilterDesign::Butterworth { order, .. } | FilterDesign::Chebyshev1 { order, .. } | FilterDesign::Bessel { order, .. } = &mut self.design {
117            *order = new_order;
118        }
119        self
120    }
121
122    pub fn with_mode(mut self, mode: FilterMode) -> Self {
123        self.mode = mode;
124        self
125    }
126
127    /// Returns the same filter starting forward passes at `start`.
128    pub fn with_start(mut self, start: FilterStart) -> Self {
129        self.start = start;
130        self
131    }
132
133    /// Designs the second-order sections for `sample_rate` Hz, validating the parameters.
134    ///
135    /// # Errors
136    ///
137    /// [`OpsError`] when the design is invalid at `sample_rate`: non-positive sample rate, zero
138    /// order, a cutoff outside `(0, Nyquist)`, band edges out of order, a bad `q` or ripple.
139    /// Explicit sections fail [`Sos::validate`].
140    pub fn design(&self, sample_rate: f64) -> Result<Sos, OpsError> {
141        match &self.design {
142            FilterDesign::Butterworth { order, band } => butterworth_sos(*order, *band, sample_rate),
143            FilterDesign::Chebyshev1 { order, ripple_db, band } => chebyshev1_sos(*order, *ripple_db, *band, sample_rate),
144            FilterDesign::Bessel { order, band } => bessel_sos(*order, *band, sample_rate),
145            FilterDesign::Notch { freq_hz, q } => notch_sos(*freq_hz, *q, sample_rate),
146            FilterDesign::Sos(sos) => {
147                sos.validate()?;
148                Ok(sos.clone())
149            }
150        }
151    }
152
153    /// `(left, right)` samples needed around a chunk for its interior to match whole-recording
154    /// filtering (impulse response decayed below [`DEFAULT_SETTLING_TOLERANCE`]).
155    ///
156    /// # Errors
157    ///
158    /// As [`Self::design`].
159    pub fn settling(&self, sample_rate: f64) -> Result<(usize, usize), OpsError> {
160        let s = self.design(sample_rate)?.settling_samples(DEFAULT_SETTLING_TOLERANCE);
161        Ok(match self.mode {
162            FilterMode::Forward => (s, 0),
163            FilterMode::ForwardBackward => (s, s),
164        })
165    }
166}
167
168pub(crate) fn check_sample_rate(fs: f64) -> Result<(), OpsError> {
169    if fs.is_finite() && fs > 0.0 { Ok(()) } else { Err(OpsError::InvalidSampleRate(fs)) }
170}
171
172pub(crate) fn check_cutoff(cutoff_hz: f64, fs: f64) -> Result<(), OpsError> {
173    let nyquist_hz = fs / 2.0;
174    if cutoff_hz.is_finite() && cutoff_hz > 0.0 && cutoff_hz < nyquist_hz {
175        Ok(())
176    } else {
177        Err(OpsError::CutoffOutOfRange { cutoff_hz, nyquist_hz })
178    }
179}
dspk_ops/filter/design/mod.rs

dspk_ops/filter/design/
mod.rs
