sound_to_patch/dsp/
lomb.rs

1//! The Lomb-Scargle periodogram for samples at irregular times.
2//!
3//! Active note frames leave gaps, and a plain FFT of a gappy track sees the
4//! note rhythm as modulation. The periodogram fits a sine and a cosine at
5//! each trial frequency to the samples that exist.
6
7use std::f64::consts::TAU;
8
9use crate::cast;
10use crate::confidence::{ConfidenceError, Estimate};
11use crate::notes::NoteList;
12use crate::units::{Decibels, Hz, Seconds};
13
14use super::frames::ACTIVE_GATE_DB;
15use super::modulation::{Pick, Track, pick_rate, rate_estimate};
16
17/// The step of the trial frequencies.
18pub const GRID_STEP_HZ: Hz = Hz::from_const(0.01);
19/// The settling time after an onset whose samples are left out.
20pub const SETTLE_S: Seconds = Seconds::from_const(0.06);
21
22/// The trial frequencies from `lo` to `hi` in 0.01 Hz steps.
23#[must_use]
24pub fn grid(lo: Hz, hi: Hz) -> Vec<Hz> {
25    let first = cast::ceil_to_usize(lo.get() / GRID_STEP_HZ.get()).unwrap_or(0);
26    let last = cast::floor_to_usize(hi.get() / GRID_STEP_HZ.get() + 1e-9).unwrap_or(0);
27    (first..=last)
28        .map(|i| Hz::clamped(cast::usize_to_f64(i) * GRID_STEP_HZ.get()))
29        .collect()
30}
31
32/// The periodogram of `values` at `times`, one power per frequency.
33///
34/// Each trial fits a mean, a sine and a cosine to the samples together:
35/// the floating-mean form of Zechmeister and Kürster, the
36/// `floating_mean=True` option of `scipy.signal.lombscargle`. The plain form subtracts the overall mean first, and on
37/// gated notes a rate at the note rate then loses to its harmonics. The
38/// power is in the legacy units `A^2 N / 4` for a sine of amplitude `A`
39/// over `N` samples, and a constant offset changes nothing.
40#[must_use]
41pub fn lomb_scargle(times: &[Seconds], values: &[f64], frequencies: &[Hz]) -> Vec<f64> {
42    let len = times.len().min(values.len());
43    if len == 0 {
44        return vec![0.0; frequencies.len()];
45    }
46    let n = cast::usize_to_f64(len);
47    let (times, values) = (&times[..len], &values[..len]);
48    let mean_y = values.iter().sum::<f64>() / n;
49    let epsilon = f64::EPSILON / 2.0;
50    frequencies
51        .iter()
52        .map(|frequency| {
53            let omega = TAU * frequency.get();
54            let (mut c, mut s, mut cc, mut cs) = (0.0, 0.0, 0.0, 0.0);
55            for time in times {
56                let (sin, cos) = (omega * time.get()).sin_cos();
57                c += cos;
58                s += sin;
59                cc += cos * cos;
60                cs += cos * sin;
61            }
62            let (c, s) = (c / n, s / n);
63            let cc_full = cc / n - c * c;
64            let ss_full = 1.0 - cc / n - s * s;
65            let cs = cs / n - c * s;
66            let tau = 0.5 * (2.0 * cs).atan2(cc_full - ss_full);
67            let (mut yc, mut ys, mut c, mut s, mut cc) = (0.0, 0.0, 0.0, 0.0, 0.0);
68            for (time, y) in times.iter().zip(values) {
69                let (sin, cos) = (omega * time.get() - tau).sin_cos();
70                yc += y * cos;
71                ys += y * sin;
72                c += cos;
73                s += sin;
74                cc += cos * cos;
75            }
76            let (c, s) = (c / n, s / n);
77            let yc = yc / n - mean_y * c;
78            let ys = ys / n - mean_y * s;
79            let ss = (1.0 - cc / n - s * s).max(epsilon);
80            let cc = (cc / n - c * c).max(epsilon);
81            2.0 * (yc * yc / cc + ys * ys / ss) * n / 4.0
82        })
83        .collect()
84}
85
86/// The samples of a level track that a modulation estimate may use: level
87/// within 40 dB of the track's peak, inside a note, and at least `settle`
88/// after the note's onset.
89#[must_use]
90pub fn gated_samples(track: &Track, notes: &NoteList, settle: Seconds) -> (Vec<Seconds>, Vec<f64>) {
91    let peak = track.values().iter().copied().fold(0.0_f64, f64::max);
92    let gate = peak
93        * Decibels::clamped(-ACTIVE_GATE_DB.get())
94            .to_amplitude_gain()
95            .unwrap_or(0.0);
96    let rate = track.rate().get();
97    track
98        .values()
99        .iter()
100        .enumerate()
101        .filter_map(|(i, &value)| {
102            let time = Seconds::clamped(cast::usize_to_f64(i) / rate);
103            let settled = notes.notes().iter().any(|note| {
104                note.sounds_at(time) && time.get() - note.onset().get() >= settle.get()
105            });
106            (value > 0.0 && value >= gate && settled).then_some((time, value))
107        })
108        .unzip()
109}
110
111/// The modulation rate of irregular samples between `lo` and `hi`.
112///
113/// The periodogram amplitude goes through the same pick as the modulation
114/// spectrum. No estimate when no line stands out or the span holds fewer
115/// than two cycles. The spread is the resolution, `1 / T` of the span.
116///
117/// # Errors
118///
119/// Never in practice; the reasons are built from counts.
120pub fn periodogram_rate(
121    times: &[Seconds],
122    values: &[f64],
123    lo: Hz,
124    hi: Hz,
125) -> Result<Estimate<Hz, Hz>, ConfidenceError> {
126    let span = match (times.first(), times.last()) {
127        (Some(first), Some(last)) => Seconds::clamped(last.get() - first.get()),
128        _ => Seconds::clamped(0.0),
129    };
130    if span.get() <= 0.0 {
131        return rate_estimate(Pick::NoLine { prominence: 0.0 }, span);
132    }
133    // The log level, so a line's amplitude is its relative depth as in the
134    // modulation spectrum; `P = A^2 N / 4` gives `A = 2 sqrt(P / N)`.
135    let logs: Vec<f64> = values
136        .iter()
137        .map(|v| v.max(f64::MIN_POSITIVE).ln())
138        .collect();
139    let n = cast::usize_to_f64(times.len().min(values.len()));
140    // One guard step outside each end, see `pick_rate`.
141    let frequencies = grid(
142        Hz::clamped(lo.get() - GRID_STEP_HZ.get()),
143        Hz::clamped(hi.get() + GRID_STEP_HZ.get()),
144    );
145    let points: Vec<(f64, f64)> = frequencies
146        .iter()
147        .zip(lomb_scargle(times, &logs, &frequencies))
148        .map(|(f, p)| (f.get(), 2.0 * (p.max(0.0) / n).sqrt()))
149        .collect();
150    rate_estimate(pick_rate(&points), span)
151}
152
153#[cfg(test)]
154mod tests {
155    use super::*;
156    use crate::notes::Note;
157    use crate::units::MidiPitch;
158
159    #[test]
160    fn the_constants_pass_the_unit_check() {
161        assert_eq!(GRID_STEP_HZ, Hz::new(0.01).unwrap());
162        assert_eq!(SETTLE_S, Seconds::new(0.06).unwrap());
163    }
164
165    #[test]
166    fn the_grid_steps_by_one_hundredth_of_a_hertz() {
167        let g = grid(Hz::new(0.2).unwrap(), Hz::new(0.25).unwrap());
168        let values: Vec<f64> = g.iter().map(|f| f.get()).collect();
169        assert_eq!(values.len(), 6);
170        assert!((values[5] - 0.25).abs() < 1e-12);
171    }
172
173    #[test]
174    fn a_sine_has_power_a_squared_n_over_4_at_its_frequency() {
175        let times: Vec<Seconds> = (0..400)
176            .map(|i| Seconds::new(cast::usize_to_f64(i) / 200.0).unwrap())
177            .collect();
178        let values: Vec<f64> = times
179            .iter()
180            .map(|t| 2.0 * (TAU * 5.0 * t.get()).cos())
181            .collect();
182        let power = lomb_scargle(&times, &values, &[Hz::new(5.0).unwrap()]);
183        assert!((power[0] - 4.0 * 400.0 / 4.0).abs() < 1e-9, "{power:?}");
184    }
185
186    #[test]
187    fn gated_samples_skip_silence_gaps_and_the_settling_time() {
188        let track = Track::new(vec![1.0; 200], Hz::new(200.0).unwrap());
189        let note = Note::new(
190            Seconds::new(0.5).unwrap(),
191            Seconds::new(0.8).unwrap(),
192            MidiPitch::new(60.0).unwrap(),
193        )
194        .unwrap();
195        let (times, _) = gated_samples(&track, &NoteList::new(vec![note]), SETTLE_S);
196        assert_eq!(times.first().map(|t| t.get()), Some(0.56));
197        assert_eq!(times.len(), 48);
198    }
199}
sound_to_patch/dsp/lomb.rs

sound_to_patch/dsp/
lomb.rs
