From b5af58581d094587b84e31f2a05d300b9e57bdaa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Robert=20J=C3=B6rdens?= <rj@quartiq.de>
Date: Tue, 19 Sep 2023 14:49:25 +0200
Subject: [PATCH] add noise test

---
 Cargo.toml      |  2 +-
 src/bin/main.rs | 11 +++++++++--
 src/psd.rs      | 17 +++++++++++++----
 src/source.rs   | 34 ++++++++++++++++++++++++++++++++--
 4 files changed, 55 insertions(+), 9 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index f57b276..99d62dc 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -27,7 +27,7 @@ anyhow = "1.0.75"
 socket2 = "0.5.3"
 idsp = "0.11.0"
 rustfft = "6.1.0"
-rand = "0.8.5"
+rand = { version = "0.8.5", features = ["small_rng"] }
 derive_builder = "0.12.0"
 
 #[build-dependencies]
diff --git a/src/bin/main.rs b/src/bin/main.rs
index 4809d36..361bbaa 100644
--- a/src/bin/main.rs
+++ b/src/bin/main.rs
@@ -30,11 +30,18 @@ pub struct Opts {
 
     #[arg(short, long, default_value_t = 4)]
     min_avg: usize,
+
+    #[arg(short, long, default_value = "mid")]
+    detrend: Detrend,
 }
 
 fn main() -> Result<()> {
     env_logger::init();
-    let Opts { source, min_avg } = Opts::parse();
+    let Opts {
+        source,
+        min_avg,
+        detrend,
+    } = Opts::parse();
 
     let (cmd_send, cmd_recv) = mpsc::channel();
     let (trace_send, trace_recv) = mpsc::sync_channel(1);
@@ -54,7 +61,7 @@ fn main() -> Result<()> {
                 dec.extend((0..4).map(|_| {
                     let mut c = PsdCascade::<{ 1 << 9 }>::default();
                     c.set_stage_depth(3);
-                    c.set_detrend(Detrend::Mid);
+                    c.set_detrend(detrend);
                     c
                 }));
                 i = 0;
diff --git a/src/psd.rs b/src/psd.rs
index af155c7..bdf4969 100644
--- a/src/psd.rs
+++ b/src/psd.rs
@@ -51,7 +51,7 @@ impl<const N: usize> Window<N> {
 }
 
 /// Detrend method
-#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord)]
+#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, clap::ValueEnum)]
 pub enum Detrend {
     /// No detrending
     None,
@@ -63,10 +63,14 @@ pub enum Detrend {
     // TODO: linear
 }
 
+impl core::fmt::Display for Detrend {
+    fn fmt(&self, f: &mut core::fmt::Formatter) -> core::fmt::Result {
+        core::fmt::Debug::fmt(self, f)
+    }
+}
+
 /// Power spectral density accumulator and decimator
 ///
-/// Note: Don't feed more than N*1e7 items without expecting loss of accuracy
-///
 /// One stage in [PsdCascade].
 #[derive(Clone)]
 pub struct Psd<const N: usize> {
@@ -162,6 +166,11 @@ pub trait PsdStage {
     /// Overlap is kept.
     /// Decimation is performed on fully processed input items.
     ///
+    /// Note: When feeding more than ~N*1e6 items expect loss of accuracy
+    /// due to rounding errors on accumulation.
+    ///
+    /// Note: Also be aware of the usual accuracy limitation of the item data type
+    ///
     /// # Args
     /// * `x`: input items
     /// * `y`: output items
@@ -435,7 +444,7 @@ impl<const N: usize> PsdCascade<N> {
 mod test {
     use super::*;
 
-    /// 36 insns per input sample: > 190 MS/s per skylake core
+    /// 36 insns per item: > 190 MS/s per skylake core
     #[test]
     #[ignore]
     fn insn() {
diff --git a/src/source.rs b/src/source.rs
index a8c6f2a..7a691a9 100644
--- a/src/source.rs
+++ b/src/source.rs
@@ -1,6 +1,7 @@
 use crate::{Frame, Loss};
 use anyhow::Result;
 use clap::Parser;
+use rand::{rngs::SmallRng, Rng, SeedableRng};
 use std::io::ErrorKind;
 use std::time::Duration;
 use std::{
@@ -34,6 +35,10 @@ pub struct SourceOpts {
     /// Single f32 raw trace in file, architecture dependent
     #[arg(short, long)]
     single: Option<String>,
+
+    /// Power law noise with psd f^noise.
+    #[arg(short, long)]
+    noise: Option<i32>,
 }
 
 #[derive(Debug)]
@@ -41,6 +46,7 @@ enum Data {
     Udp(std::net::UdpSocket),
     File(BufReader<File>),
     Single(BufReader<File>),
+    Noise((SmallRng, Vec<f64>)),
 }
 
 pub struct Source {
@@ -51,7 +57,12 @@ pub struct Source {
 
 impl Source {
     pub fn new(opts: SourceOpts) -> Result<Self> {
-        let data = if let Some(file) = &opts.file {
+        let data = if let Some(noise) = opts.noise {
+            Data::Noise((
+                SmallRng::seed_from_u64(0x7654321),
+                vec![0.0; noise.unsigned_abs() as _],
+            ))
+        } else if let Some(file) = &opts.file {
             Data::File(BufReader::with_capacity(1 << 20, File::open(file)?))
         } else if let Some(single) = &opts.single {
             Data::Single(BufReader::with_capacity(1 << 20, File::open(single)?))
@@ -70,9 +81,26 @@ impl Source {
     }
 
     pub fn get(&mut self) -> Result<Vec<Vec<f32>>> {
-        let mut buf = [0u8; 2048];
         Ok(match &mut self.data {
+            Data::Noise((rng, state)) => {
+                vec![(0..1024)
+                    .map(|_| {
+                        let mut x = (rng.gen::<f32>() - 0.5) as f64; // *6.0f64.sqrt();
+                        let diff = self.opts.noise.unwrap() > 0;
+                        for s in state.iter_mut() {
+                            if diff {
+                                (x, *s) = (x - *s, x);
+                            } else {
+                                *s += x;
+                                x = *s;
+                            }
+                        }
+                        x as _
+                    })
+                    .collect()]
+            }
             Data::File(fil) => loop {
+                let mut buf = [0u8; 2048];
                 match fil.read_exact(&mut buf[..self.opts.frame_size]) {
                     Ok(()) => {
                         let frame = Frame::from_bytes(&buf[..self.opts.frame_size])?;
@@ -86,6 +114,7 @@ impl Source {
                 }
             },
             Data::Single(fil) => loop {
+                let mut buf = [0u8; 2048];
                 match fil.read(&mut buf[..]) {
                     Ok(len) => {
                         if len == 0 && self.opts.repeat {
@@ -99,6 +128,7 @@ impl Source {
                 }
             },
             Data::Udp(socket) => {
+                let mut buf = [0u8; 2048];
                 let len = socket.recv(&mut buf[..])?;
                 let frame = Frame::from_bytes(&buf[..len])?;
                 self.loss.update(&frame);