Intermediate commit. Added some stuff related to biquads, added some stuff related to power spectra estimation

Cargo.toml

@@ -78,10 +78,17 @@ uuid = { version = "1.6.1", features = ["v4"] , optional = true}
 # Command line argument parser, for CLI apps
 clap = { version = "4.4.11", features = ["derive", "color", "help", "suggestions"] }
 
+# FFT's
+realfft = "3.3.0"
+
+[dev-dependencies]
+approx = "0.5.1"
+ndarray-rand = "0.14.0"
+
 [features]
-# default = ["f64", "cpal-api", "record"]
+default = ["f64", "cpal-api", "record"]
 # Use this for debugging extensions
-default = ["f64", "python-bindings", "record", "cpal-api"]
+# default = ["f64", "python-bindings", "record", "cpal-api"]
 
 cpal-api = ["dep:cpal"]
 record = ["dep:hdf5-sys", "dep:hdf5", "dep:chrono", "dep:uuid"]

src/config.rs

@@ -25,6 +25,7 @@ if #[cfg(feature = "python-bindings")] {
     pub use numpy::ndarray::{ArrayD, ArrayViewD, ArrayViewMutD};
     pub use numpy::ndarray::prelude::*;
     pub use numpy::{IntoPyArray,PyArray, PyArray1, PyArrayDyn, PyArrayLike1, PyReadonlyArrayDyn};
+    pub use numpy::ndarray::Zip;
     pub use pyo3::prelude::*;
     pub use pyo3::exceptions::PyValueError;
     pub use pyo3::{pymodule, types::PyModule, PyResult};
@@ -32,8 +33,8 @@ if #[cfg(feature = "python-bindings")] {
     pub use pyo3;
 
 } else {
-
     pub use ndarray::prelude::*;
+    pub use ndarray::Zip;
     pub use ndarray::{Array1, Array2, ArrayView1};
 } }
 

src/filter/biquad.rs

@@ -1,12 +1,35 @@
 use super::*;
 use crate::config::*;
-use anyhow::{Result, bail};
+use anyhow::{bail, Result};
 use num::Complex;
 
 #[cfg_attr(feature = "python-bindings", pyclass)]
 #[derive(Clone, Copy, Debug)]
 /// # A biquad is a second order recursive filter structure.
 ///
+/// This implementation only allows for normalized coefficients (a_0 = 1). It
+///  performs the following relation of output to input:
+///
+///  ```
+///  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]
+///  ```
+///
+///  The coefficients can be generated for typical standard type of biquad
+///  filters, such as low pass, high pass, bandpass (first order), low shelf,
+///  high shelf, peaking and notch filters.
+///
+///  The transfer function is:
+///
+/// ```
+///          b_0 + b_1 z^-1 + b_2 * z^-2
+///  H[z] = -----------------------------
+///          1   + a_1 z^-1 + a_2 * z^-2
+///  ```
+///
+///  And the frequency response can be found by filling in in above equation z =
+///  exp(i*omega/fs), where fs is the sampling frequency and omega is the radian
+///  frequency at which the transfer function is evaluated.
 ///
 pub struct Biquad {
     // State parameters
@@ -71,6 +94,22 @@ impl Biquad {
         }
     }
 
