From 3f58b194423d2f99b51df933ed2e7344aa0e6f59 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F." Date: Sat, 25 Nov 2023 14:58:20 +0100 Subject: [PATCH] All filter stuff with wrappers. Added some convenience functions. Added tests --- .gitignore | 2 + Cargo.toml | 18 ++- pyproject.toml | 19 +++ python/lasprs/__init__.py | 1 + python/lasprs/filter/__init__.py | 4 + src/config.rs | 2 +- src/filter.rs | 195 ++++++++++++++++++++++++++++--- src/lib.rs | 43 ++++++- test/test_filter.py | 62 ++++++++++ tools/watchdoc.sh | 12 +- 10 files changed, 330 insertions(+), 28 deletions(-) create mode 100644 pyproject.toml create mode 100644 python/lasprs/__init__.py create mode 100644 python/lasprs/filter/__init__.py create mode 100755 test/test_filter.py diff --git a/.gitignore b/.gitignore index 4fffb2f..09fa6be 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ /target /Cargo.lock +__pycache__ +python/lasprs/_lasprs* diff --git a/Cargo.toml b/Cargo.toml index 53dc9c5..21d1811 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ name = "lasprs" version = "0.1.0" edition = "2021" authors = ["J.A. de Jong "] -description = "Library for Acoustic Signal Processing (Rust edition)" +description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)" readme = "README.md" repository = "https://code.ascee.nl/ascee/lasprs" license = "MIT OR Apache-2.0" @@ -12,14 +12,26 @@ categories = [ "multimedia::audio", "scientific::dsp", "scientific::engineering"] +[lib] +name = "lasprs" +crate-type = ["cdylib", "lib"] [dependencies] anyhow = "1.0.75" -ndarray = "0.15.6" +pyo3 = { version = "0.20", features=["anyhow", "extension-module"]} + +# Optional future feature for ndarray: blas +ndarray = { version = "0.15.3", features = ["rayon"] } num = "0.4.1" rayon = "1.8.0" +numpy = { version = "0.20" } +# blas-src = { version = "0.8", features = ["openblas"] } +# openblas-src = { version = "0.10", features = ["cblas", "system"] } [features] -default = ["f64"] +# default = ["f64"] +# Use this for debugging extension +default = ["f64", "extension-module", "pyo3/extension-module"] f64 = [] f32 = [] +extension-module = ["pyo3/extension-module"] diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..dff3178 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,19 @@ +[build-system] +requires = ["maturin>=1.0,<2.0"] +build-backend = "maturin" + +[project] +name = "lasprs" +requires-python = ">=3.7" +classifiers = [ + "Programming Language :: Rust", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", +] + +dependencies = ["numpy"] + +[tool.maturin] +features = ["pyo3/extension-module", "extension-module"] +python-source = "python" +module-name = "lasprs._lasprs" diff --git a/python/lasprs/__init__.py b/python/lasprs/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/python/lasprs/__init__.py @@ -0,0 +1 @@ + diff --git a/python/lasprs/filter/__init__.py b/python/lasprs/filter/__init__.py new file mode 100644 index 0000000..61379f0 --- /dev/null +++ b/python/lasprs/filter/__init__.py @@ -0,0 +1,4 @@ +from .._lasprs import filter as _filter +Biquad = _filter.Biquad +SeriesBiquad = _filter.SeriesBiquad +BiquadBank = _filter.BiquadBank diff --git a/src/config.rs b/src/config.rs index 447e9a0..3fcc978 100644 --- a/src/config.rs +++ b/src/config.rs @@ -12,7 +12,7 @@ pub const pi: Flt = std::f64::consts::PI; use num::complex::*; pub type Cflt = Complex; -use ndarray::{Array1, Array2}; +use numpy::ndarray::{Array1, Array2}; pub type Vd = Vec; pub type Vc = Vec; diff --git a/src/filter.rs b/src/filter.rs index 5c0b495..fd0261b 100644 --- a/src/filter.rs +++ b/src/filter.rs @@ -1,9 +1,15 @@ //! # Filter implemententations for biquads, series of biquads and banks of series of biquads. //! -//! Contains `Biquad`, SeriesBiquad, and FilterBank +//! Contains [Biquad], [SeriesBiquad], and [BiquadBank]. These are all constructs that work on +//! blocks of input data, and apply filters on it. Todo: implement FIR filter. #![allow(non_snake_case)] use super::config::*; use anyhow::{bail, Result}; +use numpy::ndarray::{ArrayD, ArrayViewD, ArrayViewMutD}; +use numpy::{IntoPyArray, PyArray1, PyArrayDyn, PyArrayLike1, PyReadonlyArrayDyn}; +use pyo3::exceptions::PyValueError; +use pyo3::prelude::*; +use pyo3::{pymodule, types::PyModule, PyResult}; use rayon::prelude::*; pub trait Filter: Send { @@ -27,6 +33,7 @@ impl Clone for Box { } /// # A biquad is a second order recursive filter structure. +#[cfg_attr(feature = "extension-module", pyclass)] #[derive(Clone, Copy, Debug)] pub struct Biquad { // State parameters @@ -41,6 +48,37 @@ pub struct Biquad { a1: Flt, a2: Flt, } +#[cfg(feature = "extension-module")] +#[cfg_attr(feature = "extension-module", pymethods)] +impl Biquad { + #[new] + /// Create new biquad filter. See [Biquad::new()] + /// + pub fn new_py<'py>(coefs: PyReadonlyArrayDyn) -> PyResult { + Ok(Biquad::new(&coefs.as_slice()?)?) + } + #[pyo3(name = "unit")] + #[staticmethod] + /// See: [Biquad::unit()] + pub fn unit_py() -> Biquad { + Biquad::unit() + } + #[pyo3(name = "firstOrderHighPass")] + #[staticmethod] + /// See: [Biquad::firstOrderHighPass()] + pub fn firstOrderHighPass_py(fs: Flt, cuton_Hz: Flt) -> PyResult { + Ok(Biquad::firstOrderHighPass(fs, cuton_Hz)?) + } + #[pyo3(name = "filter")] + /// See: [Biquad::filter()] + pub fn filter_py<'py>( + &mut self, + py: Python<'py>, + input: PyArrayLike1, + ) -> PyResult<&'py PyArray1> { + Ok(self.filter(input.as_slice()?).into_pyarray(py)) + } +} impl Biquad { /// Create new biquad filter from given filter coeficients /// @@ -54,19 +92,18 @@ impl Biquad { if *a0 != 1.0 { bail!("Coefficient a0 should be equal to 1.0") } - Ok(Biquad { w1: 0., w2: 0., b0: *b0, b1: *b1, b2: *b2, a1: *a1, a2: *a2}) + Ok(Biquad { w1: 0., w2: 0., b0: *b0, b1: *b1, b2: *b2, a1: *a1, a2: *a2}) }, _ => bail!("Could not initialize biquad. Please make sure that the coefficients contain 6 terms") } } - /// Unit impulse response biquad. Input = output + /// Create unit impulse response biquad filter. Input = output fn unit() -> Biquad { let filter_coefs = &[1., 0., 0., 1., 0., 0.]; Biquad::new(filter_coefs).unwrap() } - /// Initialize biquad as first order high pass filter. /// /// @@ -110,34 +147,67 @@ impl Biquad { self.w1 = w0; inout[sample] = yn; } - println!("{:?}", inout); + // println!("{:?}", inout); } } impl Filter for Biquad { fn filter(&mut self, input: &[Flt]) -> Vec { let mut out = input.to_vec(); self.filter_inout(&mut out); - println!("{:?}", out); + // println!("{:?}", out); out } fn reset(&mut self) { self.w1 = 0.; self.w2 = 0.; } - fn clone_dyn(&self)-> Box { Box::new(self.clone()) } + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } } /// Series of biquads that filter sequentially on an input signal +/// /// # Examples /// -/// ``` -/// // Create filter coefficients for a single biquad that does output = input. +/// See [tests] /// ``` #[derive(Clone, Debug)] +#[cfg_attr(feature = "extension-module", pyclass)] pub struct SeriesBiquad { biqs: Vec, } +#[cfg(feature = "extension-module")] +#[cfg_attr(feature = "extension-module", pymethods)] +impl SeriesBiquad { + #[new] + /// Create new series filter set. See [SeriesBiquad::new()] + /// + pub fn new_py<'py>(coefs: PyReadonlyArrayDyn) -> PyResult { + Ok(SeriesBiquad::new(&coefs.as_slice()?)?) + } + #[pyo3(name = "unit")] + #[staticmethod] + /// See: [Biquad::unit()] + pub fn unit_py() -> SeriesBiquad { + SeriesBiquad::unit() + } + #[pyo3(name = "filter")] + /// See: [SeriesBiquad::filter()] + pub fn filter_py<'py>( + &mut self, + py: Python<'py>, + input: PyArrayLike1, + ) -> PyResult<&'py PyArray1> { + Ok(self.filter(input.as_slice()?).into_pyarray(py)) + } + #[pyo3(name = "reset")] + /// See: [SeriesBiquad::reset()] + pub fn reset_py(&mut self) { + self.reset(); + } +} impl SeriesBiquad { /// Create a new series biquad, having an arbitrary number of biquads. /// @@ -174,7 +244,9 @@ impl SeriesBiquad { let filter_coefs = &[1., 0., 0., 1., 0., 0.]; SeriesBiquad::new(filter_coefs).unwrap() } - fn clone_dyn(&self)-> Box { Box::new(self.clone()) } + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } } impl Filter for SeriesBiquad { @@ -191,9 +263,12 @@ impl Filter for SeriesBiquad { fn reset(&mut self) { self.biqs.iter_mut().for_each(|f| f.reset()); } - fn clone_dyn(&self)-> Box { Box::new(self.clone()) } + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } } +#[cfg_attr(feature = "extension-module", pyclass)] #[derive(Clone)] /// Multiple biquad filter that operate in parallel on a signal, and can apply a gain value to each /// of the returned values. The BiquadBank can be used to decompose a signal by running it through @@ -202,28 +277,115 @@ impl Filter for SeriesBiquad { /// /// # Detailed description /// +/// Below is an example of the signal flow is for the case of three SeriesBiquad filters, `h1`, +/// `h2` and `h3`: /// +/// ```markdown /// +/// analysis() gain application sum() +/// +/// +------+ +-----+ +------------+ +/// +-----+ h1 |---+ g1 +-------------+ | | +/// | +------+ +-----+ ++ | +------ | +/// | +-----| ++ | +/// Input | +------+ +-----+ | ++ | Output of filter() +///--------|>--+-----+ h2 |---+ g2 |--------------------| +-+ +----------------|> +/// | ---+---+ +-----+ | +-+ | +/// | +-----| ++ | +/// | +------+ +-----+ ++ | ++ | +/// +-----+ h3 |---+ g3 |-------------+ | +------ | +/// +------+ +-----+ | | +/// | +------------+ +/// | +/// | Output of analysis() method (optional) +/// +-------------|> +/// ``` pub struct BiquadBank { biqs: Vec>, gains: Vec, } +#[cfg(feature = "extension-module")] +#[cfg_attr(feature = "extension-module", pymethods)] +/// Methods to wrap it in Python impl BiquadBank { + #[new] + /// Create new biquadbank filter set. See [BiquadBank::new()] + /// + pub fn new_py<'py>(coefs: PyReadonlyArrayDyn) -> PyResult { + let mut filts = vec![]; + for col in coefs.as_array().columns() { + match col.as_slice() { + Some(colslice) => { + let new_ser = SeriesBiquad::new(colslice)?; + filts.push(new_ser.clone_dyn()); + } + None => { + return Err(PyValueError::new_err("Error generating column")); + } + } + } + Ok(BiquadBank::new(filts)) + } + #[pyo3(name = "filter")] + /// See: [BiquadBank::filter()] + pub fn filter_py<'py>( + &mut self, + py: Python<'py>, + input: PyArrayLike1, + ) -> PyResult<&'py PyArray1> { + Ok(self.filter(input.as_slice()?).into_pyarray(py)) + } + #[pyo3(name = "reset")] + /// See: [BiquadBank::reset()] + pub fn reset_py(&mut self) { + self.reset(); + } + #[pyo3(name = "set_gains")] + /// See: [BiquadBank::set_gains()] + pub fn set_gains_py<'py>(&mut self, py: Python<'py>, gains: PyArrayLike1) -> PyResult<()> { + if gains.len() != self.len() { + return Err(PyValueError::new_err("Invalid number of provided gains")); + } + self.set_gains(gains.as_slice()?); + Ok(()) + } + #[pyo3(name = "set_gains_dB")] + /// See: [BiquadBank::set_gains_dB()] + pub fn set_gains_dB_py<'py>(&mut self, py: Python<'py>, gains_dB: PyArrayLike1) -> PyResult<()> { + if gains_dB.len() != self.len() { + return Err(PyValueError::new_err("Invalid number of provided gains")); + } + self.set_gains_dB(gains_dB.as_slice()?); + Ok(()) + + } + #[pyo3(name = "len")] + /// See: [BiquadBank::len()] + pub fn len_py(&self) -> usize { + self.len() + } +} + +impl BiquadBank { /// Create new biquad bank. Initialized from given vector of series biquads. - // #[new] pub fn new(biqs: Vec>) -> BiquadBank { let gains = vec![1.0; biqs.len()]; BiquadBank { biqs, gains } } + /// Return the number of parallel filters installed. + pub fn len(&self) -> usize { + self.biqs.len() + } + /// Set the gains for each of the biquad. The gains are not used in the analyisis phase, but in /// the reconstruction phase, so when BiquadBank::filter() is run on an input signal. /// /// # Panics /// /// When gains_dB.len() != to the number of filters. - pub fn set_gains_dB(&mut self, gains_dB: Vd) { + pub fn set_gains_dB(&mut self, gains_dB: &[Flt]) { if gains_dB.len() != self.gains.len() { panic!("Invalid gains size!"); } @@ -232,6 +394,8 @@ impl BiquadBank { .zip(gains_dB) .for_each(|(g, gdB)| *g = Flt::powf(10., gdB / 20.)); } + /// Set linear gain values for each biquad. Same comments hold as for + /// [BiquadBank::set_gains_dB()]. pub fn set_gains(&mut self, gains: &[Flt]) { if gains.len() != self.gains.len() { panic!("Invalid gains size!"); @@ -273,8 +437,9 @@ impl Filter for BiquadBank { fn reset(&mut self) { self.biqs.iter_mut().for_each(|b| b.reset()); } - fn clone_dyn(&self)-> Box { Box::new(self.clone()) } - + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } } #[cfg(test)] diff --git a/src/lib.rs b/src/lib.rs index b4ef72b..7296f08 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,9 +1,40 @@ //! # Library for acoustic signal processing +//! +//! This crate contains structures and functions to perform acoustic measurements, interact with +//! data acquisition devices and apply common acoustic analysis operations on them. -#![warn(missing_docs)] -#![allow(non_snake_case)] -#![allow(non_upper_case_globals)] -#![allow(unused_imports)] -mod config; -pub mod filter; + #![warn(missing_docs)] + #![allow(non_snake_case)] + #![allow(non_upper_case_globals)] + #![allow(unused_imports)] + + mod config; + pub mod filter; + + extern crate pyo3; + #[cfg(feature = "extension-module")] + use pyo3::prelude::*; + + /// A Python module implemented in Rust. + #[cfg(feature = "extension-module")] + #[pymodule] + #[pyo3(name="_lasprs")] + fn lasprs(py: Python, m: &PyModule) -> PyResult<()> { + + pyo3_add_submodule_filter(py, &m)?; + Ok(()) + } + +/// Add filter submodule to extension +#[cfg(feature = "extension-module")] +fn pyo3_add_submodule_filter(py: Python, m: &PyModule) -> PyResult<()> { + // Add filter submodule + let filter_module = PyModule::new(py, "filter")?; + filter_module.add_class::()?; + filter_module.add_class::()?; + filter_module.add_class::()?; + m.add_submodule(filter_module)?; + + Ok(()) +} diff --git a/test/test_filter.py b/test/test_filter.py new file mode 100755 index 0000000..6e42b84 --- /dev/null +++ b/test/test_filter.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +import numpy as np +from lasprs.filter import Biquad, SeriesBiquad, BiquadBank + + +def test_biquad1(): + """ + Test if manual unit biquad matches response for manual input biquad. + """ + input_ = np.array([1.0, 0, 0, 0, 0, 0, 0]) + + b = Biquad(np.array([1., 0, 0, 1, 0, 0])) + b2 = Biquad.unit() + out1 = b.filter(input_) + out2 = b2.filter(input_) + + expected_output = input_ + assert(np.linalg.norm(out1 - out2) == 0.) + assert(np.linalg.norm(out1 - expected_output) == 0.) + +def test_seriesbiquad(): + """ + Test if manual unit biquad matches response for manual input biquad. + """ + input_ = np.array([1.0, 0, 0, 0, 0, 0, 0]) + + f1 = [1., 0, 0, 1, 0, 0] + f2 = [1., 0, 0, 1, 0, 0] + f = f1+f2 + + # Two biquads in series + b = SeriesBiquad(np.array(f)) + + # Single one + b2 = SeriesBiquad.unit() + + out1 = b.filter(input_) + out2 = b2.filter(input_) + + expected_output = input_ + assert(np.linalg.norm(out1 - out2) == 0.) + assert(np.linalg.norm(out1 - expected_output) == 0.) + +def test_biquadbank(): + """ + See if two filters + """ + input_ = np.array([1.0, 0, 0, 0, 0, 0, 0]) + f1 = [1., 0, 0, 1, 0, 0] + f2 = [1., 0, 0, 1, 0, 0] + filters = np.array([f1, f2]).T + bank = BiquadBank(filters) + bank.set_gains([0.5, 0.5]) + filtered = bank.filter(input_) + assert(bank.len() == 2) + assert(np.linalg.norm(input_ - filtered) == 0) + + +if __name__ == '__main__': + test_biquad1() + test_seriesbiquad() + test_biquadbank() diff --git a/tools/watchdoc.sh b/tools/watchdoc.sh index baaa78f..8224d0a 100755 --- a/tools/watchdoc.sh +++ b/tools/watchdoc.sh @@ -1,5 +1,11 @@ #!/bin/bash + +# Required is cargo-watch and http: # -testargs="" -testargs=-- --nocapture -cargo watch -s "cargo rustdoc --lib && cargo test ${testargs} && http target/doc" +# ```bash +# $ cargo install cargo-watch cargo-docserver` +# ``` +# +cargo watch -s "cargo rustdoc --lib && cargo test ${testargs} && cargo docserver" + +# Then open: ${browser} http://localhost:4000