idsp/
sweptsine.rs

1use core::iter::FusedIterator;
2
3use crate::{Complex, cossin::cossin};
4use num_traits::{FloatConst, real::Real};
5
6const Q: f32 = (1i64 << 32) as f32;
7
8/// Exponential sweep
9#[derive(Clone, Debug, PartialEq, PartialOrd)]
10pub struct Sweep {
11    /// Rate of exponential increase
12    pub rate: i32,
13    /// Current state
14    ///
15    /// Includes first order delta sigma modulator
16    pub state: i64,
17}
18
19/// Post-increment
20impl Iterator for Sweep {
21    type Item = i64;
22
23    #[inline]
24    fn next(&mut self) -> Option<Self::Item> {
25        const BIAS: i64 = 1 << 31;
26        let s = self.state;
27        match self
28            .state
29            .checked_add(self.rate as i64 * ((s + BIAS) >> 32))
30        {
31            Some(s) => self.state = s,
32            None => self.state = 0,
33        }
34        Some(s)
35    }
36}
37
38impl Sweep {
39    /// Create a new exponential sweep
40    #[inline]
41    pub fn new(rate: i32, state: i64) -> Self {
42        Self { rate, state }
43    }
44
45    /// Continuous time exponential sweep rate
46    #[inline]
47    pub fn rate(&self) -> f64 {
48        Real::ln_1p(self.rate as f64 / Q as f64)
49    }
50
51    /// Harmonic delay/length
52    #[inline]
53    pub fn delay(&self, harmonic: f64) -> f64 {
54        Real::ln(harmonic) / self.rate()
55    }
56
57    /// Samples per octave
58    #[inline]
59    pub fn octave(&self) -> f64 {
60        f64::LN_2() / self.rate()
61    }
62
63    /// Samples per decade
64    #[inline]
65    pub fn decade(&self) -> f64 {
66        f64::LN_10() / self.rate()
67    }
68
69    /// Current continuous time state
70    #[inline]
71    pub fn state(&self) -> f64 {
72        self.cycles() * self.rate()
73    }
74
75    /// Number of cycles per harmonic
76    #[inline]
77    pub fn cycles(&self) -> f64 {
78        self.state as f64 / (Q as f64 * self.rate as f64)
79    }
80
81    /// Evaluate integrated sweep at a given time
82    #[inline]
83    pub fn continuous(&self, t: f64) -> f64 {
84        self.cycles() * Real::exp(self.rate() * t)
85    }
86
87    /// Inverse filter
88    ///
89    /// * Stimulus `x(t)` (`AccuOsc::new(Sweep::fit(...)).collect()` plus windowing)
90    /// * Response `y(t)`
91    /// * Response FFT `Y(f)`
92    /// * Stimulus inverse filter `X'(f)`
93    /// * Transfer function `H(f) = X'(f) Y(f)'
94    /// * Impulse response `h(t)`
95    /// * Windowing each response using `delay()`
96    /// * Order responses `H_n(f)`
97    pub fn inverse_filter(&self, mut f: f32) -> Complex<f32> {
98        let rate = Real::ln_1p(self.rate as f32 / Q);
99        f /= rate;
100        let amp = 2.0 * rate * f.sqrt();
101        let inv_cycles = Q * self.rate as f32 / self.state as f32;
102        let turns = 0.125 - f * (1.0 - Real::ln(f * inv_cycles));
103        let (im, re) = Real::sin_cos(f32::TAU() * turns);
104        Complex::new(amp * re, amp * im)
105    }
106
107    /// Create new sweep
108    ///
109    /// * stop: maximum stop frequency in units of sample rate (e.g. Nyquist, 0.5)
110    /// * harmonics: number of harmonics to sweep
111    /// * cycles: number of cycles (phase wraps) per harmonic (>= 1)
112    pub fn fit(stop: f32, harmonics: f32, cycles: f32) -> Result<Self, SweepError> {
113        if !(0.0..=0.5).contains(&stop) {
114            return Err(SweepError::Stop);
115        }
116        let rate = Real::round(Q * Real::exp_m1(stop / (cycles * harmonics))) as i32;
117        let state = (rate as i64 * cycles as i64) << 32;
118        if state <= 0 {
119            return Err(SweepError::Start);
120        }
121        Ok(Self::new(rate, state))
122    }
123}
124
125/// Sweep parameter error
126#[derive(Clone, Copy, Debug, PartialEq, PartialOrd, thiserror::Error)]
127pub enum SweepError {
128    /// Start out of bounds
129    #[error("Start out of bounds")]
130    Start,
131    /// Stop out of bounds
132    #[error("Stop out of bounds")]
133    Stop,
134}
135
136/// Accumulating oscillator
137#[derive(Clone, Debug, PartialEq, PartialOrd)]
138pub struct AccuOsc<T> {
139    sweep: T,
140    /// Current phase
141    pub state: i64,
142}
143
144impl<T> AccuOsc<T> {
145    /// Create a new accumulating oscillator
146    #[inline]
147    pub fn new(sweep: T) -> Self {
148        Self { sweep, state: 0 }
149    }
150}
151
152/// Post-increment
153impl<T: Iterator<Item = i64>> Iterator for AccuOsc<T> {
154    type Item = Complex<i32>;
155
156    #[inline]
157    fn next(&mut self) -> Option<Self::Item> {
158        self.sweep.next().map(|f| {
159            let s = self.state;
160            self.state = s.wrapping_add(f);
161            let (re, im) = cossin((s >> 32) as _);
162            Complex::new(re, im)
163        })
164    }
165
166    #[inline]
167    fn size_hint(&self) -> (usize, Option<usize>) {
168        self.sweep.size_hint()
169    }
170}
171
172impl<T: FusedIterator + Iterator<Item = i64>> FusedIterator for AccuOsc<T> {}
173
174#[cfg(test)]
175mod test {
176    use super::*;
177    use crate::testing::*;
178
179    #[test]
180    fn test() {
181        let stop = 0.3;
182        let harmonics = 3000.0;
183        let cycles = 3.0;
184        let sweep = Sweep::fit(stop, harmonics, cycles).unwrap();
185        assert_eq!(sweep.rate, 0x22f40);
186        let len = sweep.delay(harmonics as _);
187        assert!(isclose(len, 240190.96, 0.0, 1e-2));
188        // Check API
189        assert!(isclose(sweep.cycles(), cycles as f64, 0.0, 1e-2));
190        assert_eq!(sweep.state(), sweep.continuous(0.0) * sweep.rate());
191        // Start/stop as desired
192        assert!((stop * 0.99..=1.01 * stop).contains(&(sweep.state() as f32 * harmonics)));
193        assert!(
194            (stop * 0.99..=1.01 * stop)
195                .contains(&((sweep.continuous(len as _) * sweep.rate()) as f32))
196        );
197        // Zero crossings and wraps
198        // Note inclusion of 0
199        for h in 0..harmonics as i32 {
200            let p = sweep.continuous(sweep.delay(h as _) as _);
201            assert!(isclose(p, h as f64 * cycles as f64, 0.0, 1e-10));
202        }
203        let sweep0 = sweep.clone();
204        for (t, p) in sweep
205            .scan(0i64, |p, f| {
206                let p0 = *p;
207                *p = p0.wrapping_add(f);
208                Some(p0 as f64 / Q.powi(2) as f64)
209            })
210            .take(len as _)
211            .enumerate()
212        {
213            // Analytic continuous time
214            let err = p - sweep0.continuous(t as _);
215            assert!(isclose(err - err.round(), 0.0, 0.0, 5e-5));
216        }
217    }
218}
idsp/sweptsine.rs

idsp/
sweptsine.rs
