Compare commits
9 Commits
eb94785a89
...
016db92eba
Author | SHA1 | Date | |
---|---|---|---|
016db92eba | |||
96187bfcf9 | |||
2efb610caa | |||
65df1c82f6 | |||
826266b8ee | |||
08ecdf6dc4 | |||
42045da5cc | |||
c843c089dd | |||
d9fbe25dc1 |
0
.gitea/workflows/workflow.yml
Normal file
0
.gitea/workflows/workflow.yml
Normal file
@ -1,6 +1,6 @@
|
||||
[package]
|
||||
name = "lasprs"
|
||||
version = "0.6.1"
|
||||
version = "0.6.2"
|
||||
edition = "2021"
|
||||
authors = ["J.A. de Jong <j.a.dejong@ascee.nl>"]
|
||||
description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)"
|
||||
|
@ -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).
|
||||
|
||||
|
||||
## 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 "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
146
examples_py/test_aps.ipynb
Normal 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
|
||||
}
|
@ -41,8 +41,8 @@
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from numpy import log10, sqrt, exp, pi\n",
|
||||
"from lasprs._lasprs import StandardFilterDescriptor, SLMSettings, FreqWeighting, TimeWeighting, SLM\n",
|
||||
"from numpy import log10\n",
|
||||
"from lasprs import StandardFilterDescriptor\n",
|
||||
"def level(a):\n",
|
||||
" return 20*np.log10(np.abs(a))"
|
||||
]
|
||||
|
@ -41,8 +41,8 @@
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from numpy import log10, sqrt, exp, pi\n",
|
||||
"from lasprs._lasprs import ZPKModel, FreqWeighting\n",
|
||||
"from numpy import log10\n",
|
||||
"from lasprs import ZPKModel, FreqWeighting\n",
|
||||
"def level(a):\n",
|
||||
" return 20*np.log10(np.abs(a))"
|
||||
]
|
||||
|
@ -17,6 +17,7 @@
|
||||
"source": [
|
||||
"# Prerequisites, uncomment below in case of errors. Also for ipympl, restart Jupyter Lab if it was not installed\n",
|
||||
"\n",
|
||||
"# Only do the following if you are in develop mode\n",
|
||||
"!cd .. && maturin develop -F python-bindings\n",
|
||||
"#!pip install ipympl scipy matplotlib"
|
||||
]
|
||||
@ -41,7 +42,7 @@
|
||||
"source": [
|
||||
"import numpy as np\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",
|
||||
"def level(a):\n",
|
||||
" return 20*np.log10(np.abs(a))"
|
||||
@ -105,27 +106,7 @@
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"t = np.linspace(0, N/fs, N, endpoint=False)\n",
|
||||
"out = slm.run(inp)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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)"
|
||||
"out = slm.run(inp, True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -141,14 +122,6 @@
|
||||
"plt.ylim(-60, 0)\n",
|
||||
"plt.legend(names)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "4303e20e-cc61-4d23-a5f0-cab683d55912",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
@ -1 +0,0 @@
|
||||
from .._lasprs import (Biquad, BiquadBank, SeriesBiquad, ZPKModel, FilterSpec)
|
@ -26,7 +26,8 @@ cfg_if::cfg_if! {
|
||||
|
||||
cfg_if::cfg_if! {
|
||||
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::exceptions::PyValueError;
|
||||
pub use pyo3::{pymodule, types::PyModule, PyResult};
|
||||
|
@ -1,5 +1,5 @@
|
||||
use crate::{config::*, ZPKModel};
|
||||
use anyhow::{anyhow, bail, Result};
|
||||
use anyhow::{anyhow, bail, Error, Result};
|
||||
use num::{traits::float, Float};
|
||||
use rayon::iter::Filter;
|
||||
use softfloat::F64;
|
||||
@ -497,6 +497,38 @@ impl StandardFilterDescriptor {
|
||||
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]
|
||||
#[pyo3(name = "genOctaveFilterSet")]
|
||||
fn genOctaveFilterSetFromFreq(fmin: Flt, fmax: Flt) -> PyResult<Vec<StandardFilterDescriptor>> {
|
||||
@ -527,8 +559,18 @@ impl StandardFilterDescriptor {
|
||||
fmax: Flt,
|
||||
append_overall: bool,
|
||||
) -> 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)?)
|
||||
}
|
||||
|
||||
#[pyo3(name="fm")]
|
||||
fn fm_py(&self) -> Option<Flt> {
|
||||
self.fm()
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
|
@ -569,7 +569,6 @@ impl ZPKModel {
|
||||
let p2 = 2. * pi * f2;
|
||||
let p3 = 2. * pi * f3;
|
||||
let p4 = 2. * pi * f4;
|
||||
println!("{b} {p1}, {p2}, {p3}, {p4}");
|
||||
|
||||
let (zeros, poles) = match wt {
|
||||
FreqWeighting::Z => {
|
||||
@ -610,11 +609,21 @@ pub enum PoleOrZero {
|
||||
|
||||
#[cfg(test)]
|
||||
mod test{
|
||||
use approx::assert_abs_diff_eq;
|
||||
use num::complex::ComplexFloat;
|
||||
|
||||
use crate::TransferFunction;
|
||||
use super::ZPKModel;
|
||||
|
||||
|
||||
#[test]
|
||||
fn test_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);
|
||||
}
|
||||
}
|
||||
|
@ -62,6 +62,12 @@ fn lasprs(m: &Bound<'_, PyModule>) -> PyResult<()> {
|
||||
m.add_class::<ps::FreqWeighting>()?;
|
||||
m.add_class::<slm::SLMSettings>()?;
|
||||
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(())
|
||||
}
|
||||
|
299
src/ps/aps.rs
299
src/ps/aps.rs
@ -1,214 +1,25 @@
|
||||
use super::timebuffer::TimeBuffer;
|
||||
use super::CrossPowerSpecra;
|
||||
use super::*;
|
||||
use crate::config::*;
|
||||
use super::{timebuffer::TimeBuffer, CrossPowerSpecra};
|
||||
use crate::{config::*, TransferFunction, ZPKModel};
|
||||
use anyhow::{bail, Error, Result};
|
||||
use derive_builder::Builder;
|
||||
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
|
||||
/// Used to compute power spectra estimations on
|
||||
/// long datasets, where nfft << length of data. This way, the variance of a
|
||||
/// single periodogram is suppressed with increasing number of averages.
|
||||
///
|
||||
/// For more information, see the book on numerical recipes.
|
||||
///
|
||||
|
||||
#[cfg(feature = "python-bindings")]
|
||||
#[cfg_attr(feature = "python-bindings", pyclass)]
|
||||
#[derive(Debug)]
|
||||
|
||||
pub struct AvPowerSpectra {
|
||||
// Power spectra estimator for single block
|
||||
ps: PowerSpectra,
|
||||
|
||||
// Settings for computing power spectra, see [ApsSettings]
|
||||
settings: ApsSettings,
|
||||
|
||||
// The number of samples to keep in the time buffer when overlapping time
|
||||
@ -221,6 +32,11 @@ pub struct AvPowerSpectra {
|
||||
/// Storage for sample data.
|
||||
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
|
||||
cur_est: CPSResult,
|
||||
}
|
||||
@ -254,8 +70,7 @@ impl AvPowerSpectra {
|
||||
/// - When providing invalid sampling frequencies
|
||||
///
|
||||
pub fn new_simple_all_averaging(fs: Flt, nfft: usize) -> AvPowerSpectra {
|
||||
let mut settings =
|
||||
ApsSettings::reasonableAcousticDefault(fs, ApsMode::AllAveraging).unwrap();
|
||||
let mut settings = ApsSettings::reasonableAcousticDefault(fs, ApsMode::default()).unwrap();
|
||||
settings.nfft = nfft;
|
||||
AvPowerSpectra::new(settings)
|
||||
}
|
||||
@ -285,19 +100,51 @@ impl AvPowerSpectra {
|
||||
|
||||
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 {
|
||||
ps,
|
||||
overlap_keep,
|
||||
settings,
|
||||
N: 0,
|
||||
freqWeighting_pwr,
|
||||
cur_est: CPSResult::default((0, 0, 0)),
|
||||
timebuf: TimeBuffer::new(),
|
||||
}
|
||||
}
|
||||
// Update result for single block
|
||||
fn update_singleblock(&mut self, timedata: ArrayView2<Flt>) {
|
||||
let Cpsnew = self.ps.compute(timedata);
|
||||
// println!("Cpsnew: {:?}", Cpsnew[[0, 0, 0]]);
|
||||
// Compute updated block of cross-spectral density
|
||||
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
|
||||
if self.cur_est.is_empty() {
|
||||
@ -310,7 +157,7 @@ impl AvPowerSpectra {
|
||||
|
||||
// Apply operation based on mode
|
||||
match self.settings.mode {
|
||||
ApsMode::AllAveraging => {
|
||||
ApsMode::AllAveraging {} => {
|
||||
let Nf = Cflt {
|
||||
re: self.N as Flt,
|
||||
im: 0.,
|
||||
@ -359,7 +206,7 @@ impl AvPowerSpectra {
|
||||
}
|
||||
}
|
||||
|
||||
ApsMode::Spectrogram => {
|
||||
ApsMode::Spectrogram {} => {
|
||||
self.cur_est = Cpsnew;
|
||||
}
|
||||
}
|
||||
@ -395,10 +242,11 @@ impl AvPowerSpectra {
|
||||
computed_single = true;
|
||||
}
|
||||
if computed_single {
|
||||
return Some(&self.cur_est);
|
||||
}
|
||||
Some(&self.cur_est)
|
||||
} else {
|
||||
None
|
||||
}
|
||||
}
|
||||
|
||||
/// Computes average (cross)power spectra, and returns all intermediate
|
||||
/// estimates that can be calculated. This is useful when plotting spectra
|
||||
@ -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)]
|
||||
mod test {
|
||||
use approx::assert_abs_diff_eq;
|
||||
@ -450,11 +321,11 @@ mod test {
|
||||
#[test]
|
||||
fn test_overlap_keep() {
|
||||
let ol = [
|
||||
Overlap::NoOverlap,
|
||||
Percentage(50.),
|
||||
Percentage(50.),
|
||||
Percentage(25.),
|
||||
Overlap::Number(10),
|
||||
Overlap::NoOverlap {},
|
||||
Percentage { pct: 50. },
|
||||
Percentage { pct: 50. },
|
||||
Percentage { pct: 25. },
|
||||
Overlap::Number { N: 10 },
|
||||
];
|
||||
let nffts = [10, 10, 1024, 10];
|
||||
let expected_keep = [0, 5, 512, 2, 10];
|
||||
@ -465,7 +336,7 @@ mod test {
|
||||
.fs(1.)
|
||||
.overlap(overlap.clone())
|
||||
.build()
|
||||
.unwrap();
|
||||
.expect("BUG: Settings cannot be build.");
|
||||
|
||||
assert_eq!(settings.get_overlap_keep(), *expected);
|
||||
}
|
||||
@ -481,15 +352,15 @@ mod test {
|
||||
let settings = ApsSettingsBuilder::default()
|
||||
.fs(fs)
|
||||
.nfft(nfft)
|
||||
.overlap(Overlap::NoOverlap)
|
||||
.overlap(Overlap::NoOverlap {})
|
||||
.mode(ApsMode::ExponentialWeighting { tau })
|
||||
.build()
|
||||
.unwrap();
|
||||
.expect("Settings cannot be empty");
|
||||
let overlap_keep = settings.get_overlap_keep();
|
||||
let mut aps = AvPowerSpectra::new(settings);
|
||||
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_zeros = Dmat::zeros((nfft, 1));
|
||||
|
||||
|
35
src/ps/apsmode.rs
Normal file
35
src/ps/apsmode.rs
Normal 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
168
src/ps/apssettings.rs
Normal 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)
|
||||
}
|
||||
}
|
@ -1,7 +1,7 @@
|
||||
//! Compute forward single sided amplitude spectra
|
||||
use crate::config::*;
|
||||
use realfft::{RealFftPlanner, RealToComplex};
|
||||
use std::sync::Arc;
|
||||
use std::{fmt::Debug, sync::Arc};
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct FFT {
|
||||
@ -14,6 +14,11 @@ pub struct FFT {
|
||||
// nfft stored as float, this is how it is required most often
|
||||
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 {
|
||||
/// Create new FFT from given nfft
|
||||
|
@ -1,12 +1,13 @@
|
||||
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)
|
||||
|
||||
// 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)]
|
||||
#[derive(Copy, Display, Debug, EnumMessage, Default, Clone, PartialEq, EnumIter)]
|
||||
pub enum FreqWeighting {
|
||||
/// A-weighting
|
||||
A,
|
||||
@ -16,6 +17,27 @@ pub enum FreqWeighting {
|
||||
#[default]
|
||||
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)]
|
||||
mod test {
|
||||
use super::*;
|
||||
|
@ -10,10 +10,16 @@ mod ps;
|
||||
mod timebuffer;
|
||||
mod window;
|
||||
mod freqweighting;
|
||||
mod apssettings;
|
||||
mod overlap;
|
||||
mod apsmode;
|
||||
use crate::config::*;
|
||||
|
||||
|
||||
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 window::{Window, WindowType};
|
||||
|
78
src/ps/overlap.rs
Normal file
78
src/ps/overlap.rs
Normal 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");
|
||||
}
|
||||
}
|
@ -85,6 +85,7 @@ impl CrossPowerSpecra for CPSResult {
|
||||
/// example the computations of spectrograms, or Welch' method of spectral
|
||||
/// estimation.
|
||||
///
|
||||
#[derive(Debug)]
|
||||
pub struct PowerSpectra {
|
||||
/// Window used in estimator. The actual Window in here is normalized with
|
||||
/// 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
|
||||
/// self.freqdata.
|
||||
fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> ArrayView2<Cflt> {
|
||||
assert!(timedata.nrows() > 0);
|
||||
let (n, nch) = timedata.dim();
|
||||
let nfft = self.nfft();
|
||||
assert!(n == nfft);
|
||||
|
||||
// Make sure enough fft engines are available
|
||||
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
|
||||
.push_column(Ccol::from_vec(vec![Cflt::new(0., 0.); nfft / 2 + 1]).view())
|
||||
.unwrap();
|
||||
|
@ -8,7 +8,7 @@ use std::collections::VecDeque;
|
||||
/// 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
|
||||
/// create, for example, overlapping windows of time data.
|
||||
#[derive(Default)]
|
||||
#[derive(Default, Debug)]
|
||||
pub struct TimeBuffer {
|
||||
data: Vec<VecDeque<Flt>>,
|
||||
}
|
||||
|
@ -1,7 +1,8 @@
|
||||
#![allow(non_snake_case)]
|
||||
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.
|
||||
fn hann(nfft: usize) -> Dcol {
|
||||
@ -72,7 +73,11 @@ fn hamming(N: usize) -> Dcol {
|
||||
/// * Blackman
|
||||
///
|
||||
/// 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 {
|
||||
/// Von Hann window
|
||||
#[default]
|
||||
@ -87,8 +92,23 @@ pub enum WindowType {
|
||||
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.
|
||||
#[derive(Clone)]
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct Window {
|
||||
/// The enum from which it is generated
|
||||
pub w: WindowType,
|
||||
|
@ -60,11 +60,11 @@
|
||||
//!
|
||||
//! ```
|
||||
//!
|
||||
mod settings;
|
||||
mod slmsettings;
|
||||
mod tw;
|
||||
mod slm;
|
||||
pub use slm::SLM;
|
||||
pub use settings::{SLMSettings, SLMSettingsBuilder};
|
||||
pub use slmsettings::{SLMSettings, SLMSettingsBuilder};
|
||||
pub use tw::TimeWeighting;
|
||||
pub use crate::ps::FreqWeighting;
|
||||
|
||||
|
@ -6,7 +6,7 @@ use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
|
||||
use rayon::prelude::*;
|
||||
use smallvec::SmallVec;
|
||||
|
||||
use super::{settings::SLMSettings, SLM_MAX_CHANNELS};
|
||||
use super::{slmsettings::SLMSettings, SLM_MAX_CHANNELS};
|
||||
use crate::{config::*, filter::Filter};
|
||||
use crate::{Biquad, Dcol, Flt, FreqWeighting, PoleOrZero, SeriesBiquad, ZPKModel};
|
||||
#[derive(Debug, Clone)]
|
||||
@ -205,7 +205,7 @@ impl SLM {
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(feature="python-bindings")]
|
||||
#[cfg(feature = "python-bindings")]
|
||||
#[cfg_attr(feature = "python-bindings", pymethods)]
|
||||
impl SLM {
|
||||
#[new]
|
||||
@ -214,8 +214,18 @@ impl SLM {
|
||||
}
|
||||
|
||||
#[pyo3(name = "run", signature=(dat, provide_output=true))]
|
||||
fn run_py(&mut self, dat: PyArrayLike1<Flt>, provide_output: bool) -> Option<Vec<Vec<Flt>>> {
|
||||
self.run(dat.as_array().as_slice()?, provide_output)
|
||||
fn run_py<'py>(
|
||||
&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")]
|
||||
|
@ -1,11 +1,14 @@
|
||||
use crate::config::*;
|
||||
|
||||
use strum::EnumMessage;
|
||||
use strum_macros::Display;
|
||||
/// 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:
|
||||
// #[cfg_attr(feature = "python-bindings", pyclass(eq))]
|
||||
// For now:
|
||||
#[cfg_attr(feature = "python-bindings", pyclass)]
|
||||
#[derive(Clone, Copy, PartialEq)]
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Display)]
|
||||
pub enum TimeWeighting {
|
||||
// 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
|
||||
@ -30,6 +33,23 @@ pub enum TimeWeighting {
|
||||
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 {
|
||||
fn default() -> Self {
|
||||
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 {});
|
||||
}
|
||||
}
|
||||
|
Loading…
Reference in New Issue
Block a user