Rolled back to older Pyo3 version to make it compatible with Rust-Numpy that is not yet ready to go to new Pyo3 API. Improved power spectra estimator by removing DC value before going into window. Later overwriting it. Improved functionality of Python bindings

Cargo.toml

@@ -12,7 +12,7 @@ categories = ["multimedia::audio", "science", "mathematics"]
 
 [lib]
 name = "lasprs"
-crate-type = ["cdylib", "rlib",]
+crate-type = ["cdylib", "rlib"]
 
 [dependencies]
 # Error handling
@@ -29,8 +29,12 @@ num = "0.4.3"
 rayon = "1.10.0"
 
 # Python bindings
-pyo3 = { version = "0.22.2", optional = true, features = ["extension-module", "anyhow"]}
+pyo3 = { version = "0.21.2", optional = true, features = [
-numpy = { git = "https://github.com/JRRudy1/rust-numpy", branch = "pyo3-0.22.0" , optional = true}
+    "extension-module",
+    "anyhow",
+] }
+# Python bindings for Numpy arrays
+numpy = { version = "0.21.0", optional = true }
 
 # White noise etc
 rand = "0.8.5"
@@ -74,12 +78,17 @@ itertools = "0.13.0"
 approx = "0.5.1"
 
 # For getting timestamps. Only useful when recording.
-chrono = {version = "0.4.38", optional = true}
+chrono = { version = "0.4.38", optional = true }
 # For getting UUIDs in recording
-uuid = { version = "1.10.0", features = ["v4"] , optional = true}
+uuid = { version = "1.10.0", features = ["v4"], optional = true }
 
 # Command line argument parser, for CLI apps
-clap = { version = "4.5.13", features = ["derive", "color", "help", "suggestions"] }
+clap = { version = "4.5.13", features = [
+    "derive",
+    "color",
+    "help",
+    "suggestions",
+] }
 
 # FFT's
 realfft = "3.3.0"
@@ -104,8 +113,21 @@ default = ["f64", "cpal-api", "record"]
 # Use this to test if everything works well in f32
 # default = ["f32", "cpal-api", "record"]
 
-# Use this for debugging extensions
+# When building Python bindings of Lasprs, you should enable the feature flag
-# default = ["f64", "python-bindings", "record", "cpal-api"]
+# "python-bindings". When debugging these, to safe the flag specification you
+# could just uncomment the line below here and comment the one above. If only
+# testing the implementations, or building the extension module, you should
+# call:
+
+# $ maturin develop -F python-bindings
+
+# Or: 
+
+# $ maturin develop --release -F python-bindings
+
+# For faster code
+
+#default = ["f64", "python-bindings", "record", "cpal-api"]
 
 pulse-api = []
 cpal-api = ["dep:cpal"]

README.md

