From 19484fd1cbe7f19e048920797271ad06d144ea80 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Wed, 10 Jul 2024 00:10:34 +0200 Subject: [PATCH] Aps working, or almost working? --- src/ps/aps.rs | 257 +++++++++++++++++++++++++++++++++++++++---- src/ps/mod.rs | 8 +- src/ps/ps.rs | 68 +++++++++++- src/ps/timebuffer.rs | 20 +++- 4 files changed, 327 insertions(+), 26 deletions(-) diff --git a/src/ps/aps.rs b/src/ps/aps.rs index b7a0c9c..b7b3f77 100644 --- a/src/ps/aps.rs +++ b/src/ps/aps.rs @@ -1,4 +1,5 @@ use super::timebuffer::TimeBuffer; +use super::CrossPowerSpecra; use super::*; use crate::config::*; use anyhow::{bail, Result}; @@ -12,6 +13,8 @@ pub enum Overlap { 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 { @@ -20,14 +23,14 @@ impl Default for Overlap { } /// Result from [AvPowerspectra.compute()] -pub enum ApsResult { +pub enum ApsResult<'a> { /// Returns all intermediate results. Useful when evolution over time should /// be visible. I.e. in case of spectrograms, but also in case the /// converging process to the average should be made visible. - AllIntermediateResults(Vec), + AllIntermediateResults(Vec), /// Give only last result back, the most recent value. - OnlyLastResult(CrossPowerSpecra), + OnlyLastResult(&'a CPS), /// No new data available. Nothing given back. None, @@ -70,6 +73,7 @@ pub struct AvPowerSpectra { // Power spectra estimator for single block ps: PowerSpectra, + // The mode with which the AvPowerSpectra is computing. mode: ApsMode, // The number of samples to keep in the time buffer when overlapping time @@ -83,7 +87,7 @@ pub struct AvPowerSpectra { timebuf: TimeBuffer, // Current estimation of the power spectra - cur_est: Cmat, + cur_est: CPS, } impl AvPowerSpectra { /// The FFT Length of estimating (cross)power spectra @@ -106,6 +110,7 @@ impl AvPowerSpectra { } // 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 { @@ -122,7 +127,7 @@ impl AvPowerSpectra { /// /// * `nfft` - FFT Length pub fn new_simple(nfft: usize) -> Result { - return AvPowerSpectra::new(nfft, None, None, None); + return AvPowerSpectra::build(nfft, None, None, None); } /// Create power spectra estimator which weighs either all data @@ -144,7 +149,7 @@ impl AvPowerSpectra { /// - `overlap` - Amount of overlap in /// - `mode` - The mode in which the /// - pub fn new( + pub fn build( nfft: usize, windowtype: Option, overlap: Option, @@ -179,10 +184,84 @@ impl AvPowerSpectra { overlap_keep, mode, N: 0, - cur_est: Cmat::default((0, 0)), + cur_est: CPS::default((0, 0, 0)), timebuf: TimeBuffer::new(), }) } + // Update result for single block + fn update_singleblock<'a, T>(&mut self, timedata: T) + where + T: AsArray<'a, Flt, Ix2>, + { + let timeblock = timedata.into(); + let Cpsnew = self.ps.compute(&timeblock); + // println!("Cpsnew: {:?}", Cpsnew[[0, 0, 0]]); + + // Initialize to zero + if self.cur_est.len() == 0 { + assert_eq!(self.N, 0); + self.cur_est = CPS::zeros(Cpsnew.raw_dim().f()); + } + + // Update the number of blocks processed + self.N += 1; + + // Apply operation based on mode + match self.mode { + ApsMode::AllAveraging => { + let Nf = Cflt { + re: self.N as Flt, + im: 0., + }; + self.cur_est = (Nf - 1.) / Nf * &self.cur_est + 1. / Nf * Cpsnew; + } + ApsMode::ExponentialWeighting { fs, tau } => { + debug_assert!(self.N > 0); + if self.N == 1 { + self.cur_est = Cpsnew; + } else { + // A sound level meter specifies a low pass filter with one + // real pole at -1/tau, for a time weighting of tau. This + // means the analogue transfer function is 1 / + // (tau*s+1). + + // Now we want to approximate this with a digital transfer + // function. The step size, or sampling period is: + // T = (nfft-overlap_keep)/fs. + + // Then, when using the matched z-transform ( + // https://en.wikipedia.org/wiki/Matched_Z-transform_method), + // an 1/(s-p) is replaced by 1/(1-exp(p*T)*z^-1). + + // So the digital transfer function will be: + // H[n] ≅ K / (1 - exp(-T/tau) * z^-1). + // , where K is a to-be-determined constant, for which we + // take the value such that the D.C. gain equals 1. To get + // the frequency response at D.C., we have to set z=1, so + // to set the D.C. to 1, we have to set K = 1-exp(-T/tau). + + // Hence: + // H[n] ≅ (1- exp(-T/tau)) / (1 - exp(-T/tau) * z^-1). + + // Or as a finite difference equation: + // + // y[n] * (1-exp(-T/tau)* z^-1) = (1-exp(-T/tau)) * x[n] + + // or, finally: + // y[n] = alpha * y[n-1] + (1-alpha) * x[n] + + // where alpha = exp(-T/tau). + let T = (self.nfft() - self.overlap_keep) as Flt / fs; + let alpha = Cflt::ONE * Flt::exp(-T / tau); + self.cur_est = alpha * &self.cur_est + (1. - alpha) * Cpsnew; + } + } + + ApsMode::Spectrogram => { + self.cur_est = Cpsnew; + } + } + } /// Computes average (cross)power spectra, and returns only the most recent /// estimate, if any can be given back. Only gives back a result when enough /// data is available. @@ -191,38 +270,178 @@ impl AvPowerSpectra { /// /// * `timedata``: New available time data. Number of columns is number of /// channels, number of rows is number of frames (samples per channel). - /// * `giveInterMediateResults` + /// * `giveInterMediateResults` - If true: returns a vector of all + /// intermediate results. Useful when plotting spectra over time. /// /// # Panics /// /// If timedata.ncols() does not match number of columns in already present /// data. - fn compute<'a, T>(&mut self, timedata: T, giveInterMediateResults: bool) -> ApsResult + pub fn compute<'a, 'b, T>( + &'a mut self, + timedata: T, + giveInterMediateResults: bool, + ) -> ApsResult<'a> where - T: AsArray<'a, Flt, Ix2>, + T: AsArray<'b, Flt, Ix2>, { // Push new data in the time buffer. self.timebuf.push(timedata); + // Storage for the result + let mut result = ApsResult::None; + + // Flag to indicate that we have obtained one result for sure. + let mut computed_single = false; + // Iterate over all blocks that can come, while let Some(timeblock) = self.timebuf.pop(self.nfft(), self.overlap_keep) { // Compute cross-power spectra for current time block - let Cpsnew = self.ps.compute(&timeblock); - // Ok, mode-determined - match self.mode { - ApsMode::AllAveraging => {} - ApsMode::ExponentialWeighting { fs, tau } => {} - ApsMode::Spectrogram => {} + self.update_singleblock(&timeblock); + + // We ha + computed_single = true; + if giveInterMediateResults { + // Initialize with empty vector + if let ApsResult::None = result { + result = ApsResult::AllIntermediateResults(Vec::new()) + } + } + + // So + if let ApsResult::AllIntermediateResults(v) = &mut result { + v.push(self.cur_est.clone()); } - // if self.cur_est.len() > 0 } - ApsResult::None + + if computed_single && !giveInterMediateResults { + // We have computed it once, but we are not interested in + // intermediate results. So we return a reference to the last data. + return ApsResult::OnlyLastResult(&self.cur_est); + } + result } /// See [AvPowerSpectra](compute()) - fn compute_last<'a, T>(&mut self, timedata: T) -> ApsResult + pub fn compute_last<'a, T>(&mut self, timedata: T) -> ApsResult where T: AsArray<'a, Flt, Ix2>, { return self.compute(timedata, false); } } + +#[cfg(test)] +mod test { + use approx::assert_abs_diff_eq; + use ndarray_rand::rand_distr::Normal; + use ndarray_rand::RandomExt; + + use super::CrossPowerSpecra; + use crate::{config::*, ps::ApsResult}; + + use super::{ApsMode, AvPowerSpectra, Overlap, WindowType, CPS}; + + #[test] + fn test_overlap_keep() { + let nfft = 10; + assert_eq!( + AvPowerSpectra::get_overlap_keep(nfft, Overlap::Percentage(50.)).unwrap(), + nfft / 2 + ); + let nfft = 11; + assert_eq!( + AvPowerSpectra::get_overlap_keep(nfft, Overlap::Percentage(50.)).unwrap(), + nfft / 2 + ); + let nfft = 1024; + assert_eq!( + AvPowerSpectra::get_overlap_keep(nfft, Overlap::Percentage(25.)).unwrap(), + nfft / 4 + ); + assert_eq!( + AvPowerSpectra::get_overlap_keep(nfft, Overlap::Number(10)).unwrap(), + 10 + ); + } + + /// When the time constant is 1.0, every second the powers approximately + /// halve. That is the subject of this test. + #[test] + fn test_expweighting() { + let nfft = 48000; + let fs = nfft as Flt; + let tau = 2.; + let mut aps = AvPowerSpectra::build( + nfft, + None, + Some(Overlap::NoOverlap), + Some(ApsMode::ExponentialWeighting { fs, tau }), + ) + .unwrap(); + assert_eq!(aps.overlap_keep, 0); + + let distr = Normal::new(1.0, 1.0).unwrap(); + let timedata_some = Dmat::random((nfft, 1), distr); + let timedata_zeros = Dmat::zeros((nfft, 1)); + + let mut first_result = CPS::zeros((0, 0, 0)); + + if let ApsResult::OnlyLastResult(v) = aps.compute_last(timedata_some.view()) { + first_result = v.clone(); + // println!("{:?}", first_result.ap(0)[0]); + } else { + assert!(false, "Should return one value"); + } + let overlap_keep = AvPowerSpectra::get_overlap_keep(nfft, Overlap::NoOverlap).unwrap(); + if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata_zeros) { + } else { + assert!(false, "Should return one value"); + } + if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata_zeros) { + let alpha = Flt::exp(-((nfft - overlap_keep) as Flt) / (fs * tau)); + // let alpha: Flt = 1.; + for i in 0..nfft / 2 + 1 { + // println!("i={i}"); + assert_abs_diff_eq!(first_result.ap(0)[i] * alpha.powi(2), v.ap(0)[i]); + } + } else { + assert!(false, "Should return one value"); + } + assert_eq!(aps.N, 3); + } + + #[test] + fn test_tf1() { + let nfft = 4800; + let distr = Normal::new(1.0, 1.0).unwrap(); + let mut timedata = Dmat::random((nfft, 1), distr); + timedata + .push_column(timedata.column(0).mapv(|a| 2. * a).view()) + .unwrap(); + + let mut aps = AvPowerSpectra::build(nfft, None, None, None).unwrap(); + if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata) { + let tf = v.tf(0, 1, None); + assert_eq!((&tf - 2.0 * Cflt::ONE).sum().abs(), 0.0); + } else { + assert!(false); + } + } + #[test] + fn test_ap() { + let nfft = 256; + let distr = Normal::new(1.0, 1.0).unwrap(); + let timedata = Dmat::random((500 * nfft, 1), distr); + let timedata_mean_square = (&timedata * &timedata).sum() / (timedata.len() as Flt); + + for wt in [Some(WindowType::Rect), Some(WindowType::Hann), Some(WindowType::Bartlett), Some(WindowType::Blackman), None] { + let mut aps = AvPowerSpectra::build(nfft, wt, None, None).unwrap(); + if let ApsResult::OnlyLastResult(v) = aps.compute_last(&timedata) { + let ap = v.ap(0); + assert_abs_diff_eq!((&ap).sum().abs(), timedata_mean_square, epsilon = 4e-3); + } else { + assert!(false); + } + } + } +} diff --git a/src/ps/mod.rs b/src/ps/mod.rs index 21d66d2..32ffdfd 100644 --- a/src/ps/mod.rs +++ b/src/ps/mod.rs @@ -3,7 +3,7 @@ //! [AvPowerSpectra] using //! Welch' method, and windows for time-windowing the data with non-rectangular //! windows (also known as 'tapers'). -//! +//! mod aps; mod fft; mod ps; @@ -11,7 +11,7 @@ mod timebuffer; mod window; use crate::config::*; -pub type CrossPowerSpecra = Array3; -pub use aps::{ApsMode, AvPowerSpectra, Overlap, ApsResult}; -pub use ps::PowerSpectra; + +pub use aps::{ApsMode, ApsResult, AvPowerSpectra, Overlap}; +pub use ps::{CrossPowerSpecra, PowerSpectra, CPS}; pub use window::{Window, WindowType}; diff --git a/src/ps/ps.rs b/src/ps/ps.rs index 7366963..3843dd3 100644 --- a/src/ps/ps.rs +++ b/src/ps/ps.rs @@ -7,12 +7,76 @@ use std::usize; use crate::Dcol; -use super::{fft::FFT, CrossPowerSpecra}; use super::window::*; +use super::fft::FFT; use std::mem::MaybeUninit; use realfft::{RealFftPlanner, RealToComplex}; +/// Cross power spectra, which is a 3D array, with the following properties: +/// +/// - The first index is the frequency index, starting at DC, ending at nfft/2. +/// - The second, and third index result in [i,j] = C_ij = p_i * conj(p_j) +/// +pub type CPS = Array3; + +/// Extra typical methods that are of use for 3D-arrays of complex numbers, that +/// are typically implemented as cross-power spectra. +pub trait CrossPowerSpecra { + /// Returns the autopower for a single channel, as a array of real values + /// (imaginary part is zero and is stripped off). + /// + /// # Args + /// + /// - `ch` - The channel number to compute autopower for. + fn ap(&self, ch: usize) -> Array1; + + /// Returns the transfer function from `chi` to `chj`, that is ~ Pj/Pi a + /// single channel, as a array of complex numbers. + /// + /// # Args + /// + /// - `chi` - The channel number of the *denominator* + /// - `chj` - The channel number of the *numerator* + /// - `chRef` - Optional, a reference channel that has the lowest noise. If + /// not given, the average of the two autopowers is used, which gives + /// always a worse result than when two a low noise reference channel is + /// used. + /// + fn tf(&self, chi: usize, chj: usize, chRef: Option) -> Array1; +} + + +impl CrossPowerSpecra for CPS { + fn ap(&self, ch: usize) -> Array1 { + // Slice out one value for all frequencies, map to only real part, and + // return. + self.slice(s![.., ch, ch]).mapv(|f| f.re) + } + // fn apsp + fn tf(&self, chi: usize, chj: usize, chRef: Option) -> Array1 { + match chRef { + None => { + let cij = self.slice(s![.., chi, chj]); + let cii = self.slice(s![.., chi, chi]); + let cjj = self.slice(s![.., chj, chj]); + Zip::from(cij) + .and(cii) + .and(cjj) + .par_map_collect(|cij, cii, cjj| 0.5 * (cij.conj() / cii + cjj / cij)) + } + Some(chr) => { + let cir = self.slice(s![.., chi, chr]); + let cjr = self.slice(s![.., chj, chr]); + + Zip::from(cir) + .and(cjr) + .par_map_collect(|cir, cjr| cjr / cir) + } + } + } +} + /// Single-sided (cross)power spectra estimator, that uses a Windowed FFT to estimate cross-power /// spectra. Window functions are documented in the `window` module. Note that /// directly using this power spectra estimator is generally not useful as it is @@ -128,7 +192,7 @@ impl PowerSpectra { /// (nfft/2+1,timedata.ncols(), timedata.ncols()). Its content is: /// [freq_index, chi, chj] = crosspower: chi*conj(chj) /// - pub fn compute<'a, T>(&mut self, tdata: T) -> CrossPowerSpecra + pub fn compute<'a, T>(&mut self, tdata: T) -> CPS where T: AsArray<'a, Flt, Ix2>, { diff --git a/src/ps/timebuffer.rs b/src/ps/timebuffer.rs index d27b104..ad4bf02 100644 --- a/src/ps/timebuffer.rs +++ b/src/ps/timebuffer.rs @@ -67,7 +67,7 @@ impl TimeBuffer { let nsamples_available = c1.len(); debug_assert!(nsamples_available == self.nsamples()); if nsamples_available < nsamples_requested { - println!("Less available than requested"); + // println!("Less available than requested"); return None; } // Number of channels @@ -159,4 +159,22 @@ mod test { assert_eq!(tres.shape(), [50, 2]); assert_eq!(t2.slice(s![..50, ..2]), tres); } + #[test] + fn test_timebuffer3() { + let mut tb = TimeBuffer::new(); + let t1 = Dmat::zeros((10,1)); + tb.push(&t1); + tb.push(&t1); + assert_eq!(tb.nsamples(), 20); + + let tres = tb.pop(10, 5).unwrap(); + assert_eq!(tres.shape(), [10, 1]); + assert_eq!(tb.nsamples(), 15); + let tres = tb.pop(10, 0).unwrap(); + assert_eq!(tb.nsamples(), 5); + let tres = tb.pop(5, 5).unwrap(); + assert_eq!(tres.shape(), [5, 1]); + assert_eq!(tb.nsamples(), 5); + // assert_eq!(t2.slice(s![..50, ..2]), tres); + } }