Aps working, or almost working?

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<CrossPowerSpecra>),
+    AllIntermediateResults(Vec<CPS>),
 
     /// 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<AvPowerSpectra> {
-        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<WindowType>,
         overlap: Option<Overlap>,
@@ -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);
+            }
+        }
+    }
+}

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<Cflt>;
-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};

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<Cflt>;
+
+/// 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<Flt>;
+
+    /// 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<usize>) -> Array1<Cflt>;
+}
+
+
+impl CrossPowerSpecra for CPS {
+    fn ap(&self, ch: usize) -> Array1<Flt> {
+        // 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<usize>) -> Array1<Cflt> {
+        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>,
     {

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);
+    }
 }