ASCEE/lasprs
5
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Bugfixing on SLM. First tests. Code seems to work, but might still have subtle bugs (as timeweighting was wrong previous to this commit).

Browse Source
This commit is contained in:
Anne de Jong 2024-08-19 23:02:13 +02:00
parent 6208e97f8a
commit d8814e6c11
13 changed files with 257 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/filter/mod.rs
View File

@ -14,7 +14,7 @@ mod zpkmodel;
mod butter;
mod octave;
pub use super::ps::FreqWeightingType;
pub use super::ps::FreqWeighting;
pub use biquad::Biquad;
pub use biquadbank::BiquadBank;
pub use octave::{StandardFilterDescriptor, G, FREQ_REF};

79
src/filter/octave.rs
View File

@ -1,6 +1,5 @@
use crate::{Flt, ZPKModel};
use anyhow::{anyhow, bail, Result};
use clap::Error;
use num::{traits::float, Float};
use rayon::iter::Filter;
use softfloat::F64;
@ -67,11 +66,11 @@ pub const FREQ_REF: Flt = 1000.;
///
/// ```rust
/// use lasprs::filter::*;
/// fn main() -> anyhow::Result<()> {
/// # fn main() -> anyhow::Result<()> {
/// let desc = StandardFilterDescriptor::Octave("16")?;
/// let filter = desc.genFilter().bilinear(48e3);
/// Ok(())
/// }
/// # Ok(())
/// # }
/// ```
///
/// ## Create a one-third octave band bandpass filter that has the frequency of 42 in its pass-band
@ -79,11 +78,27 @@ pub const FREQ_REF: Flt = 1000.;
/// ```rust
///
/// use lasprs::filter::*;
/// fn main() -> anyhow::Result<()> {
/// # fn main() -> anyhow::Result<()> {
/// let desc = StandardFilterDescriptor::filterForFreq(3, 42.)?;
/// let filter = desc.genFilter().bilinear(48e3);
/// Ok(())
/// }
/// # Ok(())
/// # }
/// ```
///
/// ## Create a set of octave band filters
///
/// ```rust
/// use lasprs::filter::*;
/// # fn main() -> anyhow::Result<()> {
/// // Generate custom set..
/// let desc = StandardFilterDescriptor::genOctaveFilterSet("16", "16k").unwrap();
/// // Or, when lazy: just generate the full set
/// let desc = StandardFilterDescriptor::fullOctaveFilterSet();
/// let filters: Vec<SeriesBiquad> = desc.iter()
/// .map(|desc| desc.genFilter().bilinear(48e3))
/// .collect();
/// # Ok(())
/// # }
/// ```
#[derive(PartialEq, Clone, Debug)]
pub struct StandardFilterDescriptor {
@ -222,38 +237,40 @@ impl StandardFilterDescriptor {
}
/// Creates a set of octave filters.
pub fn genOctaveFilterSet<T>(low_f: Option<T>, high_f: Option<T>) -> Result<Vec<Self>>
pub fn genOctaveFilterSet<T>(low_f: T, high_f: T) -> Result<Vec<Self>>
where
T: TryInto<OctaveBandDescriptor, Error = Error>,
T: TryInto<OctaveBandDescriptor, Error = anyhow::Error>,
{
let xmin = if let Some(low_f) = low_f {
low_f.try_into()?.x
} else {
-OCTAVE_NAMES_OFFSET as i32
};
let xmax = if let Some(high_f) = high_f {
high_f.try_into()?.x
} else {
(OCTAVE_NOMINAL_MIDBAND_NAMES.len() - 1) as i32
};
let xmin = low_f.try_into()?.x;
let xmax = high_f.try_into()?.x;
Ok((xmin..=xmax).map(|x| Self::Octave(x).unwrap()).collect())
}
/// Generate a full third octave bandpass filter set
pub fn fullThirdOctaveFilterSet() -> Vec<Self> {
Self::genThirdOctaveFilterSet(
THIRDOCTAVE_NOMINAL_MIDBAND_NAMES[0],
*(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.last().unwrap()),
)
.unwrap()
}
/// Generate a full octave bandpass filter set
pub fn fullOctaveFilterSet() -> Vec<Self> {
Self::genOctaveFilterSet(
OCTAVE_NOMINAL_MIDBAND_NAMES[0],
*(OCTAVE_NOMINAL_MIDBAND_NAMES.last().unwrap()),
)
.unwrap()
}
/// Creates a set of one-third octave bandpass filters.
pub fn genThirdOctaveFilterSet<T>(low_f: Option<T>, high_f: Option<T>) -> Result<Vec<Self>>
pub fn genThirdOctaveFilterSet<T>(low_f: T, high_f: T) -> Result<Vec<Self>>
where
T: TryInto<ThirdOctaveBandDescriptor, Error = Error>,
T: TryInto<ThirdOctaveBandDescriptor, Error = anyhow::Error>,
{
let xmin = if let Some(low_f) = low_f {
low_f.try_into()?.x
} else {
-THIRDOCTAVE_NAMES_OFFSET as i32
};
let xmax = if let Some(high_f) = high_f {
high_f.try_into()?.x
} else {
(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() - 1) as i32
};
let xmin = low_f.try_into()?.x;
let xmax = high_f.try_into()?.x;
Ok((xmin..=xmax)
.map(|x| Self::ThirdOctave(x).unwrap())
.collect())

1
src/filter/seriesbiquad.rs
View File

@ -103,6 +103,7 @@ impl SeriesBiquad {
SeriesBiquad::new(filter_coefs).unwrap()
}
#[allow(dead_code)]
fn clone_dyn(&self) -> Box<dyn Filter> {
Box::new(self.clone())
}

26
src/filter/zpkmodel.rs
View File

@ -1,6 +1,6 @@
use super::butter::butter_lowpass_roots;
use super::{Biquad, SeriesBiquad, TransferFunction};
use crate::{config::*, ps::FreqWeightingType};
use crate::{config::*, ps::FreqWeighting};
use itertools::{EitherOrBoth, Itertools};
use num::{zero, Complex};
use std::cmp::{max, min};
@ -73,7 +73,7 @@ pub enum FilterSpec {
/// poles should either be real, or come in complex conjugate pairs. This is
/// enforced by the way the poles and zero's are internally stored.
///
#[derive(Clone, Debug, Default)]
#[derive(Clone, Debug)]
#[cfg_attr(feature = "python-bindings", pyclass)]
pub struct ZPKModel {
// List of zeros
@ -87,6 +87,12 @@ pub struct ZPKModel {
// transform to create digital filter of this analogue one.
fwarp: Option<Flt>,
}
impl Default for ZPKModel {
fn default() -> Self {
ZPKModel{z:vec![], p:vec![], k: 1.0, fwarp: None}
}
}
impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for ZPKModel {
fn tf(&self, _fs: Flt, freq: T) -> Ccol {
let freq = freq.into();
@ -494,26 +500,26 @@ impl ZPKModel {
///
/// # Args
///
/// - `wt` - `[FreqWeightingType]` to use. i.e. A-weighting.
/// - `wt` - `[FreqWeighting]` to use. i.e. A-weighting.
///
/// # Examples
///
/// ## Get part of pulse response of digital A-filter at 48 kHz
///
/// ```
/// use lasprs::filter::{ZPKModel, FreqWeightingType, Filter};
/// use lasprs::filter::{ZPKModel, FreqWeighting, Filter};
///
/// // Sampling frequency in Hz
/// let fs = 48000.;
///
/// let mut afilter = ZPKModel::freqWeightingFilter(FreqWeightingType::A).bilinear(fs);
/// let mut afilter = ZPKModel::freqWeightingFilter(FreqWeighting::A).bilinear(fs);
/// let mut data = [0.; 1000];
/// data[0] = 1.0;
/// let out = afilter.filter(&data);
/// ```
///
pub fn freqWeightingFilter(wt: FreqWeightingType) -> ZPKModel {
if let FreqWeightingType::Z = wt {
pub fn freqWeightingFilter(wt: FreqWeighting) -> ZPKModel {
if let FreqWeighting::Z = wt {
return ZPKModel {
z: vec![],
p: vec![],
@ -548,15 +554,15 @@ impl ZPKModel {
let p4 = 2. * pi * f4;
let (zeros, poles) = match wt {
FreqWeightingType::Z => {
FreqWeighting::Z => {
unreachable!()
}
FreqWeightingType::C => {
FreqWeighting::C => {
let zeros = vec![PoleOrZero::Real2(0., 0.)];
let poles = vec![PoleOrZero::Real2(-p1, -p1), PoleOrZero::Real2(-p4, -p4)];
(zeros, poles)
}
FreqWeightingType::A => {
FreqWeighting::A => {
let poles = vec![
PoleOrZero::Real2(-p1, -p1),
PoleOrZero::Real2(-p2, -p3),

8
src/ps/aps.rs
View File

@ -4,7 +4,7 @@ use super::*;
use crate::config::*;
use anyhow::{bail, Error, Result};
use derive_builder::Builder;
use freqweighting::FreqWeightingType;
use freqweighting::FreqWeighting;
/// All settings used for computing averaged power spectra using Welch' method.
#[derive(Builder, Clone)]
@ -21,7 +21,7 @@ pub struct ApsSettings {
windowType: WindowType,
/// Kind of freqency weighting. Defaults to Z
#[builder(default)]
freqWeightingType: FreqWeightingType,
freqWeightingType: FreqWeighting,
/// FFT Length
nfft: usize,
/// Sampling frequency
@ -71,7 +71,7 @@ impl ApsSettings {
self.validate_get_overlap_keep().unwrap()
}
/// Unpack all, returns parts in tuple
pub fn get(self) -> (ApsMode, Overlap, WindowType, FreqWeightingType, usize, Flt) {
pub fn get(self) -> (ApsMode, Overlap, WindowType, FreqWeighting, usize, Flt) {
(
self.mode,
self.overlap,
@ -134,7 +134,7 @@ impl ApsSettings {
nfft,
windowType: WindowType::default(),
overlap: Overlap::default(),
freqWeightingType: FreqWeightingType::default(),
freqWeightingType: FreqWeighting::default(),
})
}

2
src/ps/fft.rs
View File

@ -11,11 +11,13 @@ pub struct FFT {
timescratch: Vec<Flt>,
// rounded down nfft/2
half_nfft_rounded: usize,
// nfft stored as float, this is how it is required most often
nfftF: Flt,
}
impl FFT {
/// Create new FFT from given nfft
#[allow(dead_code)]
pub fn newFromNFFT(nfft: usize) -> FFT {
let mut planner = RealFftPlanner::<Flt>::new();
let fft = planner.plan_fft_forward(nfft);

8
src/ps/freqweighting.rs
View File

@ -1,7 +1,7 @@
use strum_macros::{Display, EnumMessage};
/// Sound level frequency weighting type (A, C, Z)
#[derive(Display, Debug, EnumMessage, Default, Clone)]
pub enum FreqWeightingType {
pub enum FreqWeighting {
/// A-weighting
A,
/// C-weighting
@ -15,9 +15,9 @@ mod test {
use super::*;
#[test]
fn test() {
let a = FreqWeightingType::A;
let c = FreqWeightingType::C;
let z = FreqWeightingType::Z;
let a = FreqWeighting::A;
let c = FreqWeighting::C;
let z = FreqWeighting::Z;
println!("A-weighting: {a}");
println!("C-weighting: {c}");
println!("Z-weighting: {z}");

2
src/ps/mod.rs
View File

@ -13,7 +13,7 @@ mod freqweighting;
use crate::config::*;
pub use freqweighting::FreqWeightingType;
pub use freqweighting::FreqWeighting;
pub use aps::{ApsSettings, ApsSettingsBuilder,ApsMode, AvPowerSpectra, Overlap};
pub use ps::{CrossPowerSpecra, PowerSpectra, CPSResult};
pub use window::{Window, WindowType};

64
src/slm/mod.rs
View File

@ -3,9 +3,69 @@
//! Provides structs and helpers (SLMBuilder) for creating configurated Sound
//! Level Meters.
//!
//! # Usage examples
//!
//! ## Simple over-all level SLM
//!
//! This example creates a simple SLM that processes at 48kHz, uses fast time
//! weighting and A frequency weighting:
//!
//! ```
//! # use anyhow::Result;
//! use lasprs::slm::*;
//! use lasprs::filter::StandardFilterDescriptor;
//! use ndarray::Array1;
//! # fn main() -> Result<()> {
//!
//! // Generate an overall filter (no filter at all)
//! let desc = StandardFilterDescriptor::Overall().unwrap();
//!
//! // Generate settings
//! let settings = SLMSettingsBuilder::default()
//! .fs(48e3)
//! .freqWeighting(FreqWeighting::A)
//! .timeWeighting(TimeWeighting::Fast)
//! .filterDescriptors(&[desc]).build().unwrap();
//!
//! let mut slm = SLM::new(settings);
//! // Generate some data. Yes, this is not the most spectacular set
//! let mut data = Array1::zeros(48000);
//! data[0] = 1.;
//!
//! // Now apply some data. This is a kind of the SLM-s impulse response
//! let res = slm.run(&data.as_slice().unwrap()).unwrap();
//!
//! // Only one channel of result data
//! assert_eq!(res.len(), 1);
//!
//! let res = &res[0];
//! println!("Data is: {res:#?}");
//!
//! // Get the equivalent level:
//! let Leq = slm.Leq()[0];
//!
//! # Ok::<() , anyhow::Error>(())
//! # }
//! ```
//!
//! ## One-third octave band, plus overall
//!
//! ```
//! use lasprs::filter::StandardFilterDescriptor;
//! // Generate the default set of one-third octave band filters
//! let mut desc = StandardFilterDescriptor::fullThirdOctaveFilterSet();
//! desc.push(StandardFilterDescriptor::Overall().unwrap());
//!
//! // Rest of code is the same as in previous example
//!
//! ```
//!
mod settings;
mod tw;
mod slm;
pub use slm::SLM;
pub use settings::SLMSettings;
pub use tw::TimeWeightingType;
pub use settings::{SLMSettings, SLMSettingsBuilder};
pub use tw::TimeWeighting;
pub use crate::ps::FreqWeighting;
const SLM_MAX_CHANNELS: usize = 64;

47
src/slm/settings.rs
View File

@ -1,19 +1,46 @@
use super::{TimeWeighting, SLM_MAX_CHANNELS};
use crate::{filter::StandardFilterDescriptor, Flt, FreqWeighting};
use anyhow::Result;
use clap::builder;
use derive_builder::Builder;
use smallvec::SmallVec;
use crate::{Flt, FreqWeightingType, filter::StandardFilterDescriptor};
use super::TimeWeightingType;
use smallvec::{smallvec, SmallVec};
/// Settings used to create a Sound Level Meter.
#[derive(Builder, Clone)]
#[builder(setter(into))]
pub struct SLMSettings {
/// Sampling frequency in \[Hz\]
pub fs: Flt,
/// Reference level, in units of measured quantity. For sound pressure in
/// air, this is typically 20 μPa. This is also the default value.
#[builder(default = "2e-5")]
pub Lref: Flt,
pub freqWeighting: FreqWeightingType,
pub timeWeighting: TimeWeightingType,
pub filterDescriptors: SmallVec<[StandardFilterDescriptor; 64]>,
/// Frequency weightin A/C/Z applied to data. Defaults to [FreqWeighting::default()].
#[builder(default)]
pub freqWeighting: FreqWeighting,
/// For time-dependent output, the time weighthing applied (Fast / Slow)
pub timeWeighting: TimeWeighting,
/// Descriptors for the filters, maximum of 64, which is a reasonable amount
/// and - if all used - might already choke a computer.
pub filterDescriptors: Vec<StandardFilterDescriptor>,
}
impl SLMSettings {
#[cfg(test)]
mod test {
use super::*;
use anyhow::Result;
#[test]
fn test_slmsettings1() -> Result<()> {
let desc = StandardFilterDescriptor::genFilterSetInRange(1, 100., 5e3, true).unwrap();
let _ = SLMSettingsBuilder::default()
.fs(1e3)
.freqWeighting(FreqWeighting::A)
.timeWeighting(TimeWeighting::Slow)
.filterDescriptors(desc)
.build()
.unwrap();
Ok(())
}
}

57
src/slm/slm.rs
View File

@ -5,23 +5,30 @@ use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
use rayon::prelude::*;
use smallvec::SmallVec;
use super::settings::SLMSettings;
use super::{settings::SLMSettings, SLM_MAX_CHANNELS};
use crate::{config::*, filter::Filter};
use crate::{Biquad, Dcol, Flt, FreqWeightingType, PoleOrZero, SeriesBiquad, ZPKModel};
use crate::{Biquad, Dcol, Flt, FreqWeighting, PoleOrZero, SeriesBiquad, ZPKModel};
#[derive(Debug, Clone)]
struct SLMChannel {
// Statistics to store
stat: SLMStat,
// The bandpass filter for this channel
bp: SeriesBiquad,
// The rectifier filter
rect_lowpass_up: Biquad,
// The rectifier filter for decreasing values. Only for asymmetric time
// weighting (case of impulse weighting: [TimeWeighting::Impulse])
rect_lowpass_down: Option<Biquad>,
}
/// Sound Level Meter
#[derive(Debug, Clone)]
pub struct SLM {
// Number of samples processed after last run() is called.
N: usize,
Lrefsq: Flt,
prefilter: SeriesBiquad,
channels: SmallVec<[SLMChannel; 64]>,
channels: SmallVec<[SLMChannel; SLM_MAX_CHANNELS]>,
}
impl SLM {
@ -179,7 +186,7 @@ impl SLM {
}
}
#[derive(Clone, Default)]
#[derive(Debug, Clone, Default)]
/// Quantities defined as powers, i.e. square of amplitude
struct SLMStat {
// Max signal power
@ -192,3 +199,45 @@ struct SLMStat {
// Last obtained signal power, after last time run() is called.
Pt_last: Flt,
}
#[cfg(test)]
mod test {
use crate::{
siggen::Siggen,
slm::{SLMSettingsBuilder, TimeWeighting},
Flt, FreqWeighting, StandardFilterDescriptor,
};
use super::SLM;
#[test]
fn test_slm1() {
const fs: Flt = 48e3;
const N: usize = (fs/8.) as usize;
let desc = StandardFilterDescriptor::Overall().unwrap();
let settings = SLMSettingsBuilder::default()
.fs(fs)
.timeWeighting(TimeWeighting::Fast)
.freqWeighting(FreqWeighting::Z)
.filterDescriptors(&[desc])
.build()
.unwrap();
let mut siggen = Siggen::newSine(1, 1000.);
siggen.setAllMute(false);
siggen.reset(fs);
let mut data = vec![0.; N];
siggen.genSignal(&mut data);
// println!("{:#?}", data);
let mut slm = SLM::new(settings);
// println!("{slm:#?}");
let res = slm.run(&data).unwrap();
let res = &res[0];
println!("{slm:#?}");
println!("{:#?}", &res[res.len()- 100..]);
}
}

29
src/slm/tw.rs
View File

@ -1,14 +1,16 @@
use crate::Flt;
/// Time weighting to use in level detection of Sound Level Meter.
#[derive(Clone, Copy)]
pub enum TimeWeightingType {
/// Slow time weighting
pub enum TimeWeighting {
/// Slow time weighting ~ 1 s
Slow,
/// Fast time weighting
/// Fast time weighting ~ 1/8 s
Fast,
/// Impulse time weighting
/// Impulse time weighting ~ 30 ms
Impulse,
/// A custom symmetric time weighting
CustomSymmetric {
/// Custom time constant [s]
t: Flt,
},
/// A custom symmetric time weighting
@ -19,14 +21,23 @@ pub enum TimeWeightingType {
tdown: Flt,
},
}
impl TimeWeightingType {
impl Default for TimeWeighting {
fn default() -> Self {
TimeWeighting::Fast
}
}
impl TimeWeighting {
/// get the analog poles of the single pole lowpass filter required for
/// getting the 'rectified' level (detector phase of SLM).
/// getting the 'rectified' level (detection phase of SLM).
pub fn getLowpassPoles(&self) -> (Flt, Option<Flt>) {
use TimeWeightingType::*;
use TimeWeighting::*;
match self {
Slow => (-1.0, None),
Fast => (-1. / 8., None),
Fast =>
// Time constant is 1/8 s, pole is at -8 rad/s
{
(-8., None)
}
Impulse => {
// For the impulse time weighting, some source says ~ 2.9 dB/s
// drop for the decay
@ -41,7 +52,7 @@ impl TimeWeightingType {
// with 10*log10(exp(-1.0/tau)) = -10/ln(10)/tau ≅ -4.34/tau
// dB/s where ln denotes the natural logarithm. So suppose we
// have 1.5 s, we indeed get a decay rate of 2.9 dB/s
(-35e-3, Some(-1.5))
(-1. / 35e-3, Some(-1. / 1.5))
}
CustomSymmetric { t } => {
assert!(*t > 0.);

2
tools/watchdoc.sh
View File

@ -6,6 +6,6 @@
# $ cargo install cargo-watch cargo-docserver`
# ```
#
cargo watch -s "clear && cargo rustdoc -p lasprs --lib && cargo docserve"
cargo watch -s "clear && cargo doc --no-deps -p lasprs --lib && cargo docserve"
# Then open: ${browser} http://localhost:4000