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