oscillation/waveform/
sine.rs

1use num_traits::{Euclid, Float, FloatConst};
2
3use crate::Wavetable;
4
5use super::Waveform;
6
7#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, PartialOrd, Ord, serde::Serialize, serde::Deserialize)]
8pub struct Sine;
9
10impl<F> Waveform<F> for Sine
11where
12    F: Float + FloatConst + Euclid
13{
14    fn waveform(&self, theta: F) -> F
15    {
16        theta.cos()
17    }
18    fn waveform_with_dtc(&self, theta: F, duty_cycle: F) -> F
19    {
20        let half = crate::duty_cycle_default::<F>();
21
22        if duty_cycle == half
23        {
24            return self.waveform(theta)
25        }
26        let one = F::one();
27        let eps = F::from(EPS).unwrap();
28
29        let mut a = (one - duty_cycle - duty_cycle).clamp(-one, one);
30        let s = a.signum();
31        a = a.abs();
32        let p = (a + eps)/(a - one - eps);
33
34        let mut numer = (-p*s*theta.cos()).exp();
35        if !numer.is_finite()
36        {
37            numer = F::zero()
38        }
39
40        s*(p.tanh().recip() - numer/p.sinh())
41    }
42    fn wavetable<const N: usize>(&self) -> Option<Wavetable<F, N>>
43    {
44        None
45    }
46    fn wavetable_with_dtc<const N: usize>(&self, duty_cycle: F) -> Option<Wavetable<F, N>>
47    {
48        let half = crate::duty_cycle_default::<F>();
49
50        if duty_cycle == half
51        {
52            return None
53        }
54        let zero = F::zero();
55        let one = F::one();
56        let eps = F::from(EPS).unwrap();
57
58        let mut a = (one - duty_cycle - duty_cycle).clamp(-one, one);
59        let s = a.signum();
60        a = a.abs();
61        let p = (a + eps)/(a - one - eps);
62
63        #[repr(C)]
64        struct In<const N: usize>
65        {
66            i0: f64,
67            i_n: [f64; N]
68        }
69
70        let mut i_n = In {
71            i0: 0.0,
72            i_n: [0.0; N]
73        };
74        
75        {
76            let n = N.min(MAX_N);
77            let x = (-p*s).to_f64().unwrap().clamp(-MAX_I_N, MAX_I_N);
78            let i_n = unsafe {
79                core::slice::from_raw_parts_mut(&mut i_n.i0, n + 1)
80            };
81            rgsl::bessel::In_array(0, n as u32, x, i_n)
82                .expect("Bessel function failed.");
83        }
84
85        let g = p.sinh().recip();
86        let g2 = -s*(g + g);
87
88        Some(Wavetable::from_array(
89            s*(p.tanh().recip() - g*F::from(i_n.i0).unwrap()),
90            i_n.i_n.map(|i_n| (g2*F::from(i_n).unwrap(), zero))
91        ))
92    }
93}
94
95const MAX_N: usize = 64;
96const MAX_I_N: f64 = f32::MAX_EXP as f64;
97const EPS: f64 = 1.0/MAX_I_N;
98
99#[cfg(test)]
100mod test
101{
102    use core::error::Error;
103
104    use super::Sine;
105
106    #[test]
107    fn it_works() -> Result<(), Box<dyn Error>>
108    {
109        crate::tests::print_waveform(Sine)
110    }
111}
oscillation/waveform/sine.rs

oscillation/waveform/
sine.rs
