From 9e4ce617ae6e3b2eb814ddfdf8eb06bf96e2ffac Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F." Date: Sun, 27 Oct 2024 21:33:45 +0100 Subject: [PATCH] Added loose gain value to SeriesBiquad. Implemented pink noise. Bup to 0.6.5 --- Cargo.toml | 2 +- src/filter/mod.rs | 7 +- src/filter/seriesbiquad.rs | 29 ++++++--- src/siggen/mod.rs | 1 + src/siggen/noise.rs | 130 +++++++++++++++++++++++++++++++++++++ src/siggen/siggen.rs | 7 ++ src/siggen/source.rs | 33 ++++------ 7 files changed, 177 insertions(+), 32 deletions(-) create mode 100644 src/siggen/noise.rs diff --git a/Cargo.toml b/Cargo.toml index 4dd54da..26b5db8 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "lasprs" -version = "0.6.4" +version = "0.6.5" edition = "2021" authors = ["J.A. de Jong "] description = "Library for Acoustic Signal Processing (Rust edition, with optional Python bindings via pyo3)" diff --git a/src/filter/mod.rs b/src/filter/mod.rs index 83807e3..b1a004e 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -24,13 +24,14 @@ pub use dummy::DummyFilter; pub use seriesbiquad::SeriesBiquad; pub use zpkmodel::{PoleOrZero, ZPKModel, FilterSpec}; -/// Implementations of this trait are able to DSP-filter input data. +/// Implementations of this trait are able to DSP-filter input data. The filter +/// trait is implemented by, for example, [Biquad], [SeriesBiquad], and +/// [BiquadBank]. pub trait Filter: Send + Debug { - //! The filter trait is implemented by, for example, [Biquad], [SeriesBiquad], and [BiquadBank]. - /// Filter input to generate output. A vector of output floats is generated with the same /// length as input. fn filter(&mut self, input: &[Flt]) -> Vd; + /// Reset the filter state(s). In essence, this makes sure that all memory of the past is /// forgotten. fn reset(&mut self); diff --git a/src/filter/seriesbiquad.rs b/src/filter/seriesbiquad.rs index fce0149..f4ecc03 100644 --- a/src/filter/seriesbiquad.rs +++ b/src/filter/seriesbiquad.rs @@ -13,6 +13,8 @@ use anyhow::{bail, Result}; /// See (tests) /// pub struct SeriesBiquad { + // A loose gain value + gain: Flt, biqs: Vec, } @@ -23,7 +25,7 @@ impl SeriesBiquad { fn __repr__(&self) -> String { format!("{self:?}") } - + /// Create new series filter set. See [SeriesBiquad::new()] /// #[new] @@ -36,7 +38,7 @@ impl SeriesBiquad { pub fn unit_py() -> SeriesBiquad { SeriesBiquad::unit() } - + /// See: [SeriesBiquad::filter] #[pyo3(name = "filter")] pub fn filter_py<'py>( @@ -59,7 +61,12 @@ impl SeriesBiquad { /// pub fn newFromBiqs(biqs: Vec) -> SeriesBiquad { assert!(biqs.len() > 0); - SeriesBiquad { biqs } + SeriesBiquad { biqs, gain: 1.0 } + } + + /// Set an overall gain value to amplify / attenuate overall output + pub fn setGain(&mut self, g: Flt) { + self.gain = g; } /// Return reference to internally stored biquads @@ -94,7 +101,7 @@ impl SeriesBiquad { bail!("No filter coefficients given!"); } - Ok(SeriesBiquad { biqs }) + Ok(SeriesBiquad { biqs, gain: 1.0 }) } /// Unit impulse response series biquad. Input = output @@ -107,6 +114,14 @@ impl SeriesBiquad { fn clone_dyn(&self) -> Box { Box::new(self.clone()) } + + /// Filter input signal in place + pub fn filter_inout(&mut self, inout: &mut [Flt]) { + for biq in self.biqs.iter_mut() { + biq.filter_inout(inout); + } + inout.iter_mut().for_each(|io| *io *= self.gain); + } } impl Filter for SeriesBiquad { @@ -115,9 +130,7 @@ impl Filter for SeriesBiquad { //! fn filter(&mut self, input: &[Flt]) -> Vd { let mut inout = input.to_vec(); - for biq in self.biqs.iter_mut() { - biq.filter_inout(&mut inout); - } + self.filter_inout(&mut inout); inout } fn reset(&mut self) { @@ -130,7 +143,7 @@ impl Filter for SeriesBiquad { impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for SeriesBiquad { fn tf(&self, fs: Flt, freq: T) -> Ccol { let freq = freq.into(); - let mut res = self.biqs.first().unwrap().tf(fs, freq); + let mut res = self.biqs.first().unwrap().tf(fs, freq) * self.gain; for biq in self.biqs.iter().skip(1) { res = &res * biq.tf(fs, freq); } diff --git a/src/siggen/mod.rs b/src/siggen/mod.rs index d144d1a..4956212 100644 --- a/src/siggen/mod.rs +++ b/src/siggen/mod.rs @@ -21,6 +21,7 @@ //! //! ``` mod siggen; +mod noise; mod siggenchannel; mod source; mod siggencmd; diff --git a/src/siggen/noise.rs b/src/siggen/noise.rs new file mode 100644 index 0000000..2349b84 --- /dev/null +++ b/src/siggen/noise.rs @@ -0,0 +1,130 @@ +//! (Colored) noise Source implementations. Uses bilinear transform for a fixe and +//! pre-calculated set of analogue poles and zeros. +//! +use super::source::SourceImpl; +use crate::{config::*, PoleOrZero, SeriesBiquad, TransferFunction, ZPKModel}; +use rand::{rngs::SmallRng, Rng, SeedableRng}; +use rand_distr::StandardNormal; + +/// The number of poles in the filter that turns white noise into pink noise. +const PINKNOISE_ANALOG_ORDER: usize = 10; + +/// White noise source. Can be colored by applying a color filter to the source +#[derive(Clone, Debug)] +pub struct WhiteNoise { + // SmallRng is a cheap random number generator + rng: SmallRng, +} +impl WhiteNoise { + pub fn new() -> Self { + WhiteNoise { + rng: SmallRng::from_entropy(), + } + } +} +impl SourceImpl for WhiteNoise { + fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator) { + sig.for_each(|s| { + *s = self.rng.sample(StandardNormal); + }); + } + fn reset(&mut self, _fs: Flt) {} + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } +} + +#[derive(Debug, Clone)] +pub struct ColoredNoise { + // White noise generator + wn: WhiteNoise, + tmp: Vec, + analogue_blueprint: ZPKModel, + filter: SeriesBiquad, +} +impl ColoredNoise { + /// Generate a colored noise signal source that outputs pink noise (-3 dB / + /// octave ) from 20 Hz to 20 kHz. + pub fn newPinkNoise() -> Self { + let twopi = 2. * pi; + let fl = 10.; + let fu = 20e3; + // Exponentially spread poles + let fpoles_real: Vec = (0..PINKNOISE_ANALOG_ORDER) + .map(|i| fl * (fu / fl).powf((i as Flt) / (PINKNOISE_ANALOG_ORDER as Flt - 1.))) + .collect(); + + // Zeros as geometric means inbetween + let fzeros_real: Vec = (0..PINKNOISE_ANALOG_ORDER - 1) + .map(|i| (fpoles_real[i] * fpoles_real[i + 1]).sqrt()) + .collect(); + + let poles: Vec = fpoles_real + .into_iter() + .map(|f| PoleOrZero::Real1(-twopi * f)) + .collect(); + let zeros: Vec = fzeros_real + .into_iter() + .map(|f| PoleOrZero::Real1(-twopi * f)) + .collect(); + let gain = 1.; + + let analogue_blueprint = ZPKModel::new(zeros, poles, gain); + let filter = analogue_blueprint.bilinear(480000.); + Self { + wn: WhiteNoise::new(), + tmp: vec![], + analogue_blueprint, + filter, + } + } +} + +impl SourceImpl for ColoredNoise { + fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator) { + // Make temporary array of same size as output iterator + self.tmp.resize(sig.len(), 0.); + + // Generate white noise signal in place + self.wn.genSignal_unscaled(&mut self.tmp.iter_mut()); + + // Then, filter it with noise shaping filter + self.filter.filter_inout(&mut self.tmp); + + // And copy over to output iterator + sig.zip(&self.tmp).for_each(|(sig, src)| { + *sig = *src; + }); + } + fn reset(&mut self, fs: Flt) { + // Pre-computed analog poles + self.filter = self.analogue_blueprint.bilinear(fs); + let fnyq = fs / 2.; + + // We set the gain such that the signal power is the same as that for + // the white noise that enters it. We do this by calculating the sum of + // the power in each bin up to the Nyquist frequency, and calculating a + // gain that is the inverse of this value. + + // Choose a number of points. As long as its fine enough over the whole + // range (no peaks / dips are missing) At 48 kHz this means we sample + // the tf with a resolution of 10 Hz. Fine enough for slowly varying + // filter gains. + let N = 2400; + let freq = Array1::linspace(20., fnyq, 2000); + let tf_power_sum = self + .filter + .tf(fs, freq.view()) + .mapv(|t| t.abs().powi(2)) + .sum(); + // The above is sum of squares. The white noise signal is unit + // variance, meaning also unit linear level. So to compute gain, we + // have to take the square root of this ratio. + let gain = ((N as Flt) / tf_power_sum).sqrt(); + // dbg!(gain); + self.filter.setGain(gain); + } + fn clone_dyn(&self) -> Box { + Box::new(self.clone()) + } +} diff --git a/src/siggen/siggen.rs b/src/siggen/siggen.rs index ee1e8d8..e028332 100644 --- a/src/siggen/siggen.rs +++ b/src/siggen/siggen.rs @@ -21,6 +21,9 @@ use std::slice::IterMut; #[derive(Clone, Debug)] #[cfg_attr(feature = "python-bindings", pyclass)] pub struct Siggen { + // Sampling frequency in Hz + fs: Option, + // The source dynamic signal. Noise, a sine wave, sweep, etc source: Source, @@ -71,6 +74,7 @@ impl Siggen { /// - `source` - Source function pub fn new(nchannels: usize, source: Source) -> Siggen { Siggen { + fs: None, source, channels: vec![SiggenChannelConfig::new(); nchannels], source_buf: vec![], @@ -141,6 +145,8 @@ impl Siggen { match msg { SiggenCommand::ChangeSource { src } => { self.source = src; + + self.reset(self.fs.unwrap_or(48e3)); Ok(()) } SiggenCommand::ResetSiggen { fs } => { @@ -240,6 +246,7 @@ impl Siggen { /// * fs: (New) Sampling frequency \[Hz\] /// pub fn reset(&mut self, fs: Flt) { + self.fs = Some(fs); self.source.reset(fs); self.channels.iter_mut().for_each(|x| x.reset(fs)) } diff --git a/src/siggen/source.rs b/src/siggen/source.rs index dde9251..0bce2a3 100644 --- a/src/siggen/source.rs +++ b/src/siggen/source.rs @@ -1,5 +1,6 @@ //! All sources for a signal generator. Sine waves, sweeps, noise, etc. use super::sweep::{SweepParams, SweepType}; +use super::noise::{ColoredNoise, WhiteNoise}; use crate::config::*; use std::fmt::Debug; use std::ops::{Deref, DerefMut}; @@ -9,8 +10,6 @@ const twopi: Flt = 2.0 * pi; use crate::config::*; use anyhow::{bail, Result}; use rand::prelude::*; -use rand::rngs::ThreadRng; -use rand_distr::StandardNormal; /// Signal source for a signal generator. A signal source is capable of creating /// new signal data. @@ -42,7 +41,13 @@ impl Source { /// Create a white noise signal source pub fn newWhiteNoise() -> Source { Source { - src: Box::new(WhiteNoise { rng: SmallRng::from_entropy()}), + src: Box::new(WhiteNoise::new()), + } + } + /// Create a pink noise signal source + pub fn newPinkNoise() -> Source { + Source { + src: Box::new(ColoredNoise::newPinkNoise()), } } @@ -89,6 +94,11 @@ impl Source { Self::newWhiteNoise() } #[staticmethod] + #[pyo3(name = "newPinkNoise")] + fn newPinkNoise_py() -> Source { + Self::newPinkNoise() + } + #[staticmethod] #[pyo3(name = "newSweep")] fn newSweep_py( fs: Flt, @@ -119,23 +129,6 @@ impl SourceImpl for Silence { Box::new(self.clone()) } } -/// White noise source. Can be colored by applying a color filter to the source -#[derive(Clone, Debug)] -struct WhiteNoise { - // SmallRng is a cheap random number generator - rng: SmallRng -} -impl SourceImpl for WhiteNoise { - fn genSignal_unscaled(&mut self, sig: &mut dyn ExactSizeIterator) { - sig.for_each(|s| { - *s = self.rng.sample(StandardNormal); - }); - } - fn reset(&mut self, _fs: Flt) {} - fn clone_dyn(&self) -> Box { - Box::new(self.clone()) - } -} /// Sine wave, with configurable frequency #[derive(Clone, Debug)]