Compare commits

...

9 Commits

26 changed files with 705 additions and 274 deletions

View File

View File

@ -1,6 +1,6 @@
[package] [package]
name = "lasprs" name = "lasprs"
version = "0.6.1" version = "0.6.2"
edition = "2021" edition = "2021"
authors = ["J.A. de Jong <j.a.dejong@ascee.nl>"] authors = ["J.A. de Jong <j.a.dejong@ascee.nl>"]
description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)" description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)"

View File

@ -12,12 +12,17 @@ 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). Documentation is provided at [doc.rs](https://docs.rs/lasprs/latest/lasprs).
## Python bindings and examples ## Python bindings and examples
The library has Python bindings (via [pyo3](https://pyo3.rs), which can be installed via: 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 "python-bindings" $ pip install lasprs
``` ```
which pulls the library from [Pypi](https://pypi.org).
Examples of how to use the library are provided in Jupyter Notebooks, which can be found in the repository, see [lasprs/examples_py](https://code.ascee.nl/ASCEE/lasprs/src/branch/main/examples_py).
More examples will follow in the near future.

146
examples_py/test_aps.ipynb Normal file
View File

@ -0,0 +1,146 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "6763c6d0-878f-4492-840e-40d11a52f91e",
"metadata": {},
"source": [
"# Test Welch method implementation Python wrappers"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "6c06a0c3-dcf7-4d21-b1f3-eced84e31218",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"🍹 Building a mixed python/rust project\n",
"🔗 Found pyo3 bindings\n",
"🐍 Found CPython 3.10 at /home/anne/.pyenv/versions/lasprs/bin/python\n",
"Requirement already satisfied: numpy in /home/anne/.pyenv/versions/lasprs/lib/python3.10/site-packages (1.26.4)\n",
"\u001b[1m\u001b[32m Compiling\u001b[0m pyo3-build-config v0.21.2\n",
"\u001b[K\u001b[1m\u001b[32m Compiling\u001b[0m pyo3-macros-backend v0.21.2 178/189: pyo3-build-config \n",
"\u001b[K\u001b[1m\u001b[32m Compiling\u001b[0m pyo3-ffi v0.21.2=======> ] 178/189: pyo3-macros-backend, pyo3...\n",
"\u001b[1m\u001b[32m Compiling\u001b[0m pyo3 v0.21.2\n",
"\u001b[K\u001b[1m\u001b[32m Compiling\u001b[0m pyo3-macros v0.21.2=====> ] 183/189: pyo3-macros-backend, pyo3...\n",
"\u001b[K\u001b[1m\u001b[32m Compiling\u001b[0m numpy v0.21.0===========> ] 186/189: pyo3 \n",
"\u001b[K\u001b[1m\u001b[32m Compiling\u001b[0m lasprs v0.6.1 (/home/anne/wip/mycode/lasprs) \n",
"\u001b[K\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[0m\u001b[1m: type alias `Vc` is never used\u001b[0m \n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m--> \u001b[0m\u001b[0msrc/config.rs:63:10\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m63\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0mpub type Vc = Vec<Cflt>;\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m= \u001b[0m\u001b[0m\u001b[1mnote\u001b[0m\u001b[0m: `#[warn(dead_code)]` on by default\u001b[0m\n",
"\n",
"\u001b[K\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[0m\u001b[1m: type alias `Cmat` is never used\u001b[0m \n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m--> \u001b[0m\u001b[0msrc/config.rs:75:10\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m75\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0mpub type Cmat = Array2<Cflt>;\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^^^\u001b[0m\n",
"\n",
"\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[0m\u001b[1m: methods `ninchannels` and `noutchannels` are never used\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m--> \u001b[0m\u001b[0msrc/daq/api/mod.rs:25:8\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m20\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0mpub trait Stream {\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m------\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12mmethods in this trait\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m...\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m25\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m fn ninchannels(&self) -> usize;\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^^^^^^^^^^\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m...\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m28\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m fn noutchannels(&self) -> usize;\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^^^^^^^^^^^\u001b[0m\n",
"\n",
"\u001b[K\u001b[0m\u001b[1m\u001b[33mwarning\u001b[0m\u001b[0m\u001b[1m: method `setPreFilter` is never used\u001b[0m \n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m--> \u001b[0m\u001b[0msrc/siggen.rs:320:12\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m318\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0mimpl SiggenChannelConfig {\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m------------------------\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12mmethod in this implementation\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m319\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m /// Set new pre-filter that filters the source signal\u001b[0m\n",
"\u001b[0m\u001b[1m\u001b[38;5;12m320\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m pub fn setPreFilter(&mut self, pref: Option<Box<dyn Filter>>) {\u001b[0m\n",
"\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[38;5;12m|\u001b[0m\u001b[0m \u001b[0m\u001b[0m\u001b[1m\u001b[33m^^^^^^^^^^^^\u001b[0m\n",
"\n",
"\u001b[K\u001b[1m\u001b[33mwarning\u001b[0m\u001b[1m:\u001b[0m `lasprs` (lib) generated 4 warningssprs \n",
"\u001b[1m\u001b[32m Finished\u001b[0m \u001b]8;;https://doc.rust-lang.org/cargo/reference/profiles.html#default-profiles\u001b\\`dev` profile [unoptimized + debuginfo]\u001b]8;;\u001b\\ target(s) in 6.04s\n",
"📦 Built wheel for CPython 3.10 to /tmp/.tmp03vGw7/lasprs-0.6.1-cp310-cp310-linux_x86_64.whl\n",
"✏️ Setting installed package as editable\n",
"🛠 Installed lasprs-0.6.1\n"
]
}
],
"source": [
"!cd .. && maturin develop -F python-bindings"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "3e651371-e7a5-4013-a265-ab80f8b1eb41",
"metadata": {},
"outputs": [],
"source": [
"from lasprs import WindowType"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "9d6bdf3b-40f8-4c3b-85ac-4ff8b8c7dcae",
"metadata": {},
"outputs": [],
"source": [
"w = WindowType.all()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3a31899d-37de-4653-9b6e-61dcf7c683f4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[WindowType.Hann,\n",
" WindowType.Hamming,\n",
" WindowType.Rect,\n",
" WindowType.Bartlett,\n",
" WindowType.Blackman]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}

View File

@ -41,8 +41,8 @@
"source": [ "source": [
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"from numpy import log10, sqrt, exp, pi\n", "from numpy import log10\n",
"from lasprs._lasprs import StandardFilterDescriptor, SLMSettings, FreqWeighting, TimeWeighting, SLM\n", "from lasprs import StandardFilterDescriptor\n",
"def level(a):\n", "def level(a):\n",
" return 20*np.log10(np.abs(a))" " return 20*np.log10(np.abs(a))"
] ]

View File

@ -41,8 +41,8 @@
"source": [ "source": [
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"from numpy import log10, sqrt, exp, pi\n", "from numpy import log10\n",
"from lasprs._lasprs import ZPKModel, FreqWeighting\n", "from lasprs import ZPKModel, FreqWeighting\n",
"def level(a):\n", "def level(a):\n",
" return 20*np.log10(np.abs(a))" " return 20*np.log10(np.abs(a))"
] ]

View File

@ -17,6 +17,7 @@
"source": [ "source": [
"# Prerequisites, uncomment below in case of errors. Also for ipympl, restart Jupyter Lab if it was not installed\n", "# Prerequisites, uncomment below in case of errors. Also for ipympl, restart Jupyter Lab if it was not installed\n",
"\n", "\n",
"# Only do the following if you are in develop mode\n",
"!cd .. && maturin develop -F python-bindings\n", "!cd .. && maturin develop -F python-bindings\n",
"#!pip install ipympl scipy matplotlib" "#!pip install ipympl scipy matplotlib"
] ]
@ -41,7 +42,7 @@
"source": [ "source": [
"import numpy as np\n", "import numpy as np\n",
"import matplotlib.pyplot as plt\n", "import matplotlib.pyplot as plt\n",
"from numpy import log10, sqrt, exp, pi\n", "from numpy import log10\n",
"from lasprs._lasprs import StandardFilterDescriptor, SLMSettings, FreqWeighting, TimeWeighting, SLM\n", "from lasprs._lasprs import StandardFilterDescriptor, SLMSettings, FreqWeighting, TimeWeighting, SLM\n",
"def level(a):\n", "def level(a):\n",
" return 20*np.log10(np.abs(a))" " return 20*np.log10(np.abs(a))"
@ -105,27 +106,7 @@
"outputs": [], "outputs": [],
"source": [ "source": [
"t = np.linspace(0, N/fs, N, endpoint=False)\n", "t = np.linspace(0, N/fs, N, endpoint=False)\n",
"out = slm.run(inp)" "out = slm.run(inp, True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "37af9b78-73b9-4228-a758-d9a5362f79f4",
"metadata": {},
"outputs": [],
"source": [
"len(out[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4609a63f-e350-4b6b-b898-93c9e6fef804",
"metadata": {},
"outputs": [],
"source": [
"# help(plt.legend)"
] ]
}, },
{ {
@ -141,14 +122,6 @@
"plt.ylim(-60, 0)\n", "plt.ylim(-60, 0)\n",
"plt.legend(names)" "plt.legend(names)"
] ]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4303e20e-cc61-4d23-a5f0-cab683d55912",
"metadata": {},
"outputs": [],
"source": []
} }
], ],
"metadata": { "metadata": {

View File

@ -1 +0,0 @@
from .._lasprs import (Biquad, BiquadBank, SeriesBiquad, ZPKModel, FilterSpec)

View File

@ -11,7 +11,7 @@ cfg_if::cfg_if! {
/// Ratio between circumference and diameter of a circle /// Ratio between circumference and diameter of a circle
pub const pi: Flt = std::f64::consts::PI; pub const pi: Flt = std::f64::consts::PI;
} }
else if #[cfg(feature="f32")] { else if #[cfg(feature="f32")] {
/// Floating-point value, compile time option to make it either f32, or f64 /// Floating-point value, compile time option to make it either f32, or f64
@ -26,7 +26,8 @@ cfg_if::cfg_if! {
cfg_if::cfg_if! { cfg_if::cfg_if! {
if #[cfg(feature = "python-bindings")] { if #[cfg(feature = "python-bindings")] {
pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArrayDyn, PyArrayLike1, PyReadonlyArrayDyn}; pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArray2, PyArray3, PyArrayDyn, PyArrayLike1,
PyArrayLike2,PyArrayLike3,PyReadonlyArrayDyn, convert::ToPyArray};
pub use pyo3::prelude::*; pub use pyo3::prelude::*;
pub use pyo3::exceptions::PyValueError; pub use pyo3::exceptions::PyValueError;
pub use pyo3::{pymodule, types::PyModule, PyResult}; pub use pyo3::{pymodule, types::PyModule, PyResult};

View File

@ -1,5 +1,5 @@
use crate::{config::*, ZPKModel}; use crate::{config::*, ZPKModel};
use anyhow::{anyhow, bail, Result}; use anyhow::{anyhow, bail, Error, Result};
use num::{traits::float, Float}; use num::{traits::float, Float};
use rayon::iter::Filter; use rayon::iter::Filter;
use softfloat::F64; use softfloat::F64;
@ -497,6 +497,38 @@ impl StandardFilterDescriptor {
Ok(Self::genThirdOctaveFilterSet(fmin, fmax)?) Ok(Self::genThirdOctaveFilterSet(fmin, fmax)?)
} }
#[staticmethod]
#[pyo3(name = "genFilterSetInRangeWithStartStopDescriptors")]
fn genFilterSet_range(
start: StandardFilterDescriptor,
stop: StandardFilterDescriptor,
include_overall: bool,
) -> PyResult<Vec<StandardFilterDescriptor>> {
if start.b == 0 {
return Err(PyValueError::new_err(
"Overall band designator cannot be used to define range",
));
}
if start.b != stop.b {
return Err(PyValueError::new_err(
"Start and stop filter are not of same type",
));
}
if start.x > stop.x {
return Err(PyValueError::new_err(
"Start filter has higher midband frequency than stop frequency",
));
}
let mut res = (start.x..=stop.x)
.map(|x| StandardFilterDescriptor { b: start.b, x: x })
.collect::<Vec<StandardFilterDescriptor>>();
if include_overall {
res.push(StandardFilterDescriptor::Overall().expect("Overall should not err"))
}
Ok(res)
}
#[staticmethod] #[staticmethod]
#[pyo3(name = "genOctaveFilterSet")] #[pyo3(name = "genOctaveFilterSet")]
fn genOctaveFilterSetFromFreq(fmin: Flt, fmax: Flt) -> PyResult<Vec<StandardFilterDescriptor>> { fn genOctaveFilterSetFromFreq(fmin: Flt, fmax: Flt) -> PyResult<Vec<StandardFilterDescriptor>> {
@ -527,8 +559,18 @@ impl StandardFilterDescriptor {
fmax: Flt, fmax: Flt,
append_overall: bool, append_overall: bool,
) -> PyResult<Vec<StandardFilterDescriptor>> { ) -> PyResult<Vec<StandardFilterDescriptor>> {
if b == 0 {
return Err(PyValueError::new_err(
"Overall band designator cannot be used to define range",
));
}
Ok(Self::genFilterSetForRange(b, fmin, fmax, append_overall)?) Ok(Self::genFilterSetForRange(b, fmin, fmax, append_overall)?)
} }
#[pyo3(name="fm")]
fn fm_py(&self) -> Option<Flt> {
self.fm()
}
} }
#[cfg(test)] #[cfg(test)]

View File

@ -569,7 +569,6 @@ impl ZPKModel {
let p2 = 2. * pi * f2; let p2 = 2. * pi * f2;
let p3 = 2. * pi * f3; let p3 = 2. * pi * f3;
let p4 = 2. * pi * f4; let p4 = 2. * pi * f4;
println!("{b} {p1}, {p2}, {p3}, {p4}");
let (zeros, poles) = match wt { let (zeros, poles) = match wt {
FreqWeighting::Z => { FreqWeighting::Z => {
@ -610,11 +609,21 @@ pub enum PoleOrZero {
#[cfg(test)] #[cfg(test)]
mod test{ mod test{
use approx::assert_abs_diff_eq;
use num::complex::ComplexFloat;
use crate::TransferFunction;
use super::ZPKModel; use super::ZPKModel;
#[test] #[test]
fn test_A() { fn test_A() {
let Aw = ZPKModel::freqWeightingFilter(crate::FreqWeighting::A); let Aw = ZPKModel::freqWeightingFilter(crate::FreqWeighting::A);
assert_abs_diff_eq!(Aw.tf(0., &[1000.])[0].abs(), 1.0);
}
#[test]
fn test_C() {
let Cw = ZPKModel::freqWeightingFilter(crate::FreqWeighting::C);
assert_abs_diff_eq!(Cw.tf(0., &[1000.])[0].abs(), 1.0);
} }
} }

View File

@ -62,6 +62,12 @@ fn lasprs(m: &Bound<'_, PyModule>) -> PyResult<()> {
m.add_class::<ps::FreqWeighting>()?; m.add_class::<ps::FreqWeighting>()?;
m.add_class::<slm::SLMSettings>()?; m.add_class::<slm::SLMSettings>()?;
m.add_class::<slm::SLM>()?; m.add_class::<slm::SLM>()?;
m.add_class::<slm::TimeWeighting>()?;
m.add_class::<ps::WindowType>()?;
m.add_class::<ps::Overlap>()?;
m.add_class::<ps::ApsMode>()?;
m.add_class::<ps::ApsSettings>()?;
m.add_class::<ps::AvPowerSpectra>()?;
Ok(()) Ok(())
} }

View File

@ -1,214 +1,25 @@
use super::timebuffer::TimeBuffer;
use super::CrossPowerSpecra;
use super::*; use super::*;
use crate::config::*; use super::{timebuffer::TimeBuffer, CrossPowerSpecra};
use crate::{config::*, TransferFunction, ZPKModel};
use anyhow::{bail, Error, Result}; use anyhow::{bail, Error, Result};
use derive_builder::Builder;
use freqweighting::FreqWeighting; use freqweighting::FreqWeighting;
/// All settings used for computing averaged power spectra using Welch' method.
#[derive(Builder, Clone)]
#[builder(build_fn(validate = "Self::validate", error = "Error"))]
pub struct ApsSettings {
/// Mode of computation, see [ApsMode].
#[builder(default)]
mode: ApsMode,
/// Overlap in time segments. See [Overlap].
#[builder(default)]
overlap: Overlap,
/// Window applied to time segments. See [WindowType].
#[builder(default)]
windowType: WindowType,
/// Kind of freqency weighting. Defaults to Z
#[builder(default)]
freqWeightingType: FreqWeighting,
/// FFT Length
nfft: usize,
/// Sampling frequency
fs: Flt,
}
impl ApsSettingsBuilder {
fn validate(&self) -> Result<()> {
if !self.fs.is_some() {
bail!("Sampling frequency not given");
}
let fs = self.fs.unwrap();
if !fs.is_normal() {
bail!("Sampling frequency not a normal number")
}
if fs <= 0.0 {
bail!("Invalid sampling frequency given as parameter");
}
if self.nfft.is_none() {
bail!("nfft not specified")
};
let nfft = self.nfft.unwrap();
if nfft % 2 != 0 {
bail!("NFFT should be even")
}
if nfft == 0 {
bail!("Invalid NFFT, should be > 0.")
}
// Perform some checks on ApsMode
if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
if tau <= 0.0 {
bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
}
}
Ok(())
}
}
impl ApsSettings {
/// Returns nfft
pub fn nfft(&self) -> usize {
self.nfft
}
fn get_overlap_keep(&self) -> usize {
self.validate_get_overlap_keep().unwrap()
}
/// Unpack all, returns parts in tuple
pub fn get(self) -> (ApsMode, Overlap, WindowType, FreqWeighting, usize, Flt) {
(
self.mode,
self.overlap,
self.windowType,
self.freqWeightingType,
self.nfft,
self.fs,
)
}
/// Returns the amount of samples to `keep` in the time buffer when
/// overlapping time segments using [TimeBuffer].
pub fn validate_get_overlap_keep(&self) -> Result<usize> {
let nfft = self.nfft;
let overlap_keep = match self.overlap {
Overlap::Number(i) if i >= nfft => {
bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
}
// Keep 1 sample, if overlap is 1 sample etc.
Overlap::Number(i) if i < nfft => i,
// If overlap percentage is >= 100, or < 0.0 its an error
Overlap::Percentage(p) if !(0.0..100.).contains(&p) => {
bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
}
// If overlap percentage is 0, this gives
Overlap::Percentage(p) => ((p * nfft as Flt) / 100.) as usize,
Overlap::NoOverlap => 0,
_ => unreachable!(),
};
if overlap_keep >= nfft {
bail!("Computed overlap results in invalid number of overlap samples. Please make sure the FFT length is large enough, when high overlap percentages are required.");
}
Ok(overlap_keep)
}
/// Return a reasonable acoustic default with a frequency resolution around
/// ~ 10 Hz, where nfft is still an integer power of 2.
///
/// # Errors
///
/// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
///
pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
if fs < 1e3 || fs > 1e6 {
bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
}
let fs_div_10_rounded = (fs / 10.) as u32;
// 2^30 is about 1 million. We search for a two-power of an nfft that is
// the closest to fs/10. The frequency resolution is about fs/nfft.
let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
// Start wth a value that is always too large
fs as u32 * 10,
|cur, new| cur.min(new),
) as usize;
Ok(ApsSettings {
mode,
fs,
nfft,
windowType: WindowType::default(),
overlap: Overlap::default(),
freqWeightingType: FreqWeighting::default(),
})
}
/// Return sampling frequency
pub fn fs(&self) -> Flt {
self.fs
}
/// Return Nyquist frequency
pub fn fnyq(&self) -> Flt {
self.fs() / 2.
}
/// Returns a single-sided frequency array corresponding to points in Power
/// spectra computation.
pub fn getFreq(&self) -> Array1<Flt> {
let df = self.fs / self.nfft as Flt;
let K = self.nfft / 2 + 1;
Array1::linspace(0., (K - 1) as Flt * df, K)
}
}
/// Provide the overlap of blocks for computing averaged (cross) power spectra.
/// Can be provided as a percentage of the block size, or as a number of
/// samples.
#[derive(Clone, Debug)]
pub enum Overlap {
/// Overlap specified as a percentage of the total FFT length. Value should
/// be 0<=pct<100.
Percentage(Flt),
/// Number of samples to overlap
Number(usize),
/// No overlap at all, which is the same as Overlap::Number(0)
NoOverlap,
}
impl Default for Overlap {
fn default() -> Self {
Overlap::Percentage(50.)
}
}
/// The 'mode' used in computing averaged power spectra. When providing data in
/// blocks to the [AvPowerSpectra] the resulting 'current estimate' responds
/// differently, depending on the model.
#[derive(Default, Copy, Clone)]
pub enum ApsMode {
/// Averaged over all data provided. New averages can be created by calling
/// `AvPowerSpectra::reset()`
#[default]
AllAveraging,
/// In this mode, the `AvPowerSpectra` works a bit like a sound level meter,
/// where new data is weighted with old data, and old data exponentially
/// backs off. This mode only makes sense when `tau >> nfft/fs`
ExponentialWeighting {
/// Time weighting constant, follows convention of Sound Level Meters.
/// Means the data is approximately low-pass filtered with a cut-off
/// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
tau: Flt,
},
/// Spectrogram mode. Only returns the latest estimate(s).
Spectrogram,
}
/// Averaged power spectra computing engine /// Averaged power spectra computing engine
/// Used to compute power spectra estimations on /// Used to compute power spectra estimations on
/// long datasets, where nfft << length of data. This way, the variance of a /// long datasets, where nfft << length of data. This way, the variance of a
/// single periodogram is suppressed with increasing number of averages. /// single periodogram is suppressed with increasing number of averages.
/// ///
/// For more information, see the book on numerical recipes. /// For more information, see the book on numerical recipes.
///
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pyclass)]
#[derive(Debug)]
pub struct AvPowerSpectra { pub struct AvPowerSpectra {
// Power spectra estimator for single block // Power spectra estimator for single block
ps: PowerSpectra, ps: PowerSpectra,
// Settings for computing power spectra, see [ApsSettings]
settings: ApsSettings, settings: ApsSettings,
// The number of samples to keep in the time buffer when overlapping time // The number of samples to keep in the time buffer when overlapping time
@ -221,6 +32,11 @@ pub struct AvPowerSpectra {
/// Storage for sample data. /// Storage for sample data.
timebuf: TimeBuffer, timebuf: TimeBuffer,
/// Power scaling for of applied frequency weighting. Multiply each power
/// and cross-power value with these constants to apply the frequency
/// weighting.
freqWeighting_pwr: Option<Dcol>,
// Current estimation of the power spectra // Current estimation of the power spectra
cur_est: CPSResult, cur_est: CPSResult,
} }
@ -254,8 +70,7 @@ impl AvPowerSpectra {
/// - When providing invalid sampling frequencies /// - When providing invalid sampling frequencies
/// ///
pub fn new_simple_all_averaging(fs: Flt, nfft: usize) -> AvPowerSpectra { pub fn new_simple_all_averaging(fs: Flt, nfft: usize) -> AvPowerSpectra {
let mut settings = let mut settings = ApsSettings::reasonableAcousticDefault(fs, ApsMode::default()).unwrap();
ApsSettings::reasonableAcousticDefault(fs, ApsMode::AllAveraging).unwrap();
settings.nfft = nfft; settings.nfft = nfft;
AvPowerSpectra::new(settings) AvPowerSpectra::new(settings)
} }
@ -285,19 +100,51 @@ impl AvPowerSpectra {
let ps = PowerSpectra::newFromWindow(window); let ps = PowerSpectra::newFromWindow(window);
let freq = settings.getFreq();
let freqWeighting_pwr = match settings.freqWeightingType {
FreqWeighting::Z => None,
_ => {
let fw_pwr = ZPKModel::freqWeightingFilter(settings.freqWeightingType)
.tf(0., &freq)
.mapv(|a| a.abs() * a.abs());
Some(fw_pwr)
}
};
AvPowerSpectra { AvPowerSpectra {
ps, ps,
overlap_keep, overlap_keep,
settings, settings,
N: 0, N: 0,
freqWeighting_pwr,
cur_est: CPSResult::default((0, 0, 0)), cur_est: CPSResult::default((0, 0, 0)),
timebuf: TimeBuffer::new(), timebuf: TimeBuffer::new(),
} }
} }
// Update result for single block // Update result for single block
fn update_singleblock(&mut self, timedata: ArrayView2<Flt>) { fn update_singleblock(&mut self, timedata: ArrayView2<Flt>) {
let Cpsnew = self.ps.compute(timedata); // Compute updated block of cross-spectral density
// println!("Cpsnew: {:?}", Cpsnew[[0, 0, 0]]); let Cpsnew = {
let mut Cpsnew = self.ps.compute(timedata);
let dim = Cpsnew.dim();
if let Some(fw_pwr) = &self.freqWeighting_pwr {
// Frequency weighting on power is available, multiply all values with frequency weighting
// Option 1: azip with indexing
// azip!((index (i,_,_), cps in &mut Cpsnew) {
// *cps = Cflt::new(cps.re * fw_pwr[[i]], cps.im * fw_pwr[[i]]);
// });
// Option 2: broadcasting with an unwrap.
let dview_fw = fw_pwr.slice(s![..,NewAxis, NewAxis]);
azip!((c in &mut Cpsnew, f in dview_fw.broadcast(dim).expect("BUG: Cannot broadcast")) {
// Scale with frequency weighting
*c = Cflt::new(c.re * f, c.im * f);
});
}
Cpsnew
};
// Initialize to zero // Initialize to zero
if self.cur_est.is_empty() { if self.cur_est.is_empty() {
@ -310,7 +157,7 @@ impl AvPowerSpectra {
// Apply operation based on mode // Apply operation based on mode
match self.settings.mode { match self.settings.mode {
ApsMode::AllAveraging => { ApsMode::AllAveraging {} => {
let Nf = Cflt { let Nf = Cflt {
re: self.N as Flt, re: self.N as Flt,
im: 0., im: 0.,
@ -359,7 +206,7 @@ impl AvPowerSpectra {
} }
} }
ApsMode::Spectrogram => { ApsMode::Spectrogram {} => {
self.cur_est = Cpsnew; self.cur_est = Cpsnew;
} }
} }
@ -395,9 +242,10 @@ impl AvPowerSpectra {
computed_single = true; computed_single = true;
} }
if computed_single { if computed_single {
return Some(&self.cur_est); Some(&self.cur_est)
} else {
None
} }
None
} }
/// Computes average (cross)power spectra, and returns all intermediate /// Computes average (cross)power spectra, and returns all intermediate
@ -435,6 +283,29 @@ impl AvPowerSpectra {
} }
} }
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl AvPowerSpectra {
#[new]
fn new_py(s: ApsSettings) -> AvPowerSpectra {
AvPowerSpectra::new(s)
}
#[pyo3(name = "compute")]
fn compute_py<'py>(
&mut self,
py: Python<'py>,
dat: PyArrayLike2<Flt>,
) -> Bound<'py, PyArray3<Cflt>> {
let dat = dat.as_array();
if let Some(res) = self.compute_last(dat) {
let res = res.clone();
return res.to_pyarray_bound(py);
}
let res: Bound<'py, PyArray3<Cflt>> = PyArray3::zeros_bound(py, [0, 0, 0], true);
res
}
}
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use approx::assert_abs_diff_eq; use approx::assert_abs_diff_eq;
@ -450,11 +321,11 @@ mod test {
#[test] #[test]
fn test_overlap_keep() { fn test_overlap_keep() {
let ol = [ let ol = [
Overlap::NoOverlap, Overlap::NoOverlap {},
Percentage(50.), Percentage { pct: 50. },
Percentage(50.), Percentage { pct: 50. },
Percentage(25.), Percentage { pct: 25. },
Overlap::Number(10), Overlap::Number { N: 10 },
]; ];
let nffts = [10, 10, 1024, 10]; let nffts = [10, 10, 1024, 10];
let expected_keep = [0, 5, 512, 2, 10]; let expected_keep = [0, 5, 512, 2, 10];
@ -465,7 +336,7 @@ mod test {
.fs(1.) .fs(1.)
.overlap(overlap.clone()) .overlap(overlap.clone())
.build() .build()
.unwrap(); .expect("BUG: Settings cannot be build.");
assert_eq!(settings.get_overlap_keep(), *expected); assert_eq!(settings.get_overlap_keep(), *expected);
} }
@ -481,15 +352,15 @@ mod test {
let settings = ApsSettingsBuilder::default() let settings = ApsSettingsBuilder::default()
.fs(fs) .fs(fs)
.nfft(nfft) .nfft(nfft)
.overlap(Overlap::NoOverlap) .overlap(Overlap::NoOverlap {})
.mode(ApsMode::ExponentialWeighting { tau }) .mode(ApsMode::ExponentialWeighting { tau })
.build() .build()
.unwrap(); .expect("Settings cannot be empty");
let overlap_keep = settings.get_overlap_keep(); let overlap_keep = settings.get_overlap_keep();
let mut aps = AvPowerSpectra::new(settings); let mut aps = AvPowerSpectra::new(settings);
assert_eq!(aps.overlap_keep, 0); assert_eq!(aps.overlap_keep, 0);
let distr = Normal::new(1.0, 1.0).unwrap(); let distr = Normal::new(1.0, 1.0).expect("Distribution cannot be built");
let timedata_some = Dmat::random((nfft, 1), distr); let timedata_some = Dmat::random((nfft, 1), distr);
let timedata_zeros = Dmat::zeros((nfft, 1)); let timedata_zeros = Dmat::zeros((nfft, 1));

35
src/ps/apsmode.rs Normal file
View File

@ -0,0 +1,35 @@
use crate::config::*;
/// The 'mode' used in computing averaged power spectra. When providing data in
/// blocks to the [AvPowerSpectra] the resulting 'current estimate' responds
/// differently, depending on the model.
#[derive(Copy, Clone, PartialEq, Debug)]
#[cfg_attr(feature = "python-bindings", pyclass)]
pub enum ApsMode {
/// Averaged over all data provided. New averages can be created by calling
/// `AvPowerSpectra::reset()`
AllAveraging {},
/// In this mode, the `AvPowerSpectra` works a bit like a sound level meter,
/// where new data is weighted with old data, and old data exponentially
/// backs off. This mode only makes sense when `tau >> nfft/fs`
ExponentialWeighting {
/// Time weighting constant, follows convention of Sound Level Meters.
/// Means the data is approximately low-pass filtered with a cut-off
/// frequency f_c of s/tau ≅ 1 → f_c = (2 * pi * tau)^-1.
tau: Flt,
},
/// Spectrogram mode. Only returns the latest estimate(s).
Spectrogram {},
}
impl Default for ApsMode {
fn default() -> Self {
ApsMode::AllAveraging {}
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl ApsMode {
#[inline]
fn __eq__(&self, other: &Self) -> bool {
self == other
}
}

168
src/ps/apssettings.rs Normal file
View File

@ -0,0 +1,168 @@
use super::*;
use crate::config::*;
use anyhow::{bail, Error, Result};
use derive_builder::Builder;
/// All settings used for computing averaged power spectra using Welch' method.
#[derive(Builder, Clone, Debug)]
#[cfg_attr(feature = "python-bindings", pyclass)]
#[builder(build_fn(validate = "Self::validate", error = "Error"))]
pub struct ApsSettings {
/// Mode of computation, see [ApsMode].
#[builder(default)]
pub mode: ApsMode,
/// Overlap in time segments. See [Overlap].
#[builder(default)]
pub overlap: Overlap,
/// Window applied to time segments. See [WindowType].
#[builder(default)]
pub windowType: WindowType,
/// Kind of freqency weighting. Defaults to Z
#[builder(default)]
pub freqWeightingType: FreqWeighting,
/// FFT Length
pub nfft: usize,
/// Sampling frequency
pub fs: Flt,
}
impl ApsSettingsBuilder {
fn validate(&self) -> Result<()> {
if !self.fs.is_some() {
bail!("Sampling frequency not given");
}
let fs = self.fs.unwrap();
if !fs.is_normal() {
bail!("Sampling frequency not a normal number")
}
if fs <= 0.0 {
bail!("Invalid sampling frequency given as parameter");
}
if self.nfft.is_none() {
bail!("nfft not specified")
};
let nfft = self.nfft.unwrap();
if nfft % 2 != 0 {
bail!("NFFT should be even")
}
if nfft == 0 {
bail!("Invalid NFFT, should be > 0.")
}
// Perform some checks on ApsMode
if let Some(ApsMode::ExponentialWeighting { tau }) = self.mode {
if tau <= 0.0 {
bail!("Invalid time weighting constant [s]. Should be > 0 if given.");
}
}
Ok(())
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl ApsSettings {
#[new]
fn new(
mode: ApsMode,
overlap: Overlap,
windowType: WindowType,
freqWeightingType: FreqWeighting,
nfft: usize,
fs: Flt,
) -> ApsSettings {
ApsSettings {
mode,
overlap,
windowType,
freqWeightingType,
nfft,
fs,
}
}
}
impl ApsSettings {
/// Returns the amount of samples to keep in overlapping blocks of power
/// spectra.
pub fn get_overlap_keep(&self) -> usize {
self.validate_get_overlap_keep().unwrap()
}
/// Returns the amount of samples to `keep` in the time buffer when
/// overlapping time segments using [TimeBuffer].
fn validate_get_overlap_keep(&self) -> Result<usize> {
let nfft = self.nfft;
let overlap_keep = match self.overlap {
Overlap::Number { N } if N >= nfft => {
bail!("Invalid overlap number of samples. Should be < nfft, which is {nfft}.")
}
// Keep 1 sample, if overlap is 1 sample etc.
Overlap::Number { N } if N < nfft => N,
// If overlap percentage is >= 100, or < 0.0 its an error
Overlap::Percentage { pct } if !(0.0..100.).contains(&pct) => {
bail!("Invalid overlap percentage. Should be >= 0. And < 100.")
}
// If overlap percentage is 0, this gives
Overlap::Percentage { pct } => ((pct * nfft as Flt) / 100.) as usize,
Overlap::NoOverlap {} => 0,
_ => unreachable!(),
};
if overlap_keep >= nfft {
bail!("Computed overlap results in invalid number of overlap samples. Please make sure the FFT length is large enough, when high overlap percentages are required.");
}
Ok(overlap_keep)
}
/// Return a reasonable acoustic default with a frequency resolution around
/// ~ 10 Hz, where nfft is still an integer power of 2.
///
/// # Errors
///
/// If `fs` is something odd, i.e. < 1 kHz, or higher than 1 MHz.
///
pub fn reasonableAcousticDefault(fs: Flt, mode: ApsMode) -> Result<ApsSettings> {
if fs < 1e3 || fs > 1e6 {
bail!("Sampling frequency for reasonable acoustic data is >= 1 kHz and <= 1 MHz.");
}
let fs_div_10_rounded = (fs / 10.) as u32;
// 2^30 is about 1 million. We search for a two-power of an nfft that is
// the closest to fs/10. The frequency resolution is about fs/nfft.
let nfft = (0..30).map(|i| 2u32.pow(i) - fs_div_10_rounded).fold(
// Start wth a value that is always too large
fs as u32 * 10,
|cur, new| cur.min(new),
) as usize;
Ok(ApsSettings {
mode,
fs,
nfft,
windowType: WindowType::default(),
overlap: Overlap::default(),
freqWeightingType: FreqWeighting::default(),
})
}
/// Return sampling frequency
pub fn fs(&self) -> Flt {
self.fs
}
/// Return Nyquist frequency
pub fn fnyq(&self) -> Flt {
self.fs / 2.
}
/// Returns a single-sided frequency array corresponding to points in Power
/// spectra computation.
pub fn getFreq(&self) -> Array1<Flt> {
let df = self.fs / self.nfft as Flt;
let K = self.nfft / 2 + 1;
Array1::linspace(0., (K - 1) as Flt * df, K)
}
}

View File

@ -1,7 +1,7 @@
//! Compute forward single sided amplitude spectra //! Compute forward single sided amplitude spectra
use crate::config::*; use crate::config::*;
use realfft::{RealFftPlanner, RealToComplex}; use realfft::{RealFftPlanner, RealToComplex};
use std::sync::Arc; use std::{fmt::Debug, sync::Arc};
#[derive(Clone)] #[derive(Clone)]
pub struct FFT { pub struct FFT {
@ -14,6 +14,11 @@ pub struct FFT {
// nfft stored as float, this is how it is required most often // nfft stored as float, this is how it is required most often
nfftF: Flt, nfftF: Flt,
} }
impl Debug for FFT {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("Forward FFT engine, lenfth: {self.nfftF}").finish()
}
}
impl FFT { impl FFT {
/// Create new FFT from given nfft /// Create new FFT from given nfft

View File

@ -1,12 +1,13 @@
use crate::config::*; use crate::config::*;
use strum_macros::{Display, EnumMessage}; use strum::IntoEnumIterator;
use strum_macros::{Display, EnumIter, EnumMessage};
/// Sound level frequency weighting type (A, C, Z) /// Sound level frequency weighting type (A, C, Z)
// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy: // 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(eq, eq_int))]
// For now: // For now:
#[cfg_attr(feature = "python-bindings", pyclass)] #[cfg_attr(feature = "python-bindings", pyclass)]
#[derive(Display, Debug, EnumMessage, Default, Clone, PartialEq)] #[derive(Copy, Display, Debug, EnumMessage, Default, Clone, PartialEq, EnumIter)]
pub enum FreqWeighting { pub enum FreqWeighting {
/// A-weighting /// A-weighting
A, A,
@ -16,6 +17,27 @@ pub enum FreqWeighting {
#[default] #[default]
Z, Z,
} }
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl FreqWeighting {
#[staticmethod]
fn all() -> Vec<Self> {
Self::iter().collect()
}
fn __str__(&self) -> String {
format!("{self}-weighting")
}
fn letter(&self) -> String {
format!("{self}")
}
#[staticmethod]
#[pyo3(name = "default")]
fn default_py() -> Self {
Self::default()
}
}
#[cfg(test)] #[cfg(test)]
mod test { mod test {
use super::*; use super::*;

View File

@ -10,10 +10,16 @@ mod ps;
mod timebuffer; mod timebuffer;
mod window; mod window;
mod freqweighting; mod freqweighting;
mod apssettings;
mod overlap;
mod apsmode;
use crate::config::*; use crate::config::*;
pub use freqweighting::FreqWeighting; pub use freqweighting::FreqWeighting;
pub use aps::{ApsSettings, ApsSettingsBuilder,ApsMode, AvPowerSpectra, Overlap}; pub use overlap::Overlap;
pub use apssettings::{ApsSettings, ApsSettingsBuilder};
pub use apsmode::ApsMode;
pub use aps::AvPowerSpectra;
pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult}; pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult};
pub use window::{Window, WindowType}; pub use window::{Window, WindowType};

78
src/ps/overlap.rs Normal file
View File

@ -0,0 +1,78 @@
use crate::config::*;
/// Provide the overlap of blocks for computing averaged (cross) power spectra.
/// Can be provided as a percentage of the block size, or as a number of
/// samples.
#[cfg_attr(feature = "python-bindings", pyclass(frozen))]
#[derive(Clone, Debug, PartialEq)]
pub enum Overlap {
/// Overlap specified as a percentage of the total FFT length. Value should
/// be 0<=pct<100.
Percentage {
/// Percentage
pct: Flt,
},
/// Number of samples to overlap
Number {
/// N: Number of samples
N: usize,
},
/// No overlap at all, which is the same as Overlap::Number(0)
NoOverlap {},
}
impl Default for Overlap {
fn default() -> Self {
Overlap::Percentage { pct: 50. }
}
}
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl Overlap {
#[inline]
fn __eq__(&self, other: &Self) -> bool {
self == other
}
fn __str__(&self) -> String {
self.name()
}
#[staticmethod]
#[pyo3(name = "default")]
fn default_py() -> Self {
Self::default()
}
/// Export some typical settings to Python. Customs are also possible by
/// directly creating them. This is to fill up the list.
#[staticmethod]
fn some_settings() -> Vec<Overlap> {
use Overlap::*;
vec![
NoOverlap {},
Percentage { pct: 25. },
Percentage { pct: 50. },
Percentage { pct: 90. },
]
}
}
impl Overlap {
fn name(&self) -> String {
use Overlap::*;
match &self {
NoOverlap {} => "No overlap".into(),
Percentage { pct } => format! {"{} %", pct},
Number { N } => format! {"{} samples", N},
}
}
}
#[cfg(test)]
mod test {
use crate::ps::overlap::Overlap;
#[test]
fn test_overlap1() {
assert_eq!(Overlap::NoOverlap {}.name(), "No overlap");
}
}

View File

@ -85,6 +85,7 @@ impl CrossPowerSpecra for CPSResult {
/// example the computations of spectrograms, or Welch' method of spectral /// example the computations of spectrograms, or Welch' method of spectral
/// estimation. /// estimation.
/// ///
#[derive(Debug)]
pub struct PowerSpectra { pub struct PowerSpectra {
/// Window used in estimator. The actual Window in here is normalized with /// Window used in estimator. The actual Window in here is normalized with
/// the square root of the Window power. This safes one division when /// the square root of the Window power. This safes one division when
@ -140,13 +141,15 @@ impl PowerSpectra {
/// Compute FFTs of input channel data. Stores the scaled FFT data in /// Compute FFTs of input channel data. Stores the scaled FFT data in
/// self.freqdata. /// self.freqdata.
fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> ArrayView2<Cflt> { fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> ArrayView2<Cflt> {
assert!(timedata.nrows() > 0);
let (n, nch) = timedata.dim(); let (n, nch) = timedata.dim();
let nfft = self.nfft(); let nfft = self.nfft();
assert!(n == nfft); assert!(n == nfft);
// Make sure enough fft engines are available // Make sure enough fft engines are available
while nch > self.ffts.len() { while nch > self.ffts.len() {
self.ffts.push(self.ffts.last().unwrap().clone()); self.ffts
.push(self.ffts.last().expect("FFT's should not be empty").clone());
self.freqdata self.freqdata
.push_column(Ccol::from_vec(vec![Cflt::new(0., 0.); nfft / 2 + 1]).view()) .push_column(Ccol::from_vec(vec![Cflt::new(0., 0.); nfft / 2 + 1]).view())
.unwrap(); .unwrap();

View File

@ -8,7 +8,7 @@ use std::collections::VecDeque;
/// TimeBuffer, storage to add blocks of data in a ring buffer, that can be /// TimeBuffer, storage to add blocks of data in a ring buffer, that can be
/// extracted by blocks of other size. Also, we can keep samples in a buffer to /// extracted by blocks of other size. Also, we can keep samples in a buffer to
/// create, for example, overlapping windows of time data. /// create, for example, overlapping windows of time data.
#[derive(Default)] #[derive(Default, Debug)]
pub struct TimeBuffer { pub struct TimeBuffer {
data: Vec<VecDeque<Flt>>, data: Vec<VecDeque<Flt>>,
} }

View File

@ -1,7 +1,8 @@
#![allow(non_snake_case)] #![allow(non_snake_case)]
use crate::config::*; use crate::config::*;
use strum_macros::Display; use strum::IntoEnumIterator;
use strum_macros::{Display, EnumMessage, EnumIter};
/// Von Hann window, often misnamed as the 'Hanning' window. /// Von Hann window, often misnamed as the 'Hanning' window.
fn hann(nfft: usize) -> Dcol { fn hann(nfft: usize) -> Dcol {
@ -72,7 +73,11 @@ fn hamming(N: usize) -> Dcol {
/// * Blackman /// * Blackman
/// ///
/// The [WindowType::default] is [WindowType::Hann]. /// The [WindowType::default] is [WindowType::Hann].
#[derive(Display,Default, Copy, Clone, Debug)] #[derive(Display, Default, Copy, Clone, Debug, PartialEq, EnumMessage, EnumIter)]
// 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)]
pub enum WindowType { pub enum WindowType {
/// Von Hann window /// Von Hann window
#[default] #[default]
@ -87,8 +92,23 @@ pub enum WindowType {
Blackman = 4, Blackman = 4,
} }
#[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)]
impl WindowType {
#[staticmethod]
fn all() -> Vec<WindowType> {
WindowType::iter().collect()
}
#[staticmethod]
#[pyo3(name="default")]
fn default_py() -> Self {
Self::default()
}
}
/// Window (taper) computed from specified window type. /// Window (taper) computed from specified window type.
#[derive(Clone)] #[derive(Clone, Debug)]
pub struct Window { pub struct Window {
/// The enum from which it is generated /// The enum from which it is generated
pub w: WindowType, pub w: WindowType,

View File

@ -60,11 +60,11 @@
//! //!
//! ``` //! ```
//! //!
mod settings; mod slmsettings;
mod tw; mod tw;
mod slm; mod slm;
pub use slm::SLM; pub use slm::SLM;
pub use settings::{SLMSettings, SLMSettingsBuilder}; pub use slmsettings::{SLMSettings, SLMSettingsBuilder};
pub use tw::TimeWeighting; pub use tw::TimeWeighting;
pub use crate::ps::FreqWeighting; pub use crate::ps::FreqWeighting;

View File

@ -6,7 +6,7 @@ use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
use rayon::prelude::*; use rayon::prelude::*;
use smallvec::SmallVec; use smallvec::SmallVec;
use super::{settings::SLMSettings, SLM_MAX_CHANNELS}; use super::{slmsettings::SLMSettings, SLM_MAX_CHANNELS};
use crate::{config::*, filter::Filter}; use crate::{config::*, filter::Filter};
use crate::{Biquad, Dcol, Flt, FreqWeighting, PoleOrZero, SeriesBiquad, ZPKModel}; use crate::{Biquad, Dcol, Flt, FreqWeighting, PoleOrZero, SeriesBiquad, ZPKModel};
#[derive(Debug, Clone)] #[derive(Debug, Clone)]
@ -205,7 +205,7 @@ impl SLM {
} }
} }
#[cfg(feature="python-bindings")] #[cfg(feature = "python-bindings")]
#[cfg_attr(feature = "python-bindings", pymethods)] #[cfg_attr(feature = "python-bindings", pymethods)]
impl SLM { impl SLM {
#[new] #[new]
@ -214,8 +214,18 @@ impl SLM {
} }
#[pyo3(name = "run", signature=(dat, provide_output=true))] #[pyo3(name = "run", signature=(dat, provide_output=true))]
fn run_py(&mut self, dat: PyArrayLike1<Flt>, provide_output: bool) -> Option<Vec<Vec<Flt>>> { fn run_py<'py>(
self.run(dat.as_array().as_slice()?, provide_output) &mut self,
py: Python<'py>,
dat: PyArrayLike1<Flt>,
provide_output: bool,
) -> Option<Vec<Bound<'py, PyArray1<Flt>>>> {
if let Some(res) = self.run(dat.as_array().as_slice()?, provide_output) {
let vec_py_iter = res.into_iter().map(|v| PyArray1::from_vec_bound(py, v));
let vec_py = vec_py_iter.collect::<Vec<Bound<'py, PyArray1<Flt>>>>();
return Some(vec_py);
}
None
} }
#[pyo3(name = "Lmax")] #[pyo3(name = "Lmax")]

View File

@ -1,11 +1,14 @@
use crate::config::*; use crate::config::*;
use strum::EnumMessage;
use strum_macros::Display;
/// Time weighting to use in level detection of Sound Level Meter. /// Time weighting to use in level detection of Sound Level Meter.
/// ///
// Do the following when Pyo3 0.22 can finally be used combined with rust-numpy: // Do the following when Pyo3 0.22 can finally be used combined with rust-numpy:
// #[cfg_attr(feature = "python-bindings", pyclass(eq))] // #[cfg_attr(feature = "python-bindings", pyclass(eq))]
// For now: // For now:
#[cfg_attr(feature = "python-bindings", pyclass)] #[cfg_attr(feature = "python-bindings", pyclass)]
#[derive(Clone, Copy, PartialEq)] #[derive(Clone, Copy, Debug, PartialEq, Display)]
pub enum TimeWeighting { pub enum TimeWeighting {
// I know that the curly braces here are not required and add some // I know that the curly braces here are not required and add some
// boilerplate, but this is the only way Pyo3 swallows complex enums at the // boilerplate, but this is the only way Pyo3 swallows complex enums at the
@ -21,7 +24,7 @@ pub enum TimeWeighting {
/// Custom time constant [s] /// Custom time constant [s]
t: Flt, t: Flt,
}, },
/// A custom symmetric time weighting /// A custom symmetric time weighting
CustomAsymmetric { CustomAsymmetric {
/// Time weighting when level is increasing /// Time weighting when level is increasing
@ -30,6 +33,23 @@ pub enum TimeWeighting {
tdown: Flt, tdown: Flt,
}, },
} }
#[cfg_attr(feature = "python-bindings", pymethods)]
impl TimeWeighting {
// This method is still required in Pyo3 0.21, not anymore in 0.22
fn __eq__(&self, other: &Self) -> bool {
self == other
}
fn __str__(&self) -> String {
format!("{self}")
}
#[staticmethod]
fn all_standards() -> Vec<TimeWeighting> {
use TimeWeighting::*;
vec![Slow {}, Fast {}, Impulse {}]
}
}
impl Default for TimeWeighting { impl Default for TimeWeighting {
fn default() -> Self { fn default() -> Self {
TimeWeighting::Fast {} TimeWeighting::Fast {}
@ -75,3 +95,15 @@ impl TimeWeighting {
} }
} }
} }
#[cfg(test)]
mod test {
use crate::slm::TimeWeighting;
#[test]
fn test_tw() {
println!("Impulse: {}", TimeWeighting::Impulse {});
println!("Fast : {}", TimeWeighting::Fast {});
println!("Slow : {}", TimeWeighting::Slow {});
}
}