harmos_storage/signal/spectrum/
taper.rs

1/* ***
2 *
3 * Harmos
4 *
5 * SPDX-License-Identifier: Apache-2.0
6 *
7 * ***
8 */
9
10/* === Modules === */
11
12/* === Re-Exports === */
13
14/* === Imports === */
15
16// Standard Library Imports
17use std::f64::consts::PI;
18// Third Party Imports
19
20// Local Imports
21
22/* === Constants === */
23
24/* === Macros === */
25
26/* === Type Aliases === */
27
28/* === Structs === */
29
30/* === Enums === */
31
32#[derive(Clone, Copy, Debug, PartialEq, Eq)]
33pub enum Taper {
34    Rectangular,
35    Hamming,
36    Hann,
37    Blackman,
38}
39
40#[derive(Clone, Copy, Debug, PartialEq)]
41pub enum Band {
42    Rectangular { start: f64, end: f64 },
43    RaisedCosine { start: f64, end: f64, transition: f64 },
44    Kaiser { start: f64, end: f64, beta: f64 },
45    Gaussian { start: f64, end: f64, sigma: f64 },
46}
47
48/* === Traits === */
49
50/* === Implementations === */
51
52impl Taper {
53    #[must_use]
54    pub fn coefficient(&self, index: usize, length: usize) -> f64 {
55        let length = length as f64;
56        let index = index as f64;
57
58        match self {
59            Self::Rectangular => 1.0,
60            Self::Hamming => 0.54 - 0.46 * (2.0 * PI * index / length).cos(),
61            Self::Hann => 0.5 * (1.0 - (2.0 * PI * index / length).cos()),
62            Self::Blackman => 0.42 - 0.5 * (2.0 * PI * index / length).cos() + 0.08 * (4.0 * PI * index / length).cos(),
63        }
64    }
65
66    #[must_use]
67    pub fn coefficients(&self, length: usize) -> Vec<f64> {
68        (0..length).map(|index| self.coefficient(index, length)).collect()
69    }
70}
71
72impl Band {
73    #[must_use]
74    pub fn coefficient(&self, frequency: f64) -> f64 {
75        let (start, end) = self.bounds();
76        if frequency < start || frequency > end {
77            return 0.0;
78        }
79
80        match self {
81            Self::Rectangular { .. } => 1.0,
82            Self::RaisedCosine { transition, .. } => {
83                if frequency < start + transition {
84                    0.5 * (1.0 - (PI * (frequency - start) / transition).cos())
85                } else if frequency > end - transition {
86                    0.5 * (1.0 - (PI * (end - frequency) / transition).cos())
87                } else {
88                    1.0
89                }
90            }
91            Self::Kaiser { beta, .. } => {
92                let half = (end - start) / 2.0;
93                let offset = (frequency - (start + end) / 2.0) / half;
94                bessel_i0(beta * (1.0 - offset * offset).sqrt()) / bessel_i0(*beta)
95            }
96            Self::Gaussian { sigma, .. } => {
97                let offset = (frequency - (start + end) / 2.0) / sigma;
98                (-0.5 * offset * offset).exp()
99            }
100        }
101    }
102
103    #[must_use]
104    pub fn bounds(&self) -> (f64, f64) {
105        match self {
106            Self::Rectangular { start, end }
107            | Self::RaisedCosine { start, end, .. }
108            | Self::Kaiser { start, end, .. }
109            | Self::Gaussian { start, end, .. } => (*start, *end),
110        }
111    }
112}
113
114/* === Public Functions === */
115
116/* === Private Functions === */
117
118fn bessel_i0(argument: f64) -> f64 {
119    let magnitude = argument.abs();
120
121    if magnitude < 3.75 {
122        let square = (argument / 3.75).powi(2);
123        1.0 + square
124            * (3.515_622_9
125                + square * (3.089_942_4 + square * (1.206_749_2 + square * (0.265_973_2 + square * (0.036_076_8 + square * 0.004_581_3)))))
126    } else {
127        let ratio = 3.75 / magnitude;
128        (magnitude.exp() / magnitude.sqrt())
129            * (0.398_942_28
130                + ratio
131                    * (0.013_285_92
132                        + ratio
133                            * (0.002_253_19
134                                + ratio
135                                    * (-0.001_575_65
136                                        + ratio
137                                            * (0.009_162_81
138                                                + ratio
139                                                    * (-0.020_577_06
140                                                        + ratio * (0.026_355_37 + ratio * (-0.016_476_33 + ratio * 0.003_923_77))))))))
141    }
142}
143
144/* === Tests === */
harmos_storage/signal/spectrum/taper.rs

harmos_storage/signal/spectrum/
taper.rs
