diff --git a/Cargo.toml b/Cargo.toml index b36a913..264dd68 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -16,7 +16,7 @@ crate-type = ["cdylib", "rlib",] [dependencies] # Error handling -anyhow = "1.0.75" +anyhow = "1.0.86" # Numerics # Optional future feature for ndarray: blas @@ -90,6 +90,12 @@ parking_lot = "0.12.3" # Builder code derive_builder = "0.20.0" +# Stack-allocated vectors +smallvec = "1.13.2" + +# Compile time constant floating point operations +softfloat = "1.0.0" + [dev-dependencies] ndarray-rand = "0.14.0" diff --git a/src/config.rs b/src/config.rs index 749bd22..fee5277 100644 --- a/src/config.rs +++ b/src/config.rs @@ -36,7 +36,7 @@ if #[cfg(feature = "python-bindings")] { } } pub use ndarray::prelude::*; -pub use ndarray::{Array1, Array2, ArrayView1}; +pub use ndarray::{Array1, Array2, ArrayView1, ArrayViewMut1}; pub use ndarray::Zip; use num::complex::Complex; diff --git a/src/filter/biquad.rs b/src/filter/biquad.rs index c7816de..de0b5d6 100644 --- a/src/filter/biquad.rs +++ b/src/filter/biquad.rs @@ -30,7 +30,11 @@ use num::Complex; /// And the frequency response can be found by filling in in above equation z = /// exp(i*omega/fs), where fs is the sampling frequency and omega is the radian /// frequency at which the transfer function is evaluated. +/// +/// ## Implementation details /// +/// The implementaion is so-called "Direct-form 2", see +/// [https://en.wikipedia.org/wiki/Digital_biquad_filter]. pub struct Biquad { // State parameters w1: Flt, @@ -124,6 +128,29 @@ impl Biquad { } } + /// Re-initialize state. *This is an advanced function. You should know what + /// you are doing!*. If not, please use any other function like + /// [Biquad::reset]. + pub fn setNextOutputX0(&mut self, out: Flt) { + let (b0, b1, b2, a1, a2) = (self.b0, self.b1, self.b2, self.a1, self.a2); + let w = out / (b1 + b2 - b0 * (a1 + a2)); + self.w1 = w; + self.w2 = w; + } + + /// Change the gain value such that it matches `val` at frequency `freq`. + /// Does not change the phase at the given frequency. + pub fn setGainAt(mut self, freq: Flt, required_gain: Flt) -> Biquad { + assert!(required_gain > 0.); + let freq = [freq]; + let cur_gain_at_freq = self.tf(-1.0, &freq)[0].abs(); + let gain_fac = required_gain / cur_gain_at_freq; + self.b0 *= gain_fac; + self.b1 *= gain_fac; + self.b2 *= gain_fac; + self + } + /// Construct a Biquad with 0 initial state from coefficients given as /// arguments. /// @@ -202,17 +229,22 @@ impl Biquad { Ok(Biquad::fromCoefs(b0, b1, 0., a1, 0.)) } + #[inline] + /// Filter single sample, outputs by overwriting input sample. + pub fn filter_inout_single(&mut self, sample: &mut Flt) { + let w0 = *sample - self.a1 * self.w1 - self.a2 * self.w2; + let yn = self.b0 * w0 + self.b1 * self.w1 + self.b2 * self.w2; + self.w2 = self.w1; + self.w1 = w0; + *sample = yn; + } + /// Filter input signal, output by overwriting input slice. #[inline] pub fn filter_inout(&mut self, inout: &mut [Flt]) { for sample in inout.iter_mut() { - let w0 = *sample - self.a1 * self.w1 - self.a2 * self.w2; - let yn = self.b0 * w0 + self.b1 * self.w1 + self.b2 * self.w2; - self.w2 = self.w1; - self.w1 = w0; - *sample = yn; + self.filter_inout_single(sample); } - // println!("{:?}", inout); } /// Create new biquad using bilinear transform. Optionally pre-warps the @@ -470,4 +502,21 @@ mod test { // let freq = &[0., 10.,100.,1000., 2000.]; // println!("{:?}", b3.tf(fs, freq)); } + + #[test] + fn test_setOutput1() { + let mut f = Biquad::firstOrderHighPass(10., 1.).unwrap(); + f.setNextOutputX0(1.0); + let mut sample = 0.; + f.filter_inout_single(&mut sample); + assert_abs_diff_eq!(sample, 1.0); + } + #[test] + fn test_setOutput2() { + let mut f = Biquad::bilinear_zpk(1.0, None, Some(PoleOrZero::Real1(-1.)), None, None); + f.setNextOutputX0(4.2); + let mut sample = 0.; + f.filter_inout_single(&mut sample); + assert_abs_diff_eq!(sample, 4.2, epsilon = 1e-6); + } } diff --git a/src/filter/mod.rs b/src/filter/mod.rs index 6ca5c7e..f41b2cc 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -12,10 +12,12 @@ mod dummy; mod seriesbiquad; mod zpkmodel; mod butter; +mod octave; pub use super::ps::FreqWeightingType; pub use biquad::Biquad; pub use biquadbank::BiquadBank; +pub use octave::{StandardFilterDescriptor, G, FREQ_REF}; pub use dummy::DummyFilter; pub use seriesbiquad::SeriesBiquad; pub use zpkmodel::{PoleOrZero, ZPKModel, FilterSpec}; diff --git a/src/filter/octave.rs b/src/filter/octave.rs new file mode 100644 index 0000000..4f1f781 --- /dev/null +++ b/src/filter/octave.rs @@ -0,0 +1,542 @@ +use crate::{Flt, ZPKModel}; +use anyhow::{anyhow, bail, Result}; +use clap::Error; +use num::{traits::float, Float}; +use rayon::iter::Filter; +use softfloat::F64; +use std::{borrow::Cow, cmp::Ordering}; + +/// Names of standard octave filters +const OCTAVE_NOMINAL_MIDBAND_NAMES: [&str; 12] = [ + "8", "16", "31.5", "63", "125", "250", "500", "1k", "2k", "4k", "8k", "16k", +]; +const OCTAVE_NAMES_OFFSET: i32 = 7; + +const MIN_MIDBAND_FREQ: Flt = 8.; +const MAX_MIDBAND_FREQ: Flt = 20e3; + +const THIRDOCTAVE_NOMINAL_MIDBAND_NAMES: [&'static str; 33] = [ + "12.5", "16", "20", "25", "31.5", "40", "50", "63", "80", "100", "125", "160", "200", "250", + "315", "400", "500", "630", "800", "1k", "1.25k", "1.6k", "2k", "2.5k", "3.15k", "4k", "5k", + "6.3k", "8k", "10k", "12.5k", "16k", "20k", +]; +const THIRDOCTAVE_NAMES_OFFSET: i32 = 19; + +/// Return the num x-value for a certain 'name', like '16', or '1k' +fn nominal_octave_designator(name: &str) -> Result { + debug_assert!(OCTAVE_NOMINAL_MIDBAND_NAMES[OCTAVE_NAMES_OFFSET as usize] == "1k"); + Ok(OCTAVE_NOMINAL_MIDBAND_NAMES + .iter() + .position(|i| *i == name) + .ok_or(anyhow!( + "Cannot find name in list of OCTAVE_NOMINAL_MIDBAND_NAMES" + ))? as i32 + - OCTAVE_NAMES_OFFSET) +} + +fn nominal_thirdoctave_designator(name: &str) -> Result { + debug_assert!(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES[THIRDOCTAVE_NAMES_OFFSET as usize] == "1k"); + Ok(THIRDOCTAVE_NOMINAL_MIDBAND_NAMES + .iter() + .position(|i| *i == name) + .ok_or(anyhow!( + "Cannot find name in list of THIRDOCTAVE_NOMINAL_MIDBAND_NAMES" + ))? as i32 + - THIRDOCTAVE_NAMES_OFFSET) +} +// Raise a^b. In const-mode. +const fn powf(a: Flt, b: Flt) -> Flt { + let a = a as f64; + let b = b as f64; + let a = softfloat::F64::from_native_f64(a); + let b = softfloat::F64::from_native_f64(b); + softfloat::F64::exp(b.mul(a.ln())).to_native_f64() as Flt +} + +/// Octave ratio. We use G_10, which is 10^(3/10) ≅ 1.995 +pub const G: Flt = powf(10., 0.3); +/// Reference freuqency, 1kHz +pub const FREQ_REF: Flt = 1000.; + +/// Standard filter descriptor. Used to generate bandpass filters that are +/// compliant with IEC 61260 (1995). +/// +/// # Examples +/// +/// ## Create a 16 Hz octave band digital filter running at 48kHz. +/// +/// ```rust +/// use lasprs::filter::*; +/// fn main() -> anyhow::Result<()> { +/// let desc = StandardFilterDescriptor::Octave("16")?; +/// let filter = desc.genFilter().bilinear(48e3); +/// Ok(()) +/// } +/// ``` +/// +/// ## Create a one-third octave band bandpass filter that has the frequency of 42 in its pass-band +/// +/// ```rust +/// +/// use lasprs::filter::*; +/// fn main() -> anyhow::Result<()> { +/// let desc = StandardFilterDescriptor::filterForFreq(3, 42.)?; +/// let filter = desc.genFilter().bilinear(48e3); +/// Ok(()) +/// } +/// ``` +#[derive(PartialEq, Clone, Debug)] +pub struct StandardFilterDescriptor { + /// b and x. Bandwidth and offset w.r.t. reference frequency. + /// + /// Band width fraction of an octave. 1 means full octave of bandwidth. 3 + /// means 1/3th octave, 6 means 1/6th, and so on. + /// + /// If bx is None, it means we do not filter at all (an overall channel) + b: u32, + x: i32, +} + +impl StandardFilterDescriptor { + /// Create analog filter specification from descriptor + pub fn genFilter(&self) -> ZPKModel { + let order = 5; + if let Some((fl, fu)) = self.fl_fh() { + ZPKModel::butter(crate::FilterSpec::Bandpass { fl, fu, order }) + } else { + ZPKModel::default() + } + } + // Check whether a certain midband frequency of created + // StandardFilterDescriptor is in the allowed range. This is a helper + // function that is used to check wheter created StandardFilterDescriptors + // are valid. + fn check_fmid_in_range(&self) -> Result<()> { + if let Some(fm) = self.fm() { + if fm < MIN_MIDBAND_FREQ { + bail!( + "Invalid x. Computed filter center frequency is {} Hz, which is too low. Lowest allowed is {} Hz", + fm, MIN_MIDBAND_FREQ + ) + } else if fm > 20e3 { + bail!( + "Invalid x. Computed filter center frequency is {} Hz, which is too high. Highest allowed is {} Hz", + fm, MAX_MIDBAND_FREQ + ) + } + } + + Ok(()) + } + /// Create new standard filter descriptor `b` from given relative bandwidth + /// and band designator `x`. If not sure what `x` and `b` are, see + /// documentation on [StandardFilterDescriptor::genFilterSetByDesignator]. + pub fn build(b: u32, x: i32) -> Result { + let desc = StandardFilterDescriptor { b, x }; + desc.check_fmid_in_range()?; + match b { + 0 => Ok(desc), + 1 => Ok(desc), + 3 => Ok(desc), + 6 => Ok(desc), + 12 => Ok(desc), + 24 => Ok(desc), + _ => bail!( + "Bandwidth {} is invalid. Please choose a value from 0, 1, 3, 6, 12 or 24", + b + ), + } + } + /// Generate filter descriptor. Practically applies no filtering at all. + pub fn Overall() -> Result { + Ok(StandardFilterDescriptor { b: 0, x: 0 }) + } + /// Generate filter descriptor for octave band. + /// + /// # Args + /// + /// - `band_descr` - band designator. Can be '1k', or 0. + pub fn Octave(band_descr: T) -> Result + where + T: TryInto, + { + let x = band_descr.try_into()?.x; + Ok(StandardFilterDescriptor { b: 1, x }) + } + /// Generate filter descriptor for one-third octave band. + /// + /// # Args + /// + /// - `x` - band designator + pub fn ThirdOctave(band_descr: T) -> Result + where + T: TryInto, + { + let x = band_descr.try_into()?.x; + Ok(StandardFilterDescriptor { b: 3, x }) + } + + /// Searches for a filter with `1/b` relative bandwidth w.r.t one octave + /// that has frequency `f` in its pass-band. + /// + pub fn filterForFreq(b: u32, f: Flt) -> Result { + if f < MIN_MIDBAND_FREQ || f > MAX_MIDBAND_FREQ { + bail!("Invalid frequency. Please use search frequency between 8 Hz and 20 kHz") + } + match b { + 0 => Self::Overall(), + 1 | 3 | 6 | 12 | 24 => { + let mut desc = StandardFilterDescriptor { b, x: 0 }; + + let f_in_range = |desc: &StandardFilterDescriptor| { + let (fl, fh) = desc.fl_fh().unwrap(); + // println!("fl: {fl}, fh: {fh}"); + if f < fl { + Ordering::Less + } else if f > fh { + Ordering::Greater + } else { + Ordering::Equal + } + }; + + loop { + let fm = desc.fm().unwrap(); + // eprintln!("Fmid: {fm:.2e}"); + // eprintln!("desc: {desc:#?}"); + let ord = f_in_range(&desc); + // Bands for midband frequencies are a bit wider here + if fm < MIN_MIDBAND_FREQ - 3. || fm > MAX_MIDBAND_FREQ * 1.1 { + bail!("Frequency not in range"); + } + match ord { + Ordering::Equal => break, + Ordering::Less => desc = StandardFilterDescriptor { b, x: desc.x - 1 }, + Ordering::Greater => desc = StandardFilterDescriptor { b, x: desc.x + 1 }, + } + } + Ok(desc) + } + _ => Self::build(b, 0), + } + } + + /// Creates a set of octave filters. + pub fn genOctaveFilterSet(low_f: Option, high_f: Option) -> Result> + where + T: TryInto, + { + 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 + }; + Ok((xmin..=xmax).map(|x| Self::Octave(x).unwrap()).collect()) + } + + /// Creates a set of one-third octave bandpass filters. + pub fn genThirdOctaveFilterSet(low_f: Option, high_f: Option) -> Result> + where + T: TryInto, + { + 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 + }; + Ok((xmin..=xmax) + .map(|x| Self::ThirdOctave(x).unwrap()) + .collect()) + } + + /// Generate a filter set using designators + /// + /// # Args + /// + /// - `b` - Inverse of the relative bandwidth w.r.t. one octave. `b=0` means + /// overall, `b=1` is one octave, `b=3`` is one-third, etc. + /// - `xmin` - Band designator of lowest band. Midband frequency can be computed as [FREQ_REF]*[G]^(`xmin/b`) + /// - `xmax` - Band designator of lowest band. Midband frequency can be computed as [FREQ_REF]*[G]^(`xmax/b`) + /// - `include_overall` - If `true`, adds an overall filter (a no-op) as the last designator in the list + pub fn genFilterSetByDesignator( + b: u32, + xmin: i32, + xmax: i32, + include_overall: bool, + ) -> Result> { + if xmin > xmax { + bail!("xmin should be <= xmax"); + } + let cap = (xmax - xmin) as usize + if include_overall { 1 } else { 0 }; + let mut res = Vec::with_capacity(cap); + + for x in xmin..=xmax { + res.push(StandardFilterDescriptor::build(b, x)?); + } + + if include_overall { + res.push(StandardFilterDescriptor::Overall()?) + } + + Ok(res) + } + + /// Convenience function for creating a filter bank. Creates a set of + /// standard filters with relative bandwidth `b`, that has `fl` in the + /// lowest bandpass filter and `fu` in the highest. + /// + /// # Other args + /// + /// - `include_overall` - If `true`, adds an overall filter (a no-op) as the + pub fn genFilterSetInRange( + b: u32, + fl: Flt, + fu: Flt, + include_overall: bool, + ) -> Result> { + let xmin = StandardFilterDescriptor::filterForFreq(b, fl)?.x; + let xmax = StandardFilterDescriptor::filterForFreq(b, fu)?.x; + StandardFilterDescriptor::genFilterSetByDesignator(b, xmin, xmax, include_overall) + } + + /// Returns the midband frequency in \[Hz\] + pub fn fm(&self) -> Option { + if self.b == 0 { + None + } else { + let b = self.b as Flt; + let x = self.x as Flt; + Some(FREQ_REF * G.powf(x / b)) + } + } + + /// Cuton frequency and cut-off frequency, in \[Hz\]. + /// Returns none if it does not apply, for [FilterDescriptor::Overall]. + pub fn fl_fh(&self) -> Option<(Flt, Flt)> { + match self.b { + 0 => None, + b => { + let fm = self.fm().unwrap(); + let b = b as Flt; + let fl = fm * G.powf(-1. / (2. * b)); + let fu = fm * G.powf(1. / (2. * b)); + Some((fl, fu)) + } + } + } + /// Give a common name to filter, specifically the filters are named after + /// the midband frequency. + pub fn name(&self) -> Cow<'static, str> { + let x = self.x; + match self.b { + 0 => Cow::Borrowed("Overall"), + 1 => OctaveBandDescriptor { x }.name(), + 3 => ThirdOctaveBandDescriptor { x }.name(), + 6 => { + if x % 2 == 0 { + ThirdOctaveBandDescriptor { x: x / 2 }.name() + } else { + Default::default() + } + } + 12 => { + if x % 2 == 0 { + StandardFilterDescriptor { + b: self.b / 2, + x: self.x / 2, + } + .name() + } else { + Default::default() + } + } + _ => unreachable!(), + } + } +} + +/// A valid descriptor for a standard one-third octave band +pub struct ThirdOctaveBandDescriptor { + x: i32, +} +/// A valid descriptor for a standard octave band +pub struct OctaveBandDescriptor { + x: i32, +} + +impl TryFrom for OctaveBandDescriptor { + type Error = anyhow::Error; + fn try_from(x: i32) -> Result { + if x + OCTAVE_NAMES_OFFSET < 0 + || x + OCTAVE_NAMES_OFFSET >= OCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 + { + bail!( + "Invalid filter designator x. Should be >= -{OCTAVE_NAMES_OFFSET} and < {}", + OCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 - OCTAVE_NAMES_OFFSET + ); + } + Ok(Self { x }) + } +} +impl TryFrom<&str> for OctaveBandDescriptor { + type Error = anyhow::Error; + fn try_from(name: &str) -> Result { + Ok(Self { + x: nominal_octave_designator(name)?, + }) + } +} +impl TryFrom for ThirdOctaveBandDescriptor { + type Error = anyhow::Error; + fn try_from(x: i32) -> Result { + if x + THIRDOCTAVE_NAMES_OFFSET < 0 + || x + THIRDOCTAVE_NAMES_OFFSET >= THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 + { + bail!( + "Invalid filter designator x. Should be >= -{THIRDOCTAVE_NAMES_OFFSET} and < {}", + THIRDOCTAVE_NOMINAL_MIDBAND_NAMES.len() as i32 - THIRDOCTAVE_NAMES_OFFSET + ); + } + Ok(Self { x }) + } +} +impl TryFrom<&str> for ThirdOctaveBandDescriptor { + type Error = anyhow::Error; + fn try_from(name: &str) -> Result { + Ok(Self { + x: nominal_thirdoctave_designator(name)?, + }) + } +} + +trait BandDescriptor { + fn name(&self) -> Cow<'static, str>; +} +impl BandDescriptor for OctaveBandDescriptor { + fn name(&self) -> Cow<'static, str> { + Cow::Borrowed( + OCTAVE_NOMINAL_MIDBAND_NAMES + .get((self.x + OCTAVE_NAMES_OFFSET) as usize) + .map(|s| *s) + .unwrap_or_default(), + ) + } +} +impl BandDescriptor for ThirdOctaveBandDescriptor { + fn name(&self) -> Cow<'static, str> { + Cow::Borrowed( + THIRDOCTAVE_NOMINAL_MIDBAND_NAMES + .get((self.x + THIRDOCTAVE_NAMES_OFFSET) as usize) + .map(|s| *s) + .unwrap_or_default(), + ) + } +} + +// impl Into for &str { +// fn into(self) -> OctaveBandDescriptor { +// OctaveBandDescriptor{x: self} +// } +// } +#[cfg(test)] +mod test { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn test_finder() { + // assert_eq!( + // StandardFilterDescriptor::filterForFreq(0, 1000.).unwrap(), + // StandardFilterDescriptor::Overall().unwrap() + // ); + // assert_eq!( + // StandardFilterDescriptor::filterForFreq(1, 1e3).unwrap(), + // StandardFilterDescriptor::Octave(0).unwrap() + // ); + assert_eq!( + StandardFilterDescriptor::filterForFreq(1, 8.).unwrap(), + StandardFilterDescriptor::Octave(-OCTAVE_NAMES_OFFSET).unwrap() + ); + // assert_eq!( + // StandardFilterDescriptor::filterForFreq(3, 1000.).unwrap(), + // StandardFilterDescriptor::ThirdOctave(0).unwrap() + // ); + // assert_eq!( + // StandardFilterDescriptor::filterForFreq(3, 12.).unwrap(), + // StandardFilterDescriptor::ThirdOctave(-THIRDOCTAVE_NAMES_OFFSET).unwrap() + // ); + } + + #[test] + fn test_builders() { + assert_eq!( + StandardFilterDescriptor::Octave("8").unwrap(), + StandardFilterDescriptor::Octave(-OCTAVE_NAMES_OFFSET).unwrap() + ); + assert_eq!( + StandardFilterDescriptor::Octave("2k").unwrap(), + StandardFilterDescriptor::Octave(1).unwrap() + ); + assert_eq!( + StandardFilterDescriptor::ThirdOctave("12.5").unwrap(), + StandardFilterDescriptor::ThirdOctave(-THIRDOCTAVE_NAMES_OFFSET).unwrap() + ); + assert_eq!( + StandardFilterDescriptor::ThirdOctave("2k").unwrap(), + StandardFilterDescriptor::ThirdOctave(3).unwrap() + ); + } + #[test] + #[should_panic] + fn out_range_octave1() { + StandardFilterDescriptor::Octave("4").unwrap(); + } + #[test] + #[should_panic] + fn out_range_octave2() { + StandardFilterDescriptor::Octave("7").unwrap(); + } + + #[test] + fn test_name_and_approx() { + assert_eq!( + StandardFilterDescriptor::filterForFreq(1, 16e3).unwrap(), + StandardFilterDescriptor::Octave("16k").unwrap(), + ); + } + + #[test] + fn test_names() { + assert_eq!(StandardFilterDescriptor::Octave(1).unwrap().name(), "2k"); + assert_eq!( + StandardFilterDescriptor::ThirdOctave(1).unwrap().name(), + "1.25k" + ); + } + + #[test] + fn test_octave() { + assert_eq!(nominal_octave_designator("1k").unwrap(), 0); + assert_eq!(nominal_octave_designator("2k").unwrap(), 1); + } + #[test] + fn test_thirdoctave() { + assert_eq!(nominal_thirdoctave_designator("1k").unwrap(), 0); + assert_eq!(nominal_thirdoctave_designator("2k").unwrap(), 3); + } + + #[test] + fn test_G() { + assert_abs_diff_eq!(G, (10 as Flt).powf(3. / 10.)); + } +} diff --git a/src/filter/seriesbiquad.rs b/src/filter/seriesbiquad.rs index 10b6515..b899f22 100644 --- a/src/filter/seriesbiquad.rs +++ b/src/filter/seriesbiquad.rs @@ -62,6 +62,11 @@ impl SeriesBiquad { SeriesBiquad { biqs } } + /// Return reference to internally stored biquads + pub fn getBiquads(&self) -> &Vec { + &self.biqs + } + /// Create a new series biquad, having an arbitrary number of biquads. /// /// # Arguments diff --git a/src/filter/zpkmodel.rs b/src/filter/zpkmodel.rs index 29e29fc..0c146bc 100644 --- a/src/filter/zpkmodel.rs +++ b/src/filter/zpkmodel.rs @@ -176,9 +176,10 @@ impl ZPKModel { /// - `poles` - list like struct of poles. Can be a `Vec` or an /// `&[ZeroOrPole]`. /// - `k` - linear gain. - pub fn new(zeros: T, poles: T, k: Flt) -> ZPKModel + pub fn new(zeros: T, poles: U, k: Flt) -> ZPKModel where T: Into>, + U: Into>, { let z = zeros.into(); let p = poles.into(); diff --git a/src/lib.rs b/src/lib.rs index b7f34a3..9656be9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -41,6 +41,7 @@ pub mod ps; pub mod siggen; use filter::*; pub mod rt; +pub mod slm; /// A Python module implemented in Rust. #[cfg(feature = "python-bindings")] diff --git a/src/slm/mod.rs b/src/slm/mod.rs index a6583d6..fd2370d 100644 --- a/src/slm/mod.rs +++ b/src/slm/mod.rs @@ -1,15 +1,11 @@ //! Sound Level Meter (SLM) module. -//! +//! //! Provides structs and helpers (SLMBuilder) for creating configurated Sound //! Level Meters. -//! - - -/// Sound Level Meter -struct SLM { - -} - -impl SLM { - -} \ No newline at end of file +//! +mod settings; +mod tw; +mod slm; +pub use slm::SLM; +pub use settings::SLMSettings; +pub use tw::TimeWeightingType; \ No newline at end of file diff --git a/src/slm/settings.rs b/src/slm/settings.rs new file mode 100644 index 0000000..64aa6f3 --- /dev/null +++ b/src/slm/settings.rs @@ -0,0 +1,19 @@ +use derive_builder::Builder; +use smallvec::SmallVec; +use crate::{Flt, FreqWeightingType, filter::StandardFilterDescriptor}; +use super::TimeWeightingType; + + +#[derive(Builder, Clone)] +pub struct SLMSettings { + pub fs: Flt, + pub Lref: Flt, + pub freqWeighting: FreqWeightingType, + pub timeWeighting: TimeWeightingType, + pub filterDescriptors: SmallVec<[StandardFilterDescriptor; 64]>, + +} + +impl SLMSettings { + +} \ No newline at end of file diff --git a/src/slm/slm.rs b/src/slm/slm.rs new file mode 100644 index 0000000..6045053 --- /dev/null +++ b/src/slm/slm.rs @@ -0,0 +1,194 @@ +use derive_builder::Builder; +use itertools::Itertools; +use ndarray::ArrayView1; +use rayon::iter::{IntoParallelRefIterator, ParallelIterator}; +use rayon::prelude::*; +use smallvec::SmallVec; + +use super::settings::SLMSettings; +use crate::{config::*, filter::Filter}; +use crate::{Biquad, Dcol, Flt, FreqWeightingType, PoleOrZero, SeriesBiquad, ZPKModel}; +struct SLMChannel { + stat: SLMStat, + bp: SeriesBiquad, + rect_lowpass_up: Biquad, + rect_lowpass_down: Option, +} + +/// Sound Level Meter +pub struct SLM { + // Number of samples processed after last run() is called. + N: usize, + Lrefsq: Flt, + prefilter: SeriesBiquad, + channels: SmallVec<[SLMChannel; 64]>, +} + +impl SLM { + // Create simple first order lowpass filter with unit D.C. gain and given + // real pole. + fn lpfilter_from_pole(fs: Flt, p: PoleOrZero) -> Biquad { + Biquad::bilinear_zpk(fs, None, Some(p), Some(1.0), None).setGainAt(0., 1.) + } + /// Create new Sound Level Meter from given settings + pub fn new(settings: SLMSettings) -> Self { + let fs = settings.fs; + let prefilter = ZPKModel::freqWeightingFilter(settings.freqWeighting).bilinear(fs); + let channels = settings + .filterDescriptors + .iter() + .map(|descriptor| { + // Generate bandpass filter + let bp = descriptor.genFilter().bilinear(fs); + // Initalize statistics with defaults + let stat = SLMStat::default(); + + // Generate rectifier filter for upwards + let poles = settings.timeWeighting.getLowpassPoles(); + + let rect_lowpass_up = Self::lpfilter_from_pole(fs, PoleOrZero::Real1(poles.0)); + + let rect_lowpass_down = if let Some(p) = poles.1 { + Some(Self::lpfilter_from_pole(fs, PoleOrZero::Real1(p))) + } else { + None + }; + SLMChannel { + stat, + bp, + rect_lowpass_up, + rect_lowpass_down, + } + }) + .collect(); + SLM { + prefilter, + channels, + Lrefsq: settings.Lref.powi(2), + N: 0, + } + } + /// Push new time data through sound level meter. Returns L(t) data for each + /// channel. + /// + /// # Args + /// + /// - `td`: Time data + pub fn run(&mut self, td: &[Flt]) -> Option>> { + if td.len() == 0 { + return None; + } + let prefiltered = self.prefilter.filter(td); + let Lt_iter = self.channels.par_iter_mut().map(|ch| { + let mut tmp = ch.bp.filter(&prefiltered); + let mut N = self.N; + + // Filtered squared + let mut filtered_squared = { + let mut tmp_view = ArrayViewMut1::from(&mut tmp); + tmp_view.mapv_inplace(|a| a * a); + tmp_view + }; + + // Update Lpk, Leq + filtered_squared.for_each(|sample_pwr| { + let new_pk = sample_pwr.abs(); + if new_pk > ch.stat.Ppk { + ch.stat.Ppk = new_pk; + } + // Update equivalent level + ch.stat.Peq = (ch.stat.Peq * N as Flt + sample_pwr) / (N as Flt + 1.); + N += 1; + }); + + // Run filtered_squared signal throug rectifier + if let Some(rectifier_down) = &mut ch.rect_lowpass_down { + filtered_squared.mapv_inplace(|sample_sq| { + let mut fup = sample_sq; + let mut fdown = sample_sq; + + // Filter in up-filter + let rectifier_up = &mut ch.rect_lowpass_up; + rectifier_up.filter_inout_single(&mut fup); + // Filter in down-filter + rectifier_down.filter_inout_single(&mut fdown); + + // Check who wins + if fup > fdown { + rectifier_down.setNextOutputX0(fup); + fup + } else { + rectifier_up.setNextOutputX0(fdown); + fdown + } + }); + } else { + // Filter in place + let rectifier = &mut ch.rect_lowpass_up; + rectifier.filter_inout(filtered_squared.as_slice_mut().unwrap()); + } + + // Update max signal power gotten so far + let rectified = &mut filtered_squared; + rectified.for_each(|val| { + if *val > ch.stat.Pmax { + ch.stat.Pmax = *val; + } + }); + // Update last signal power coming from SLM + ch.stat.Pt_last = *filtered_squared.last().unwrap(); + tmp + }); + let Lt: Vec<_> = Lt_iter.collect(); + self.N += td.len(); + Some(Lt) + } + + /// Number of channels in SLM + pub fn nch(&self) -> usize { + self.channels.len() + } + + fn levels_from(&self, stat_returner: T) -> Dcol + where + T: Fn(&SLMChannel) -> Flt, + { + Dcol::from_iter( + self.channels + .iter() + .map(|ch| 20. * Flt::log10(stat_returner(ch) / self.Lrefsq)), + ) + } + + /// Get max levels for each channel + pub fn Lmax(&self) -> Dcol { + self.levels_from(|ch| ch.stat.Pmax) + } + /// Get peak levels for each channel + pub fn Lpk(&self) -> Dcol { + self.levels_from(|ch| ch.stat.Ppk) + } + + /// Get equivalent levels for each channel + pub fn Leq(&self) -> Dcol { + self.levels_from(|ch| ch.stat.Peq) + } + /// Get last value of level vs time + pub fn Ltlast(&self) -> Dcol { + self.levels_from(|ch| ch.stat.Pt_last) + } +} + +#[derive(Clone, Default)] +/// Quantities defined as powers, i.e. square of amplitude +struct SLMStat { + // Max signal power + Pmax: Flt, + // Peak signal power + Ppk: Flt, + // Equivalent signal power + Peq: Flt, + + // Last obtained signal power, after last time run() is called. + Pt_last: Flt, +} diff --git a/src/slm/tw.rs b/src/slm/tw.rs new file mode 100644 index 0000000..6540ae6 --- /dev/null +++ b/src/slm/tw.rs @@ -0,0 +1,57 @@ +use crate::Flt; +#[derive(Clone, Copy)] +pub enum TimeWeightingType { + /// Slow time weighting + Slow, + /// Fast time weighting + Fast, + /// Impulse time weighting + Impulse, + /// A custom symmetric time weighting + CustomSymmetric { + t: Flt, + }, + /// A custom symmetric time weighting + CustomAsymmetric { + /// Time weighting when level is increasing + tup: Flt, + /// Time weighting when level is decreasing + tdown: Flt, + }, +} +impl TimeWeightingType { + /// get the analog poles of the single pole lowpass filter required for + /// getting the 'rectified' level (detector phase of SLM). + pub fn getLowpassPoles(&self) -> (Flt, Option) { + use TimeWeightingType::*; + match self { + Slow => (-1.0, None), + Fast => (-1. / 8., None), + Impulse => { + // For the impulse time weighting, some source says ~ 2.9 dB/s + // drop for the decay + // [https://www.nti-audio.com/en/support/know-how/fast-slow-impulse-time-weighting-what-do-they-mean]. + // + // Other source + // [https://support.dewesoft.com/en/support/solutions/articles/14000139949-exponential-averaging-fast-f-slow-s-impulse-i-] + // say a time constant of 1.5 s. Are they compatible? + + // Compute decay rate in dB/s from the filter time constant. An + // initial value drops as exp(-t/tau). So in 1 s the level drops + // 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)) + } + CustomSymmetric { t } => { + assert!(*t > 0.); + (-*t, None) + } + CustomAsymmetric { tup, tdown } => { + assert!(*tup > 0.); + assert!(*tdown > 0.); + (-*tup, Some(-*tdown)) + } + } + } +}