All filter stuff with wrappers. Added some convenience functions. Added tests

.gitignore

@@ -1,2 +1,4 @@
 /target
 /Cargo.lock
+__pycache__
+python/lasprs/_lasprs*

Cargo.toml

@@ -3,7 +3,7 @@ name = "lasprs"
 version = "0.1.0"
 edition = "2021"
 authors = ["J.A. de Jong <j.a.dejong@ascee.nl>"]
-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"]

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"

python/lasprs/__init__.py

@@ -0,0 +1 @@
+

python/lasprs/filter/__init__.py

@@ -0,0 +1,4 @@
+from .._lasprs import filter as _filter
+Biquad = _filter.Biquad
+SeriesBiquad = _filter.SeriesBiquad
+BiquadBank = _filter.BiquadBank

src/config.rs

@@ -12,7 +12,7 @@ pub const pi: Flt = std::f64::consts::PI;
 use num::complex::*;
 pub type Cflt = Complex<Flt>;
 
-use ndarray::{Array1, Array2};
+use numpy::ndarray::{Array1, Array2};
 pub type Vd = Vec<Flt>;
 pub type Vc = Vec<Cflt>;
 

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<dyn Filter> {
 }
 
 /// # 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<Flt>) -> PyResult<Self> {
+        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<Biquad> {
+        Ok(Biquad::firstOrderHighPass(fs, cuton_Hz)?)
+    }
+    #[pyo3(name = "filter")]
+    /// See: [Biquad::filter()]
+    pub fn filter_py<'py>(
+        &mut self,
+        py: Python<'py>,
+        input: PyArrayLike1<Flt>,
+    ) -> PyResult<&'py PyArray1<Flt>> {
+        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<Flt> {
         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<dyn Filter> { Box::new(self.clone()) }
+    fn clone_dyn(&self) -> Box<dyn Filter> {
+        Box::new(self.clone())
+    }
 }
 
 /// Series of biquads that filter sequentially on an input signal
+///
 /// # Examples
 ///
-/// ```
+/// See [tests]
-/// // Create filter coefficients for a single biquad that does output = input.
 /// ```
 #[derive(Clone, Debug)]
+#[cfg_attr(feature = "extension-module", pyclass)]
 pub struct SeriesBiquad {
     biqs: Vec<Biquad>,
 }
 
+#[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<Flt>) -> PyResult<Self> {
+        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<Flt>,
+    ) -> PyResult<&'py PyArray1<Flt>> {
+        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<dyn Filter> { Box::new(self.clone()) }
+    fn clone_dyn(&self) -> Box<dyn Filter> {
+        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<dyn Filter> { Box::new(self.clone()) }
+    fn clone_dyn(&self) -> Box<dyn Filter> {
+        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<Box<dyn Filter>>,
     gains: Vec<Flt>,
 }
 
+#[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<Flt>) -> PyResult<Self> {
+        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<Flt>,
+    ) -> PyResult<&'py PyArray1<Flt>> {
+        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<Flt>) -> 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<Flt>) -> 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<Box<dyn Filter>>) -> 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<dyn Filter> { Box::new(self.clone()) }
+    fn clone_dyn(&self) -> Box<dyn Filter> {
-
+        Box::new(self.clone())
+    }
 }
 
 #[cfg(test)]

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)]
+ #![warn(missing_docs)]
-#![allow(non_snake_case)]
+ #![allow(non_snake_case)]
-#![allow(non_upper_case_globals)]
+ #![allow(non_upper_case_globals)]
-#![allow(unused_imports)]
+ #![allow(unused_imports)]
-mod config;
+
-pub mod filter;
+ 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::Biquad>()?;
+    filter_module.add_class::<filter::SeriesBiquad>()?;
+    filter_module.add_class::<filter::BiquadBank>()?;
+    m.add_submodule(filter_module)?;
+
+    Ok(())
+}
 

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()

tools/watchdoc.sh

@@ -1,5 +1,11 @@
 #!/bin/bash
+
+# Required is cargo-watch and http:
 #
-testargs=""
+# ```bash
-testargs=-- --nocapture
+# $ cargo install cargo-watch cargo-docserver`
-cargo watch -s "cargo rustdoc --lib && cargo test ${testargs} && http target/doc"
+# ```
+#
+cargo watch -s "cargo rustdoc --lib && cargo test ${testargs} && cargo docserver"
+
+# Then open: ${browser} http://localhost:4000