+    /// Construct a Biquad with 0 initial state from coefficients given as
+    /// arguments.
+    ///
+    /// *CAREFUL*: No checks are don  on validity / stability of the created filter!
+    fn fromCoefs(b0: Flt, b1: Flt, b2: Flt, a1: Flt, a2: Flt) -> Biquad {
+        Biquad {
+            w1: 0.,
+            w2: 0.,
+            b0,
+            b1,
+            b2,
+            a1,
+            a2,
+        }
+    }
+
     /// Create unit impulse response biquad filter. Input = output
     fn unit() -> Biquad {
         let filter_coefs = &[1., 0., 0., 1., 0., 0.];
@@ -101,19 +140,37 @@ impl Biquad {
         let facnum = 2. * fs * tau / (1. + 2. * fs * tau);
         let facden = (1. - 2. * fs * tau) / (1. + 2. * fs * tau);
 
-        let coefs = [
+        Ok(Biquad::fromCoefs(
             facnum,  // b0
             -facnum, // b1
-            0.,      // b2
+            0.,      // b2,
-            1.,      // a0
             facden,  // a1
             0.,      // a2
-        ];
+        ))
+    }
 
-        Ok(Biquad::new(&coefs).unwrap())
+    /// First order low pass filter (one pole in the real axis). No pre-warping
+    /// correction done.
+    pub fn firstOrderLowPass(fs: Flt, fc: Flt) -> Result<Biquad> {
+        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"),
+            _ => (),
+        }
+        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.;
+
+        Ok(Biquad::fromCoefs(b0, b1, b2, a1, a2))
     }
 
     /// Filter input signal, output by overwriting input slice.
+    #[inline]
     pub fn filter_inout(&mut self, inout: &mut [Flt]) {
         for sample in inout.iter_mut() {
             let w0 = *sample - self.a1 * self.w1 - self.a2 * self.w2;
@@ -125,6 +182,12 @@ impl Biquad {
         // println!("{:?}", inout);
     }
 }
+impl Default for Biquad {
+    /// Unit impulse (does not transform signal whatsoever)
+    fn default() -> Self {
+        Biquad::unit()
+    }
+}
 
 impl Filter for Biquad {
     fn filter(&mut self, input: &[Flt]) -> Vec<Flt> {
@@ -150,11 +213,13 @@ impl TransferFunction for Biquad {
             num / den
         });
         res
-        // re
     }
 }
 #[cfg(test)]
 mod test {
+    use approx::assert_abs_diff_eq;
+    use num::complex::ComplexFloat;
+
     use super::*;
 
     #[test]
@@ -164,4 +229,18 @@ mod test {
         let filtered = ser.filter(&inp);
         assert_eq!(&filtered, &inp);
     }
+
+    #[test]
+    fn test_firstOrderLowpass() {
+        let fs = 10.;
+        let fc = 1.;
+        let b = Biquad::firstOrderLowPass(fs, fc).unwrap();
+        let mut freq = Dcol::from_elem((3), 0.);
+        freq[1] = fc;
+        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.));
+    }
 }

src/lib.rs

@@ -13,11 +13,11 @@ mod config;
 use config::*;
 
 pub use config::Flt;
-// pub mod window;
-// pub mod ps;
 pub mod filter;
 pub mod daq;
+pub mod ps;
 pub mod siggen;
+mod math;
 
 use filter::*;
 use daq::*;

src/ps/aps.rs

@@ -0,0 +1 @@
+//! Averaged power spectra module

src/ps/fft.rs

@@ -0,0 +1,50 @@
+//! Compute forward single sided amplitude spectra
+use crate::config::*;
+use realfft::{RealFftPlanner, RealToComplex};
+use std::sync::Arc;
+
+#[derive(Clone)]
+pub struct FFT {
+    // The fft engine
+    fft: Arc<dyn RealToComplex<Flt>>,
+    // Copy over time data, as it is used as scratch data in the fft engine
+    timescratch: Vec<Flt>,
+    // rounded down nfft/2
+    half_nfft_rounded: usize,
+    nfftF: Flt,
+}
+
+impl FFT {
+    /// Create new FFT from given nfft
+    pub fn newFromNFFT(nfft: usize) -> FFT {
+        let mut planner = RealFftPlanner::<Flt>::new();
+        let fft = planner.plan_fft_forward(nfft);
+        Self::new(fft)
+    }
+    /// Create new fft engine from given fft engine
+    pub fn new(fft: Arc<dyn RealToComplex<Flt>>) -> FFT {
+        let nfft = fft.len();
+        let timescratch = vec![0.; nfft];
+        FFT {
+            fft,
+            timescratch,
+            half_nfft_rounded: nfft / 2,
+            nfftF: nfft as Flt,
+        }
+    }
+    pub fn process<'a, T>(&mut self, time: &[Flt], freq: T)
+    where
+        T: Into<ArrayViewMut<'a, Cflt, Ix1>>,
+    {
+        let mut freq = freq.into();
+        self.timescratch.copy_from_slice(time);
+        let _ = self
+            .fft
+            .process(&mut self.timescratch, freq.as_slice_mut().unwrap());
+
+        freq[0] /= self.nfftF;
+        freq[self.half_nfft_rounded] /= self.nfftF;
+        freq.slice_mut(s![1..self.half_nfft_rounded])
+            .par_mapv_inplace(|x| 2. * x / self.nfftF);
+    }
+}

