From 378daa2228c17429fb02e85a2e25dcb5fd4a1344 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 16 Oct 2024 15:00:44 +0800 Subject: [PATCH 01/27] Add kaiser_atten --- sci-rs/src/signal/filter/design/kaiser.rs | 46 +++++++++++++++++++++++ sci-rs/src/signal/filter/design/mod.rs | 2 + 2 files changed, 48 insertions(+) create mode 100644 sci-rs/src/signal/filter/design/kaiser.rs diff --git a/sci-rs/src/signal/filter/design/kaiser.rs b/sci-rs/src/signal/filter/design/kaiser.rs new file mode 100644 index 00000000..eb532beb --- /dev/null +++ b/sci-rs/src/signal/filter/design/kaiser.rs @@ -0,0 +1,46 @@ +use core::f64::consts::PI; +use core::ops::Mul; +use num_traits::{real::Real, MulAdd, Pow, ToPrimitive}; + +/// Compute the attenuation of a Kaiser FIR filter. +/// +/// Given the number of taps `N` and the transition width `width`, compute the +/// attenuation `a` in dB, given by Kaiser's formula: +/// ```custom +/// a = 2.285 * (N - 1) * pi * width + 7.95 +/// ``` +/// +/// # Parameters +/// * `numtaps`: int +/// The number of taps in the FIR filter. +/// * `width`: float +/// The desired width of the transition region between passband and +/// stopband (or, in general, at any discontinuity) for the filter, +/// expressed as a fraction of the Nyquist frequency. +/// +/// # Returns +/// `a`: The attenuation of the ripple, in dB. +/// +/// # Examples +/// Suppose we want to design a FIR filter using the Kaiser window method +/// that will have 211 taps and a transition width of 9 Hz for a signal that +/// is sampled at 480 Hz. Expressed as a fraction of the Nyquist frequency, +/// the width is 9/(0.5*480) = 0.0375. The approximate attenuation (in dB) +/// is computed as follows: +/// ``` +/// use sci_rs::signal::filter::design::kaiser_atten; +/// assert_eq!(64.48099630593983 , kaiser_atten(211, 0.0375)); +/// ``` +/// +/// # See Also +/// [kaiserord], [kaiser_beta] +pub fn kaiser_atten(numtaps: u32, width: F) -> F +where + F: Real + MulAdd, +{ + MulAdd::mul_add( + width, + F::from(numtaps - 1).unwrap() * F::from(2.285 * PI).unwrap(), + F::from(7.95).unwrap(), + ) +} diff --git a/sci-rs/src/signal/filter/design/mod.rs b/sci-rs/src/signal/filter/design/mod.rs index 9079bec7..5cdfebd5 100644 --- a/sci-rs/src/signal/filter/design/mod.rs +++ b/sci-rs/src/signal/filter/design/mod.rs @@ -4,6 +4,7 @@ mod cplx; mod filter_output; mod filter_type; mod iirfilter; +mod kaiser; mod lp2bp_zpk; mod lp2bs_zpk; mod lp2hp_zpk; @@ -19,6 +20,7 @@ use cplx::*; pub use filter_output::*; pub use filter_type::*; pub use iirfilter::*; +pub use kaiser::*; pub use lp2bp_zpk::*; pub use lp2bs_zpk::*; pub use lp2hp_zpk::*; From d2c78ec2109cb3a3e780db8601758ad67cfc01eb Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Thu, 17 Oct 2024 15:38:37 +0800 Subject: [PATCH 02/27] Add kaiser_beta --- sci-rs/src/signal/filter/design/kaiser.rs | 46 +++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/sci-rs/src/signal/filter/design/kaiser.rs b/sci-rs/src/signal/filter/design/kaiser.rs index eb532beb..d4c6b10f 100644 --- a/sci-rs/src/signal/filter/design/kaiser.rs +++ b/sci-rs/src/signal/filter/design/kaiser.rs @@ -2,6 +2,52 @@ use core::f64::consts::PI; use core::ops::Mul; use num_traits::{real::Real, MulAdd, Pow, ToPrimitive}; +/// Compute the Kaiser parameter `beta`, given the attenuation `a`. +/// +/// # Parameters +/// * `a`: float +/// The desired attenuation in the stopband and maximum ripple in the passband, in dB. This +/// should be a *positive* number. +/// +/// # Returns +/// * `beta`: float +/// The `beta` parameter to be used in the formula for a Kaiser window. +/// +/// # References +/// Oppenheim, Schafer, "Discrete-Time Signal Processing", p.475-476. +/// +/// Examples +/// -------- +/// Suppose we want to design a lowpass filter, with 65 dB attenuation +/// in the stop band. The Kaiser window parameter to be used in the +/// window method is computed by ``kaiser_beta(65.)``: +/// +/// ``` +/// use sci_rs::signal::filter::design::kaiser_beta; +/// assert_eq!(6.20426, kaiser_beta(65.)); +/// ``` +pub fn kaiser_beta(a: F) -> F +where + F: Real + MulAdd + Pow, + >::Output: MulAdd, +{ + if a > F::from(50).unwrap() { + F::from(0.1102).unwrap() * (a - F::from(8.7).unwrap()) + } else if a > F::from(21).unwrap() { + let a = a - F::from(21).unwrap(); + + MulAdd::mul_add( + a.pow(F::from(0.4).unwrap()), + F::from(0.5842).unwrap(), + F::from(0.07886).unwrap() * a, + ) + } else if a > F::from(0).unwrap() { + F::from(0).unwrap() + } else { + panic!("Expected a positive input.") + } +} + /// Compute the attenuation of a Kaiser FIR filter. /// /// Given the number of taps `N` and the transition width `width`, compute the From 7fdef5cc77db493c86c1cdf72e229da8d7fe593b Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Fri, 18 Oct 2024 17:43:47 +0800 Subject: [PATCH 03/27] Add kaiserord --- sci-rs/src/signal/filter/design/kaiser.rs | 151 ++++++++++++++++++++++ 1 file changed, 151 insertions(+) diff --git a/sci-rs/src/signal/filter/design/kaiser.rs b/sci-rs/src/signal/filter/design/kaiser.rs index d4c6b10f..67be1e45 100644 --- a/sci-rs/src/signal/filter/design/kaiser.rs +++ b/sci-rs/src/signal/filter/design/kaiser.rs @@ -90,3 +90,154 @@ where F::from(7.95).unwrap(), ) } + +/// Determine the filter window parameters for the Kaiser window method. +/// +/// The parameters returned by this function are generally used to create +/// a finite impulse response filter using the window method, with either +/// `firwin` or `firwin2`. +/// +/// # Parameters +/// * `ripple`: float +/// Upper bound for the deviation (in dB) of the magnitude of the +/// filter's frequency response from that of the desired filter (not +/// including frequencies in any transition intervals). That is, if w +/// is the frequency expressed as a fraction of the Nyquist frequency, +/// A(w) is the actual frequency response of the filter and D(w) is the +/// desired frequency response, the design requirement is that: +/// ```abs(A(w) - D(w))) < 10**(-ripple/20)``` +/// for 0 <= w <= 1 and w not in a transition interval. +/// * `width`: float +/// Width of transition region, normalized so that 1 corresponds to pi +/// radians / sample. That is, the frequency is expressed as a fraction +/// of the Nyquist frequency. +/// +/// # Returns +/// * `numtaps`: int +/// The length of the Kaiser window. +/// * `beta`: float +/// The beta parameter for the Kaiser window. +/// +/// # Notes +/// There are several ways to obtain the Kaiser window: +/// +/// - ``signal.windows.kaiser(numtaps, beta, sym=True)`` +/// - ``signal.get_window(beta, numtaps)`` +/// - ``signal.get_window(('kaiser', beta), numtaps)`` +/// +/// The empirical equations discovered by Kaiser are used. +/// +/// # References +/// Oppenheim, Schafer, "Discrete-Time Signal Processing", pp.475-476. +/// +/// # Scipy Example +/// We will use the Kaiser window method to design a lowpass FIR filter +/// for a signal that is sampled at 1000 Hz. +/// +/// We want at least 65 dB rejection in the stop band, and in the pass +/// band the gain should vary no more than 0.5%. +/// +/// We want a cutoff frequency of 175 Hz, with a transition between the +/// pass band and the stop band of 24 Hz. That is, in the band [0, 163], +/// the gain varies no more than 0.5%, and in the band [187, 500], the +/// signal is attenuated by at least 65 dB. +/// +/// ```custom,{class=language-python} +/// >>> import numpy as np +/// >>> from scipy.signal import kaiserord, firwin, freqz +/// >>> import matplotlib.pyplot as plt +/// >>> fs = 1000.0 +/// >>> cutoff = 175 +/// >>> width = 24 +/// ``` +/// +/// The Kaiser method accepts just a single parameter to control the pass +/// band ripple and the stop band rejection, so we use the more restrictive +/// of the two. In this case, the pass band ripple is 0.005, or 46.02 dB, +/// so we will use 65 dB as the design parameter. +/// +/// Use `kaiserord` to determine the length of the filter and the +/// parameter for the Kaiser window. +/// +/// ```custom,{class=language-python} +/// >>> numtaps, beta = kaiserord(65, width/(0.5*fs)) +/// >>> numtaps +/// 167 +/// >>> beta +/// 6.20426 +/// ``` +/// +/// Use `firwin` to create the FIR filter. +/// +/// ```custom,{class=language-python} +/// >>> taps = firwin(numtaps, cutoff, window=('kaiser', beta), +/// ... scale=False, fs=fs) +/// ``` +/// +/// Compute the frequency response of the filter. ``w`` is the array of +/// frequencies, and ``h`` is the corresponding complex array of frequency +/// responses. +/// +/// ```custom,{class=language-python} +/// >>> w, h = freqz(taps, worN=8000) +/// >>> w *= 0.5*fs/np.pi # Convert w to Hz. +/// ``` +/// +/// Compute the deviation of the magnitude of the filter's response from +/// that of the ideal lowpass filter. Values in the transition region are +/// set to ``nan``, so they won't appear in the plot. +/// +/// ```custom,{class=language-python} +/// >>> ideal = w < cutoff # The "ideal" frequency response. +/// >>> deviation = np.abs(np.abs(h) - ideal) +/// >>> deviation[(w > cutoff - 0.5*width) & (w < cutoff + 0.5*width)] = np.nan +/// ``` +/// +/// Plot the deviation. A close look at the left end of the stop band shows +/// that the requirement for 65 dB attenuation is violated in the first lobe +/// by about 0.125 dB. This is not unusual for the Kaiser window method. +/// +/// ```custom,{class=language-python} +/// >>> plt.plot(w, 20*np.log10(np.abs(deviation))) +/// >>> plt.xlim(0, 0.5*fs) +/// >>> plt.ylim(-90, -60) +/// >>> plt.grid(alpha=0.25) +/// >>> plt.axhline(-65, color='r', ls='--', alpha=0.3) +/// >>> plt.xlabel('Frequency (Hz)') +/// >>> plt.ylabel('Deviation from ideal (dB)') +/// >>> plt.title('Lowpass Filter Frequency Response') +/// >>> plt.show() +/// ``` +/// +/// See Also +/// -------- +/// [kaiser_beta], [kaiser_atten] +/// +pub fn kaiserord(ripple: F, width: F) -> (F, F) +where + F: Real + MulAdd + Pow, +{ + let a = ripple.abs(); + if a < F::from(8).unwrap() { + panic!("Requested maximum ripple attenuation is too small for the Kaiser formula."); + } + let beta = kaiser_beta(a); + let numtaps = + F::from(1).unwrap() + (a - F::from(7.95).unwrap()) / (F::from(2.285 * PI).unwrap() * width); + (numtaps.ceil(), beta) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn kaiserord_test_a() { + let fs = 1000.; + let cutoff = 175.; + let width = 24.; + let ripple = 65.; + + assert_eq!((167., 6.20426), kaiserord(ripple, width / (0.5 * fs))); + } +} From 0025daf2e38787c702464ab74ca49e36ef302e96 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 22 Oct 2024 22:23:05 +0800 Subject: [PATCH 04/27] Add signal::windows namespace and helper functions To implement firwin, there is a need to get_windows on all window types defined in the Scipy signal.windows namespace, most of which are implemented in a _windows.py file. This commit keeps all the window variants that needs to be implemented as comments in an Enum for future implementation, along with some functions that are used across the implementation of all window types. The intended strategy to combat the variadic function type associated to the Scipy's get_windows function is to enclose all function arguments in a struct that represents the given function type, while ensuring all structs representing the Windows type implement a common trait. --- sci-rs/src/signal/mod.rs | 10 ++++ sci-rs/src/signal/windows/mod.rs | 92 ++++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+) create mode 100644 sci-rs/src/signal/windows/mod.rs diff --git a/sci-rs/src/signal/mod.rs b/sci-rs/src/signal/mod.rs index 8aede037..ce9571cd 100644 --- a/sci-rs/src/signal/mod.rs +++ b/sci-rs/src/signal/mod.rs @@ -8,6 +8,16 @@ pub mod wave; #[cfg(feature = "std")] pub mod convolve; +/// Window functions +/// This contains all window functions in the +/// [scipy.signal.windows](https://docs.scipy.org/doc/scipy/reference/signal.windows.html#module-scipy.signal.windows) +/// namespace. +/// The convenience function +/// [`get_windows`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window) +/// in the [scipy.signal](https://docs.scipy.org/doc/scipy/reference/signal.html#window-functions) +/// namespace is located here. +pub mod windows; + /// Signal Resampling #[cfg(feature = "std")] pub mod resample; diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs new file mode 100644 index 00000000..d1452091 --- /dev/null +++ b/sci-rs/src/signal/windows/mod.rs @@ -0,0 +1,92 @@ +#[cfg(feature = "alloc")] +use alloc::vec::Vec; +use num_traits::{real::Real, Float}; + +/// Corresponding window representation for tuple-structs of [Window] variants. +#[cfg(feature = "alloc")] +pub trait GetWindow +where + W: Real, +{ + /// Returns a window of given length and type. + /// + /// # Parameters + /// `self`: + /// The type of window to construct, typically consists of at least the following + /// arguments: + /// * `Nx`: usize + /// The number of samples in the window + /// * `fftbins/~sym`: bool + /// If fftbins=true/sym=false, create a “periodic” window, ready to use with ifftshift and + /// be multiplied by the result of an FFT. This is the default behaviour in scipy. + /// If fftbins=false/sym=true, create a "symmetric" window, for use in filter design. + /// * `*args`: + /// Other arguments relevant to the window type. + /// + /// # Reference + /// + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec; +} + +/// Private function for windows implementing [GetWindow] +/// Handle small or incorrect window lengths. +#[inline(always)] +fn len_guard(m: usize) -> bool { + m <= 1 +} + +/// Private function for windows implementing [GetWindow] +/// Extend window by 1 sample if needed for DFT-even symmetry. +#[inline(always)] +fn extend(m: usize, sym: bool) -> (usize, bool) { + if !sym { + (m + 1, true) + } else { + (m, false) + } +} + +/// Private function for windows implementing [GetWindow] +/// Truncate window by 1 sample if needed for DFT-even symmetry. +#[inline(always)] +fn truncate(mut w: Vec, needed: bool) -> Vec { + if needed { + w.pop(); + } + w +} + +/// todo +// Ordering is as in accordance with +// https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html. +#[derive(Debug, Clone, PartialEq)] +pub enum Window { + // Boxcar, + // Triangle, + // Blackman, + // Hamming, + // Hann, + // Bartlett, + // Flattop, + // Parzen, + // Bohman, + // BlackmanHarris, + // Nuttall, + // BartHann, + // Cosine, + // Exponential, + // Tukey, + // Taylor, + // Lanczos, + // Kaiser Window. + // Kaiser, // Needs Beta + // KaiserBesselDerived, // Needs Beta + // Gaussian, // Needs Standard Deviation + // Generic weighted sum of cosine term windows. + // GenericCosine(GenericCosine), // Needs Weighting Coefficients + // GeneralGaussian, // Needs Power, Width + // GeneralHamming, // Needs Window Coefficients. + // Dpss, // Needs Normalized Half-Bandwidth. + // Chebwin, // Needs Attenuation. +} From 1a2feb6b51aa1675d86f12a6f45ccf3884645b8d Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 23 Oct 2024 18:54:36 +0800 Subject: [PATCH 05/27] Add Boxcar window type --- sci-rs/src/signal/windows/boxcar.rs | 83 +++++++++++++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 9 +++- 2 files changed, 91 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/boxcar.rs diff --git a/sci-rs/src/signal/windows/boxcar.rs b/sci-rs/src/signal/windows/boxcar.rs new file mode 100644 index 00000000..0410fc90 --- /dev/null +++ b/sci-rs/src/signal/windows/boxcar.rs @@ -0,0 +1,83 @@ +use super::{extend, len_guard}; +use num_traits::real::Real; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `Boxcar` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct Boxcar { + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Whether the window is symmetric. (Has no effect for boxcar.) + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Boxcar { + /// Returns a Boxcar struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, sym: bool) -> Self { + Boxcar { m, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for Boxcar +where + W: Real, +{ + /// Return a window of type: boxcar or rectangular window. + /// + /// Also known as a rectangular window or Dirichlet window, this is equivalent to no window at + /// all. + /// + /// # Parameters + /// `self`: [Boxcar] + /// + /// # Example + /// ``` + /// use sci_rs::signal::windows::{Boxcar, GetWindow}; + /// + /// let nx = 5; + /// let boxcar = Boxcar::new(nx, false); + /// let window: Vec = boxcar.get_window(); + /// assert_eq!(vec![1.; nx], window); + /// ``` + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + if len_guard(self.m) { + return Vec::::new(); + } + let (m, needs_trunc) = extend(self.m, self.sym); + + if !needs_trunc { + vec![W::from(1).unwrap(); m] + } else { + vec![W::from(1).unwrap(); m - 1] + } + } +} + +#[cfg(test)] +mod test { + use super::*; + + #[test] + fn case_a() { + let nx = 5; + let boxcar = Boxcar::new(nx, false); + let window: Vec = boxcar.get_window(); + assert_eq!(vec![1.; nx], window); + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index d1452091..3983130f 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -57,12 +57,19 @@ fn truncate(mut w: Vec, needed: bool) -> Vec { w } +mod boxcar; +pub use boxcar::Boxcar; + /// todo // Ordering is as in accordance with // https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html. #[derive(Debug, Clone, PartialEq)] +// Is it possible for the enums which wraps the structs to only require the generic that the struct +// implements GetWindow? pub enum Window { - // Boxcar, + /// [Boxcar] window, also known as a rectangular window or Dirichlet window; This is equivalent + /// to no window at all. + Boxcar(Boxcar), // Triangle, // Blackman, // Hamming, From e860413a869ddaf0257651515a732e44a33e8f02 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Fri, 25 Oct 2024 18:09:50 +0800 Subject: [PATCH 06/27] Add triangle window --- sci-rs/src/signal/windows/mod.rs | 5 +- sci-rs/src/signal/windows/triangle.rs | 202 ++++++++++++++++++++++++++ 2 files changed, 206 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/triangle.rs diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 3983130f..fc1393dc 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -58,7 +58,9 @@ fn truncate(mut w: Vec, needed: bool) -> Vec { } mod boxcar; +mod triangle; pub use boxcar::Boxcar; +pub use triangle::Triangle; /// todo // Ordering is as in accordance with @@ -70,7 +72,8 @@ pub enum Window { /// [Boxcar] window, also known as a rectangular window or Dirichlet window; This is equivalent /// to no window at all. Boxcar(Boxcar), - // Triangle, + /// [Triangle] window. + Triangle(Triangle), // Blackman, // Hamming, // Hann, diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs new file mode 100644 index 00000000..a832fb72 --- /dev/null +++ b/sci-rs/src/signal/windows/triangle.rs @@ -0,0 +1,202 @@ +use super::{extend, len_guard, truncate}; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `Triangle` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct Triangle { + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Triangle { + /// Returns a Triangle struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, sym: bool) -> Self { + Triangle { m, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for Triangle +where + W: Real, +{ + /// Return a window of type: Triangle. + /// + /// This is the not the same as the Bartlett window. + /// + /// # Parameters + /// `self`: [Triangle] + /// + /// # Returns + /// `w`: `vec` + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # Example + /// ``` + /// use sci_rs::signal::windows::{Triangle, GetWindow}; + /// + /// let nx = 8; + /// let tri = Triangle::new(nx, true); + /// assert_eq!(vec![0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375, 0.125], tri.get_window()); + /// + /// let nx = 9; + /// let tri = Triangle::new(nx, false); + /// assert_eq!(vec![0.1, 0.3, 0.5, 0.7, 0.9, 0.9, 0.7, 0.5, 0.3], tri.get_window()); + /// ``` + /// + /// # References + /// + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + if len_guard(self.m) { + return Vec::::new(); + } + let (m, needs_trunc) = extend(self.m, self.sym); + + let mut n = (1..=((m + 1) / 2)).map(|x| W::from(x).unwrap()); + let m_f: W = W::from(m).unwrap(); + let w: Vec = match m % 2 { + 0 => { + let mut w: Vec = n + .map(|n| (W::from(2).unwrap() * n - W::from(1).unwrap()) / m_f) + .collect(); + w.extend(w.clone().iter().rev()); + w + } + 1 => { + let mut w: Vec = n + .map(|n| W::from(2).unwrap() * n / (m_f + W::from(1).unwrap())) + .collect(); + w.extend(w.clone().iter().rev().skip(1)); + w + } + _ => panic!(), + }; + + truncate(w, needs_trunc) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_relative_eq; + + #[test] + #[cfg(feature = "alloc")] + fn case_even_true() { + // from scipy.signal.windows import triangle + // triangle(n) + + let upper = 1_000; + for i in 0..upper { + let nx = 2 * i; + let tri = Triangle::new(nx, true); + let expected: Vec = (0..nx) + .into_iter() + .chain((0..nx).rev()) + .filter(|n| n % 2 == 1) + .map(|n| n as f64 / nx as f64) + .collect(); + let result: Vec = tri.get_window(); + + assert_eq!(expected, result); + // for (&e, t) in expected.iter().zip(tri.get_window::()) { + // assert_relative_eq!(e, t, max_relative = 1e-7) + // } + } + } + + #[test] + #[cfg(feature = "alloc")] + fn case_even_false() { + // from scipy.signal import get_window + // get_window('triangle', 4) + + let upper = 1_000; + for i in 0..upper { + let nx = 2 * i; + let tri = Triangle::new(nx, false); + let expected: Vec = (1..=(nx / 2) + 1) + .into_iter() + .chain((1..(nx / 2) + 1).rev()) + .take(nx) + .map(|n| n as f64 / ((nx / 2) + 1) as f64) + .collect(); + + let result: Vec = tri.get_window(); + assert_eq!(expected, result); + // for (&e, t) in expected.iter().zip(tri.get_window::()) { + // assert_relative_eq!(e, t, max_relative = 1e-7) + // } + } + } + + #[test] + #[cfg(feature = "alloc")] + fn case_odd_true() { + // from scipy.signal.windows import triangle + // triangle(nx) + + let upper = 1_000; + for i in 1..upper { + let nx = 2 * i + 1; + let tri = Triangle::new(nx, true); + let expected: Vec<_> = (1..=(nx + 1) / 2) + .into_iter() + .chain((1..=(nx + 1) / 2).rev().skip(1)) + .map(|n| (n as f64) / ((nx + 1) as f64 / 2.)) + .collect(); + + let result: Vec = tri.get_window(); + assert_eq!(expected, result); + // for (&e, t) in expected.iter().zip(tri.get_window::()) { + // assert_relative_eq!(e, t); + // } + } + } + + #[test] + #[cfg(feature = "alloc")] + fn case_odd_false() { + // from scipy.signal import get_window + // get_window('triangle', 5) + + let upper = 1_000; + for i in 1..upper { + let nx = 2 * i + 1; + let tri = Triangle::new(nx, false); + let expected: Vec<_> = (0..=nx) + .into_iter() + .chain((0..=nx).into_iter().rev()) + .filter(|n| n % 2 == 1) + .take(nx) + .map(|n| n as f64 / (nx as f64 + 1.)) + .collect(); + + let result: Vec = tri.get_window(); + assert_eq!(expected, result); + // for (&e, t) in expected.iter().zip(tri.get_window::()) { + // assert_relative_eq!(e, t); + // } + } + } +} From 87384f4692e065afd7e4b2a5ff9a1483f484c55b Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 19 Nov 2024 21:18:32 +0800 Subject: [PATCH 07/27] Add general cosine window --- sci-rs/src/signal/windows/general_cosine.rs | 213 ++++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 17 +- 2 files changed, 226 insertions(+), 4 deletions(-) create mode 100644 sci-rs/src/signal/windows/general_cosine.rs diff --git a/sci-rs/src/signal/windows/general_cosine.rs b/sci-rs/src/signal/windows/general_cosine.rs new file mode 100644 index 00000000..3ac09dcc --- /dev/null +++ b/sci-rs/src/signal/windows/general_cosine.rs @@ -0,0 +1,213 @@ +use super::{extend, len_guard, truncate}; +use nalgebra::RealField; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `GeneralCosine` for use in [GetWindow]. +#[cfg(feature = "alloc")] +#[derive(Debug, Clone, PartialEq)] +pub struct GeneralCosine +where + F: Real, +{ + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Sequence of weighting coefficients. + pub a: Vec, // Is there a better type, such as impl Into? + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl GeneralCosine +where + F: Real, +{ + /// Returns a GeneralCosine struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `a`: `Vec` + /// Sequence of weighting coefficients. This uses the convention of being centered on the + /// origin, so these will typically all be positive numbers, not alternating sign. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, a: Vec, sym: bool) -> Self { + GeneralCosine { m, a, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for GeneralCosine +where + F: Real, + W: Real + Float + RealField, +{ + /// Return a window of type: GeneralCosine. + /// + /// The general cosine window is a generic weighted sum of cosine terms. + /// + /// # Parameters + /// `self`: [GeneralCosine] + /// + /// # Returns + /// `w`: `vec` + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # Example + /// We can create a `flat-top` window named "HFT90D" as following: + /// ``` + /// use sci_rs::signal::windows::{GeneralCosine, GetWindow}; + /// + /// let hfd90 = [1., 1.942604, 1.340318, 0.440811, 0.043097].into(); + /// let window: Vec = GeneralCosine::new(30, hfd90, true).get_window(); + /// ``` + /// + /// This is equivalent to the Python code: + /// ```custom,{class=language-python} + /// from scipy.signal.windows import general_cosine + /// HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097] + /// window = general_cosine(30, HFT90D, sym=False) + /// ``` + /// + /// # References + /// + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + if len_guard(self.m) { + return Vec::::new(); + } + let (m, needs_trunc) = extend(self.m, self.sym); + + let linspace = (0..m).map(|i| W::from(i).unwrap()); + let fac = linspace.map(|i| W::two_pi() * i / W::from(m - 1).unwrap() - W::pi()); + let w: Vec<_> = self + .a + .iter() + .enumerate() + .map(|(k, a)| { + fac.clone() + .map(move |f| Float::cos(f * W::from(k).unwrap()) * W::from(*a).unwrap()) + }) + .fold( + vec![W::from(0).unwrap(); fac.clone().count()], + // Would this have made more sense if using ndarray::Array1? No, Array1 has no pop. + |acc, x| acc.iter().zip(x).map(|(&a, b)| a + b).collect(), + ) + .into_iter() + .collect(); + + truncate(w, needs_trunc) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn general_cosine_scipy_eg() { + // Created with + // >>> from scipy.signal.windows import general_cosine + // >>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097] + // >>> window = general_cosine(30, HFT90D, sym=False) + let expected = vec![ + -1.04083409e-16, + -3.49808964e-03, + -1.85322176e-02, + -5.60667246e-02, + -1.25488810e-01, + -2.22198500e-01, + -3.14696393e-01, + -3.38497088e-01, + -2.04818447e-01, + 1.72651725e-01, + 8.38783500e-01, + 1.76097559e+00, + 2.81469639e+00, + 3.80321808e+00, + 4.51005597e+00, + 4.76683000e+00, + 4.51005597e+00, + 3.80321808e+00, + 2.81469639e+00, + 1.76097559e+00, + 8.38783500e-01, + 1.72651725e-01, + -2.04818447e-01, + -3.38497088e-01, + -3.14696393e-01, + -2.22198500e-01, + -1.25488810e-01, + -5.60667246e-02, + -1.85322176e-02, + -3.49808964e-03, + ]; + + let hfd90 = [1., 1.942604, 1.340318, 0.440811, 0.043097].into(); + let gc = GeneralCosine::new(30, hfd90, false); + assert_vec_eq(expected, gc.get_window()); + } + + #[test] + fn constant() { + let x = 0.42; + let n = 6; + let expected = vec![x; n]; + + let actual: Vec = GeneralCosine::new(n, vec![x], true).get_window(); + assert_eq!(expected, actual); + } + + #[test] + fn case_b() { + // Created with + // >>> from scipy.signal.windows import general_cosine + // >>> n = 20 + // >>> a = [0.42, 0.50] + // >>> window = general_cosine(30, a, sym=False) + let n = 20; + let a = vec![0.42, 0.50]; + let expected = vec![ + -0.08, + -0.05552826, + 0.0154915, + 0.12610737, + 0.2654915, + 0.42, + 0.5745085, + 0.71389263, + 0.8245085, + 0.89552826, + 0.92, + 0.89552826, + 0.8245085, + 0.71389263, + 0.5745085, + 0.42, + 0.2654915, + 0.12610737, + 0.0154915, + -0.05552826, + ]; + + assert_vec_eq(expected, GeneralCosine::new(n, a, false).get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index fc1393dc..354940cb 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -1,5 +1,6 @@ #[cfg(feature = "alloc")] use alloc::vec::Vec; +use nalgebra::RealField; use num_traits::{real::Real, Float}; /// Corresponding window representation for tuple-structs of [Window] variants. @@ -58,17 +59,24 @@ fn truncate(mut w: Vec, needed: bool) -> Vec { } mod boxcar; +mod general_cosine; mod triangle; pub use boxcar::Boxcar; +pub use general_cosine::GeneralCosine; pub use triangle::Triangle; -/// todo +/// This collects all structs that implement the [GetWindow] trait. +/// This allows for running `.get_window()` on the struct, which can then be, for example, used in +/// Firwin. // Ordering is as in accordance with // https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html. #[derive(Debug, Clone, PartialEq)] // Is it possible for the enums which wraps the structs to only require the generic that the struct // implements GetWindow? -pub enum Window { +pub enum Window +where + F: Real, +{ /// [Boxcar] window, also known as a rectangular window or Dirichlet window; This is equivalent /// to no window at all. Boxcar(Boxcar), @@ -93,8 +101,9 @@ pub enum Window { // Kaiser, // Needs Beta // KaiserBesselDerived, // Needs Beta // Gaussian, // Needs Standard Deviation - // Generic weighted sum of cosine term windows. - // GenericCosine(GenericCosine), // Needs Weighting Coefficients + /// [GeneralCosine] window, a generic weighted sum of cosine term windows. + // Needs Weighting Coefficients + GeneralCosine(GeneralCosine), // GeneralGaussian, // Needs Power, Width // GeneralHamming, // Needs Window Coefficients. // Dpss, // Needs Normalized Half-Bandwidth. From b3b4ff9f0dc4a89b18686f8181cb9d1649426aa5 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Fri, 22 Nov 2024 18:31:06 +0800 Subject: [PATCH 08/27] Add Blackman window --- sci-rs/src/signal/windows/blackman.rs | 120 ++++++++++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 5 +- 2 files changed, 124 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/blackman.rs diff --git a/sci-rs/src/signal/windows/blackman.rs b/sci-rs/src/signal/windows/blackman.rs new file mode 100644 index 00000000..5f6f2718 --- /dev/null +++ b/sci-rs/src/signal/windows/blackman.rs @@ -0,0 +1,120 @@ +use super::GeneralCosine; +use nalgebra::RealField; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `Blackman` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct Blackman { + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Blackman { + /// Returns a Blackman struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, sym: bool) -> Self { + Blackman { m, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for Blackman +where + W: Real + Float + RealField, +{ + /// Return a window of type: Blackman. + /// + /// The Blackman window is a taper formed by using the first three terms of a summation of + /// cosines. It was designed to have close to the minimal leakage possible. It is close to + /// optimal, only slightly worse than a Kaiser window. + /// + /// # Parameters + /// `self`: [Blackman] + /// + /// # Returns + /// `w`: `vec` + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # Example + /// ``` + /// use sci_rs::signal::windows::{Blackman, GetWindow}; + /// + /// let nx = 8; + /// let b = Blackman::new(nx, true); + /// ``` + /// + /// # References + /// + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + GeneralCosine::new( + self.m, + [0.42, 0.50, 0.08] + .into_iter() + .map(|n| W::from(n).unwrap()) + .collect(), + self.sym, + ) + .get_window() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn blackman_20() { + // Created with + // >>> from scipy.signal.windows import blackman + // >>> blackman(20) + let expected = vec![ + -1.38777878e-17, + 1.02226199e-02, + 4.50685843e-02, + 1.14390287e-01, + 2.26899356e-01, + 3.82380768e-01, + 5.66665187e-01, + 7.52034438e-01, + 9.03492728e-01, + 9.88846031e-01, + 9.88846031e-01, + 9.03492728e-01, + 7.52034438e-01, + 5.66665187e-01, + 3.82380768e-01, + 2.26899356e-01, + 1.14390287e-01, + 4.50685843e-02, + 1.02226199e-02, + -1.38777878e-17, + ]; + assert_vec_eq(expected, Blackman::new(20, true).get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 354940cb..bd1e1181 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -58,9 +58,11 @@ fn truncate(mut w: Vec, needed: bool) -> Vec { w } +mod blackman; mod boxcar; mod general_cosine; mod triangle; +pub use blackman::Blackman; pub use boxcar::Boxcar; pub use general_cosine::GeneralCosine; pub use triangle::Triangle; @@ -82,7 +84,8 @@ where Boxcar(Boxcar), /// [Triangle] window. Triangle(Triangle), - // Blackman, + /// [Blackman] window. + Blackman(Blackman), // Hamming, // Hann, // Bartlett, From e5033da8a4970053a8de902828caf4368b0cb740 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sat, 23 Nov 2024 14:00:45 +0800 Subject: [PATCH 09/27] Add more documentation to signal module --- sci-rs/src/signal/filter/design/kaiser.rs | 2 ++ sci-rs/src/signal/filter/mod.rs | 4 +++- sci-rs/src/signal/mod.rs | 19 ++++++++++++++----- sci-rs/src/signal/wave/mod.rs | 7 ++++--- 4 files changed, 23 insertions(+), 9 deletions(-) diff --git a/sci-rs/src/signal/filter/design/kaiser.rs b/sci-rs/src/signal/filter/design/kaiser.rs index 67be1e45..68b0e85b 100644 --- a/sci-rs/src/signal/filter/design/kaiser.rs +++ b/sci-rs/src/signal/filter/design/kaiser.rs @@ -26,6 +26,8 @@ use num_traits::{real::Real, MulAdd, Pow, ToPrimitive}; /// use sci_rs::signal::filter::design::kaiser_beta; /// assert_eq!(6.20426, kaiser_beta(65.)); /// ``` +/// # See Also +/// [kaiser_atten], [kaiserord] pub fn kaiser_beta(a: F) -> F where F: Real + MulAdd + Pow, diff --git a/sci-rs/src/signal/filter/mod.rs b/sci-rs/src/signal/filter/mod.rs index 805c65b7..ae7bd936 100644 --- a/sci-rs/src/signal/filter/mod.rs +++ b/sci-rs/src/signal/filter/mod.rs @@ -13,7 +13,9 @@ pub use kalmanfilt::kalman::kalman_filter; /// pub use gaussfilt as gaussian_filter; -/// Digital IIR/FIR filter design +/// Digital IIR/FIR filter design +/// Functions located in the [`Filter design` section of +/// `scipy.signal`](https://docs.scipy.org/doc/scipy/reference/signal.html#filter-design). pub mod design; mod ext; diff --git a/sci-rs/src/signal/mod.rs b/sci-rs/src/signal/mod.rs index ce9571cd..5dd9a505 100644 --- a/sci-rs/src/signal/mod.rs +++ b/sci-rs/src/signal/mod.rs @@ -1,16 +1,22 @@ -/// Digital Filtering +/// Digital Filtering +/// Contains functions from [Filtering section of +/// `scipy.signal`](https://docs.scipy.org/doc/scipy/reference/signal.html#filtering). pub mod filter; -/// Signal Generation +/// Signal Generation +/// Contains functions from the [Waveforms section of +/// `scipy.signal`](). pub mod wave; -/// Convolution +/// Convolution +/// Contains functions from the [Convolution section of +/// `scipy.signal`](). #[cfg(feature = "std")] pub mod convolve; /// Window functions /// This contains all window functions in the -/// [scipy.signal.windows](https://docs.scipy.org/doc/scipy/reference/signal.windows.html#module-scipy.signal.windows) +/// [`scipy.signal.windows`](https://docs.scipy.org/doc/scipy/reference/signal.windows.html#module-scipy.signal.windows) /// namespace. /// The convenience function /// [`get_windows`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html#scipy.signal.get_window) @@ -18,6 +24,9 @@ pub mod convolve; /// namespace is located here. pub mod windows; -/// Signal Resampling +/// Signal Resampling +/// This contains only the +/// [`resample`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html#scipy.signal.resample) +/// function from `scipy.signal`. #[cfg(feature = "std")] pub mod resample; diff --git a/sci-rs/src/signal/wave/mod.rs b/sci-rs/src/signal/wave/mod.rs index 87372fcc..90bd9ea0 100644 --- a/sci-rs/src/signal/wave/mod.rs +++ b/sci-rs/src/signal/wave/mod.rs @@ -1,7 +1,6 @@ use nalgebra::RealField; use ndarray::{Array, ArrayBase, Data, Dimension, RawData}; -/// """ /// Return a periodic square-wave waveform. /// /// The square wave has a period ``2*pi``, has value +1 from 0 to @@ -30,15 +29,18 @@ use ndarray::{Array, ArrayBase, Data, Dimension, RawData}; /// -------- /// A 5 Hz waveform sampled at 500 Hz for 1 second: /// +/// ```custom,{class=language-python} /// >>> import numpy as np /// >>> from scipy import signal /// >>> import matplotlib.pyplot as plt /// >>> t = np.linspace(0, 1, 500, endpoint=False) /// >>> plt.plot(t, signal.square(2 * np.pi * 5 * t)) /// >>> plt.ylim(-2, 2) +/// ``` /// /// A pulse-width modulated sine wave: /// +/// ```custom,{class=language-python} /// >>> plt.figure() /// >>> sig = np.sin(2 * np.pi * t) /// >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2) @@ -47,8 +49,7 @@ use ndarray::{Array, ArrayBase, Data, Dimension, RawData}; /// >>> plt.subplot(2, 1, 2) /// >>> plt.plot(t, pwm) /// >>> plt.ylim(-1.5, 1.5) -/// -/// """ +/// ``` pub fn square(t: &ArrayBase, duty: F) -> Array where F: RealField, From 35aad91993ecf696a96813cb4dddbc8613010431 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Sun, 24 Nov 2024 12:52:45 +0800 Subject: [PATCH 10/27] Add general gaussian window --- sci-rs/src/signal/windows/general_gaussian.rs | 148 ++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 6 +- 2 files changed, 153 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/general_gaussian.rs diff --git a/sci-rs/src/signal/windows/general_gaussian.rs b/sci-rs/src/signal/windows/general_gaussian.rs new file mode 100644 index 00000000..9f19801a --- /dev/null +++ b/sci-rs/src/signal/windows/general_gaussian.rs @@ -0,0 +1,148 @@ +use super::{extend, len_guard, truncate}; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `GeneralCosine` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct GeneralGaussian +where + F: Real, +{ + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Shape parameter. + pub p: F, + /// The standard deviation, σ. + pub sigma: F, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl GeneralGaussian +where + F: Real, +{ + /// Returns a GeneralGaussian struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `p` : float + /// Shape parameter. p = 1 is identical to `gaussian`, p = 0.5 is + /// the same shape as the Laplace distribution. + /// * `sig` : float + /// The standard deviation, σ. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, p: F, sigma: F, sym: bool) -> Self { + GeneralGaussian { m, p, sigma, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for GeneralGaussian +where + F: Real, + W: Real, +{ + /// Return a window with a generalized gaussian shape. + /// + /// # Parameters + /// `self`: [GeneralGaussian] + /// + /// # Returns + /// `w`: `vec` + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # Notes + /// The generalized Gaussian window is defined as + /// $$w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} }$$ + /// the half-power point is at + /// $$(2 \log(2))^{1/(2 p)} \sigma$$ + /// + /// # Example + /// ``` + /// use sci_rs::signal::windows::{GeneralGaussian, GetWindow}; + /// let window: Vec = GeneralGaussian::new(51, 1.5, 7., true).get_window(); + /// ``` + /// + /// This is equivalent to the Python code: + /// ```custom,{class=language-python} + /// from scipy import signal + /// window = signal.windows.general_gaussian(51, p=1.5, sig=7) + /// ``` + /// + /// # References + /// + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + if len_guard(self.m) { + return Vec::::new(); + } + let (m, needs_trunc) = extend(self.m, self.sym); + + let n = (0..m).map(|v| { + W::from(v).unwrap() - (W::from(m).unwrap() - W::from(1).unwrap()) / W::from(2).unwrap() + }); + let sig = W::from(self.sigma).unwrap(); + let two_p = W::from(self.p).unwrap() * W::from(2).unwrap(); + let w = n + .into_iter() + .map(|v| { + (v / sig) + .abs() + .powf(two_p) + .mul(W::from(-0.5).unwrap()) + .exp() + }) + .collect(); + + truncate(w, needs_trunc) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn general_gaussian_case_a() { + let gc = GeneralGaussian::new(20, 0.8, 4.1, true); + let expected = vec![ + 0.14688425, 0.2008071, 0.2687279, 0.35164027, 0.44932387, 0.55972797, 0.67829536, + 0.79725628, 0.90478285, 0.98289486, 0.98289486, 0.90478285, 0.79725628, 0.67829536, + 0.55972797, 0.44932387, 0.35164027, 0.2687279, 0.2008071, 0.14688425, + ]; + + assert_vec_eq(expected, gc.get_window()); + } + + #[test] + fn general_gaussian_case_b() { + let gc = GeneralGaussian::new(20, 0.8, 4.1, false); + let expected = vec![ + 0.12465962, 0.17219048, 0.23293383, 0.30828863, 0.3987132, 0.5031538, 0.61839303, + 0.73835825, 0.85338105, 0.94904249, 1., 0.94904249, 0.85338105, 0.73835825, 0.61839303, + 0.5031538, 0.3987132, 0.30828863, 0.23293383, 0.17219048, + ]; + + assert_vec_eq(expected, gc.get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index bd1e1181..942659f7 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -61,10 +61,12 @@ fn truncate(mut w: Vec, needed: bool) -> Vec { mod blackman; mod boxcar; mod general_cosine; +mod general_gaussian; mod triangle; pub use blackman::Blackman; pub use boxcar::Boxcar; pub use general_cosine::GeneralCosine; +pub use general_gaussian::GeneralGaussian; pub use triangle::Triangle; /// This collects all structs that implement the [GetWindow] trait. @@ -107,7 +109,9 @@ where /// [GeneralCosine] window, a generic weighted sum of cosine term windows. // Needs Weighting Coefficients GeneralCosine(GeneralCosine), - // GeneralGaussian, // Needs Power, Width + /// [GeneralGaussian] window. + // Needs Power, Width + GeneralGaussian(GeneralGaussian), // GeneralHamming, // Needs Window Coefficients. // Dpss, // Needs Normalized Half-Bandwidth. // Chebwin, // Needs Attenuation. From 575e428f091d0650d834d652faeefedd5d1bc83a Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Mon, 25 Nov 2024 17:09:08 +0800 Subject: [PATCH 11/27] Add General Hamming --- sci-rs/src/signal/windows/general_hamming.rs | 196 +++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 6 +- 2 files changed, 201 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/general_hamming.rs diff --git a/sci-rs/src/signal/windows/general_hamming.rs b/sci-rs/src/signal/windows/general_hamming.rs new file mode 100644 index 00000000..e3c93c9d --- /dev/null +++ b/sci-rs/src/signal/windows/general_hamming.rs @@ -0,0 +1,196 @@ +use super::GeneralCosine; +use nalgebra::RealField; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::vec::Vec; + +/// Collection of arguments for window `GeneralHamming` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct GeneralHamming +where + F: Real, +{ + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// The window coefficient, α. + pub alpha: F, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl GeneralHamming +where + F: Real, +{ + /// Returns a GeneralHamming struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `alpha` : float + /// The window coefficient, α. + /// * `sym`: bool + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, alpha: F, sym: bool) -> Self { + GeneralHamming { m, alpha, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for GeneralHamming +where + F: Real, + W: Real + Float + RealField, +{ + /// Return a generalized Hamming window. + /// + /// The generalized Hamming window is constructed by multiplying a rectangular + /// window by one period of a cosine function [[1][1]]. + /// + /// # Parameters + /// * `M` : usize + /// Number of points in the output window. If zero, an empty array is returned. An + /// exception is thrown when it is negative. + /// * `alpha` : float + /// The window coefficient, α. + /// * `sym` : bool, optional + /// When True (default), generates a symmetric window, for use in filter design. + /// When False, generates a periodic window, for use in spectral analysis. + /// + /// # Returns + /// `w` : ndarray + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # See Also + /// hamming, hann + /// + /// # Notes + /// The generalized Hamming window is defined as + /// $$w(n) = \alpha - \left(1 - \alpha\right) + /// \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$$ + /// Both the common Hamming window and Hann window are special cases of the generalized Hamming + /// window with :math:`\alpha` = 0.54 and :math:`\alpha` = 0.5, respectively [[2][2]]. + /// + /// # References + /// + /// [[1]] DSPRelated, "Generalized Hamming Window Family", + /// + /// [[2]] Wikipedia, "Window function", + /// [[3]] Riccardo Piantanida ESA, "Sentinel-1 Level 1 Detailed Algorithm Definition", + /// + /// [[4]] Matthieu Bourbigot ESA, "Sentinel-1 Product Definition", + /// + /// [[5]] Scipy, + /// + /// + /// # Examples + /// The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming + /// windows in the processing of spaceborne Synthetic Aperture Radar (SAR) + /// data [[3][3]]. The facility uses various values for the :math:`\alpha` + /// parameter based on operating mode of the SAR instrument. Some common + /// :math:`\alpha` values include 0.75, 0.7 and 0.52 [[4][4]]. As an example, we + /// plot these different windows. + /// + /// ```custom,{class=language-python} + /// >>> import numpy as np + /// >>> from scipy.signal.windows import general_hamming + /// >>> from scipy.fft import fft, fftshift + /// >>> import matplotlib.pyplot as plt + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> fig1, spatial_plot = plt.subplots() + /// >>> spatial_plot.set_title("Generalized Hamming Windows") + /// >>> spatial_plot.set_ylabel("Amplitude") + /// >>> spatial_plot.set_xlabel("Sample") + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> fig2, freq_plot = plt.subplots() + /// >>> freq_plot.set_title("Frequency Responses") + /// >>> freq_plot.set_ylabel("Normalized magnitude [dB]") + /// >>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]") + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> for alpha in [0.75, 0.7, 0.52]: + /// ... window = general_hamming(41, alpha) + /// ... spatial_plot.plot(window, label="{:.2f}".format(alpha)) + /// ... A = fft(window, 2048) / (len(window)/2.0) + /// ... freq = np.linspace(-0.5, 0.5, len(A)) + /// ... response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) + /// ... freq_plot.plot(freq, response, label="{:.2f}".format(alpha)) + /// >>> freq_plot.legend(loc="upper right") + /// >>> spatial_plot.legend(loc="upper right") + /// ``` + /// + /// The equivalent is: + /// ``` + /// use sci_rs::signal::windows::{GetWindow, GeneralHamming}; + /// let window = GeneralHamming::new(41, 0.75, true); + /// ``` + /// + /// [1]: #references + /// [2]: #references + /// [3]: #references + /// [4]: #references + /// [5]: #references + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + GeneralCosine::new( + self.m, + [self.alpha, F::from(1).unwrap() - self.alpha].into(), + self.sym, + ) + .get_window() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn general_hamming_case_a() { + // from scipy.signal.windows import general_hamming + // general_hamming(15, 0.65) + let gh = GeneralHamming::new(15, 0.65, true); + + let expected = [ + 0.3, 0.3346609, 0.43177857, 0.57211767, 0.72788233, 0.86822143, 0.9653391, 1., + 0.9653391, 0.86822143, 0.72788233, 0.57211767, 0.43177857, 0.3346609, 0.3, + ] + .into(); + assert_vec_eq(expected, gh.get_window()); + } + + #[test] + fn general_hamming_case_b() { + // from scipy.signal.windows import general_hamming + // general_hamming(15, 0.65) + let gh = GeneralHamming::new(15, 0.65, false); + + let expected = [ + 0.3, 0.33025909, 0.41580429, 0.54184405, 0.68658496, 0.825, 0.93315595, 0.99235166, + 0.99235166, 0.93315595, 0.825, 0.68658496, 0.54184405, 0.41580429, 0.33025909, + ] + .into(); + assert_vec_eq(expected, gh.get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 942659f7..d48a6285 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -62,11 +62,13 @@ mod blackman; mod boxcar; mod general_cosine; mod general_gaussian; +mod general_hamming; mod triangle; pub use blackman::Blackman; pub use boxcar::Boxcar; pub use general_cosine::GeneralCosine; pub use general_gaussian::GeneralGaussian; +pub use general_hamming::GeneralHamming; pub use triangle::Triangle; /// This collects all structs that implement the [GetWindow] trait. @@ -112,7 +114,9 @@ where /// [GeneralGaussian] window. // Needs Power, Width GeneralGaussian(GeneralGaussian), - // GeneralHamming, // Needs Window Coefficients. + /// [GeneralHamming] window. + // Needs Window Coefficients. + GeneralHamming(GeneralHamming), // Dpss, // Needs Normalized Half-Bandwidth. // Chebwin, // Needs Attenuation. } From a052c208e2c55419473df9eefc48da6c98ab68fa Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 26 Nov 2024 16:23:13 +0800 Subject: [PATCH 12/27] Add Hamming window --- sci-rs/src/signal/windows/hamming.rs | 171 +++++++++++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 5 +- 2 files changed, 175 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/hamming.rs diff --git a/sci-rs/src/signal/windows/hamming.rs b/sci-rs/src/signal/windows/hamming.rs new file mode 100644 index 00000000..90016f5c --- /dev/null +++ b/sci-rs/src/signal/windows/hamming.rs @@ -0,0 +1,171 @@ +use super::GeneralHamming; +use nalgebra::RealField; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::vec::Vec; + +/// Collection of arguments for window `Hamming` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct Hamming { + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Hamming { + /// Returns a Hamming struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `sym`: bool + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, sym: bool) -> Self { + Hamming { m, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for Hamming +where + W: Real + Float + RealField, +{ + /// Return a Hamming window. + /// + /// The Hamming window is a taper formed by using a raised cosine with non-zero endpoints, + /// optimized to minimize the nearest side lobe. + /// + /// # Parameters + /// * `M` : int + /// Number of points in the output window. If zero, an empty array is returned. An + /// exception is thrown when it is negative. + /// * `sym` : bool, optional + /// When True (default), generates a symmetric window, for use in filter + /// design. + /// When False, generates a periodic window, for use in spectral analysis. + /// + /// # Returns + /// `w` : ndarray + /// The window, with the maximum value normalized to 1 (though the value 1 + /// does not appear if `M` is even and `sym` is True). + /// + /// # Notes + /// The Hamming window is defined as + /// + /// $$w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1$$ + /// + /// The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and is described in + /// Blackman and Tukey. It was recommended for smoothing the truncated autocovariance function + /// in the time domain. Most references to the Hamming window come from the signal processing + /// literature, where it is used as one of many windowing functions for smoothing values. It + /// is also known as an apodization (which means "removing the foot", i.e. smoothing + /// discontinuities at the beginning and end of the sampled signal) or tapering function. + /// + /// # References + /// [[1]] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover + /// Publications, New York. + /// [[2]] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The University of Alberta + /// Press, 1975, pp. 109-110. + /// [[3]] Wikipedia, "Window function", + /// [[4]] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, "Numerical Recipes", + /// Cambridge University Press, 1986, page 425. + /// [[5]] Scipy, + /// + /// + /// Examples + /// -------- + /// Plot the window and its frequency response: + /// + /// ```custom,{class=language-python} + /// >>> import numpy as np + /// >>> from scipy import signal + /// >>> from scipy.fft import fft, fftshift + /// >>> import matplotlib.pyplot as plt + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> window = signal.windows.hamming(51) + /// >>> plt.plot(window) + /// >>> plt.title("Hamming window") + /// >>> plt.ylabel("Amplitude") + /// >>> plt.xlabel("Sample") + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> plt.figure() + /// >>> A = fft(window, 2048) / (len(window)/2.0) + /// >>> freq = np.linspace(-0.5, 0.5, len(A)) + /// >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) + /// >>> plt.plot(freq, response) + /// >>> plt.axis([-0.5, 0.5, -120, 0]) + /// >>> plt.title("Frequency response of the Hamming window") + /// >>> plt.ylabel("Normalized magnitude [dB]") + /// >>> plt.xlabel("Normalized frequency [cycles per sample]") + /// ``` + /// + /// The equivalent is: + /// ``` + /// use sci_rs::signal::windows::{GetWindow, Hamming}; + /// let window: Vec = Hamming::new(51, true).get_window(); + /// ``` + /// + /// [1]: #references + /// [2]: #references + /// [3]: #references + /// [4]: #references + /// [5]: #references + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + GeneralHamming::::new(self.m, W::from(0.54).unwrap(), self.sym).get_window() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn hamming_case_a() { + // from scipy.signal.windows import hamming + // hamming(17) + let h = Hamming::new(17, true); + let expected = [ + 0.08, 0.11501542, 0.21473088, 0.36396562, 0.54, 0.71603438, 0.86526912, 0.96498458, 1., + 0.96498458, 0.86526912, 0.71603438, 0.54, 0.36396562, 0.21473088, 0.11501542, 0.08, + ] + .into(); + + assert_vec_eq(expected, h.get_window()); + } + + #[test] + fn hamming_case_b() { + // from scipy.signal.windows import hamming + // hamming(17, false) + let h = Hamming::new(17, false); + let expected = [ + 0.08, 0.11106277, 0.2000559, 0.33496036, 0.49755655, 0.66588498, 0.81721193, + 0.93109988, 0.99216763, 0.99216763, 0.93109988, 0.81721193, 0.66588498, 0.49755655, + 0.33496036, 0.2000559, 0.11106277, + ] + .into(); + + assert_vec_eq(expected, h.get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index d48a6285..096969e2 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -63,12 +63,14 @@ mod boxcar; mod general_cosine; mod general_gaussian; mod general_hamming; +mod hamming; mod triangle; pub use blackman::Blackman; pub use boxcar::Boxcar; pub use general_cosine::GeneralCosine; pub use general_gaussian::GeneralGaussian; pub use general_hamming::GeneralHamming; +pub use hamming::Hamming; pub use triangle::Triangle; /// This collects all structs that implement the [GetWindow] trait. @@ -90,7 +92,8 @@ where Triangle(Triangle), /// [Blackman] window. Blackman(Blackman), - // Hamming, + /// [Hamming] window. + Hamming(Hamming), // Hann, // Bartlett, // Flattop, From 4696efcf1f238b7a5a9a2e8288f3949d1ccfb217 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Thu, 28 Nov 2024 18:12:36 +0800 Subject: [PATCH 13/27] Add kaiser window --- sci-rs/src/signal/windows/kaiser.rs | 228 ++++++++++++++++++++++++++++ sci-rs/src/signal/windows/mod.rs | 7 +- 2 files changed, 233 insertions(+), 2 deletions(-) create mode 100644 sci-rs/src/signal/windows/kaiser.rs diff --git a/sci-rs/src/signal/windows/kaiser.rs b/sci-rs/src/signal/windows/kaiser.rs new file mode 100644 index 00000000..2029f07e --- /dev/null +++ b/sci-rs/src/signal/windows/kaiser.rs @@ -0,0 +1,228 @@ +use super::{extend, len_guard, truncate}; +use crate::special::Bessel; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; + +#[cfg(feature = "alloc")] +use alloc::{vec, vec::Vec}; + +/// Collection of arguments for window `Kaiser` for use in [GetWindow]. +#[cfg(feature = "alloc")] +#[derive(Debug, Clone, PartialEq)] +pub struct Kaiser +where + F: Real, +{ + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Shape parameter. + pub beta: F, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Kaiser +where + F: Real, +{ + /// Returns a Kaiser struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `beta` : float + /// Shape parameter, determines trade-off between main-lobe width and side lobe level. As + /// beta gets large, the window narrows. + /// * `sym`: + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, beta: F, sym: bool) -> Self { + Kaiser { m, beta, sym } + } +} + +impl GetWindow for Kaiser +where + F: Real, + W: Real + Bessel, +{ + /// Return a [Kaiser] window. + /// + /// The Kaiser window is a taper formed by using a Bessel function. + /// + /// Parameters + /// ---------- + /// * `M` : int + /// Number of points in the output window. If zero, an empty array is returned. An + /// exception is thrown when it is negative. + /// * `beta` : float + /// Shape parameter, determines trade-off between main-lobe width and side lobe level. As + /// beta gets large, the window narrows. + /// * `sym` : bool, optional + /// When True (default), generates a symmetric window, for use in filter design. + /// When False, generates a periodic window, for use in spectral analysis. + /// + /// Returns + /// ------- + /// `w` : ndarray + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// Notes + /// ----- + /// The Kaiser window is defined as + /// + /// $$w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}} \right)/I_0(\beta)$$ + /// + /// with + /// + /// $$\quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2}$$, + /// + /// where $I_0$ is the modified zeroth-order Bessel function. + /// + /// The Kaiser was named for Jim Kaiser, who discovered a simple approximation to the DPSS + /// window based on Bessel functions. + /// The Kaiser window is a very good approximation to the Digital Prolate Spheroidal Sequence, + /// or Slepian window, which is the transform which maximizes the energy in the main lobe of + /// the window relative to total energy. + /// + /// The Kaiser can approximate other windows by varying the beta parameter. [(Some literature + /// uses alpha = beta/pi.)]([4]) + /// + /// ==== ======================= + /// beta Window shape + /// ==== ======================= + /// 0 Rectangular + /// 5 Similar to a Hamming + /// 6 Similar to a Hann + /// 8.6 Similar to a Blackman + /// ==== ======================= + /// + /// A beta value of 14 is probably a good starting point. Note that as beta gets large, the + /// window narrows, and so the number of samples needs to be large enough to sample the + /// increasingly narrow spike, otherwise NaNs will be returned. + /// + /// Most references to the Kaiser window come from the signal processing literature, where it + /// is used as one of many windowing functions for smoothing values. It is also known as an + /// apodization (which means "removing the foot", i.e. smoothing discontinuities at the + /// beginning and end of the sampled signal) or tapering function. + /// + /// # References + /// [[1]] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by digital computer", + /// Editors: F.F. Kuo and J.F. Kaiser, p 218-285. John Wiley and Sons, New York, (1966). + /// [[2]] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The University of Alberta + /// Press, 1975, pp. 177-178. + /// [[3]] Wikipedia, "Window function", + /// [[4]] F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier + /// transform," Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978. + /// :doi:`10.1109/PROC.1978.10837`. + /// [[5]] + /// [Scipy]() + /// + /// # Examples + /// Plot the window and its frequency response: + /// + /// ```custom,{class=language-python} + /// >>> import numpy as np + /// >>> from scipy import signal + /// >>> from scipy.fft import fft, fftshift + /// >>> import matplotlib.pyplot as plt + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> window = signal.windows.kaiser(51, beta=14) + /// >>> plt.plot(window) + /// >>> plt.title(r"Kaiser window ($\beta$=14)") + /// >>> plt.ylabel("Amplitude") + /// >>> plt.xlabel("Sample") + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> plt.figure() + /// >>> A = fft(window, 2048) / (len(window)/2.0) + /// >>> freq = np.linspace(-0.5, 0.5, len(A)) + /// >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) + /// >>> plt.plot(freq, response) + /// >>> plt.axis([-0.5, 0.5, -120, 0]) + /// >>> plt.title(r"Frequency response of the Kaiser window ($\beta$=14)") + /// >>> plt.ylabel("Normalized magnitude [dB]") + /// >>> plt.xlabel("Normalized frequency [cycles per sample]") + /// ``` + /// + /// The equivalent is: + /// ``` + /// use sci_rs::signal::windows::{GetWindow, Kaiser}; + /// let window: Vec = Kaiser::new(51, 14., true).get_window(); + /// ``` + /// + /// Due to current implementation limitations, note that the following might not possible: + /// ``` + /// use sci_rs::signal::windows::{GetWindow, Kaiser}; + /// let window: Vec = Kaiser::new(51, 14., true).get_window(); + /// println!("window = {:?}", window); + /// ``` + /// + /// [1]: #references + /// [2]: #references + /// [3]: #references + /// [4]: #references + /// [5]: #references + fn get_window(&self) -> Vec { + if len_guard(self.m) { + return Vec::::new(); + } + let (m, needs_trunc) = extend(self.m, self.sym); + let n = (0..m); + let alpha = W::from(self.m - 1).unwrap() / W::from(2).unwrap(); + let beta = W::from(self.beta).unwrap(); + // let w: Vec = n + // .map(|ni| W::from(ni).unwrap() - alpha) + // .map(|ni| (ni / alpha).powf(W::from(2).unwrap())) + // .map(|ni| (W::one() - ni).sqrt()) + // .map(|ni| i0(beta * ni) / i0(beta)) + // .collect(); + let w: Vec = n + .map(|ni| { + (beta + * (W::one() + - ((W::from(ni).unwrap() - alpha) / alpha).powf(W::from(2).unwrap())) + .sqrt()) + .i0() + / (beta.i0()) + }) + .collect(); + truncate(w, needs_trunc) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn kaiser_17_8_true() { + // from scipy.signal.windows import kaiser + // kaiser(17, beta = 0.8) + let expected = vec![ + 0.85725436, 0.88970403, 0.9183205, 0.94289618, 0.96325245, 0.97924114, 0.99074569, + 0.99768219, 1., 0.99768219, 0.99074569, 0.97924114, 0.96325245, 0.94289618, 0.9183205, + 0.88970403, 0.85725436, + ]; + let k = Kaiser::new(17, 0.8, true); + + assert_vec_eq(expected, k.get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 096969e2..4c9c3929 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -64,6 +64,7 @@ mod general_cosine; mod general_gaussian; mod general_hamming; mod hamming; +mod kaiser; mod triangle; pub use blackman::Blackman; pub use boxcar::Boxcar; @@ -71,6 +72,7 @@ pub use general_cosine::GeneralCosine; pub use general_gaussian::GeneralGaussian; pub use general_hamming::GeneralHamming; pub use hamming::Hamming; +pub use kaiser::Kaiser; pub use triangle::Triangle; /// This collects all structs that implement the [GetWindow] trait. @@ -107,8 +109,9 @@ where // Tukey, // Taylor, // Lanczos, - // Kaiser Window. - // Kaiser, // Needs Beta + /// [Kaiser] window. + // Needs Beta + Kaiser(Kaiser), // KaiserBesselDerived, // Needs Beta // Gaussian, // Needs Standard Deviation /// [GeneralCosine] window, a generic weighted sum of cosine term windows. From fe309d32ce3dc720e446af9883c77958e9013cb6 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Fri, 29 Nov 2024 15:51:29 +0800 Subject: [PATCH 14/27] Add library level errors The error variants here should guide the user as to the wrong use of functions, instead of taking down the entire program with a panic, which is especially necessary in an embedded environment. --- sci-rs/src/error/mod.rs | 27 +++++++++++++++++++++++++++ sci-rs/src/lib.rs | 3 +++ 2 files changed, 30 insertions(+) create mode 100644 sci-rs/src/error/mod.rs diff --git a/sci-rs/src/error/mod.rs b/sci-rs/src/error/mod.rs new file mode 100644 index 00000000..c4fe9ade --- /dev/null +++ b/sci-rs/src/error/mod.rs @@ -0,0 +1,27 @@ +use core::{error, fmt}; + +/// Errors raised whilst running sci-rs. +#[derive(Debug, PartialEq, Eq)] +#[cfg(feature = "alloc")] +pub enum Error { + /// Argument parsed into function were invalid. + InvalidArg { + /// The invalid arg + arg: alloc::string::String, + /// Explaining why arg is invalid. + reason: alloc::string::String, + }, + /// Two or more optional arguments passed into functions conflict. + ConflictArg { + /// Explaining what arg is invalid. + reason: alloc::string::String, + }, +} + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + todo!() + } +} + +impl error::Error for Error {} diff --git a/sci-rs/src/lib.rs b/sci-rs/src/lib.rs index 7e47de09..c19eec65 100644 --- a/sci-rs/src/lib.rs +++ b/sci-rs/src/lib.rs @@ -41,3 +41,6 @@ pub mod special; /// Debug plotting #[cfg(feature = "plot")] pub mod plot; + +/// Errors +pub mod error; From ca2a0645660f4314b264fdf1e084f08356b1d73e Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 3 Dec 2024 18:07:19 +0800 Subject: [PATCH 15/27] Add Nuttall window --- sci-rs/src/signal/windows/mod.rs | 5 +- sci-rs/src/signal/windows/nuttall.rs | 192 +++++++++++++++++++++++++++ 2 files changed, 196 insertions(+), 1 deletion(-) create mode 100644 sci-rs/src/signal/windows/nuttall.rs diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 4c9c3929..ab6c22d9 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -65,6 +65,7 @@ mod general_gaussian; mod general_hamming; mod hamming; mod kaiser; +mod nuttall; mod triangle; pub use blackman::Blackman; pub use boxcar::Boxcar; @@ -73,6 +74,7 @@ pub use general_gaussian::GeneralGaussian; pub use general_hamming::GeneralHamming; pub use hamming::Hamming; pub use kaiser::Kaiser; +pub use nuttall::Nuttall; pub use triangle::Triangle; /// This collects all structs that implement the [GetWindow] trait. @@ -102,7 +104,8 @@ where // Parzen, // Bohman, // BlackmanHarris, - // Nuttall, + /// [Nuttall] window. + Nuttall(Nuttall), // BartHann, // Cosine, // Exponential, diff --git a/sci-rs/src/signal/windows/nuttall.rs b/sci-rs/src/signal/windows/nuttall.rs new file mode 100644 index 00000000..d3d37f5b --- /dev/null +++ b/sci-rs/src/signal/windows/nuttall.rs @@ -0,0 +1,192 @@ +use super::GeneralCosine; +use nalgebra::RealField; +use num_traits::{real::Real, Float}; + +#[cfg(feature = "alloc")] +use super::GetWindow; +#[cfg(feature = "alloc")] +use alloc::vec::Vec; + +/// Collection of arguments for window `Nuttall` for use in [GetWindow]. +#[derive(Debug, Clone, PartialEq)] +pub struct Nuttall { + /// Number of points in the output window. If zero, an empty array is returned in [GetWindow]. + pub m: usize, + /// Whether the window is symmetric. + /// + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub sym: bool, +} + +impl Nuttall { + /// Returns a Nuttall struct. + /// + /// # Parameters + /// * `m`: + /// Number of points in the output window. If zero, an empty array is returned. + /// * `sym`: bool + /// When true, generates a symmetric window, for use in filter design. + /// When false, generates a periodic window, for use in spectral analysis. + pub fn new(m: usize, sym: bool) -> Self { + Nuttall { m, sym } + } +} + +#[cfg(feature = "alloc")] +impl GetWindow for Nuttall +where + W: Real + Float + RealField, +{ + /// Return a minimum 4-term Blackman-Harris window according to Nuttall. + /// + /// This variation is called "Nuttall4c" by Heinzel. [[2]] + /// + /// # Parameters + /// * `M` : int + /// Number of points in the output window. If zero, an empty array is returned. An + /// exception is thrown when it is negative. + /// * `sym` : bool, optional + /// When True (default), generates a symmetric window, for use in filter + /// design. + /// When False, generates a periodic window, for use in spectral analysis. + /// + /// # Returns + /// `w` : ndarray + /// The window, with the maximum value normalized to 1 (though the value 1 does not appear + /// if `M` is even and `sym` is True). + /// + /// # References + /// [[1]] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE Transactions on + /// Acoustics, Speech, and Signal Processing, vol. 29, no. 1, pp. 84-91, Feb 1981. + /// :doi:`10.1109/TASSP.1981.1163506`. + /// [[2]] Heinzel G. et al., "Spectrum and spectral density estimation by the Discrete Fourier + /// transform (DFT), including a comprehensive list of window functions and some new flat-top + /// windows", February 15, 2002 + /// [[3]] Scipy, + /// + /// + /// Examples + /// -------- + /// Plot the window and its frequency response: + /// + /// ```custom,{class=language-python} + /// >>> import numpy as np + /// >>> from scipy import signal + /// >>> from scipy.fft import fft, fftshift + /// >>> import matplotlib.pyplot as plt + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> window = signal.windows.nuttall(51) + /// >>> plt.plot(window) + /// >>> plt.title("Nuttall window") + /// >>> plt.ylabel("Amplitude") + /// >>> plt.xlabel("Sample") + /// ``` + /// + /// ```custom,{class=language-python} + /// >>> plt.figure() + /// >>> A = fft(window, 2048) / (len(window)/2.0) + /// >>> freq = np.linspace(-0.5, 0.5, len(A)) + /// >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) + /// >>> plt.plot(freq, response) + /// >>> plt.axis([-0.5, 0.5, -120, 0]) + /// >>> plt.title("Frequency response of the Hamming window") + /// >>> plt.ylabel("Normalized magnitude [dB]") + /// >>> plt.xlabel("Normalized frequency [cycles per sample]") + /// ``` + /// + /// The equivalent is: + /// ``` + /// use sci_rs::signal::windows::{GetWindow, Nuttall}; + /// let window: Vec = Nuttall::new(51, true).get_window(); + /// ``` + /// + /// [1]: #references + /// [2]: #references + /// [3]: #references + #[cfg(feature = "alloc")] + fn get_window(&self) -> Vec { + GeneralCosine::::new( + self.m, + [0.3635819, 0.4891775, 0.1365995, 0.0106411] + .map(|f| W::from(f).unwrap()) + .into_iter() + .collect(), + self.sym, + ) + .get_window() + } +} + +#[cfg(test)] +mod tests { + use super::*; + use approx::assert_abs_diff_eq; + + #[test] + fn nuttall_case_a() { + // from scipy.signal.windows import nuttall + // nuttall(17) + let h = Nuttall::new(17, true); + let expected = [ + 3.62800000e-04, + 4.15908007e-03, + 2.52055665e-02, + 8.96224370e-02, + 2.26982400e-01, + 4.44360497e-01, + 7.01958233e-01, + 9.16185585e-01, + 1.00000000e+00, + 9.16185585e-01, + 7.01958233e-01, + 4.44360497e-01, + 2.26982400e-01, + 8.96224370e-02, + 2.52055665e-02, + 4.15908007e-03, + 3.62800000e-04, + ] + .into(); + + assert_vec_eq(expected, h.get_window()); + } + + #[test] + fn nuttall_case_b() { + // from scipy.signal.windows import nuttall + // nuttall(17, false) + let h = Nuttall::new(17, false); + let expected = [ + 3.62800000e-04, + 3.64256817e-03, + 2.10918726e-02, + 7.36770505e-02, + 1.87084736e-01, + 3.73448574e-01, + 6.11072447e-01, + 8.39394793e-01, + 9.80852709e-01, + 9.80852709e-01, + 8.39394793e-01, + 6.11072447e-01, + 3.73448574e-01, + 1.87084736e-01, + 7.36770505e-02, + 2.10918726e-02, + 3.64256817e-03, + ] + .into(); + + assert_vec_eq(expected, h.get_window()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} From cd9b27d1014de504b248b058e7c66d8f0634ec0f Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 4 Dec 2024 23:54:06 +0800 Subject: [PATCH 16/27] Add convenience GetWindow impl for Window variants --- sci-rs/src/signal/windows/mod.rs | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index ab6c22d9..d14691b5 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -1,8 +1,10 @@ -#[cfg(feature = "alloc")] -use alloc::vec::Vec; +use crate::special; use nalgebra::RealField; use num_traits::{real::Real, Float}; +#[cfg(feature = "alloc")] +use alloc::vec::Vec; + /// Corresponding window representation for tuple-structs of [Window] variants. #[cfg(feature = "alloc")] pub trait GetWindow @@ -129,3 +131,23 @@ where // Dpss, // Needs Normalized Half-Bandwidth. // Chebwin, // Needs Attenuation. } + +impl GetWindow for Window +where + F: Real, + W: Real + Float + RealField + special::Bessel, +{ + fn get_window(&self) -> Vec { + match &self { + Window::Boxcar(x) => x.get_window(), + Window::Triangle(x) => x.get_window(), + Window::Blackman(x) => x.get_window(), + Window::Hamming(x) => x.get_window(), + Window::Nuttall(x) => x.get_window(), + Window::Kaiser(x) => x.get_window(), + Window::GeneralCosine(x) => x.get_window(), + Window::GeneralGaussian(x) => x.get_window(), + Window::GeneralHamming(x) => x.get_window(), + } + } +} From df21ea3208549daa40089959d1c7290c432f209b Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 4 Dec 2024 17:32:28 +0800 Subject: [PATCH 17/27] Add scipy.signal.windows.get_window function We also introduce the convenience enum to go along with the function. --- sci-rs/src/signal/windows/mod.rs | 159 +++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index d14691b5..e8b19ba8 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -151,3 +151,162 @@ where } } } + +/// This provides a set of enum variants that for use in [get_window]. +#[derive(Debug, Clone, PartialEq)] // Derive eq? +pub enum GetWindowBuilder<'a, F> +where + F: Real, +{ + /// [Boxcar] window, also known as a rectangular window or Dirichlet window; This is equivalent + /// to no window at all. + Boxcar, + /// [Triangle] window. + Triangle, + /// [Blackman] window. + Blackman, + /// [Hamming] window. + Hamming, + // Hann, + // Bartlett, + // Flattop, + // Parzen, + // Bohman, + // BlackmanHarris, + /// [Nuttall] window. + Nuttall, + // BartHann, + // Cosine, + // Exponential, + // Tukey, + // Taylor, + // Lanczos, + /// [Kaiser] window. + Kaiser { + /// Shape parameter `β`, please refer to [Kaiser]. + beta: F, + }, + // KaiserBesselDerived, // Needs Beta + // Gaussian, // Needs Standard Deviation + /// [GeneralCosine] window: Generic weighted sum of cosine term windows. + GeneralCosine { + /// Weighting Coefficients `a`, please refer to [GeneralCosine]. + weights: &'a [F], + }, + /// [GeneralGaussian] window: Generalized Gaussian Shape. + GeneralGaussian { + /// Shape parameter. + p: F, + /// The standard deviation, σ. + width: F, + }, + /// [GeneralHamming] window. + // Needs Window Coefficients. + GeneralHamming { + /// Window coefficient, ɑ + coefficient: F, + }, + // Dpss, // Needs Normalized Half-Bandwidth. + // Chebwin, // Needs Attenuation. +} + +/// Return a window of a given length and type. +/// +/// Parameters +/// ---------- +/// * `window`: [GetWindowBuilder] +/// The type of window to create. See below for more details. +/// * `Nx`: usize +/// The number of samples in the window. +/// * `fftbins`: bool, optional +/// If True (default), create a "periodic" window, ready to use with `ifftshift` and be +/// multiplied by the result of an FFT (see also :func:`~scipy.fft.fftfreq`). +/// If False, create a "symmetric" window, for use in filter design. +/// +/// Returns +/// ------- +/// * `get_window` : ndarray +/// Returns a window of length `Nx` and type `window` +/// +/// Notes +/// ----- +/// Window types: +/// * [Boxcar] +/// * [Triangle] +/// * [Blackman] +/// * [Hamming] +// Hann, +// Bartlett, +// Flattop, +// Parzen, +// Bohman, +// BlackmanHarris, +/// * [Nuttall] +// BartHann, +// Cosine, +// Exponential, +// Tukey, +// Taylor, +// Lanczos, +/// * [Kaiser] // Needs Beta +// KaiserBesselDerived, // Needs Beta +// Gaussian, // Needs Standard Deviation +/// * [GeneralCosine] +/// * [GeneralGaussian] // Needs Power, Width +/// * [GeneralHamming] // Needs Window Coefficients. +// Dpss, // Needs Normalized Half-Bandwidth. +// Chebwin, // Needs Attenuation. +/// # References +/// +pub fn get_window(window: GetWindowBuilder<'_, F>, nx: usize, fftbins: Option) -> Window +where + F: Real, +{ + match window { + GetWindowBuilder::Boxcar => Window::Boxcar(Boxcar { + m: nx, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::Triangle => Window::Triangle(Triangle { + m: nx, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::Blackman => Window::Blackman(Blackman { + m: nx, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::Hamming => Window::Hamming(Hamming { + m: nx, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::Nuttall => Window::Nuttall(Nuttall { + m: nx, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::Kaiser { beta } => Window::Kaiser(Kaiser { + m: nx, + beta, + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::GeneralCosine { weights } => Window::GeneralCosine(GeneralCosine { + m: nx, + a: weights.into(), + sym: !fftbins.unwrap_or(true), + }), + GetWindowBuilder::GeneralGaussian { p, width } => { + Window::GeneralGaussian(GeneralGaussian { + m: nx, + p, + sigma: width, + sym: !fftbins.unwrap_or(true), + }) + } + GetWindowBuilder::GeneralHamming { coefficient } => { + Window::GeneralHamming(GeneralHamming { + m: nx, + alpha: coefficient, + sym: !fftbins.unwrap_or(true), + }) + } + } +} From fe8c79ad255c5b276efb9717c15fd2c363c69998 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Fri, 6 Dec 2024 22:27:07 +0800 Subject: [PATCH 18/27] Add firwin There still is a need for double specification of window size across both firwin and the Window struct that uses it, unless we have a string to enum converter I suppose...; Either way, its suboptimal currently. --- sci-rs/src/signal/filter/design/firwin.rs | 857 ++++++++++++++++++++++ sci-rs/src/signal/filter/design/mod.rs | 2 + sci-rs/src/signal/windows/mod.rs | 22 + 3 files changed, 881 insertions(+) create mode 100644 sci-rs/src/signal/filter/design/firwin.rs diff --git a/sci-rs/src/signal/filter/design/firwin.rs b/sci-rs/src/signal/filter/design/firwin.rs new file mode 100644 index 00000000..1b490aa7 --- /dev/null +++ b/sci-rs/src/signal/filter/design/firwin.rs @@ -0,0 +1,857 @@ +use super::filter_type::FilterBandType; +use super::iirfilter_dyn; +use super::{kaiser_atten, kaiser_beta}; +use crate::signal::windows::get_window; +use crate::{error, error::Error, special}; +use core::cmp::Ord; +use nalgebra::RealField; +use num_traits::{real::Real, Float, MulAdd, Pow}; + +#[cfg(feature = "alloc")] +use crate::signal::{ + windows, + windows::{GetWindow, GetWindowBuilder}, + windows::{Hamming, Kaiser}, +}; +#[cfg(feature = "alloc")] +use alloc::{string::String, vec, vec::Vec}; + +/// Validation of [firwin_dyn] input. +fn firwin_dyn_validate( + numtaps: &usize, + cutoff: &[F], + pass_zero: &FilterBandType, + width: &Option, + window: &Option<&impl GetWindow>, +) -> Result<(), Error> { + if cutoff.is_empty() { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "At least one cutoff frequency must be given.".into(), + }); + } + + // Whilst it may be faster to write + // if *cutoff.iter().min().unwrap() <= F::zero() || *cutoff.iter().max().unwrap() >= F::one() { + // + // vec.min() requires the Ord trait, which is very difficult to impose onto f64 + let minimal = cutoff.iter().fold(F::max_value(), |acc, &m| acc.min(m)); + let maximal = cutoff.iter().fold(F::min_value(), |acc, &m| acc.max(m)); + if minimal <= F::zero() || maximal >= F::one() { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: + "Invalid cutoff frequency: frequencies must be greater than 0 and less than fs/2." + .into(), + }); + } + if cutoff.windows(2).any(|x| x[0] >= x[1]) { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "Invalid cutoff frequencies: the frequencies must be strictly increasing." + .into(), + }); + } + + match pass_zero { + FilterBandType::Lowpass => { + if cutoff.len() != 1 { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have one element if pass_zero is Lowpass.".into(), + }); + } + } + FilterBandType::Bandstop => { + if cutoff.len() < 2 { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have at least two elements if pass_zero is Bandstop." + .into(), + }); + } + } + FilterBandType::Highpass => { + if cutoff.len() != 1 { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have one element if pass_zero is Highpass.".into(), + }); + } + } + FilterBandType::Bandpass => { + if cutoff.len() < 2 { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have at least two elements if pass_zero is Bandpass." + .into(), + }); + } + } + } + + // While this was silently ignored in Scipy, we make this explicit here. + // Cannot use != on &impl GetWindow. + if window.is_some() && width.is_some() { + return Err(Error::InvalidArg { + arg: "cutoff".into(), + reason: "Setting both window and width to something is silently ignored only in Scipy." + .into(), + }); + } + + // This is here only because impl GetWindow on GetWindowBuilder is non-trivial to fetch numtaps + // outside the struct/enum-variants. + // if let Some(w) = *window { + // if w.m != numtaps { + // return Err(Error::ConfictArg { + // arg: "window".into(), + // reason: "Window has m value differing from numtaps", + // }); + // } + // } + Ok(()) +} + +/// FIR filter design using the window method. +/// +/// This function computes the coefficients of a finite impulse response filter. The filter will +/// have linear phase; it will be Type I if `numtaps` is odd and Type II if `numtaps` is even. +/// +/// Type II filters always have zero response at the Nyquist frequency, so a [error::Error] is +/// returned if firwin is called with `numtaps` even and having a passband whose right end is at +/// the Nyquist frequency. +/// +/// # Parameters +/// * `numtaps`: usize +/// Length of the filter (number of coefficients, i.e. the filter order + 1). `numtaps` must +/// be odd if a passband includes the Nyquist frequency. +/// * `cutoff`: 1-D array_like +/// Cutoff frequency of filter (expressed in the same units as `fs`) OR an array of cutoff +/// frequencies (that is, band edges). In the latter case, the frequencies in `cutoff` should +/// be positive and monotonically increasing between 0 and `fs/2`. The values 0 and `fs/2` must +/// not be included in `cutoff`. +/// * `width`: float or None, optional +/// If `width` is not None, then assume it is the approximate width of the transition region +/// (expressed in the same units as `fs`) for use in Kaiser FIR filter design. In this case, +/// the `window` argument is ignored. +/// * `window` : string or tuple of string and parameter values, optional +/// Desired window to use. See [GetWindow] for a list of windows and required parameters. +/// Defaults to Hamming. +/// Please set `sym=True` if you wish to have similar behaviour to scipy. +/// * `pass_zero` : [FilterBandType], optional +/// If True, the gain at the frequency 0 (i.e., the "DC gain") is 1. +/// If False, the DC gain is +/// 0. Can also be a string argument for the desired filter type (equivalent to ``btype`` in +/// [IIR design functions](iirfilter_dyn)). +/// +/// * `scale` : bool, optional +/// Set to True to scale the coefficients so that the frequency response is exactly unity at a +/// certain frequency. That frequency is either: +/// +/// * 0 (DC) if the first passband starts at 0 (i.e. pass_zero +/// is True) +/// * `fs/2` (the Nyquist frequency) if the first passband ends at +/// `fs/2` (i.e the filter is a single band highpass filter); +/// center of first passband otherwise +/// +/// * fs : float, optional +/// The sampling frequency of the signal. Each frequency in `cutoff` +/// must be between 0 and ``fs/2``. Default is 2. +/// +/// Returns +/// ------- +/// `h` : Result<`W`, _> +/// Coefficients of length `numtaps` FIR filter. +/// +/// Raises +/// ------ +/// Error : +/// If any value in `cutoff` is less than or equal to 0 or greater than or equal to ``fs/2``, +/// if the values in `cutoff` are not strictly monotonically increasing, or if `numtaps` is +/// even but a passband includes the Nyquist frequency. +/// +/// See Also +/// -------- +/// firwin2 +/// firls +/// minimum_phase +/// remez +/// +/// Examples +/// -------- +/// * Low-pass from 0 to f: +/// +/// ```custom,{class=language-python} +/// >>> from scipy import signal +/// >>> numtaps = 3 +/// >>> f = 0.1 +/// >>> signal.firwin(numtaps, f) +/// array([ 0.06799017, 0.86401967, 0.06799017]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx:: assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Lowpass, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![0.06799017, 0.86401967, 0.06799017]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// Use a specific window function: +/// +/// ```custom,{class=language-python} +/// >>> signal.firwin(numtaps, f, window='nuttall') +/// array([ 3.56607041e-04, 9.99286786e-01, 3.56607041e-04]) +/// ``` +/// Sci-rs: +// TODO: This still needs work so that its not Some(&Nuttall::new(...)) which might be mistakenly +// different from what's above it. +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Nuttall; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64], +/// None, +/// Some(&Nuttall::new(3, true)), +/// &FilterBandType::Lowpass, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![3.56607041e-04, 9.99286786e-01, 3.56607041e-04]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// High-pass ('stop' from 0 to f): +/// +/// ```custom,{class=language-python} +/// >>> signal.firwin(numtaps, f, pass_zero=False) +/// array([-0.00859313, 0.98281375, -0.00859313]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Highpass, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![-0.00859313, 0.98281375, -0.00859313]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// Band-pass: +/// +/// ```custom,{class=language-python} +/// >>> f1, f2 = 0.1, 0.2 +/// >>> signal.firwin(numtaps, [f1, f2], pass_zero=False) +/// array([ 0.06301614, 0.88770441, 0.06301614]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64, 0.2f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Bandpass, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![ 0.06301614, 0.88770441, 0.06301614]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// Band-stop: +/// +/// ```custom,{class=language-python} +/// >>> signal.firwin(numtaps, [f1, f2]) +/// array([-0.00801395, 1.0160279 , -0.00801395]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64, 0.2f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Bandstop, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![-0.00801395, 1.0160279 , -0.00801395]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// Multi-band (passbands are `[0, f1]`, `[f2, f3]` and `[f4, 1]`): +/// +/// ```custom,{class=language-python} +/// >>> f3, f4 = 0.3, 0.4 +/// >>> signal.firwin(numtaps, [f1, f2, f3, f4]) +/// array([-0.01376344, 1.02752689, -0.01376344]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64, 0.2f64, 0.3f64, 0.4f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Bandstop, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![-0.01376344, 1.02752689, -0.01376344]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// Multi-band (passbands are `[f1, f2]` and `[f3,f4]`): +/// +/// ```custom,{class=language-python} +/// >>> signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False) +/// array([ 0.04890915, 0.91284326, 0.04890915]) +/// ``` +/// Sci-rs: +/// ``` +/// use approx::assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::Hamming; +/// +/// let window: Vec = firwin_dyn( +/// 3, +/// &[0.1f64, 0.2f64, 0.3f64, 0.4f64], +/// None, +/// None::<&Hamming>, +/// &FilterBandType::Bandpass, +/// None, +/// None, +/// ) +/// .unwrap(); +/// let expected = vec![0.04890915, 0.91284326, 0.04890915]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +// In accordance with https://github.com/scipy/scipy/pull/16315, nyq as an argument should not be +// provided. +#[cfg(feature = "alloc")] +pub fn firwin_dyn( + numtaps: usize, + cutoff: &[F], + width: Option, + window: Option<&impl GetWindow>, + pass_zero: &FilterBandType, // Union with bools to follow scipy? + scale: Option, + fs: Option, +) -> Result, Error> +where + W: Real + Float + RealField + special::Bessel, + F: Real + PartialOrd + Float + RealField + MulAdd + Pow, +{ + let fs = match fs { + None => F::from(2).unwrap(), + Some(x) => x, + }; + let nyq = fs / F::from(2).unwrap(); + let cutoff: Vec<_> = cutoff.iter().map(|&c| c / nyq).collect(); + + firwin_dyn_validate(&numtaps, &cutoff, pass_zero, &width, &window)?; + + // Get and apply the window function. + let win: Vec = if let Some(x) = window { + // ?warn if window.sym != true + x.get_window() + } else if let Some(width) = width { + let atten = kaiser_atten(numtaps.try_into().unwrap(), width / nyq); + let beta = kaiser_beta(atten); + let k = get_window(GetWindowBuilder::Kaiser { beta }, numtaps, Some(false)); + k.get_window() + } else { + let h = get_window(GetWindowBuilder::::Hamming, numtaps, Some(false)); + h.get_window() + }; + + let pass_zero: bool = match pass_zero { + FilterBandType::Lowpass => true, + FilterBandType::Bandstop => true, + FilterBandType::Highpass => false, + FilterBandType::Bandpass => false, + }; + let pass_nyquist = (cutoff.len() % 2 == 1) ^ pass_zero; + if pass_nyquist && numtaps.is_multiple_of(2) { + return Err(Error::InvalidArg { + arg: "numtaps".into(), + reason: "A filter with an even number of coefficients must have zero response at the Nyquist frequency." + .into(), + }); + } + + let cutoff: Vec = { + let mut tmp = Vec::::new(); + if pass_zero { + tmp.push(F::zero()); + } + tmp.extend_from_slice(&cutoff); + if pass_nyquist { + tmp.push(F::one()); + } + tmp + }; + if !cutoff.len().is_multiple_of(2) { + unreachable!(); + // return Err(Error::InvalidArg { + // arg: "cutoff".into(), + // reason: "Parity of cutoff given for type of Filter is incorrect.".into(), + // }); + } + let bands: Vec<_> = cutoff.chunks_exact(2).collect(); + let scale_frequency = { + // for use in scale branch later + let left = bands[0][0]; + let right = bands[0][1]; + if left == F::zero() { + F::zero() + } else if right == F::one() { + F::one() + } else { + F::from(0.5).unwrap() * (left + right) + } + }; + + // Build up the coefficients. + // alpha: scalar = 0.5 * (numtaps - 1) + let alpha = F::from(0.5).unwrap() * F::from(numtaps - 1).unwrap(); + // m: 1Darray[numtaps] = 0..numtaps - alpha // lifetimes + // h: 1Darray[numtaps] = sum([right * sinc(right *m) - left * sinc(left * m)) for (left, right) + // in bands]) + let h: Vec = bands + .into_iter() + .map(|b| { + let left = b[0]; + let right = b[1]; + let m = (0..numtaps).map(|mi| F::from(mi).unwrap() - alpha); + let h: Vec = m + .map(|mi| { + if !mi.is_zero() { + right * (right * mi * F::pi()).sinc() - left * (left * mi * F::pi()).sinc() + } else { + right - left + } + }) + .collect(); + h + }) + .fold(vec![W::zero(); numtaps], |mut acc, x| { + acc.iter_mut() + .zip(x) + .map(|(a, xi)| *a += W::from(xi).unwrap()) + .collect::<()>(); + acc + }); + let mut h: Vec = h.into_iter().zip(win).map(|(hi, wi)| hi * wi).collect(); + + let scale = scale.unwrap_or(true); + if scale { + // m: 1Darray[numtaps] = 0..numtaps - alpha // lifetimes + let m = (0..numtaps).map(|mi| F::from(mi).unwrap() - alpha); + // c: 1Darray[numtaps] = np.cos(np.pi * m * scale_frequency) + let c = m + .map(|mi| F::pi() * mi * scale_frequency) + .map(|x| W::from(Real::cos(x)).unwrap()); + // s: scalar = np.sum(h*c) + let s: W = h + .iter() + .zip(c) + .map(|(&hi, ci)| hi * ci) + .fold(W::zero(), |acc, x| acc + x); + // h: 1Darray[numtaps] = h/s + h.iter_mut().map(|hi| *hi /= s).collect::<()>(); + } + + Ok(h) +} + +#[cfg(test)] +mod test { + use super::*; + use crate::signal::windows; + use approx::assert_abs_diff_eq; + + #[test] + fn invalid_args() { + type E = crate::error::Error; + { + let cutoff = Vec::::new(); + let empty_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Lowpass, + None, + None, + ); + assert_eq!( + empty_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "At least one cutoff frequency must be given.".into() + } + ); + } + { + let cutoff: Vec = vec![-0.1, 3.]; + let invalid_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandpass, + None, + None, + ); + assert_eq!( + invalid_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "Invalid cutoff frequency: frequencies must be greater than 0 and less than fs/2.".into() + } + ); + } + { + let cutoff: Vec = vec![0.2, 0.2]; + let decreasing_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandpass, + None, + None, + ); + assert_eq!( + decreasing_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: + "Invalid cutoff frequencies: the frequencies must be strictly increasing." + .into() + } + ); + } + { + // Lowpass + let cutoff: Vec = vec![0.2, 0.7]; + let lowpass_invalid_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Lowpass, + None, + None, + ); + assert_eq!( + lowpass_invalid_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have one element if pass_zero is Lowpass.".into() + } + ); + } + { + // Bandstop + let cutoff: Vec = vec![0.7]; + let bandstop_invalid_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandstop, + None, + None, + ); + assert_eq!( + bandstop_invalid_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have at least two elements if pass_zero is Bandstop." + .into() + } + ); + } + { + // Highpass + let cutoff: Vec = vec![0.2, 0.7]; + let highpass_invalid_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Highpass, + None, + None, + ); + assert_eq!( + highpass_invalid_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have one element if pass_zero is Highpass.".into() + } + ); + } + { + // Bandpass + let cutoff: Vec = vec![0.2]; + let bandpass_invalid_cutoff: Result, E> = firwin_dyn( + 3, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandpass, + None, + None, + ); + assert_eq!( + bandpass_invalid_cutoff.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "cutoff must have at least two elements if pass_zero is Bandpass." + .into() + } + ); + } + { + let cutoff: Vec = vec![0.2, 0.7]; + let invalid_numtaps: Result, E> = firwin_dyn( + 4, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandstop, + None, + None, + ); + assert_eq!( + invalid_numtaps.unwrap_err(), + E::InvalidArg { + arg: "numtaps".into(), + reason: "A filter with an even number of coefficients must have zero response at the Nyquist frequency.".into() + } + ); + } + } + + #[test] + fn conflicting_args() { + type E = crate::error::Error; + { + let cutoff: Vec = vec![0.2, 0.5, 0.7]; + let window = Hamming::new(12, true); + let conflicting: Result, E> = firwin_dyn( + 5, + &cutoff, + Some(0.3), + Some(&window), + &FilterBandType::Bandpass, + None, + None, + ); + assert_eq!( + conflicting.unwrap_err(), + E::InvalidArg { + arg: "cutoff".into(), + reason: "Setting both window and with to something is silently ignored only in Scipy." + .into() + } + ); + } + } + + #[test] + fn bandstop() { + // from scipy.signal import firwin + // firwin(numtaps=5, cutoff=[0.2, 0.7], width=None, # window=None, + // pass_zero='bandstop', scale=True, fs=None) + let cutoff: Vec = vec![0.2, 0.7]; + let window: Result, _> = firwin_dyn( + 5, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandstop, + None, + None, + ); + let expected = vec![0.05126868, -0.08050021, 1.05846306, -0.08050021, 0.05126868]; + assert_vec_eq(expected, window.unwrap()); + } + + #[test] + fn bandpass() { + // from scipy.signal import firwin + // firwin(numtaps=5, cutoff=[0.2, 0.7], width=None, # window=None, + // pass_zero='bandpass', scale=True, fs=None) + let cutoff: Vec = vec![0.2, 0.7]; + let window: Result, _> = firwin_dyn( + 5, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandpass, + None, + None, + ); + let expected = vec![-0.04340507, 0.06815306, 0.89611567, 0.06815306, -0.04340507]; + assert_vec_eq(expected, window.unwrap()); + } + + #[test] + fn lowpass() { + // from scipy.signal import firwin + // firwin(numtaps=5, cutoff=[0.2], width=None, # window=None, + // pass_zero='lowpass', scale=True, fs=None) + let cutoff: Vec = vec![0.2]; + let window: Result, _> = firwin_dyn( + 5, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Lowpass, + None, + None, + ); + let expected = vec![0.02840647, 0.23700821, 0.46917063, 0.23700821, 0.02840647]; + assert_vec_eq(expected, window.unwrap()); + } + + #[test] + fn highpass() { + // from scipy.signal import firwin + // firwin(numtaps=5, cutoff=[0.2], width=None, # window=None, + // pass_zero='highpass', scale=True, fs=None) + let cutoff: Vec = vec![0.2]; + let window: Result, _> = firwin_dyn( + 5, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Highpass, + None, + None, + ); + let expected = vec![-0.01238356, -0.1033217, 0.81812371, -0.1033217, -0.01238356]; + assert_vec_eq(expected, window.unwrap()); + } + + #[test] + fn different_fs() { + // from scipy.signal import firwin + // firwin(numtaps=5, cutoff=[2], width=None, # window=None, + // pass_zero='lowpass', scale=True, fs=10) + let cutoff: Vec = vec![2.]; + let window: Result, _> = firwin_dyn( + 5, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Lowpass, + None, + Some(10.), + ); + let expected = vec![0.01008727, 0.22034079, 0.53914388, 0.22034079, 0.01008727]; + assert_vec_eq(expected, window.unwrap()); + } + + #[track_caller] + fn assert_vec_eq(a: Vec, b: Vec) { + for (a, b) in a.into_iter().zip(b) { + assert_abs_diff_eq!(a, b, epsilon = 1e-6); + } + } +} diff --git a/sci-rs/src/signal/filter/design/mod.rs b/sci-rs/src/signal/filter/design/mod.rs index 5cdfebd5..192bba47 100644 --- a/sci-rs/src/signal/filter/design/mod.rs +++ b/sci-rs/src/signal/filter/design/mod.rs @@ -3,6 +3,7 @@ mod butter; mod cplx; mod filter_output; mod filter_type; +mod firwin; mod iirfilter; mod kaiser; mod lp2bp_zpk; @@ -19,6 +20,7 @@ pub use butter::*; use cplx::*; pub use filter_output::*; pub use filter_type::*; +pub use firwin::*; pub use iirfilter::*; pub use kaiser::*; pub use lp2bp_zpk::*; diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index e8b19ba8..5183d043 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -256,6 +256,28 @@ where /// * [GeneralHamming] // Needs Window Coefficients. // Dpss, // Needs Normalized Half-Bandwidth. // Chebwin, // Needs Attenuation. +/// +/// Examples +/// ----- +/// ``` +/// use approx:: assert_abs_diff_eq; +/// use sci_rs::signal::filter::design::{firwin_dyn, FilterBandType}; +/// use sci_rs::signal::windows::{get_window, GetWindow, GetWindowBuilder}; +/// +/// let window_struct = get_window(GetWindowBuilder::::Hamming, 3, None); +/// let window: Vec = window_struct.get_window(); +/// let expected = vec![0.08, 0.77, 0.77]; +/// +/// fn assert_vec_eq(a: Vec, b: Vec) { +/// for (a, b) in a.into_iter().zip(b) { +/// assert_abs_diff_eq!(a, b, epsilon = 1e-6); +/// } +/// } +/// +/// assert_vec_eq(window, expected); +/// ``` +/// +/// /// # References /// pub fn get_window(window: GetWindowBuilder<'_, F>, nx: usize, fftbins: Option) -> Window From bff6c9dc1dd0c081b7c95d62b9cb6bb0bfbb9c4c Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 10 Dec 2024 06:49:06 +0800 Subject: [PATCH 19/27] Fix lack of use of num_traits::identities Older commits didn't utilise some properities afforded by Real and RealField. Also added a useless(?) cfg(attr = "alloc") compiler attribute for consistency. --- sci-rs/src/signal/filter/design/kaiser.rs | 19 +++++++++---------- sci-rs/src/signal/windows/boxcar.rs | 4 ++-- sci-rs/src/signal/windows/general_hamming.rs | 8 ++------ sci-rs/src/signal/windows/kaiser.rs | 1 + sci-rs/src/signal/windows/triangle.rs | 4 ++-- 5 files changed, 16 insertions(+), 20 deletions(-) diff --git a/sci-rs/src/signal/filter/design/kaiser.rs b/sci-rs/src/signal/filter/design/kaiser.rs index 68b0e85b..c2c1fb43 100644 --- a/sci-rs/src/signal/filter/design/kaiser.rs +++ b/sci-rs/src/signal/filter/design/kaiser.rs @@ -1,5 +1,4 @@ -use core::f64::consts::PI; -use core::ops::Mul; +use nalgebra::RealField; use num_traits::{real::Real, MulAdd, Pow, ToPrimitive}; /// Compute the Kaiser parameter `beta`, given the attenuation `a`. @@ -43,8 +42,8 @@ where F::from(0.5842).unwrap(), F::from(0.07886).unwrap() * a, ) - } else if a > F::from(0).unwrap() { - F::from(0).unwrap() + } else if a > F::zero() { + F::zero() } else { panic!("Expected a positive input.") } @@ -84,11 +83,11 @@ where /// [kaiserord], [kaiser_beta] pub fn kaiser_atten(numtaps: u32, width: F) -> F where - F: Real + MulAdd, + F: Real + MulAdd + RealField, { MulAdd::mul_add( width, - F::from(numtaps - 1).unwrap() * F::from(2.285 * PI).unwrap(), + F::from(numtaps - 1).unwrap() * F::from(2.285).unwrap() * F::pi(), F::from(7.95).unwrap(), ) } @@ -217,16 +216,16 @@ where /// pub fn kaiserord(ripple: F, width: F) -> (F, F) where - F: Real + MulAdd + Pow, + F: Real + MulAdd + Pow + RealField, { - let a = ripple.abs(); + let a = Real::abs(ripple); if a < F::from(8).unwrap() { panic!("Requested maximum ripple attenuation is too small for the Kaiser formula."); } let beta = kaiser_beta(a); let numtaps = - F::from(1).unwrap() + (a - F::from(7.95).unwrap()) / (F::from(2.285 * PI).unwrap() * width); - (numtaps.ceil(), beta) + F::one() + (a - F::from(7.95).unwrap()) / (F::from(2.285).unwrap() * F::pi() * width); + (Real::ceil(numtaps), beta) } #[cfg(test)] diff --git a/sci-rs/src/signal/windows/boxcar.rs b/sci-rs/src/signal/windows/boxcar.rs index 0410fc90..5c5faaaa 100644 --- a/sci-rs/src/signal/windows/boxcar.rs +++ b/sci-rs/src/signal/windows/boxcar.rs @@ -62,9 +62,9 @@ where let (m, needs_trunc) = extend(self.m, self.sym); if !needs_trunc { - vec![W::from(1).unwrap(); m] + vec![W::one(); m] } else { - vec![W::from(1).unwrap(); m - 1] + vec![W::one(); m - 1] } } } diff --git a/sci-rs/src/signal/windows/general_hamming.rs b/sci-rs/src/signal/windows/general_hamming.rs index e3c93c9d..0ad2b356 100644 --- a/sci-rs/src/signal/windows/general_hamming.rs +++ b/sci-rs/src/signal/windows/general_hamming.rs @@ -145,12 +145,8 @@ where /// [5]: #references #[cfg(feature = "alloc")] fn get_window(&self) -> Vec { - GeneralCosine::new( - self.m, - [self.alpha, F::from(1).unwrap() - self.alpha].into(), - self.sym, - ) - .get_window() + GeneralCosine::new(self.m, [self.alpha, F::one() - self.alpha].into(), self.sym) + .get_window() } } diff --git a/sci-rs/src/signal/windows/kaiser.rs b/sci-rs/src/signal/windows/kaiser.rs index 2029f07e..e468fdcd 100644 --- a/sci-rs/src/signal/windows/kaiser.rs +++ b/sci-rs/src/signal/windows/kaiser.rs @@ -172,6 +172,7 @@ where /// [3]: #references /// [4]: #references /// [5]: #references + #[cfg(feature = "alloc")] fn get_window(&self) -> Vec { if len_guard(self.m) { return Vec::::new(); diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index a832fb72..894ad6d1 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -76,14 +76,14 @@ where let w: Vec = match m % 2 { 0 => { let mut w: Vec = n - .map(|n| (W::from(2).unwrap() * n - W::from(1).unwrap()) / m_f) + .map(|n| (W::from(2).unwrap() * n - W::one()) / m_f) .collect(); w.extend(w.clone().iter().rev()); w } 1 => { let mut w: Vec = n - .map(|n| W::from(2).unwrap() * n / (m_f + W::from(1).unwrap())) + .map(|n| W::from(2).unwrap() * n / (m_f + W::one())) .collect(); w.extend(w.clone().iter().rev().skip(1)); w From 6a456def1b11e41000561ffc7867b043a3d34811 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 15 Jul 2025 19:22:27 +0800 Subject: [PATCH 20/27] Fix clippy errors in Triangle --- sci-rs/src/signal/windows/triangle.rs | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index 894ad6d1..7a5c080b 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -71,7 +71,7 @@ where } let (m, needs_trunc) = extend(self.m, self.sym); - let mut n = (1..=((m + 1) / 2)).map(|x| W::from(x).unwrap()); + let mut n = (1..=(m.div_ceil(2))).map(|x| W::from(x).unwrap()); let m_f: W = W::from(m).unwrap(); let w: Vec = match m % 2 { 0 => { @@ -111,7 +111,6 @@ mod tests { let nx = 2 * i; let tri = Triangle::new(nx, true); let expected: Vec = (0..nx) - .into_iter() .chain((0..nx).rev()) .filter(|n| n % 2 == 1) .map(|n| n as f64 / nx as f64) @@ -136,7 +135,6 @@ mod tests { let nx = 2 * i; let tri = Triangle::new(nx, false); let expected: Vec = (1..=(nx / 2) + 1) - .into_iter() .chain((1..(nx / 2) + 1).rev()) .take(nx) .map(|n| n as f64 / ((nx / 2) + 1) as f64) @@ -160,9 +158,8 @@ mod tests { for i in 1..upper { let nx = 2 * i + 1; let tri = Triangle::new(nx, true); - let expected: Vec<_> = (1..=(nx + 1) / 2) - .into_iter() - .chain((1..=(nx + 1) / 2).rev().skip(1)) + let expected: Vec<_> = (1..=nx.div_ceil(2)) + .chain((1..=nx.div_ceil(2)).rev().skip(1)) .map(|n| (n as f64) / ((nx + 1) as f64 / 2.)) .collect(); @@ -185,8 +182,7 @@ mod tests { let nx = 2 * i + 1; let tri = Triangle::new(nx, false); let expected: Vec<_> = (0..=nx) - .into_iter() - .chain((0..=nx).into_iter().rev()) + .chain((0..=nx).rev()) .filter(|n| n % 2 == 1) .take(nx) .map(|n| n as f64 / (nx as f64 + 1.)) From 776075426410ff02851e5f9eb0c64f45d25030f4 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 23 Sep 2025 09:32:47 +0800 Subject: [PATCH 21/27] fix: Add zero numtap check in Firwin validation Scipy returns a value error in this case. --- sci-rs/src/signal/filter/design/firwin.rs | 28 ++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/sci-rs/src/signal/filter/design/firwin.rs b/sci-rs/src/signal/filter/design/firwin.rs index 1b490aa7..b9ccfe72 100644 --- a/sci-rs/src/signal/filter/design/firwin.rs +++ b/sci-rs/src/signal/filter/design/firwin.rs @@ -30,6 +30,12 @@ fn firwin_dyn_validate( reason: "At least one cutoff frequency must be given.".into(), }); } + if *numtaps == 0 { + return Err(Error::InvalidArg { + arg: "numtaps".into(), + reason: "Invalid numtaps: Nonzero-numtaps is expected!".into(), + }); + } // Whilst it may be faster to write // if *cutoff.iter().min().unwrap() <= F::zero() || *cutoff.iter().max().unwrap() >= F::one() { @@ -624,6 +630,26 @@ mod test { } ); } + { + // numtaps = 0 has ValueError in Python + let cutoff: Vec = vec![0.2, 0.2]; + let decreasing_cutoff: Result, E> = firwin_dyn( + 0, + &cutoff, + None, + None::<&Hamming>, + &FilterBandType::Bandpass, + None, + None, + ); + assert_eq!( + decreasing_cutoff.unwrap_err(), + E::InvalidArg { + arg: "numtaps".into(), + reason: "Invalid numtaps: Nonzero-numtaps is expected!".into() + } + ); + } { // Lowpass let cutoff: Vec = vec![0.2, 0.7]; @@ -746,7 +772,7 @@ mod test { conflicting.unwrap_err(), E::InvalidArg { arg: "cutoff".into(), - reason: "Setting both window and with to something is silently ignored only in Scipy." + reason: "Setting both window and width to something is silently ignored only in Scipy." .into() } ); From c300bc837028ba357c31f80f4dafd1407d888d2c Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Mon, 17 Nov 2025 22:05:10 +0800 Subject: [PATCH 22/27] Fix improper error message from window/width both being some Wrong argument was being signalled in the message. --- sci-rs/src/signal/filter/design/firwin.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sci-rs/src/signal/filter/design/firwin.rs b/sci-rs/src/signal/filter/design/firwin.rs index b9ccfe72..125f6820 100644 --- a/sci-rs/src/signal/filter/design/firwin.rs +++ b/sci-rs/src/signal/filter/design/firwin.rs @@ -100,7 +100,7 @@ fn firwin_dyn_validate( // Cannot use != on &impl GetWindow. if window.is_some() && width.is_some() { return Err(Error::InvalidArg { - arg: "cutoff".into(), + arg: "window/width".into(), reason: "Setting both window and width to something is silently ignored only in Scipy." .into(), }); @@ -771,7 +771,7 @@ mod test { assert_eq!( conflicting.unwrap_err(), E::InvalidArg { - arg: "cutoff".into(), + arg: "window/width".into(), reason: "Setting both window and width to something is silently ignored only in Scipy." .into() } From 50b02ee2734b2c1fa8dbd850052610c6f432d947 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 24 Sep 2025 10:52:45 +0800 Subject: [PATCH 23/27] feat: Add get_window macro Currently, it is still ordered really messily... Does not handle proper imports... Also does not scope into the appropriate namespace due to the way macros just are.... --- sci-rs/src/signal/windows/mod.rs | 141 +++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 5183d043..77918017 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -332,3 +332,144 @@ where } } } + +// Convenience macro for down-stream users +/// This provides the same convenience function as in the scipy.signal namespace. Returns a window +/// of given length and type. +/// +/// # Parameters: +/// - window: The type of window to create. +/// - m: Number of samples in the window. +/// - fftbins: If true (default), creates a "periodic" window, ready to use with ifftshift and be +/// multiplied by the result of an FFT. If False, create a "symmetric" window, for use in filter design. +// #[cfg(feature = "std")] +#[macro_export] +macro_rules! get_window { + // Parameterized windows (tuple form with sym) + ( ("kaiser", $beta:expr), $m:expr, $sym:expr ) => { + $crate::signal::windows::Kaiser::new($m, $beta, !$sym).get_window() + }; + // ( ("kaiser_bessel_derived", $beta:expr), $m:expr, $sym:expr ) => { + // $crate::signal::windows::KaiserBesselDerived::new($m, $beta, !$sym).get_window() + // }; + ( ("gaussian", $std:expr), $m:expr, $sym:expr ) => { + $crate::signal::windows::Gaussian::new($m, $std, $sym).get_window() + }; + ( ("general_cosine", $($coeff:expr),+), $m:expr, $sym:expr ) => { + $crate::signal::windows::GeneralCosine::new($m, &vec![$($coeff),+], $sym).get_window() + }; + ( ("general_gaussian", $power:expr, $width:expr), $m:expr, $sym:expr ) => { + $crate::signal::windows::GeneralGaussian::new($m, $power, $width, $sym).get_window() + }; + ( ("general_hamming", $alpha:expr), $m:expr, $sym:expr ) => { + $crate::signal::windows::GeneralHamming::new($m, $alpha, $sym).get_window() + }; + // ( ("dpss", $bandwidth:expr), $m:expr, $sym:expr ) => { + // $crate::signal::windows::Dpss::new($m, $bandwidth, $sym).get_window() + // }; + ( ("chebwin", $attenuation:expr), $m:expr, $sym:expr ) => { + $crate::signal::windows::Chebwin::new($m, $attenuation, $sym).get_window() + }; + + // Parameterized windows (tuple form without sym - defaults to true) + ( ("kaiser", $beta:expr), $m:expr ) => { + get_window!(("kaiser", $beta), $m, true) + }; + // ( ("kaiser_bessel_derived", $beta:expr), $m:expr ) => { + // get_window!(("kaiser_bessel_derived", $beta), $m, true) + // }; + // ( ("gaussian", $std:expr), $m:expr ) => { + // get_window!(("gaussian", $std), $m, true) + // }; + ( ("general_cosine", $($coeff:expr),+), $m:expr ) => { + get_window!(("general_cosine", $($coeff),+), $m, true) + }; + ( ("general_gaussian", $power:expr, $width:expr), $m:expr ) => { + get_window!(("general_gaussian", $power, $width), $m, true) + }; + ( ("general_hamming", $alpha:expr), $m:expr ) => { + get_window!(("general_hamming", $alpha), $m, true) + }; + // ( ("dpss", $bandwidth:expr), $m:expr ) => { + // get_window!(("dpss", $bandwidth), $m, true) + // }; + // ( ("chebwin", $attenuation:expr), $m:expr ) => { + // get_window!(("chebwin", $attenuation), $m, true) + // }; + + // Simple windows (string form with sym) + ("boxcar", $m:expr, $sym:expr) => { + $crate::signal::windows::Boxcar::new($m, $sym).get_window() + }; + ("triang", $m:expr, $sym:expr) => { + get_window!("triangle", $m, $sym) + }; + ("triangle", $m:expr, $sym:expr) => { + $crate::signal::windows::Triangle::new($m, $sym).get_window() + }; + ("blackman", $m:expr, $sym:expr) => { + $crate::signal::windows::Blackman::new($m, $sym).get_window() + }; + ("hamming", $m:expr, $sym:expr) => { + $crate::signal::windows::Hamming::new($m, $sym).get_window() + }; + // ("hann", $m:expr, $sym:expr) => { + // Hann::new($m, $sym).get_window() + // }; + // ("bartlett", $m:expr, $sym:expr) => { + // Bartlett::new($m, $sym).get_window() + // }; + // ("flattop", $m:expr, $sym:expr) => { + // FlatTop::new($m, $sym).get_window() + // }; + // ("parzen", $m:expr, $sym:expr) => { + // Parzen::new($m, $sym).get_window() + // }; + // ("bohman", $m:expr, $sym:expr) => { + // Bohman::new($m, $sym).get_window() + // }; + // ("blackmanharris", $m:expr, $sym:expr) => { + // BlackmanHarris::new($m, $sym).get_window() + // }; + ("nuttall", $m:expr, $sym:expr) => { + $crate::signal::windows::Nuttall::new($m, $sym).get_window() + }; + // ("barthann", $m:expr, $sym:expr) => { + // BartHann::new($m, $sym).get_window() + // }; + // ("cosine", $m:expr, $sym:expr) => { + // Cosine::new($m, $sym).get_window() + // }; + // ("exponential", $m:expr, $sym:expr) => { + // Exponential::new($m, $sym).get_window() + // }; + // ("tukey", $m:expr, $sym:expr) => { + // Tukey::new($m, $sym).get_window() + // }; + // ("taylor", $m:expr, $sym:expr) => { + // Taylor::new($m, $sym).get_window() + // }; + // ("lanczos", $m:expr, $sym:expr) => { + // Lanczos::new($m, $sym).get_window() + // }; + + // Simple windows (string form without sym - defaults to true) + ("boxcar", $m:expr) => { get_window!("boxcar", $m, true) }; + ("triang", $m:expr) => { get_window!("triang", $m, true) }; + ("triangle", $m:expr) => { get_window!("triang", $m, true) }; + ("blackman", $m:expr) => { get_window!("blackman", $m, true) }; + ("hamming", $m:expr) => { get_window!("hamming", $m, true) }; + // ("hann", $m:expr) => { get_window!("hann", $m, true) }; + // ("bartlett", $m:expr) => { get_window!("bartlett", $m, true) }; + // ("flattop", $m:expr) => { get_window!("flattop", $m, true) }; + // ("parzen", $m:expr) => { get_window!("parzen", $m, true) }; + // ("bohman", $m:expr) => { get_window!("bohman", $m, true) }; + // ("blackmanharris", $m:expr) => { get_window!("blackmanharris", $m, true) }; + ("nuttall", $m:expr) => { get_window!("nuttall", $m, true) }; + // ("barthann", $m:expr) => { get_window!("barthann", $m, true) }; + // ("cosine", $m:expr) => { get_window!("cosine", $m, true) }; + // ("exponential", $m:expr) => { get_window!("exponential", $m, true) }; + // ("tukey", $m:expr) => { get_window!("tukey", $m, true) }; + // ("taylor", $m:expr) => { get_window!("taylor", $m, true) }; + // ("lanczos", $m:expr) => { get_window!("lanczos", $m, true) }; +} From 09c93a4fdb4c73d0d5432b9659ac900ca929af49 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Thu, 25 Sep 2025 20:24:00 +0800 Subject: [PATCH 24/27] docs: Add macro use in Triangle Scoped behind standard as a doc test to demonstrate how to utilise the macro. Mandatory imports should be removed in the future as the macro is made more explicit. --- sci-rs/src/signal/windows/triangle.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index 7a5c080b..79a28afb 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -62,6 +62,16 @@ where /// assert_eq!(vec![0.1, 0.3, 0.5, 0.7, 0.9, 0.9, 0.7, 0.5, 0.3], tri.get_window()); /// ``` /// + /// Alternatively, this can be run with the convenience macro from [crate::signal] namespace. + /// (`std` feature is required for macros.) + /// #[cfg(feature = "std")] + /// ``` + /// use sci_rs::get_window; + /// use sci_rs::signal::windows::GetWindow; + /// let actual = get_window!("triangle", 8); + /// assert_eq!(vec![0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375, 0.125], actual); + /// ``` + /// /// # References /// #[cfg(feature = "alloc")] From ab8e89ff0e877016f74f94f9adbdd30c262bd4ac Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Thu, 25 Sep 2025 20:24:56 +0800 Subject: [PATCH 25/27] test: Add get_window macro test in Triangle This is to make sure the macro scoping works. --- sci-rs/src/signal/windows/triangle.rs | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index 79a28afb..82275653 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -205,4 +205,16 @@ mod tests { // } } } + + #[test] + #[cfg(feature = "std")] + fn from_macro() { + use crate::get_window; + + let actual = get_window!("triangle", 8); + assert_eq!( + vec![0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375, 0.125], + actual + ); + } } From 16152779737d7246dab82bbada7c4af7aba8f417 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Tue, 2 Dec 2025 16:33:05 +0800 Subject: [PATCH 26/27] Hack macro-scoping into crate::signal for `get_window` Utilise `#[doc(inline)]` and `#[doc(never)]` to hide name from crate root and show up in `::signal` namespace instead. --- sci-rs/src/signal/mod.rs | 12 ++++++++++++ sci-rs/src/signal/windows/mod.rs | 14 ++++---------- sci-rs/src/signal/windows/triangle.rs | 4 ++-- 3 files changed, 18 insertions(+), 12 deletions(-) diff --git a/sci-rs/src/signal/mod.rs b/sci-rs/src/signal/mod.rs index 5dd9a505..b9860ccb 100644 --- a/sci-rs/src/signal/mod.rs +++ b/sci-rs/src/signal/mod.rs @@ -24,6 +24,18 @@ pub mod convolve; /// namespace is located here. pub mod windows; +// Convenience macro for down-stream users +/// This provides the same convenience function as in the scipy.signal namespace. Returns a window +/// of given length and type. +/// +/// # Parameters: +/// - `window`: The type of window to create. +/// - `m`: Number of samples in the window. +/// - `fftbins`: If true (default), creates a "periodic" window, ready to use with ifftshift and be +/// multiplied by the result of an FFT. If False, create a "symmetric" window, for use in filter design. +#[doc(inline)] +pub use crate::_signal_windows_getWindow as get_window; + /// Signal Resampling /// This contains only the /// [`resample`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html#scipy.signal.resample) diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 77918017..53325c03 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -333,18 +333,12 @@ where } } -// Convenience macro for down-stream users -/// This provides the same convenience function as in the scipy.signal namespace. Returns a window -/// of given length and type. -/// -/// # Parameters: -/// - window: The type of window to create. -/// - m: Number of samples in the window. -/// - fftbins: If true (default), creates a "periodic" window, ready to use with ifftshift and be -/// multiplied by the result of an FFT. If False, create a "symmetric" window, for use in filter design. +/// Documentation is provided only at the public point of the interface. It is mangled here to +/// provide scoping. // #[cfg(feature = "std")] +#[doc(hidden)] #[macro_export] -macro_rules! get_window { +macro_rules! _signal_windows_getWindow { // Parameterized windows (tuple form with sym) ( ("kaiser", $beta:expr), $m:expr, $sym:expr ) => { $crate::signal::windows::Kaiser::new($m, $beta, !$sym).get_window() diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index 82275653..8dd778ce 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -66,7 +66,7 @@ where /// (`std` feature is required for macros.) /// #[cfg(feature = "std")] /// ``` - /// use sci_rs::get_window; + /// use sci_rs::signal::get_window; /// use sci_rs::signal::windows::GetWindow; /// let actual = get_window!("triangle", 8); /// assert_eq!(vec![0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375, 0.125], actual); @@ -209,7 +209,7 @@ mod tests { #[test] #[cfg(feature = "std")] fn from_macro() { - use crate::get_window; + use crate::signal::get_window; let actual = get_window!("triangle", 8); assert_eq!( From 94a951f9b2ad82f1cef57d6af38f42d7d061e004 Mon Sep 17 00:00:00 2001 From: SpookyYomo <48710653+SpookyYomo@users.noreply.github.com> Date: Wed, 3 Dec 2025 15:18:26 +0800 Subject: [PATCH 27/27] Remove need to import GetWindow trait for users of `get_window` Scope the import within the macro. --- sci-rs/src/signal/windows/mod.rs | 120 +++++++++++++++----------- sci-rs/src/signal/windows/triangle.rs | 1 - 2 files changed, 72 insertions(+), 49 deletions(-) diff --git a/sci-rs/src/signal/windows/mod.rs b/sci-rs/src/signal/windows/mod.rs index 53325c03..a0e6cb8e 100644 --- a/sci-rs/src/signal/windows/mod.rs +++ b/sci-rs/src/signal/windows/mod.rs @@ -340,30 +340,37 @@ where #[macro_export] macro_rules! _signal_windows_getWindow { // Parameterized windows (tuple form with sym) - ( ("kaiser", $beta:expr), $m:expr, $sym:expr ) => { + ( ("kaiser", $beta:expr), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Kaiser::new($m, $beta, !$sym).get_window() - }; + }}; // ( ("kaiser_bessel_derived", $beta:expr), $m:expr, $sym:expr ) => { // $crate::signal::windows::KaiserBesselDerived::new($m, $beta, !$sym).get_window() // }; - ( ("gaussian", $std:expr), $m:expr, $sym:expr ) => { + ( ("gaussian", $std:expr), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Gaussian::new($m, $std, $sym).get_window() - }; - ( ("general_cosine", $($coeff:expr),+), $m:expr, $sym:expr ) => { + }}; + ( ("general_cosine", $($coeff:expr),+), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::GeneralCosine::new($m, &vec![$($coeff),+], $sym).get_window() - }; - ( ("general_gaussian", $power:expr, $width:expr), $m:expr, $sym:expr ) => { + }}; + ( ("general_gaussian", $power:expr, $width:expr), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::GeneralGaussian::new($m, $power, $width, $sym).get_window() - }; - ( ("general_hamming", $alpha:expr), $m:expr, $sym:expr ) => { + }}; + ( ("general_hamming", $alpha:expr), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::GeneralHamming::new($m, $alpha, $sym).get_window() - }; - // ( ("dpss", $bandwidth:expr), $m:expr, $sym:expr ) => { + }}; + // ( ("dpss", $bandwidth:expr), $m:expr, $sym:expr ) => {{ + // use $crate::signal::windows::GetWindow; // $crate::signal::windows::Dpss::new($m, $bandwidth, $sym).get_window() - // }; - ( ("chebwin", $attenuation:expr), $m:expr, $sym:expr ) => { + // }}; + ( ("chebwin", $attenuation:expr), $m:expr, $sym:expr ) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Chebwin::new($m, $attenuation, $sym).get_window() - }; + }}; // Parameterized windows (tuple form without sym - defaults to true) ( ("kaiser", $beta:expr), $m:expr ) => { @@ -392,60 +399,77 @@ macro_rules! _signal_windows_getWindow { // }; // Simple windows (string form with sym) - ("boxcar", $m:expr, $sym:expr) => { + ("boxcar", $m:expr, $sym:expr) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Boxcar::new($m, $sym).get_window() - }; + }}; ("triang", $m:expr, $sym:expr) => { get_window!("triangle", $m, $sym) }; - ("triangle", $m:expr, $sym:expr) => { + ("triangle", $m:expr, $sym:expr) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Triangle::new($m, $sym).get_window() - }; - ("blackman", $m:expr, $sym:expr) => { + }}; + ("blackman", $m:expr, $sym:expr) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Blackman::new($m, $sym).get_window() - }; - ("hamming", $m:expr, $sym:expr) => { + }}; + ("hamming", $m:expr, $sym:expr) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Hamming::new($m, $sym).get_window() - }; - // ("hann", $m:expr, $sym:expr) => { + }}; + // ("hann", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Hann::new($m, $sym).get_window() - // }; - // ("bartlett", $m:expr, $sym:expr) => { + // }}; + // ("bartlett", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Bartlett::new($m, $sym).get_window() - // }; - // ("flattop", $m:expr, $sym:expr) => { + // }}; + // ("flattop", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // FlatTop::new($m, $sym).get_window() - // }; - // ("parzen", $m:expr, $sym:expr) => { + // }}; + // ("parzen", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Parzen::new($m, $sym).get_window() - // }; - // ("bohman", $m:expr, $sym:expr) => { + // }}; + // ("bohman", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Bohman::new($m, $sym).get_window() - // }; - // ("blackmanharris", $m:expr, $sym:expr) => { + // }}; + // ("blackmanharris", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // BlackmanHarris::new($m, $sym).get_window() - // }; - ("nuttall", $m:expr, $sym:expr) => { + // }}; + ("nuttall", $m:expr, $sym:expr) => {{ + use $crate::signal::windows::GetWindow; $crate::signal::windows::Nuttall::new($m, $sym).get_window() - }; - // ("barthann", $m:expr, $sym:expr) => { + }}; + // ("barthann", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // BartHann::new($m, $sym).get_window() - // }; - // ("cosine", $m:expr, $sym:expr) => { + // }}; + // ("cosine", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Cosine::new($m, $sym).get_window() - // }; - // ("exponential", $m:expr, $sym:expr) => { + // }}; + // ("exponential", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Exponential::new($m, $sym).get_window() - // }; - // ("tukey", $m:expr, $sym:expr) => { + // }}; + // ("tukey", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Tukey::new($m, $sym).get_window() - // }; - // ("taylor", $m:expr, $sym:expr) => { + // }}; + // ("taylor", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Taylor::new($m, $sym).get_window() - // }; - // ("lanczos", $m:expr, $sym:expr) => { + // }}; + // ("lanczos", $m:expr, $sym:expr) => {{ + // use $crate::signal::windows::GetWindow; // Lanczos::new($m, $sym).get_window() - // }; + // }}; // Simple windows (string form without sym - defaults to true) ("boxcar", $m:expr) => { get_window!("boxcar", $m, true) }; diff --git a/sci-rs/src/signal/windows/triangle.rs b/sci-rs/src/signal/windows/triangle.rs index 8dd778ce..e1d48037 100644 --- a/sci-rs/src/signal/windows/triangle.rs +++ b/sci-rs/src/signal/windows/triangle.rs @@ -67,7 +67,6 @@ where /// #[cfg(feature = "std")] /// ``` /// use sci_rs::signal::get_window; - /// use sci_rs::signal::windows::GetWindow; /// let actual = get_window!("triangle", 8); /// assert_eq!(vec![0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375, 0.125], actual); /// ```