Skip to main content

sva_analysis/stable/
trajectory.rs

1// Concern: per-frame trajectory and a tolerance-based monotonicity verdict | Non-concern: computing level or width (envelope.rs/stereo.rs) | IO: (frames, samples) -> Trajectory
2
3use sva_samples::{EnvelopeFrame, StereoImage};
4
5use crate::decibels::to_db;
6use crate::frame::spectral_frames;
7
8/// Loudness's just-noticeable difference sits near 1 dB; half of that absorbs windowing and
9/// quantization jitter in an RMS trace without hiding a real, audible decay.
10pub const LEVEL_TOLERANCE_DB: f64 = 0.5;
11/// The same slack for a linear ratio (width) rather than a level.
12pub const RATIO_TOLERANCE: f64 = 0.05;
13
14#[derive(Clone, Copy, Debug, PartialEq)]
15pub struct TrajectoryFrame {
16    pub t_secs: f64,
17    pub rms_db: f64,
18    pub peak_db: f64,
19    pub width: Option<f64>,
20    pub centroid_hz: f64,
21    pub flatness: f64,
22}
23
24#[derive(Clone, Copy, Debug, PartialEq)]
25pub enum Direction {
26    Rising,
27    Falling,
28    Flat,
29}
30
31#[derive(Clone, Debug, PartialEq)]
32pub struct Verdict {
33    pub monotonic: bool,
34    pub direction: Direction,
35    pub violations: Vec<f64>,
36}
37
38#[derive(Clone, Debug, PartialEq)]
39pub struct Trajectory {
40    pub frames: Vec<TrajectoryFrame>,
41    pub level: Option<Verdict>,
42    pub width: Option<Verdict>,
43    pub centroid: Option<Verdict>,
44}
45
46pub fn analyze(
47    envelope: &[EnvelopeFrame],
48    stereo: Option<&StereoImage>,
49    samples: &[f32],
50    sample_rate: f64,
51    frame_secs: f64,
52) -> Trajectory {
53    let start_secs = envelope.first().map_or(0.0, |f| f.t_secs);
54    let spectral = spectral_frames(samples, sample_rate, start_secs, frame_secs, frame_secs);
55    let frames: Vec<TrajectoryFrame> = envelope
56        .iter()
57        .enumerate()
58        .map(|(i, e)| {
59            let (centroid_hz, flatness) = spectral
60                .get(i)
61                .map(|f| spectral_measures(&f.mags, f.bin_hz))
62                .unwrap_or((0.0, 0.0));
63            TrajectoryFrame {
64                t_secs: e.t_secs,
65                rms_db: to_db(e.rms),
66                peak_db: to_db(e.peak),
67                width: stereo.and_then(|s| s.frames.get(i)).map(|f| f.width),
68                centroid_hz,
69                flatness,
70            }
71        })
72        .collect();
73
74    let level = verdict(
75        &frames
76            .iter()
77            .map(|f| (f.t_secs, f.rms_db))
78            .collect::<Vec<_>>(),
79        Tolerance::Absolute(LEVEL_TOLERANCE_DB),
80    );
81    let width = stereo.and_then(|_| {
82        verdict(
83            &frames
84                .iter()
85                .filter_map(|f| f.width.map(|w| (f.t_secs, w)))
86                .collect::<Vec<_>>(),
87            Tolerance::Relative(RATIO_TOLERANCE),
88        )
89    });
90    let centroid = verdict(
91        &frames
92            .iter()
93            .map(|f| (f.t_secs, f.centroid_hz))
94            .collect::<Vec<_>>(),
95        Tolerance::Relative(RATIO_TOLERANCE),
96    );
97
98    Trajectory {
99        frames,
100        level,
101        width,
102        centroid,
103    }
104}
105
106/// Centroid as `spectrum::analyze` computes it; flatness is the Wiener entropy — the power
107/// spectrum's geometric mean over its arithmetic mean, 1.0 for white noise and near 0 for a
108/// single tone.
109fn spectral_measures(mags: &[f64], bin_hz: f64) -> (f64, f64) {
110    let power: Vec<f64> = mags.iter().map(|m| m * m).collect();
111    let total: f64 = power.iter().sum();
112    let centroid_hz = if total > 0.0 {
113        power
114            .iter()
115            .enumerate()
116            .map(|(k, p)| k as f64 * bin_hz * p)
117            .sum::<f64>()
118            / total
119    } else {
120        0.0
121    };
122    let floor = 1e-12;
123    let n = power.len().max(1) as f64;
124    let log_mean = power.iter().map(|p| (p + floor).ln()).sum::<f64>() / n;
125    let arithmetic_mean = (total + floor) / n;
126    let flatness = (log_mean.exp() / arithmetic_mean).clamp(0.0, 1.0);
127    (centroid_hz, flatness)
128}
129
130#[derive(Clone, Copy)]
131enum Tolerance {
132    Absolute(f64),
133    Relative(f64),
134}
135
136/// Tracks a running extreme in the claimed direction; a frame outside a slack band around it
137/// (additive in dB, multiplicative for a ratio) is a violation, and resets the extreme there so
138/// one real step does not keep re-triggering. A step of `-inf` to `-inf` — a window silent
139/// throughout, or none at all — is the one head-to-tail difference with no sign to read.
140fn verdict(points: &[(f64, f64)], tolerance: Tolerance) -> Option<Verdict> {
141    let span = (points.len() / 5).max(1).min(points.len());
142    let head = mean(&points[..span]);
143    let tail = mean(&points[points.len() - span..]);
144    if (tail - head).is_nan() {
145        return None;
146    }
147    let flat_band = match tolerance {
148        Tolerance::Absolute(t) => t,
149        Tolerance::Relative(t) => head.abs() * t,
150    };
151    let direction = if (tail - head).abs() <= flat_band {
152        Direction::Flat
153    } else if tail > head {
154        Direction::Rising
155    } else {
156        Direction::Falling
157    };
158
159    let rising = matches!(direction, Direction::Rising);
160    let mut extreme = points[0].1;
161    let mut violations = Vec::new();
162    for &(t, v) in &points[1..] {
163        let ok = match (direction, tolerance) {
164            (Direction::Flat, Tolerance::Absolute(tol)) => (v - points[0].1).abs() <= tol,
165            (Direction::Flat, Tolerance::Relative(tol)) => {
166                (v - points[0].1).abs() <= points[0].1.abs() * tol
167            }
168            (_, Tolerance::Absolute(tol)) if rising => v >= extreme - tol,
169            (_, Tolerance::Absolute(tol)) => v <= extreme + tol,
170            (_, Tolerance::Relative(tol)) if rising => v >= extreme * (1.0 - tol),
171            (_, Tolerance::Relative(tol)) => v <= extreme * (1.0 + tol),
172        };
173        if ok {
174            extreme = if rising {
175                extreme.max(v)
176            } else {
177                extreme.min(v)
178            };
179        } else {
180            violations.push(t);
181            extreme = v;
182        }
183    }
184    Some(Verdict {
185        monotonic: violations.is_empty(),
186        direction,
187        violations,
188    })
189}
190
191fn mean(points: &[(f64, f64)]) -> f64 {
192    points.iter().map(|&(_, v)| v).sum::<f64>() / points.len() as f64
193}
194
195#[cfg(test)]
196mod tests {
197    use super::*;
198
199    fn envelope(rms: &[f64]) -> Vec<EnvelopeFrame> {
200        rms.iter()
201            .enumerate()
202            .map(|(i, &r)| EnvelopeFrame {
203                t_secs: i as f64 * 0.01,
204                rms: r,
205                peak: r,
206            })
207            .collect()
208    }
209
210    #[test]
211    fn a_clean_decay_reads_monotonic_falling() {
212        let e = envelope(&[1.0, 0.5, 0.25, 0.125, 0.0625, 0.03125]);
213        let t = analyze(&e, None, &vec![0.0f32; 600], 44100.0, 0.01);
214        let level = t.level.expect("a sounding window has a level verdict");
215        assert_eq!(level.direction, Direction::Falling);
216        assert!(level.monotonic, "{:?}", level.violations);
217    }
218
219    /// The whole point of a tolerance band: real audio jitters by a fraction of a dB frame to
220    /// frame even while genuinely decaying, and a strict `<=` would flag every one of those.
221    #[test]
222    fn a_decay_with_natural_micro_variation_still_reads_monotonic() {
223        let mut rms = vec![1.0];
224        for i in 1..40 {
225            let trend = 1.0 * 0.9f64.powi(i);
226            let jitter = if i % 2 == 0 { 1.02 } else { 0.99 };
227            rms.push(trend * jitter);
228        }
229        let e = envelope(&rms);
230        let t = analyze(&e, None, &vec![0.0f32; 4000], 44100.0, 0.01);
231        let level = t.level.expect("a sounding window has a level verdict");
232        assert!(level.monotonic, "{:?}", level.violations);
233        assert_eq!(level.direction, Direction::Falling);
234    }
235
236    #[test]
237    fn a_real_swell_in_the_middle_of_a_decay_is_flagged() {
238        let mut rms = vec![1.0, 0.8, 0.6, 0.4];
239        rms.extend([0.9, 0.85]);
240        rms.extend([0.3, 0.2, 0.1, 0.05]);
241        let e = envelope(&rms);
242        let t = analyze(&e, None, &vec![0.0f32; 1000], 44100.0, 0.01);
243        let level = t.level.expect("a sounding window has a level verdict");
244        assert!(!level.monotonic);
245        assert!(!level.violations.is_empty());
246    }
247
248    #[test]
249    fn a_steady_level_reads_flat_not_falling() {
250        let e = envelope(&[0.5; 20]);
251        let t = analyze(&e, None, &vec![0.0f32; 2000], 44100.0, 0.01);
252        let level = t.level.expect("a sounding window has a level verdict");
253        assert_eq!(level.direction, Direction::Flat);
254        assert!(level.monotonic);
255    }
256
257    #[test]
258    fn width_is_absent_without_a_stereo_image() {
259        let e = envelope(&[0.5, 0.4]);
260        let t = analyze(&e, None, &vec![0.0f32; 200], 44100.0, 0.01);
261        assert!(t.width.is_none());
262    }
263}