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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

use serde::{Deserialize, Serialize};
use std::f32::consts::PI;

#[derive(Clone, PartialEq, Serialize, Deserialize)]
#[must_use]
pub struct WindowSinc {
    m: usize,
    fc: f32,
    bw: f32,
    taps: Vec<f32>,
    latency: usize,
}

impl WindowSinc {
    /// Creates a new [`WindowSinc`] instance.
    ///
    /// # Panics
    ///
    /// Panics if `cutoff` or `bandwidth` ratio to `sample_rate` is greater than `0.5`.
    pub fn new(sample_rate: f32, cutoff: f32, bandwidth: f32) -> Self {
        let fc = cutoff / sample_rate;
        let bw = bandwidth / sample_rate;
        assert!(
            (0.0..=0.5).contains(&fc),
            "cutoff frequency can not be greater than 1/2 the sampling rate: {cutoff} / {sample_rate}",
        );
        assert!(
            (0.0..=0.5).contains(&bw),
            "transition bandwidth can not be greater than 1/2 the sampling rate: {bandwidth} / {sample_rate}",
        );

        let m = (4.0 / bw) as usize; // Approximation
        let latency = m / 2; // Middle sample of FIR

        let mut h = Self::blackman_window(m);

        // Apply window sinc filter
        let p = 2.0 * PI * fc;
        for (i, h) in h.iter_mut().enumerate() {
            let i = i as f32 - latency as f32;
            *h *= if i == 0.0 { p } else { (p * i).sin() / i };
        }

        // Normalize
        let sum_inv = 1.0 / h.iter().sum::<f32>();
        for h in &mut h {
            *h *= sum_inv;
        }

        Self {
            m,
            fc,
            bw,
            taps: h,
            latency,
        }
    }

    fn blackman_window(m: usize) -> Vec<f32> {
        let p1 = 2.0 * PI / m as f32;
        let p2 = 4.0 * PI / m as f32;

        // Force N to be symmetrical
        let n = if m % 2 == 0 { m + 1 } else { m };
        let mut h = vec![0.0; n];

        for (i, h) in h.iter_mut().enumerate() {
            let i = i as f32;
            *h = 0.42 - 0.5 * (p1 * i).cos() + 0.8 * (p2 * i).cos();
        }

        h
    }

    #[inline]
    #[must_use]
    pub const fn taps(&self) -> &Vec<f32> {
        &self.taps
    }

    pub fn spectral_invert(&mut self) {
        let mut i = 1.0;
        for h in &mut self.taps {
            i *= -1.0;
            *h *= i;
        }
        self.taps[self.latency] += 1.0;
    }

    #[inline]
    #[must_use]
    pub fn len(&self) -> usize {
        self.taps.len()
    }

    #[inline]
    #[must_use]
    pub fn is_empty(&self) -> bool {
        self.taps.is_empty()
    }

    #[inline]
    #[must_use]
    pub const fn latency(&self) -> usize {
        self.latency
    }
}

impl std::fmt::Debug for WindowSinc {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("WindowSinc")
            .field("m", &self.m)
            .field("fc", &self.fc)
            .field("bw", &self.bw)
            .field("taps_len", &self.taps.len())
            .field("latency", &self.latency)
            .finish()
    }
}