neco-stft 0.1.0

Backend-agnostic real FFT facade, windows, and STFT
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

use std::f64::consts::PI;

pub fn hann(n: usize) -> Vec<f64> {
    (0..n)
        .map(|i| 0.5 * (1.0 - (2.0 * PI * i as f64 / n as f64).cos()))
        .collect()
}

pub fn kaiser_bessel_derived(n: usize, alpha: f64) -> Vec<f64> {
    let half = n / 2;
    let kaiser = kaiser(half + 1, alpha);

    let mut cumsum = vec![0.0; half + 1];
    cumsum[0] = kaiser[0];
    for i in 1..=half {
        cumsum[i] = cumsum[i - 1] + kaiser[i];
    }
    let total = cumsum[half];

    let mut window = vec![0.0; n];
    for i in 0..half {
        window[i] = (cumsum[i] / total).sqrt();
    }
    for i in 0..half {
        window[n - 1 - i] = window[i];
    }

    window
}

fn kaiser(n: usize, alpha: f64) -> Vec<f64> {
    let m = (n - 1) as f64;
    let denom = bessel_i0(PI * alpha);
    (0..n)
        .map(|i| {
            let t = 2.0 * i as f64 / m - 1.0;
            bessel_i0(PI * alpha * (1.0 - t * t).sqrt()) / denom
        })
        .collect()
}

fn bessel_i0(x: f64) -> f64 {
    let mut sum = 1.0;
    let mut term = 1.0;
    let x_half = x / 2.0;
    for k in 1..=30 {
        term *= (x_half / k as f64) * (x_half / k as f64);
        sum += term;
        if term < 1e-16 * sum {
            break;
        }
    }
    sum
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn hann_endpoints_are_near_zero() {
        let window = hann(1024);
        assert!(window[0].abs() < 1e-15);
        assert!(window[1023] < 1e-3);
    }

    #[test]
    fn hann_squared_cola_for_quarter_hop() {
        let n = 1024;
        let hop = n / 4;
        let window = hann(n);
        let test_len = n * 4;
        let mut sum = vec![0.0; test_len];
        let mut pos = 0;
        while pos + n <= test_len {
            for i in 0..n {
                sum[pos + i] += window[i] * window[i];
            }
            pos += hop;
        }
        let steady_start = n;
        let steady_end = test_len - n;
        let ref_val = sum[steady_start];
        for value in &sum[steady_start..steady_end] {
            assert!((*value - ref_val).abs() < 1e-10);
        }
    }

    #[test]
    fn kbd_is_symmetric() {
        let window = kaiser_bessel_derived(2048, 4.0);
        for i in 0..1024 {
            assert!((window[i] - window[2047 - i]).abs() < 1e-12);
        }
    }
}