From 7e9cf734d0f3453d78751f00beb09721535ba97f Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Wed, 3 Jul 2024 20:01:12 +0200 Subject: [PATCH] Biguad, biquadbank, dummy filters, etc. A bit on all topics --- src/filter/biquad.rs | 76 +++++++++++------ src/filter/biquadbank.rs | 4 +- src/filter/dummy.rs | 24 ++++++ src/filter/mod.rs | 8 +- src/filter/seriesbiquad.rs | 4 +- src/ps/aps.rs | 6 +- src/ps/ps.rs | 31 +++++-- src/ps/window.rs | 167 ++++++++++++++++--------------------- src/timebuffer.rs | 8 ++ 9 files changed, 188 insertions(+), 140 deletions(-) create mode 100644 src/filter/dummy.rs diff --git a/src/filter/biquad.rs b/src/filter/biquad.rs index 00ce934..ee550fe 100644 --- a/src/filter/biquad.rs +++ b/src/filter/biquad.rs @@ -10,7 +10,7 @@ use num::Complex; /// This implementation only allows for normalized coefficients (a_0 = 1). It /// performs the following relation of output to input: /// -/// ``` +/// ```math /// y[n] = - a_1 * y[n-1] - a_2 * y[n-2] /// + b_0 * x[n] + b_1 * x[n-1] + b_2 * x[n-2] /// ``` @@ -21,7 +21,7 @@ use num::Complex; /// /// The transfer function is: /// -/// ``` +/// ```math /// b_0 + b_1 z^-1 + b_2 * z^-2 /// H[z] = ----------------------------- /// 1 + a_1 z^-1 + a_2 * z^-2 @@ -62,11 +62,32 @@ impl Biquad { #[pyo3(name = "firstOrderHighPass")] #[staticmethod] /// See: [Biquad::firstOrderHighPass()] - pub fn firstOrderHighPass_py(fs: Flt, cuton_Hz: Flt) -> PyResult { - Ok(Biquad::firstOrderHighPass(fs, cuton_Hz)?) + pub fn firstOrderHighPass_py(fs: Flt, fc: Flt) -> PyResult { + Ok(Biquad::firstOrderHighPass(fs, fc)?) } - #[pyo3(name = "filter")] + + /// See: [Biquad::firstOrderLowPass()] + #[pyo3(name = "firstOrderLowPass")] + #[staticmethod] + pub fn firstOrderLowPass_py(fs: Flt, fc: Flt) -> PyResult { + Ok(Biquad::firstOrderLowPass(fs, fc)?) + } + + /// See: [Biquad::tf()] + #[pyo3(name = "tf")] + pub fn tf_py<'py>( + &self, + py: Python<'py>, + fs: Flt, + freq: PyArrayLike1, + ) -> PyResult> { + let freq = freq.as_array(); + let res = PyArray1::from_array_bound(py, &self.tf(fs, freq)); + Ok(res) + } + /// See: [Biquad::filter()] + #[pyo3(name = "filter")] pub fn filter_py<'py>( &mut self, py: Python<'py>, @@ -151,22 +172,21 @@ impl Biquad { /// First order low pass filter (one pole in the real axis). No pre-warping /// correction done. + /// + /// * `fs` - Sampling frequency [Hz] + /// * `fc` - Cut-off frequency (-3 dB point) [Hz] pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result { - match fc { - x if fc <= 0. => bail!("Cuton frequency should be > 0"), - x if fc >= fs / 2. => bail!("Cuton frequency should be smaller than Nyquist frequency"), - _ => (), + if fc <= 0. { + bail!("Cuton frequency, given: should be > 0") } - let w0: Flt = 2. * pi * fc / fs; - let cw = Flt::cos(w0); - let b0: Flt = 2. * pi * fc * (cw + 1.) / (2. * pi * fc * cw + 2. * pi * fc - cw + 1.); - let b1: Flt = 2. * pi * fc * (cw + 1.) / (2. * pi * fc * cw + 2. * pi * fc - cw + 1.); - let b2: Flt = 0.; - let a1: Flt = (2. * pi * fc * cw + 2. * pi * fc + cw - 1.) - / (2. * pi * fc * cw + 2. * pi * fc - cw + 1.); - let a2: Flt = 0.; + if fc >= fs / 2. { + bail!("Cuton frequency should be smaller than Nyquist frequency") + } + let b0 = pi*fc/(pi*fc+fs); + let b1 = b0; + let a1 = (pi*fc-fs)/(pi*fc+fs); - Ok(Biquad::fromCoefs(b0, b1, b2, a1, a2)) + Ok(Biquad::fromCoefs(b0, b1, 0., a1, 0.)) } /// Filter input signal, output by overwriting input slice. @@ -204,8 +224,9 @@ impl Filter for Biquad { Box::new(*self) } } -impl TransferFunction for Biquad { - fn tf(&self, fs: Flt, freq: VdView) -> Ccol { +impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for Biquad { + fn tf(&self, fs: Flt, freq: T) -> Ccol { + let freq = freq.into(); let res = freq.mapv(|f| { let z = Complex::exp(I * 2. * pi * f / fs); let num = self.b0 + self.b1 / z + self.b2 / z / z; @@ -232,15 +253,16 @@ mod test { #[test] fn test_firstOrderLowpass() { - let fs = 10.; - let fc = 1.; + let fs = 1e5; + let fc = 10.; let b = Biquad::firstOrderLowPass(fs, fc).unwrap(); - let mut freq = Dcol::from_elem((3), 0.); + let mut freq = Dcol::from_elem(5, 0.); freq[1] = fc; - freq[2] = fs/2.; + freq[2] = fs / 2.; let tf = b.tf(fs, freq.view()); - assert_abs_diff_eq!(tf[0].re,1.); - assert_abs_diff_eq!(tf[0].im,0.); - assert_abs_diff_eq!(tf[1].abs(),1./Flt::sqrt(2.)); + // println!("{:?}", tf); + assert_abs_diff_eq!(tf[0].re, 1., epsilon = 1e-6); + assert_abs_diff_eq!(tf[0].im, 0.); + assert_abs_diff_eq!(tf[1].abs(), 1. / Flt::sqrt(2.), epsilon = 1e-6); } } diff --git a/src/filter/biquadbank.rs b/src/filter/biquadbank.rs index 8b77914..c2fa116 100644 --- a/src/filter/biquadbank.rs +++ b/src/filter/biquadbank.rs @@ -66,8 +66,8 @@ impl BiquadBank { &mut self, py: Python<'py>, input: PyArrayLike1, - ) -> PyResult<&'py PyArray1> { - Ok(self.filter(input.as_slice()?).into_pyarray(py)) + ) -> PyResult> { + Ok(self.filter(input.as_slice()?).into_pyarray_bound(py)) } #[pyo3(name = "reset")] /// See: [BiquadBank::reset()] diff --git a/src/filter/dummy.rs b/src/filter/dummy.rs new file mode 100644 index 0000000..afea697 --- /dev/null +++ b/src/filter/dummy.rs @@ -0,0 +1,24 @@ +use itertools::Itertools; + +use super::*; + +/// A Dummy fillter just does 'nothing', its input equals its output +#[derive(Clone, Copy, Debug)] +pub struct DummyFilter; + +impl Filter for DummyFilter { + fn filter(&mut self, input: &[Flt]) -> Vec { + // Just returns an allocated copy + input.to_vec() + } + fn reset(&mut self) { } + fn clone_dyn(&self) -> Box { + Box::new(*self) + } +} +impl<'a, T: AsArray<'a, Flt>> TransferFunction<'a, T> for DummyFilter { + fn tf(&self, _fs: Flt, freq: T) -> Ccol { + let freq = freq.into(); + Ccol::ones(freq.len()) + } +} \ No newline at end of file diff --git a/src/filter/mod.rs b/src/filter/mod.rs index f59ebc2..2c87412 100644 --- a/src/filter/mod.rs +++ b/src/filter/mod.rs @@ -7,6 +7,7 @@ pub use super::config::*; mod biquad; mod biquadbank; +mod dummy; mod seriesbiquad; pub use biquad::Biquad; @@ -31,11 +32,14 @@ pub trait Filter: Send { /// Implementations are able to generate transfer functions of itself -pub trait TransferFunction: Send { +pub trait TransferFunction<'a, T>: Send +where + T: AsArray<'a, Flt>, +{ /// Compute frequency response (i.e. transfer function from input to output) /// /// Args - fn tf(&self, fs: Flt, freq: VdView) -> Ccol; + fn tf(&self, fs: Flt, freq: T) -> Ccol; } impl Clone for Box { diff --git a/src/filter/seriesbiquad.rs b/src/filter/seriesbiquad.rs index 2197527..ee482a0 100644 --- a/src/filter/seriesbiquad.rs +++ b/src/filter/seriesbiquad.rs @@ -22,11 +22,9 @@ pub struct SeriesBiquad { #[cfg(feature = "python-bindings")] #[cfg_attr(feature = "python-bindings", pymethods)] impl SeriesBiquad { - - - #[new] /// Create new series filter set. See [SeriesBiquad::new()] /// + #[new] pub fn new_py<'py>(coefs: PyReadonlyArrayDyn) -> PyResult { Ok(SeriesBiquad::new(coefs.as_slice()?)?) } diff --git a/src/ps/aps.rs b/src/ps/aps.rs index 0fb847c..3e8fe76 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -1 +1,5 @@ -//! Averaged power spectra module \ No newline at end of file +//! Averaged power spectra module. 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, \ No newline at end of file diff --git a/src/ps/ps.rs b/src/ps/ps.rs index c4f6401..45fa83c 100644 --- a/src/ps/ps.rs +++ b/src/ps/ps.rs @@ -1,3 +1,5 @@ +//! Power spectra estimator, that uses a Windowed FFT to estimate cross-power +//! spectra. Window functions are documented in the `window` module. use crate::config::*; use ndarray::parallel::prelude::*; use num::pow::Pow; @@ -14,7 +16,10 @@ use std::mem::MaybeUninit; use realfft::{RealFftPlanner, RealToComplex}; /// Singlesided cross-Power spectra computation engine. -struct PowerSpectra { +/// +/// Computes the signal(s) auto power and cross-power spectrum in each frequency +/// bin. +pub struct PowerSpectra { // Window used in estimator pub window: Window, // The window power, is corrected for in power spectra estimants @@ -100,6 +105,11 @@ impl PowerSpectra { /// Compute cross power spectra from input time data. First axis is /// frequency, second axis is channel i, third axis is channel j. + /// + /// # Argument + /// + /// * `tdata` - Input time data. This is a 2D array, where the first axis is + /// time and the second axis is the channel number. pub fn compute<'a, T>(&mut self, tdata: T) -> Array3 where T: Into>, @@ -130,9 +140,10 @@ impl PowerSpectra { Zip::from(&mut out) .and(chi) .and(chj) - .for_each(|out, chi, chjc| + .for_each(|out, chi, chjc|{ // Loop over frequency components - *out = 0.5 * chi * chjc + *out = 0.5 * chi * chjc; + } ); // The DC component has no 0.5 correction, as it only @@ -275,12 +286,14 @@ mod test { #[test] fn test_parseval_with_window() { - const nfft: usize = 48000; - let rect = Window::new(WindowType::Hann, nfft); - let mut ps = PowerSpectra::newFromWindow(rect); + // A sufficiently high value is required here, to show that it works. + const nfft: usize = 2usize.pow(20); + let window = Window::new(WindowType::Hann, nfft); + // let window = Window::new(WindowType::Rect, nfft); + let mut ps = PowerSpectra::newFromWindow(window); // Start with a time signal - let t: Dmat = Dmat::random((nfft, 1), StandardNormal); + let t: Dmat = 2.*Dmat::random((nfft, 1), StandardNormal); let tavg = t.sum()/(nfft as Flt); let t_dc_power = tavg.powi(2); @@ -295,7 +308,9 @@ mod test { let fpower = power.sum().abs(); assert_ulps_eq!(t_dc_power, power[(0,0,0)].abs(), epsilon = Flt::EPSILON * (nfft as Flt).powi(2)); - assert_ulps_eq!(signal_pwr, fpower, epsilon = Flt::EPSILON * (nfft as Flt).powi(2)); + + // This one fails when nfft is too short. + assert_ulps_eq!(signal_pwr, fpower, epsilon = 1e-2); } diff --git a/src/ps/window.rs b/src/ps/window.rs index 59216c1..7d19091 100644 --- a/src/ps/window.rs +++ b/src/ps/window.rs @@ -1,56 +1,70 @@ -//! Window functions designed for Welch' method. Implementations are given for common window -//! functions, as well as optimal 'jump' values `R`, that result in a certain overlap. +//! Window functions designed for Welch' method. Implementations are given for +//! the 5 classical window functions: +//! +//! * Rect - rectangular +//! * Hann - Von Hann window (sometimes wrongly called "Hanning") +//! * Bartlett +//! * Hamming +//! * Blackman +//! #![allow(non_snake_case)] use crate::config::*; -#[macro_use] -use strum_macros::{Display}; +use strum_macros::Display; -fn linspace(nfft: usize) -> Dcol { - Dcol::linspace(0., nfft as Flt, nfft) -} /// Von Hann window, often misnamed as the 'Hanning' window. -fn hann(nfft: usize) -> (Dcol, usize) { +fn hann(nfft: usize) -> Dcol { let nfftF = nfft as Flt; - ( - // The Window - linspace(nfft).mapv(|i| (pi * i / (nfftF+1.)).sin().powi(2)), - // The hopp size - (nfft) / 2, - ) + Dcol::from_iter((0..nfft).map(|i| { + // 0.5*(1-cos(2*pi*i/(n-1))) + 0.5 * (1. - Flt::cos(2. * pi * i as Flt / (nfftF - 1.))) + })) } -fn rect(nfft: usize) -> (Dcol, usize) { - (Dcol::ones(nfft), nfft) +/// Rectangular window +fn rect(nfft: usize) -> Dcol { + Dcol::ones(nfft) } -fn blackman(N: usize) -> (Dcol, usize) { - let a0 = 7938. / 18608.; - let a1 = 9240. / 18608.; - let a2 = 1430. / 18608.; +// Blackman window +fn blackman(N: usize) -> Dcol { + // Exact a0 coefficient, approximate value is 0.42 + const a0: Flt = 7938. / 18608.; + // Exact a1 coefficient, approximate value is 0.50 + const a1: Flt = 9240. / 18608.; + // Exact a2 coefficient, approximate value is 0.08 + const a2: Flt = 1430. / 18608.; + let Nf = N as Flt; - let lin = linspace(N); - ( - a0 - a1 * ((2. * pi / Nf) * lin.clone()).mapv(|x| x.cos()) - + a2 * ((4. * pi / Nf) * lin).mapv(|x| x.cos()), - // The hop size - N / 3, - ) + + Dcol::from_iter((0..N).map(|i| { + let iF = i as Flt; + a0 - a1 * Flt::cos(2. * pi * iF / (Nf - 1.)) + a2 * Flt::cos(4. * pi * iF / (Nf - 1.)) + })) } -fn bartlett(nfft: usize) -> (Dcol, usize) { - let Nf = nfft as Flt; - ( - (1. - (2. * (linspace(nfft) - (Nf - 1.) / 2.) / Nf)).mapv(|x| x.abs()), - // The hop size - nfft / 2, - ) +fn bartlett(N: usize) -> Dcol { + let Nf = N as Flt; + Dcol::from_iter((0..N).map(|i| { + let iF = i as Flt; + if i <= (N - 1) / 2 { + 2. * iF / (Nf - 1.) + } else { + 2. - 2. * iF / (Nf - 1.) + } + })) } -fn hamming(nfft: usize) -> (Dcol, usize) { - let alpha = 25.0 / 46.0; - let beta = (1. - alpha) / 2.; - let Nf = nfft as Flt; - ( - alpha + 2. * beta * ((2. * pi / (Nf - 0.)) * linspace(nfft)).mapv(|x| x.sin()), - // The hop size - (nfft) / 2, +fn hamming(N: usize) -> Dcol { + let Nf = N as Flt; + // Approx 0.54 + const a0: Flt = 25.0 / 46.0; + // Approx 0.46 + const a1: Flt = 1. - a0; + + Dcol::from_iter((0..N).map(|i| + // Alessio et al. + a0 - a1 * Flt::cos(2. * pi * i as Flt / (Nf - 1.)) + + // end of map + ) + // end from iter ) } @@ -90,13 +104,14 @@ impl Window { if nfft % 2 != 0 { panic!("Requires even nfft"); } - let (win, R) = match w { + let win = match w { WindowType::Hann => hann(nfft), WindowType::Hamming => hamming(nfft), WindowType::Rect => rect(nfft), WindowType::Bartlett => bartlett(nfft), WindowType::Blackman => blackman(nfft), }; + let R = nfft/2; Window { w, win, R } } /// Convenience function that returns the length of the window. @@ -107,62 +122,20 @@ impl Window { // #[cfg(test)] mod test { + use approx::assert_abs_diff_eq; + use super::*; - #[test] - fn test_linspace() { - assert!(linspace(2)[0] == 0.); - // println!("{:?}", linspace(3)); - assert!(linspace(3)[1] == 1.); - assert!(linspace(4).len() == 4); - } #[test] - fn test_cola_hann() { - let nfft = 66; - let hann = Window::new(WindowType::Hann, nfft); - println!("{:?}", hann.win); - - let mut hanntot = Dcol::zeros(hann.len() * 4); - - assert!(2 * hann.R == nfft); - hanntot.slice_mut(s![0..nfft]).assign(&hann.win); - hanntot - .slice_mut(s![hann.R..nfft + hann.R]) - .scaled_add(1.0, &hann.win); - - hanntot - .slice_mut(s![nfft..2 * nfft]) - .scaled_add(1.0, &hann.win); - hanntot - .slice_mut(s![nfft + hann.R..2 * nfft + hann.R]) - .scaled_add(1.0, &hann.win); - - hanntot - .slice_mut(s![2 * nfft..3 * nfft]) - .scaled_add(1.0, &hann.win); - hanntot - .slice_mut(s![2 * nfft + hann.R..3 * nfft + hann.R]) - .scaled_add(1.0, &hann.win); - - println!("{:?}", hanntot); - } - - #[test] - fn tets_cola_hamming() { - let nfft = 25; - let ham = Window::new(WindowType::Hamming, nfft); - let mut hamtot = Dcol::zeros(ham.len() * 3); - - assert!(2 * ham.R == nfft); - hamtot.slice_mut(s![0..nfft]).scaled_add(1.0, &ham.win); - // println!("{:?}", hamtot); - hamtot - .slice_mut(s![ham.R..nfft + ham.R]) - .scaled_add(1.0, &ham.win); - hamtot - .slice_mut(s![nfft..2 * nfft]) - .scaled_add(1.0, &ham.win); - // println!("{:?}", hamtot); - // hantot.slice_mut(s![1+2*han1.R..nfft+1+2*han1.R]).scaled_add(1.0, &han2.win); + fn test_windows(){ + let Hann = hann(11); + let Hamming = hamming(11); + let Bartlett = bartlett(11); + let Blackmann = bartlett(11); + // let h = hann(11); + assert_eq!(Hann[5] , 1.); + assert_eq!(Hamming[5] , 1.); + assert_eq!(Bartlett[5] , 1.); + assert_eq!(Blackmann[5] , 1.); } } diff --git a/src/timebuffer.rs b/src/timebuffer.rs index e69de29..90e4f1b 100644 --- a/src/timebuffer.rs +++ b/src/timebuffer.rs @@ -0,0 +1,8 @@ +use super::*; +/// A time buffer is used as intermediate storage for samples, to spare up +/// enough values to do 'something' with. Typical application: getting enough +/// samples to perform an FFT. The use case is computing average power spectra. +/// +struct TimeBuffer { + +} \ No newline at end of file