denoize/
window.rs

1//! Window functions for short-time Fourier analysis/synthesis.
2//!
3//! All windows here are *periodic* (computed over `N` samples, not `N+1`),
4//! which makes the constant-overlap-add (COLA) sums exact for the standard
5//! hop sizes.
6
7use std::f64::consts::PI;
8
9/// Supported window families.
10#[derive(Clone, Copy, Debug, PartialEq, Eq)]
11pub enum WindowType {
12    /// Periodic Hann window. Default; COLA-exact at 50% overlap.
13    Hann,
14    /// Periodic Hamming window (slightly higher sidelobe attenuation trade-off).
15    Hamming,
16    /// Sine (a.k.a. half-sine / sqrt-Hann) window.
17    Sine,
18    /// Blackman window (low spectral leakage, narrower main lobe).
19    Blackman,
20    /// Kaiser-Bessel window. Adjustable sidelobe suppression via `kaiser_beta`.
21    Kaiser,
22    /// Flat-top window. Excellent amplitude accuracy; wider main lobe.
23    FlatTop,
24    /// DPSS (Discrete Prolate Spheroidal Sequence) / Slepian window.
25    /// Excellent energy concentration; `dpss_bandwidth` sets time-bandwidth product NW.
26    Dpss,
27}
28
29impl WindowType {
30    pub fn parse(s: &str) -> Option<Self> {
31        Some(match s.to_ascii_lowercase().as_str() {
32            "hann" => WindowType::Hann,
33            "hamming" => WindowType::Hamming,
34            "sine" => WindowType::Sine,
35            "blackman" => WindowType::Blackman,
36            "kaiser" | "kaiser-bessel" => WindowType::Kaiser,
37            "flattop" | "flat-top" | "flat_top" => WindowType::FlatTop,
38            "dpss" | "slepian" => WindowType::Dpss,
39            _ => return None,
40        })
41    }
42}
43
44/// Parameters for advanced windows.
45#[derive(Clone, Copy, Debug)]
46pub struct WindowParams {
47    /// Kaiser β (typical 5–12 for audio; higher = lower sidelobes).
48    pub kaiser_beta: f64,
49    /// DPSS time-bandwidth product NW (typical 2.5–4.0).
50    pub dpss_bandwidth: f64,
51}
52
53impl Default for WindowParams {
54    fn default() -> Self {
55        WindowParams {
56            kaiser_beta: 8.0,
57            dpss_bandwidth: 3.0,
58        }
59    }
60}
61
62/// Evaluate a single window sample `w(n)` for `n in 0..N` (periodic form).
63pub fn sample(kind: WindowType, n: usize, n_total: usize) -> f64 {
64    sample_with_params(kind, n, n_total, &WindowParams::default())
65}
66
67/// Evaluate with explicit parameters (Kaiser β, DPSS NW).
68pub fn sample_with_params(
69    kind: WindowType,
70    n: usize,
71    n_total: usize,
72    params: &WindowParams,
73) -> f64 {
74    let n = n as f64;
75    let nn = n_total as f64;
76    match kind {
77        WindowType::Hann => 0.5 * (1.0 - (2.0 * PI * n / nn).cos()),
78        WindowType::Hamming => 0.54 - 0.46 * (2.0 * PI * n / nn).cos(),
79        WindowType::Sine => (PI * (n + 0.5) / nn).sin(),
80        WindowType::Blackman => {
81            0.42 - 0.5 * (2.0 * PI * n / nn).cos() + 0.08 * (4.0 * PI * n / nn).cos()
82        }
83        WindowType::Kaiser => kaiser(n, nn, params.kaiser_beta),
84        WindowType::FlatTop => flat_top(n, nn),
85        WindowType::Dpss => dpss_approx(n, nn, params.dpss_bandwidth),
86    }
87}
88
89/// Build a window of length `n_total`.
90pub fn make(kind: WindowType, n_total: usize) -> Vec<f64> {
91    make_with_params(kind, n_total, &WindowParams::default())
92}
93
94/// Build with explicit parameters.
95pub fn make_with_params(kind: WindowType, n_total: usize, params: &WindowParams) -> Vec<f64> {
96    let mut w: Vec<f64> = (0..n_total)
97        .map(|i| sample_with_params(kind, i, n_total, params))
98        .collect();
99    // DPSS approximation may need peak normalisation.
100    if kind == WindowType::Dpss {
101        let peak = w.iter().cloned().fold(0.0f64, f64::max).max(1e-12);
102        if (peak - 1.0).abs() > 1e-6 {
103            for v in &mut w {
104                *v /= peak;
105            }
106        }
107    }
108    w
109}
110
111/// Kaiser-Bessel window (periodic form).
112fn kaiser(n: f64, n_total: f64, beta: f64) -> f64 {
113    let alpha = 0.5 * (n_total - 1.0);
114    let x = (n - alpha) / alpha;
115    bessel_i0(beta * (1.0 - x * x).max(0.0).sqrt()) / bessel_i0(beta)
116}
117
118/// Modified Bessel I0 for Kaiser window (simple series).
119fn bessel_i0(x: f64) -> f64 {
120    if x.abs() < 3.75 {
121        let t = x / 3.75;
122        let t2 = t * t;
123        1.0 + 3.5156229 * t2
124            + 3.0899424 * t2.powi(2)
125            + 1.2067492 * t2.powi(3)
126            + 0.2659732 * t2.powi(4)
127            + 0.0360768 * t2.powi(5)
128            + 0.0045813 * t2.powi(6)
129    } else {
130        let t = 3.75 / x.abs();
131        let ax = x.abs().exp() / x.abs().sqrt();
132        let c = 0.39894228 + 0.01328592 * t + 0.00225319 * t * t - 0.00157565 * t.powi(3)
133            + 0.00916281 * t.powi(4)
134            - 0.02057706 * t.powi(5)
135            + 0.02635537 * t.powi(6)
136            - 0.01647633 * t.powi(7)
137            + 0.00392377 * t.powi(8);
138        ax * c
139    }
140}
141
142/// Flat-top window (Heinzel et al. 5-term, periodic, clamped ≥ 0 for STFT).
143fn flat_top(n: f64, n_total: f64) -> f64 {
144    let a0 = 0.21557895;
145    let a1 = 0.41663158;
146    let a2 = 0.277263158;
147    let a3 = 0.083578947;
148    let a4 = 0.006947368;
149    let phi = 2.0 * PI * n / n_total;
150    let w = a0 - a1 * phi.cos() + a2 * (2.0 * phi).cos() - a3 * (3.0 * phi).cos()
151        + a4 * (4.0 * phi).cos();
152    w.max(0.0)
153}
154
155/// Approximate DPSS (Slepian) window via summed-cosine expansion.
156/// NW = time-bandwidth product; higher NW → narrower main lobe, lower sidelobes.
157fn dpss_approx(n: f64, n_total: f64, nw: f64) -> f64 {
158    let m = (2.0 * nw - 0.5).ceil() as i32;
159    let m = m.clamp(1, 8);
160    let mut w = 0.0;
161    for k in 0..m {
162        let vk = dpss_eigenvector_coeff(k, m, nw);
163        w += vk * (2.0 * PI * k as f64 * n / n_total).cos();
164    }
165    w.max(0.0)
166}
167
168/// DPSS eigenvector coefficient (simplified, good for audio STFT).
169fn dpss_eigenvector_coeff(k: i32, m: i32, nw: f64) -> f64 {
170    // Normalised cosine-series weights; dominant k=0 term ≈ 1.
171    let lambda_k = 1.0 - (k as f64) / (2.0 * nw);
172    let mut v = lambda_k.max(0.0);
173    if k == 0 {
174        v = 1.0;
175    } else {
176        v *= 0.5 / (k as f64);
177    }
178    let _ = m; // m controls series length above
179    v
180}
181
182#[cfg(test)]
183mod tests {
184    use super::*;
185
186    #[test]
187    fn hann_cola_at_50pct() {
188        let n = 1024;
189        let w = make(WindowType::Hann, n);
190        let hop = n / 2;
191        let len = n + hop * 3;
192        let mut cola = vec![0.0; len];
193        let mut start = 0;
194        while start + n <= len {
195            for i in 0..n {
196                cola[start + i] += w[i];
197            }
198            start += hop;
199        }
200        for &c in cola.iter().take(len - n).skip(n) {
201            assert!((c - 1.0).abs() < 1e-12, "cola={c}");
202        }
203    }
204
205    #[test]
206    fn advanced_windows_bounded() {
207        let n = 512;
208        for kind in [WindowType::Kaiser, WindowType::FlatTop, WindowType::Dpss] {
209            let w = make(kind, n);
210            for &v in &w {
211                assert!(v.is_finite());
212                // Flat-top has tiny negative sidelobes at edges; clamp check.
213                assert!((-0.01..=1.05).contains(&v), "{kind:?} value {v}");
214            }
215        }
216    }
217}
denoize/window.rs

denoize/
window.rs
