Add general gaussian window

sci-rs/src/signal/windows/general_gaussian.rs

@@ -0,0 +1,152 @@
+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 GeneralGaussian<A>
+where
+    A: Float,
+{
+    /// Number of points in the output window. If zero, an empty array is returned in [GetWindow].
+    pub m: usize,
+    /// Shape parameter.
+    pub p: A,
+    /// The standard deviation, σ.
+    pub sigma: A,
+    /// 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<A> GeneralGaussian<A>
+where
+    A: Float,
+{
+    /// 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: A, sigma: A, sym: bool) -> Self {
+        GeneralGaussian { m, p, sigma, sym }
+    }
+}
+
+#[cfg(feature = "alloc")]
+impl<A> GetWindow for GeneralGaussian<A>
+where
+    A: Real + Float,
+{
+    /// Return a window with a generalized gaussian shape.
+    ///
+    /// # Parameters
+    /// `self`: [GeneralGaussian]
+    ///
+    /// # Returns
+    /// `w`: `vec<F>`  
+    ///     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<f64> = 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
+    /// <https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.general_cosine.html>
+    #[cfg(feature = "alloc")]
+    fn get_window<F>(&self) -> Vec<F>
+    where
+        F: Real,
+    {
+        if len_guard(self.m) {
+            return Vec::<F>::new();
+        }
+        let (m, needs_trunc) = extend(self.m, self.sym);
+
+        let n = (0..m).map(|v| {
+            F::from(v).unwrap() - (F::from(m).unwrap() - F::from(1).unwrap()) / F::from(2).unwrap()
+        });
+        let sig = F::from(self.sigma).unwrap();
+        let two_p = F::from(self.p).unwrap() * F::from(2).unwrap();
+        let w = n
+            .into_iter()
+            .map(|v| {
+                (v / sig)
+                    .abs()
+                    .powf(two_p)
+                    .mul(F::from(-0.5).unwrap())
+                    .exp()
+            })
+            .collect();
+
+        return 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<f64>, b: Vec<f64>) {
+        for (a, b) in a.into_iter().zip(b) {
+            assert_abs_diff_eq!(a, b, epsilon = 1e-6);
+        }
+    }
+}

sci-rs/src/signal/windows/mod.rs

@@ -60,10 +60,12 @@ fn truncate<'a, W>(mut w: Vec<W>, needed: bool) -> Vec<W> {
 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.  
@@ -106,7 +108,9 @@ where
     /// Generic weighted sum of cosine term windows.
     // Needs Weighting Coefficients
     GeneralCosine(GeneralCosine<A>),
-    // GeneralGaussian, // Needs Power, Width
+    /// GeneralGaussian.
+    // Needs Power, Width
+    GeneralGaussian(GeneralGaussian<A>),
     // GeneralHamming, // Needs Window Coefficients.
     // Dpss, // Needs Normalized Half-Bandwidth.
     // Chebwin, // Needs Attenuation.