src/ps/mod.rs

@@ -0,0 +1,4 @@
+//! Power spectra, averaged power spectra, etc. This module contains several 
+mod window;
+mod ps;
+mod fft;

src/ps/ps.rs

@@ -0,0 +1,302 @@
+use crate::config::*;
+use ndarray::parallel::prelude::*;
+use num::pow::Pow;
+use reinterpret::reinterpret_slice;
+use std::sync::Arc;
+use std::usize;
+
+use crate::Dcol;
+
+use super::fft::FFT;
+use super::window::*;
+use std::mem::MaybeUninit;
+
+use realfft::{RealFftPlanner, RealToComplex};
+
+/// Singlesided cross-Power spectra computation engine.
+struct PowerSpectra {
+    // Window used in estimator
+    pub window: Window,
+    // The window power, is corrected for in power spectra estimants
+    pub sqrt_win_pwr: Flt,
+
+    ffts: Vec<FFT>,
+
+    // Time-data buffer used for multiplying signals with Window
+    timedata: Array2<Flt>,
+    // Frequency domain buffer used for storage of signal FFt's in inbetween stage
+    freqdata: Array2<Cflt>,
+}
+
+impl PowerSpectra {
+    /// Return the FFT length used in power spectra computations
+    pub fn nfft(&self) -> usize {
+        self.window.win.len()
+    }
+    /// Create new power spectra estimator. Uses FFT size from window length
+    ///
+    /// # Panics
+    ///
+    /// - If win.len() != nfft
+    /// - if nfft == 0
+    pub fn newFromWindow(window: Window) -> PowerSpectra {
+        let nfft = window.win.len();
+        let win_pwr = window.win.mapv(|w| w.powi(2)).sum()/(nfft as Flt);
+        assert!(nfft > 0);
+        assert!(nfft % 2 == 0);
+
+        let mut planner = RealFftPlanner::<Flt>::new();
+        let fft = planner.plan_fft_forward(nfft);
+
+        let Fft = FFT::new(fft);
+
+        PowerSpectra {
+            window,
+            sqrt_win_pwr: Flt::sqrt(win_pwr),
+            ffts: vec![Fft],
+            timedata: Array2::zeros((nfft, 1)),
+            freqdata: Array2::zeros((nfft / 2 + 1, 1)),
+        }
+    }
+
+    // Compute FFTs of input channel data.
+    fn compute_ffts(&mut self, timedata: ArrayView2<Flt>) -> &Array2<Cflt> {
+        let (n, nch) = timedata.dim();
+        let nfft = self.nfft();
+        assert!(n == nfft);
+
+        // Make sure enough fft engines are available
+        while nch > self.ffts.len() {
+            self.ffts.push(self.ffts.last().unwrap().clone());
+            self.freqdata
+                .push_column(Ccol::from_vec(vec![Cflt::new(0., 0.); nfft / 2 + 1]).view())
+                .unwrap();
+            self.timedata
+                .push_column(Dcol::zeros(nfft).view())
+                .unwrap();
+        }
+
+        assert!(n == self.nfft());
+        assert!(n == self.window.win.len());
+        let sqrt_win_pwr = self.sqrt_win_pwr;
+
+        // Multiply signals with window function, and compute fft's for each channel
+        Zip::from(timedata.axis_iter(Axis(1)))
+            .and(self.timedata.axis_iter_mut(Axis(1)))
+            .and(&mut self.ffts)
+            .and(self.freqdata.axis_iter_mut(Axis(1)))
+            .par_for_each(|time_in,mut time, fft, mut freq| {
+
+                // Multiply with window and copy over to local time data buffer
+                azip!((t in &mut time, &tin in time_in, &win in &self.window.win) *t=tin*win/sqrt_win_pwr);
+
+                let tslice = time.as_slice().unwrap();
+                let fslice = freq.as_slice_mut().unwrap();
+                fft.process(tslice, fslice);
+            });
+
+        &self.freqdata
+    }
+
+    /// Compute cross power spectra from input time data. First axis is
+    /// frequency, second axis is channel i, third axis is channel j.
+    pub fn compute<'a, T>(&mut self, tdata: T) -> Array3<Cflt>
+    where
+        T: Into<ArrayView<'a, Flt, Ix2>>,
+    {
+        let tdata = tdata.into();
+        let clen = self.nfft() / 2 + 1;
+        let nchannel = tdata.ncols();
+        let win_pwr = self.sqrt_win_pwr;
+
+        // Compute fft of input data, and store in self.freqdata
+        let fd = self.compute_ffts(tdata);
+        let fdconj = fd.mapv(|c| c.conj());
+
+        let result = Array3::uninit((clen, nchannel, nchannel));
+        let mut result: Array3<Cflt> = unsafe { result.assume_init() };
+
+        // Loop over result axis one and channel i IN PARALLEL
+        Zip::from(result.axis_iter_mut(Axis(1)))
+            .and(fd.axis_iter(Axis(1)))
+            .par_for_each(|mut out, chi| {
+                // out: channel i of output 3D array, channel j all
+                // chi: channel i 
+                Zip::from(out.axis_iter_mut(Axis(1)))
+                    .and(fdconj.axis_iter(Axis(1)))
+                    .for_each(|mut out, chj| {
+                        // out: channel i, j
+                        // chj: channel j conjugated
+                        Zip::from(&mut out)
+                            .and(chi)
+                            .and(chj)
+                            .for_each(|out, chi, chjc| 
+                                // Loop over frequency components
+                                *out = 0.5 * chi * chjc
+                            );
+                        
+                        // The DC component has no 0.5 correction, as it only
+                        // occurs ones in a (double-sided) power spectrum. So
+                        // here we undo the 0.5 of 4 lines above here.
+                        out[0] *= 2.;
+                        out[clen-1] *= 2.;
+
+                    });
+            });
+        result
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use approx::{abs_diff_eq, assert_relative_eq, assert_ulps_eq, ulps_eq};
+    // For absolute value
+    use num::complex::ComplexFloat;
+    use rand_distr::StandardNormal;
+
+    /// Generate a sine wave at the order i
+    fn generate_sinewave(nfft: usize,order: usize) -> Dcol {
+        Dcol::from_iter((0..nfft).map(|i| {
+            Flt::sin(i as Flt/(nfft) as Flt * order as Flt * 2.*pi)
+        }))
+    }
+    /// Generate a sine wave at the order i
+    fn generate_cosinewave(nfft: usize,order: usize) -> Dcol {
+        Dcol::from_iter((0..nfft).map(|i| {
+            Flt::cos(i as Flt/(nfft) as Flt * order as Flt * 2.*pi)
+        }))
+    }
+
+    use super::*;
+    #[test]
+    /// Test whether DC part of single-sided FFT has right properties
+    fn test_fft_DC() {
+        const nfft: usize = 10;
+        let rect = Window::new(WindowType::Rect, nfft);
+        let mut ps = PowerSpectra::newFromWindow(rect);
+
+        let td = Dmat::ones((nfft, 1));
+
+        let fd = ps.compute_ffts(td.view());
+        // println!("{:?}", fd);
+        assert_relative_eq!(fd[(0, 0)].re, 1.);
+        assert_relative_eq!(fd[(0, 0)].im, 0.);
+        let abs_fneq0 = fd.slice(s![1.., 0]).sum();
+        assert_relative_eq!(abs_fneq0.re, 0.);
+        assert_relative_eq!(abs_fneq0.im, 0.);
+    }
+
+    /// Test whether AC part of single-sided FFT has right properties
+    #[test]
+    fn test_fft_AC() {
+        const nfft: usize = 256;
+        let rect = Window::new(WindowType::Rect, nfft);
+        let mut ps = PowerSpectra::newFromWindow(rect);
+
+        // Start with a time signal
+        let mut t: Dmat = Dmat::default((nfft, 0));
+        t.push_column(generate_sinewave(nfft,1).view())
+            .unwrap();
+        // println!("{:?}", t);
+
+        let fd = ps.compute_ffts(t.view());
+        // println!("{:?}", fd);
+        assert_relative_eq!(fd[(0, 0)].re, 0., epsilon = Flt::EPSILON * nfft as Flt);
+        assert_relative_eq!(fd[(0, 0)].im, 0., epsilon = Flt::EPSILON * nfft as Flt);
+
+        assert_relative_eq!(fd[(1, 0)].re, 0., epsilon = Flt::EPSILON * nfft as Flt);
+        assert_ulps_eq!(fd[(1, 0)].im, -1., epsilon = Flt::EPSILON * nfft as Flt);
+
+        // Sum of all terms at frequency index 2 to ...
+        let sum_higher_freqs_abs = Cflt::abs(fd.slice(s![2.., 0]).sum());
+        assert_ulps_eq!(
+            sum_higher_freqs_abs,
+            0.,
+            epsilon = Flt::EPSILON * nfft as Flt
+        );
+    }
+
+
+    /// Thest whether power spectra scale properly. Signals with amplitude of 1
+    /// should come back with a power of 0.5. DC offsets should come in as
+    /// value^2 at frequency index 0.
+    #[test]
+    fn test_ps_scale() {
+
+        const nfft: usize = 124;
+        let rect = Window::new(WindowType::Rect, nfft);
+        let mut ps = PowerSpectra::newFromWindow(rect);
+
+        // Start with a time signal
+        let mut t: Dmat = Dmat::default((nfft, 0));
+        t.push_column(generate_cosinewave(nfft,1).view())
+            .unwrap();
+        let dc_component = 0.25;
+        let dc_power = dc_component.pow(2);
+        t.mapv_inplace(|t| t + dc_component);
+
+        let power = ps.compute(t.view());
+        assert_relative_eq!(power[(0, 0,0)].re, dc_power, epsilon = Flt::EPSILON * nfft as Flt);
+        assert_relative_eq!(power[(1, 0,0)].re, 0.5, epsilon = Flt::EPSILON * nfft as Flt);
+        assert_relative_eq!(power[(1, 0,0)].im, 0.0, epsilon = Flt::EPSILON * nfft as Flt);
+
+    }
+
+    use ndarray_rand::RandomExt;
+    // Test parseval's theorem for some random data
+    #[test]
+    fn test_parseval() {
+
+        const nfft: usize = 512;
+        let rect = Window::new(WindowType::Rect, nfft);
+        let mut ps = PowerSpectra::newFromWindow(rect);
+
+        // Start with a time signal
+        let t: Dmat = Dmat::random((nfft, 1), StandardNormal);
+
+        let tavg = t.sum()/(nfft as Flt);
+        let t_dc_power = tavg.powi(2);
+        // println!("dc power in time domain: {:?}", t_dc_power);
+
+        let signal_pwr = t.mapv(|t| t.powi(2)).sum()/(nfft as Flt);
+        // println!("Total signal power in time domain: {:?} ", signal_pwr);
+
+        let power = ps.compute(t.view());
+        // println!("freq domain power: {:?}", power);
+
+        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));
+
+    }
+
+    // Test parseval's theorem for some random data
+    #[test]
+    fn test_parseval_with_window() {
+
+        const nfft: usize = 48000;
+        let rect = Window::new(WindowType::Hann, nfft);
+        let mut ps = PowerSpectra::newFromWindow(rect);
+
+        // Start with a time signal
+        let t: Dmat = Dmat::random((nfft, 1), StandardNormal);
+
+        let tavg = t.sum()/(nfft as Flt);
+        let t_dc_power = tavg.powi(2);
+        // println!("dc power in time domain: {:?}", t_dc_power);
+
+        let signal_pwr = t.mapv(|t| t.powi(2)).sum()/(nfft as Flt);
+        // println!("Total signal power in time domain: {:?} ", signal_pwr);
+
+        let power = ps.compute(t.view());
+        // println!("freq domain power: {:?}", power);
+
+        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));
+
+    }
+
+}

