embedded_dsp/
lut.rs

1//! Lookup-table trigonometry.
2//!
3//! Alternative to [`crate::fast_math`]'s polynomial approximations: a
4//! 256-entry sin table (1KB), computed at compile time via a Taylor-series
5//! `const fn` so no runtime initialization is needed. Trades accuracy
6//! (~2% error) for speed on cores without hardware sin/cos or an FPU.
7
8use core::f32::consts::PI;
9
10/// 256-entry sine lookup table covering 0 to 2*pi.
11/// Each entry is scaled by 32768 for fixed-point arithmetic.
12const SIN_TABLE: [i16; 256] = generate_sin_table();
13
14/// Generate the sine lookup table at compile time.
15const fn generate_sin_table() -> [i16; 256] {
16    let mut table = [0i16; 256];
17    let mut i = 0;
18    while i < 256 {
19        // sin() isn't const, so approximate with a Taylor series:
20        // sin(x) = x - x^3/3! + x^5/5! - x^7/7! + x^9/9! - x^11/11!
21        // The truncated series only converges well for |x| <= pi, so
22        // reduce the table angle into (-pi, pi] first -- without this,
23        // entries in the upper quarter of the table (close to 2*pi) blow
24        // up well past +/-1.0 and silently saturate when scaled by 32768.
25        let angle_rad = (i as f32 * 2.0 * PI) / 256.0;
26        let x = if angle_rad > PI {
27            angle_rad - 2.0 * PI
28        } else {
29            angle_rad
30        };
31        let x2 = x * x;
32        let x3 = x2 * x;
33        let x5 = x3 * x2;
34        let x7 = x5 * x2;
35        let x9 = x7 * x2;
36        let x11 = x9 * x2;
37
38        let sin_val = x - x3 / 6.0 + x5 / 120.0 - x7 / 5040.0 + x9 / 362880.0 - x11 / 39916800.0;
39
40        table[i] = (sin_val * 32768.0) as i16;
41        i += 1;
42    }
43    table
44}
45
46/// Fast sine using the lookup table.
47///
48/// Input: angle in radians (any range).
49/// Output: sine value scaled by 32768 (-32768 to 32767).
50///
51/// To convert to float: `result as f32 / 32768.0`.
52#[inline(always)]
53pub fn fast_sin_i16(angle: f32) -> i16 {
54    // Normalize angle to 0..2*pi and map to 0..255. `%` keeps the sign of
55    // the dividend, so negative angles land in (-1.0, 0.0) here; nudge
56    // those into [0.0, 1.0) before the (saturating) cast to usize, or
57    // they'd all collapse to index 0 instead of wrapping around the table.
58    let mut normalized = (angle % (2.0 * PI)) / (2.0 * PI);
59    if normalized < 0.0 {
60        normalized += 1.0;
61    }
62    let index = ((normalized * 256.0) as usize) & 0xFF;
63    SIN_TABLE[index]
64}
65
66/// Fast cosine using the lookup table (cos(x) = sin(x + pi/2)).
67#[inline(always)]
68pub fn fast_cos_i16(angle: f32) -> i16 {
69    fast_sin_i16(angle + PI / 2.0)
70}
71
72#[cfg(test)]
73mod tests {
74    extern crate std;
75    use super::*;
76
77    fn ref_sin(a: f32) -> f32 {
78        a.sin()
79    }
80    fn ref_cos(a: f32) -> f32 {
81        a.cos()
82    }
83
84    #[test]
85    fn sin_accuracy() {
86        let test_angles = [
87            0.0,
88            PI / 6.0,
89            PI / 4.0,
90            PI / 3.0,
91            PI / 2.0,
92            PI,
93            3.0 * PI / 2.0,
94        ];
95        for angle in test_angles.iter() {
96            let fast = fast_sin_i16(*angle) as f32 / 32768.0;
97            let accurate = ref_sin(*angle);
98            let error = (fast - accurate).abs();
99            assert!(error < 0.02, "error too large at angle {angle}: {error}");
100        }
101    }
102
103    #[test]
104    fn cos_accuracy() {
105        let test_angles = [0.0, PI / 6.0, PI / 4.0, PI / 3.0, PI / 2.0, PI];
106        for angle in test_angles.iter() {
107            let fast = fast_cos_i16(*angle) as f32 / 32768.0;
108            let accurate = ref_cos(*angle);
109            let error = (fast - accurate).abs();
110            assert!(error < 0.02, "error too large at angle {angle}: {error}");
111        }
112    }
113
114    #[test]
115    fn periodicity() {
116        let angle = PI / 4.0;
117        let sin1 = fast_sin_i16(angle);
118        let sin2 = fast_sin_i16(angle + 2.0 * PI);
119        let sin3 = fast_sin_i16(angle + 4.0 * PI);
120        assert!(
121            (sin1 - sin2).abs() <= 1000,
122            "periodicity error: {sin1} vs {sin2}"
123        );
124        assert!(
125            (sin1 - sin3).abs() <= 1000,
126            "periodicity error: {sin1} vs {sin3}"
127        );
128    }
129
130    #[test]
131    fn pythagorean_identity() {
132        let angles = [0.0, 0.5, 1.0, PI / 4.0, PI / 3.0, 2.0, 3.14];
133        for &a in &angles {
134            let s = fast_sin_i16(a) as f32 / 32768.0;
135            let c = fast_cos_i16(a) as f32 / 32768.0;
136            let pyth = s * s + c * c;
137            assert!(
138                (pyth - 1.0).abs() < 0.01,
139                "pythagorean identity failed at {a}: sin^2+cos^2 = {pyth}"
140            );
141        }
142    }
143
144    #[test]
145    fn trig_symmetry() {
146        let angles = [0.1, 0.785, 1.5, 2.3];
147        for &a in &angles {
148            let sin_pos = fast_sin_i16(a);
149            let sin_neg = fast_sin_i16(-a);
150            assert!(
151                (sin_pos + sin_neg).abs() <= 1000,
152                "sin odd symmetry failed at {a}: sin(a)={sin_pos}, sin(-a)={sin_neg}"
153            );
154
155            let cos_pos = fast_cos_i16(a);
156            let cos_neg = fast_cos_i16(-a);
157            assert!(
158                (cos_pos - cos_neg).abs() <= 1000,
159                "cos even symmetry failed at {a}: cos(a)={cos_pos}, cos(-a)={cos_neg}"
160            );
161        }
162    }
163}
embedded_dsp/lut.rs

embedded_dsp/
lut.rs
