diff --git a/Cargo.toml b/Cargo.toml index 329b26e..86cd6f7 100644 --- a/Cargo.toml +++ b/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"] diff --git a/src/config.rs b/src/config.rs index fe3dad4..13f7f94 100644 --- a/src/config.rs +++ b/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}; } } diff --git a/src/filter/biquad.rs b/src/filter/biquad.rs index 7805014..00ce934 100644 --- a/src/filter/biquad.rs +++ b/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 - 1., // a0 + 0., // b2, 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 { + 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 { @@ -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.)); + } } diff --git a/src/lib.rs b/src/lib.rs index 8cbc788..da2a130 100644 --- a/src/lib.rs +++ b/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::*; diff --git a/src/ps/aps.rs b/src/ps/aps.rs new file mode 100644 index 0000000..0fb847c --- /dev/null +++ b/src/ps/aps.rs @@ -0,0 +1 @@ +//! Averaged power spectra module \ No newline at end of file diff --git a/src/ps/fft.rs b/src/ps/fft.rs new file mode 100644 index 0000000..4bda503 --- /dev/null +++ b/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>, + // Copy over time data, as it is used as scratch data in the fft engine + timescratch: Vec, + // 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::::new(); + let fft = planner.plan_fft_forward(nfft); + Self::new(fft) + } + /// Create new fft engine from given fft engine + pub fn new(fft: Arc>) -> 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>, + { + 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); + } +} \ No newline at end of file diff --git a/src/ps/mod.rs b/src/ps/mod.rs new file mode 100644 index 0000000..f352071 --- /dev/null +++ b/src/ps/mod.rs @@ -0,0 +1,4 @@ +//! Power spectra, averaged power spectra, etc. This module contains several +mod window; +mod ps; +mod fft; \ No newline at end of file diff --git a/src/ps/ps.rs b/src/ps/ps.rs new file mode 100644 index 0000000..c4f6401 --- /dev/null +++ b/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, + + // Time-data buffer used for multiplying signals with Window + timedata: Array2, + // Frequency domain buffer used for storage of signal FFt's in inbetween stage + freqdata: Array2, +} + +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::::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) -> &Array2 { + 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 + where + T: Into>, + { + 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 = 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)); + + } + +} diff --git a/src/ps/window.rs b/src/ps/window.rs new file mode 100644 index 0000000..59216c1 --- /dev/null +++ b/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); + } +} diff --git a/src/siggen.rs b/src/siggen.rs index 92d2b8c..bef9ce1 100644 --- a/src/siggen.rs +++ b/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>, } +#[cfg(feature="python-bindings")] #[cfg_attr(feature = "python-bindings", pymethods)] impl Siggen { #[pyo3(name = "newWhiteNoise")] diff --git a/src/timebuffer.rs b/src/timebuffer.rs new file mode 100644 index 0000000..e69de29