src/ps/window.rs

@@ -0,0 +1,168 @@
+//! 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.
+#![allow(non_snake_case)]
+use crate::config::*;
+
+#[macro_use]
+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) {
+    let nfftF = nfft as Flt;
+    (
+        // The Window
+        linspace(nfft).mapv(|i| (pi * i / (nfftF+1.)).sin().powi(2)),
+        // The hopp size
+        (nfft) / 2,
+    )
+}
+fn rect(nfft: usize) -> (Dcol, usize) {
+    (Dcol::ones(nfft), nfft)
+}
+fn blackman(N: usize) -> (Dcol, usize) {
+    let a0 = 7938. / 18608.;
+    let a1 = 9240. / 18608.;
+    let a2 = 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,
+    )
+}
+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 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,
+    )
+}
+
+/// Window type descriptors. Used for storage
+#[derive(Display, Clone, Debug)]
+pub enum WindowType {
+    /// Von Hann window
+    Hann = 0,
+    /// Hamming window
+    Hamming = 1,
+    /// Boxcar / rectangular window
+    Rect = 2,
+    /// Bartlett window
+    Bartlett = 3,
+    /// Blackman window
+    Blackman = 4,
+}
+
+/// Window (taper) computed from specified window type.
+#[derive(Clone)]
+pub struct Window {
+    /// The enum from which it is generated
+    pub w: WindowType,
+    /// The actual window computed from specified nfft
+    pub win: Dcol,
+    /// The 'optimal' number of samples of shift per window average (hop size).
+    pub R: usize,
+}
+impl Window {
+    /// Create new window based on type and fft length. FFT length should be even. The (Window)
+    /// struct contains type and generated window in the `win` member.
+    ///
+    /// # Panics
+    ///
+    /// If nfft %2 != 0
+    pub fn new(w: WindowType, nfft: usize) -> Window {
+        if nfft % 2 != 0 {
+            panic!("Requires even nfft");
+        }
+        let (win, R) = match w {
+            WindowType::Hann => hann(nfft),
+            WindowType::Hamming => hamming(nfft),
+            WindowType::Rect => rect(nfft),
+            WindowType::Bartlett => bartlett(nfft),
+            WindowType::Blackman => blackman(nfft),
+        };
+        Window { w, win, R }
+    }
+    /// Convenience function that returns the length of the window.
+    pub fn len(&self) -> usize {
+        self.win.len()
+    }
+}
+//
+#[cfg(test)]
+mod test {
+    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);
+    }
+}

src/siggen.rs

@@ -27,6 +27,7 @@ use rand::prelude::*;
 use rand::rngs::ThreadRng;
 use rand_distr::StandardNormal;
 
+/// Ratio between circumference and radius of a circle
 const twopi: Flt = 2.0 * pi;
 
 /// Source for the signal generator. Implementations are sine waves, sweeps, noise.
@@ -148,6 +149,7 @@ pub struct Siggen {
     // Output buffers (for filtered source signal)
     chout_buf: Vec<Vec<Flt>>,
 }
+#[cfg(feature="python-bindings")]
 #[cfg_attr(feature = "python-bindings", pymethods)]
 impl Siggen {
     #[pyo3(name = "newWhiteNoise")]