scirs2_fft/sparse_fft/
windowing.rs

1//! Window function implementations for Sparse FFT
2//!
3//! This module provides various window functions that can be applied to signals
4//! before performing sparse FFT operations to reduce spectral leakage.
5
6use crate::error::{FFTError, FFTResult};
7use scirs2_core::numeric::Complex64;
8use scirs2_core::numeric::NumCast;
9use std::f64::consts::PI;
10use std::fmt::Debug;
11
12use super::config::WindowFunction;
13
14/// Apply a window function to the signal
15#[allow(dead_code)]
16pub fn apply_window<T>(
17    signal: &[T],
18    window_function: WindowFunction,
19    kaiser_beta: f64,
20) -> FFTResult<Vec<Complex64>>
21where
22    T: NumCast + Copy + Debug + 'static,
23{
24    // Convert input to complex
25    let signal_complex: Vec<Complex64> = signal
26        .iter()
27        .map(|&val| {
28            let val_f64 = NumCast::from(val)
29                .ok_or_else(|| FFTError::ValueError(format!("Could not convert {val:?} to f64")))?;
30            Ok(Complex64::new(val_f64, 0.0))
31        })
32        .collect::<FFTResult<Vec<_>>>()?;
33
34    let n = signal_complex.len();
35
36    // If no windowing is required, return the original signal
37    if window_function == WindowFunction::None {
38        return Ok(signal_complex);
39    }
40
41    // Apply the selected window _function
42    let windowed_signal = match window_function {
43        WindowFunction::None => signal_complex, // Already handled above, but included for completeness
44
45        WindowFunction::Hann => {
46            let mut windowed = signal_complex;
47            for (i, sample) in windowed.iter_mut().enumerate() {
48                let window_val = 0.5 * (1.0 - (2.0 * PI * i as f64 / (n - 1) as f64).cos());
49                *sample *= window_val;
50            }
51            windowed
52        }
53
54        WindowFunction::Hamming => {
55            let mut windowed = signal_complex;
56            for (i, sample) in windowed.iter_mut().enumerate() {
57                let window_val = 0.54 - 0.46 * (2.0 * PI * i as f64 / (n - 1) as f64).cos();
58                *sample *= window_val;
59            }
60            windowed
61        }
62
63        WindowFunction::Blackman => {
64            let mut windowed = signal_complex;
65            for (i, sample) in windowed.iter_mut().enumerate() {
66                let angle = 2.0 * PI * i as f64 / (n - 1) as f64;
67                let window_val = 0.42 - 0.5 * angle.cos() + 0.08 * (2.0 * angle).cos();
68                *sample *= window_val;
69            }
70            windowed
71        }
72
73        WindowFunction::FlatTop => {
74            let mut windowed = signal_complex;
75            for (i, sample) in windowed.iter_mut().enumerate() {
76                let angle = 2.0 * PI * i as f64 / (n - 1) as f64;
77                let window_val = 0.21557895 - 0.41663158 * angle.cos()
78                    + 0.277263158 * (2.0 * angle).cos()
79                    - 0.083578947 * (3.0 * angle).cos()
80                    + 0.006947368 * (4.0 * angle).cos();
81                *sample *= window_val;
82            }
83            windowed
84        }
85
86        WindowFunction::Kaiser => {
87            let mut windowed = signal_complex;
88            let _beta = kaiser_beta;
89            let alpha = (n - 1) as f64 / 2.0;
90            let i0_beta = modified_bessel_i0(_beta);
91
92            for (i, sample) in windowed.iter_mut().enumerate() {
93                let x = _beta * (1.0 - ((i as f64 - alpha) / alpha).powi(2)).sqrt();
94                let window_val = modified_bessel_i0(x) / i0_beta;
95                *sample *= window_val;
96            }
97            windowed
98        }
99    };
100
101    Ok(windowed_signal)
102}
103
104/// Modified Bessel function of the first kind, order 0
105/// Used for Kaiser window computation
106#[allow(dead_code)]
107fn modified_bessel_i0(x: f64) -> f64 {
108    let mut sum = 1.0;
109    let mut term = 1.0;
110    let half_x = x / 2.0;
111
112    for k in 1..=50 {
113        term *= (half_x / k as f64).powi(2);
114        sum += term;
115        if term < 1e-12 * sum {
116            break;
117        }
118    }
119
120    sum
121}
122
123#[cfg(test)]
124mod tests {
125    use super::*;
126
127    #[test]
128    fn test_apply_window_none() {
129        let signal = vec![1.0, 2.0, 3.0, 4.0];
130        let result = apply_window(&signal, WindowFunction::None, 14.0).unwrap();
131
132        assert_eq!(result.len(), 4);
133        assert_eq!(result[0], Complex64::new(1.0, 0.0));
134        assert_eq!(result[1], Complex64::new(2.0, 0.0));
135    }
136
137    #[test]
138    fn test_apply_window_hann() {
139        let signal = vec![1.0; 4];
140        let result = apply_window(&signal, WindowFunction::Hann, 14.0).unwrap();
141
142        assert_eq!(result.len(), 4);
143        // First and last samples should be zero for Hann window
144        assert!((result[0].re).abs() < 1e-10);
145        assert!((result[3].re).abs() < 1e-10);
146    }
147
148    #[test]
149    fn test_apply_window_hamming() {
150        let signal = vec![1.0; 4];
151        let result = apply_window(&signal, WindowFunction::Hamming, 14.0).unwrap();
152
153        assert_eq!(result.len(), 4);
154        // Hamming window should not be zero at endpoints
155        assert!(result[0].re > 0.0);
156        assert!(result[3].re > 0.0);
157    }
158
159    #[test]
160    fn test_modified_bessel_i0() {
161        // Test known values
162        let result = modified_bessel_i0(0.0);
163        assert!((result - 1.0).abs() < 1e-10);
164
165        let result = modified_bessel_i0(1.0);
166        assert!((result - 1.2660658777520084).abs() < 1e-10);
167    }
168}
scirs2_fft/sparse_fft/windowing.rs

scirs2_fft/sparse_fft/
windowing.rs
