From 476479b693df0e87a7607faeddfc2b4e5f846063 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Wed, 28 Aug 2024 14:59:10 +0200 Subject: [PATCH] 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 | 38 ++++++++++++++++++++++++++-------- README.md | 5 +++-- src/daq/api/mod.rs | 5 ++++- src/daq/datatype.rs | 5 ++++- src/daq/mod.rs | 5 ++++- src/daq/qty.rs | 4 +++- src/daq/streamerror.rs | 7 ++++++- src/filter/octave.rs | 3 ++- src/ps/fft.rs | 10 +++++---- src/ps/freqweighting.rs | 6 +++++- src/ps/ps.rs | 45 ++++++++++++++++++++++------------------- src/slm/mod.rs | 4 ++-- src/slm/settings.rs | 1 + src/slm/slm.rs | 1 + src/slm/tw.rs | 6 +++++- 15 files changed, 100 insertions(+), 45 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index d7f06ec..1e5206d 100644 --- a/Cargo.toml +++ b/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"]} -numpy = { git = "https://github.com/JRRudy1/rust-numpy", branch = "pyo3-0.22.0" , optional = true} +pyo3 = { version = "0.21.2", optional = true, features = [ + "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 -# default = ["f64", "python-bindings", "record", "cpal-api"] +# When building Python bindings of Lasprs, you should enable the feature flag +# "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"] diff --git a/README.md b/README.md index f022c61..67327f0 100644 --- a/README.md +++ b/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" ``` diff --git a/src/daq/api/mod.rs b/src/daq/api/mod.rs index 29f352b..487f7da 100644 --- a/src/daq/api/mod.rs +++ b/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 { diff --git a/src/daq/datatype.rs b/src/daq/datatype.rs index ed53a63..52679a9 100644 --- a/src/daq/datatype.rs +++ b/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")] diff --git a/src/daq/mod.rs b/src/daq/mod.rs index 408463a..9e9a3c7 100644 --- a/src/daq/mod.rs +++ b/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 diff --git a/src/daq/qty.rs b/src/daq/qty.rs index c3049b7..e5849d2 100644 --- a/src/daq/qty.rs +++ b/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 diff --git a/src/daq/streamerror.rs b/src/daq/streamerror.rs index b3f4539..5bc4a3c 100644 --- a/src/daq/streamerror.rs +++ b/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( diff --git a/src/filter/octave.rs b/src/filter/octave.rs index 603661e..d388784 100644 --- a/src/filter/octave.rs +++ b/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 { diff --git a/src/ps/fft.rs b/src/ps/fft.rs index e4585b7..a7b708b 100644 --- a/src/ps/fft.rs +++ b/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>, + T: Into>, + U: Into>, { 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); } -} \ No newline at end of file +} diff --git a/src/ps/freqweighting.rs b/src/ps/freqweighting.rs index 1369102..d13d7bd 100644 --- a/src/ps/freqweighting.rs +++ b/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 diff --git a/src/ps/ps.rs b/src/ps/ps.rs index 46bbde5..959a18d 100644 --- a/src/ps/ps.rs +++ b/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; /// 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) -> Array1; } - impl CrossPowerSpecra for CPSResult { fn ap(&self, ch: usize) -> Array1 { // 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 - pub window: Window, - /// The window power, is corrected for in power spectra estimants - pub sqrt_win_pwr: Flt, + /// Window used in estimator. The actual Window in here is normalized with + /// the square root of the Window power. This safes one division when + /// processing time data. + pub window_normalized: Window, ffts: Vec, @@ -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, - sqrt_win_pwr: Flt::sqrt(win_pwr), + window_normalized: window, 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) -> &Array2 { + fn compute_ffts(&mut self, timedata: ArrayView2) -> ArrayView2 { 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()); - let sqrt_win_pwr = self.sqrt_win_pwr; + assert!(n == self.window_normalized.win.len()); // 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(); - let fslice = freq.as_slice_mut().unwrap(); - fft.process(tslice, fslice); + fft.process(&time_tmp_storage, &mut freq); + freq[0] = DC + 0. * I; }); - &self.freqdata + self.freqdata.view() } /// Compute cross power spectra from input time data. First axis is diff --git a/src/slm/mod.rs b/src/slm/mod.rs index 58f39a0..597e0e3 100644 --- a/src/slm/mod.rs +++ b/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); diff --git a/src/slm/settings.rs b/src/slm/settings.rs index 78b7d08..674f090 100644 --- a/src/slm/settings.rs +++ b/src/slm/settings.rs @@ -28,6 +28,7 @@ pub struct SLMSettings { pub filterDescriptors: Vec, } +#[cfg(feature="python-bindings")] #[cfg_attr(feature = "python-bindings", pymethods)] impl SLMSettings { #[new] diff --git a/src/slm/slm.rs b/src/slm/slm.rs index 8246a44..14c3b8f 100644 --- a/src/slm/slm.rs +++ b/src/slm/slm.rs @@ -205,6 +205,7 @@ impl SLM { } } +#[cfg(feature="python-bindings")] #[cfg_attr(feature = "python-bindings", pymethods)] impl SLM { #[new] diff --git a/src/slm/tw.rs b/src/slm/tw.rs index 9b3f3d1..7d1b594 100644 --- a/src/slm/tw.rs +++ b/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