@@ -12,11 +12,12 @@ and processing of (multi) sensor data in real time on a PC and output results.
 
 Documentation is provided at [doc.rs](https://docs.rs/lasprs/latest/lasprs).
 
-## Python bindings
+
+## Python bindings and examples
 
 The library has Python bindings (via [pyo3](https://pyo3.rs), which can be installed via:
 
 ```
-$ pip install git+https://code.ascee.nl/ascee/lasprs --install-option "--all-features"
+$ pip install git+https://code.ascee.nl/ascee/lasprs --install-option "python-bindings"
 ```
 

src/daq/api/mod.rs

@@ -32,7 +32,10 @@ pub trait Stream {
 }
 
 /// Stream API descriptor: type and corresponding text
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+//#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(strum_macros::EnumMessage, Debug, Clone, PartialEq, Serialize, Deserialize, strum_macros::Display)]
 #[allow(dead_code)]
 pub enum StreamApiDescr {

src/daq/datatype.rs

@@ -7,8 +7,11 @@ use crate::config::*;
 
 /// Data type description for samples coming from a stream
 #[derive(strum_macros::EnumMessage, PartialEq, Copy, Debug, Clone, Serialize, Deserialize)]
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+//#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[allow(dead_code)]
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
 pub enum DataType {
     /// 32-bit floats
     #[strum(message = "F32", detailed_message = "32-bits floating points")]

src/daq/mod.rs

@@ -50,7 +50,10 @@ use api::*;
 use crate::config::*;
 /// Stream types that can be started
 ///
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+// #[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(PartialEq, Clone, Copy)]
 pub enum StreamType {
     /// Input-only stream

src/daq/qty.rs

@@ -7,8 +7,10 @@ use strum_macros;
 use serde::{Serialize, Deserialize};
 
 /// Physical quantities that are I/O of a Daq device.
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+// #[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(PartialEq, Serialize, Deserialize, strum_macros::EnumMessage, Debug, Clone, Copy)]
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
 #[allow(dead_code)]
 pub enum Qty {
     /// Number

src/daq/streamerror.rs

@@ -3,7 +3,12 @@ use crate::config::*;
 
 /// Errors that happen in a stream
 #[derive(strum_macros::EnumMessage, PartialEq, Debug, Clone, Display, Copy)]
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+//#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
+
 pub enum StreamError {
     /// Input overrun
     #[strum(

src/filter/octave.rs

@@ -478,6 +478,7 @@ impl BandDescriptor for ThirdOctaveBandDescriptor {
     }
 }
 
+#[cfg(feature = "python-bindings")]
 #[cfg_attr(feature = "python-bindings", pymethods)]
 impl StandardFilterDescriptor {
     #[pyo3(name = "genFilter")]
@@ -511,7 +512,7 @@ impl StandardFilterDescriptor {
     }
 
     fn __repr__(&self) -> String {
-        format!{"{:#?}", self}
+        format! {"{:#?}", self}
     }
 
     fn __str__(&self) -> String {

src/ps/fft.rs

@@ -34,12 +34,14 @@ impl FFT {
             nfftF: nfft as Flt,
         }
     }
-    pub fn process<'a, T>(&mut self, time: &[Flt], freq: T)
+    pub fn process<'a, T, U>(&mut self, time: T, freq: U)
     where
-        T: Into<ArrayViewMut<'a, Cflt, Ix1>>,
+        T: Into<ArrayView<'a, Flt, Ix1>>,
+        U: Into<ArrayViewMut<'a, Cflt, Ix1>>,
     {
         let mut freq = freq.into();
-        self.timescratch.copy_from_slice(time);
+        let time = time.into();
+        self.timescratch.copy_from_slice(time.as_slice().unwrap());
         let _ = self
             .fft
             .process(&mut self.timescratch, freq.as_slice_mut().unwrap());
@@ -49,4 +51,4 @@ impl FFT {
         freq.slice_mut(s![1..self.half_nfft_rounded])
             .par_mapv_inplace(|x| 2. * x / self.nfftF);
     }
-}
+}

src/ps/freqweighting.rs

@@ -1,7 +1,11 @@
 use crate::config::*;
 use strum_macros::{Display, EnumMessage};
 /// Sound level frequency weighting type (A, C, Z)
-#[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+// #[cfg_attr(feature = "python-bindings", pyclass(eq, eq_int))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(Display, Debug, EnumMessage, Default, Clone, PartialEq)]
 pub enum FreqWeighting {
     /// A-weighting

src/ps/ps.rs

@@ -7,17 +7,17 @@ use std::usize;
 
 use crate::Dcol;
 
-use super::window::*;
 use super::fft::FFT;
+use super::window::*;
 use std::mem::MaybeUninit;
 
 use realfft::{RealFftPlanner, RealToComplex};
 
 /// Cross power spectra, which is a 3D array, with the following properties:
-/// 
+///
 /// - The first index is the frequency index, starting at DC, ending at nfft/2.
 /// - The second, and third index result in `[i,j]` = C_ij = p_i  * conj(p_j)
-/// 
+///
 pub type CPSResult = Array3<Cflt>;
 
 /// Extra typical methods that are of use for 3D-arrays of complex numbers, that
@@ -46,7 +46,6 @@ pub trait CrossPowerSpecra {
     fn tf(&self, chi: usize, chj: usize, chRef: Option<usize>) -> Array1<Cflt>;
 }
 
-
 impl CrossPowerSpecra for CPSResult {
     fn ap(&self, ch: usize) -> Array1<Flt> {
         // Slice out one value for all frequencies, map to only real part, and
@@ -87,10 +86,10 @@ impl CrossPowerSpecra for CPSResult {
 /// estimation.
 ///
 pub struct PowerSpectra {
-    /// Window used in estimator
+    /// Window used in estimator. The actual Window in here is normalized with
-    pub window: Window,
+    /// the square root of the Window power. This safes one division when
-    /// The window power, is corrected for in power spectra estimants
+    /// processing time data.
-    pub sqrt_win_pwr: Flt,
+    pub window_normalized: Window,
 
     ffts: Vec<FFT>,
 
@@ -103,7 +102,7 @@ pub struct PowerSpectra {
 impl PowerSpectra {
     /// Returns the FFT length used in power spectra computations
     pub fn nfft(&self) -> usize {
-        self.window.win.len()
+        self.window_normalized.win.len()
     }
     /// Create new power spectra estimator. Uses FFT size from window length
     ///
@@ -116,9 +115,12 @@ impl PowerSpectra {
     ///
     /// - `window` - A `Window` struct, from which NFFT is also used.
     ///
-    pub fn newFromWindow(window: Window) -> PowerSpectra {
+    pub fn newFromWindow(mut window: Window) -> PowerSpectra {
         let nfft = window.win.len();
         let win_pwr = window.win.mapv(|w| w.powi(2)).sum() / (nfft as Flt);
+        let sqrt_win_pwr = Flt::sqrt(win_pwr);
+        window.win.mapv_inplace(|v| v / sqrt_win_pwr);
+
         assert!(nfft > 0);
         assert!(nfft % 2 == 0);
 
@@ -128,8 +130,7 @@ impl PowerSpectra {
         let Fft = FFT::new(fft);
 
         PowerSpectra {
-            window,
+            window_normalized: window,
-            sqrt_win_pwr: Flt::sqrt(win_pwr),
             ffts: vec![Fft],
             timedata: Array2::zeros((nfft, 1)),
             freqdata: Array2::zeros((nfft / 2 + 1, 1)),
@@ -138,7 +139,7 @@ impl PowerSpectra {
 
     /// Compute FFTs of input channel data. Stores the scaled FFT data in
     /// self.freqdata.
-    fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> &Array2<Cflt> {
+    fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> ArrayView2<Cflt> {
         let (n, nch) = timedata.dim();
         let nfft = self.nfft();
         assert!(n == nfft);
@@ -153,25 +154,27 @@ impl PowerSpectra {
         }
 
         assert!(n == self.nfft());
-        assert!(n == self.window.win.len());
+        assert!(n == self.window_normalized.win.len());
-        let sqrt_win_pwr = self.sqrt_win_pwr;
 
         // Multiply signals with window function, and compute fft's for each channel
         Zip::from(timedata.axis_iter(Axis(1)))
             .and(self.timedata.axis_iter_mut(Axis(1)))
             .and(&mut self.ffts)
             .and(self.freqdata.axis_iter_mut(Axis(1)))
-            .par_for_each(|time_in,mut time, fft, mut freq| {
+            .par_for_each(|time_in, mut time_tmp_storage, fft, mut freq| {
+                let DC = time_in.mean().unwrap();
 
+                azip!((t in &mut time_tmp_storage, &tin in time_in, &win in &self.window_normalized.win) {
+                // Substract DC value from time data, as this leaks into
+                // positive frequencies due to windowing.
                 // Multiply with window and copy over to local time data buffer
-                azip!((t in &mut time, &tin in time_in, &win in &self.window.win) *t=tin*win/sqrt_win_pwr);
+                    *t=(tin-DC)*win});
 
-                let tslice = time.as_slice().unwrap();
+                fft.process(&time_tmp_storage, &mut freq);
-                let fslice = freq.as_slice_mut().unwrap();
+                freq[0] = DC + 0. * I;
-                fft.process(tslice, fslice);
             });
 
-        &self.freqdata
+        self.freqdata.view()
     }
 
     /// Compute cross power spectra from input time data. First axis is

src/slm/mod.rs

@@ -24,7 +24,7 @@
 //! let settings = SLMSettingsBuilder::default()
 //!     .fs(48e3)
 //!     .freqWeighting(FreqWeighting::A)
-//!     .timeWeighting(TimeWeighting::Fast)
+//!     .timeWeighting(TimeWeighting::Fast{})
 //!     .filterDescriptors(&[desc]).build().unwrap();
 //! 
 //! let mut slm = SLM::new(settings);
@@ -33,7 +33,7 @@
 //! data[0] = 1.;
 //! 
 //! // Now apply some data. This is a kind of the SLM-s impulse response
-//! let res = slm.run(&data.as_slice().unwrap()).unwrap();
+//! let res = slm.run(data.as_slice().unwrap(), true).unwrap();
 //! 
 //! // Only one channel of result data
 //! assert_eq!(res.len(), 1);

src/slm/settings.rs

@@ -28,6 +28,7 @@ pub struct SLMSettings {
     pub filterDescriptors: Vec<StandardFilterDescriptor>,
 }
 
+#[cfg(feature="python-bindings")]
 #[cfg_attr(feature = "python-bindings", pymethods)]
 impl SLMSettings {
     #[new]

src/slm/slm.rs

@@ -205,6 +205,7 @@ impl SLM {
     }
 }
 
+#[cfg(feature="python-bindings")]
 #[cfg_attr(feature = "python-bindings", pymethods)]
 impl SLM {
     #[new]

src/slm/tw.rs

@@ -1,6 +1,10 @@
 use crate::config::*;
 /// Time weighting to use in level detection of Sound Level Meter.
-#[cfg_attr(feature = "python-bindings", pyclass(eq))]
+/// 
+// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
+// #[cfg_attr(feature = "python-bindings", pyclass(eq))]
+// For now:
+#[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(Clone, Copy, PartialEq)]
 pub enum TimeWeighting {
     // I know that the curly braces here are not required and add some