sva_samples/measure/
loudness.rs

1// Concern: measures BS.1770 K-weighted loudness across a node's components | Non-concern: per-component level (envelope.rs), designing a biquad (biquad.rs) | IO: (planes, sample rate) -> Loudness
2
3use crate::biquad::{Coeffs, State};
4
5/// K-weighting as the analog prototypes BS.1770 bilinear-transforms, not the cookbook shapes:
6/// stating them as a design is what makes every rate but 48 kHz right too.
7const SHELF_HZ: f64 = 1681.97;
8const SHELF_Q: f64 = 0.70718;
9const SHELF_GAIN_DB: f64 = 3.99984;
10const SHELF_MID: f64 = 0.499_666_774_154_541_6;
11const HIGHPASS_HZ: f64 = 38.1355;
12const HIGHPASS_Q: f64 = 0.50033;
13
14fn shelf(sr: f64) -> Coeffs {
15    let k = (std::f64::consts::PI * SHELF_HZ / sr).tan();
16    let high = 10f64.powf(SHELF_GAIN_DB / 20.0);
17    let mid = high.powf(SHELF_MID);
18    let a0 = 1.0 + k / SHELF_Q + k * k;
19    Coeffs {
20        b0: (high + mid * k / SHELF_Q + k * k) / a0,
21        b1: 2.0 * (k * k - high) / a0,
22        b2: (high - mid * k / SHELF_Q + k * k) / a0,
23        a1: 2.0 * (k * k - 1.0) / a0,
24        a2: (1.0 - k / SHELF_Q + k * k) / a0,
25    }
26}
27
28fn highpass(sr: f64) -> Coeffs {
29    let k = (std::f64::consts::PI * HIGHPASS_HZ / sr).tan();
30    let a0 = 1.0 + k / HIGHPASS_Q + k * k;
31    Coeffs {
32        b0: 1.0,
33        b1: -2.0,
34        b2: 1.0,
35        a1: 2.0 * (k * k - 1.0) / a0,
36        a2: (1.0 - k / HIGHPASS_Q + k * k) / a0,
37    }
38}
39
40const OFFSET_DB: f64 = -0.691;
41const ABSOLUTE_GATE: f64 = -70.0;
42const INTEGRATED_GATE_LU: f64 = -10.0;
43
44const RANGE_GATE_LU: f64 = -20.0;
45
46const MOMENTARY_SECS: f64 = 0.4;
47const SHORT_TERM_SECS: f64 = 3.0;
48
49const STEP_SECS: f64 = 0.1;
50
51const PEAK_NOTE: &str = "sample peak, not true peak: no oversampling exists here, so an inter-sample peak above \
52     this figure is not measured";
53
54#[derive(Clone, Copy, Debug, PartialEq)]
55pub struct LoudnessFrame {
56    /// The block's START, the global seconds every other representation reports.
57    pub t: f64,
58    pub lufs: f64,
59}
60
61#[derive(Clone, Debug, PartialEq)]
62pub struct Loudness {
63    pub integrated_lufs: Option<f64>,
64    pub range_lu: Option<f64>,
65    pub momentary_max_lufs: Option<f64>,
66    pub short_term_max_lufs: Option<f64>,
67    pub sample_peak: f64,
68    pub sample_peak_dbfs: Option<f64>,
69    pub peak_note: &'static str,
70    pub momentary: Vec<LoudnessFrame>,
71    pub short_term: Vec<LoudnessFrame>,
72}
73
74/// Every component weighs 1.0: BS.1770's 1.41 surround lift needs a channel identity an
75/// indexed component does not carry.
76pub fn analyze(planes: &[&[f64]], sr: f64, start_secs: f64) -> Loudness {
77    let sample_peak = planes
78        .iter()
79        .flat_map(|p| p.iter())
80        .fold(0.0f64, |acc, v| acc.max(v.abs()));
81    let scale = headroom(sample_peak);
82    let squares: Vec<Vec<f64>> = planes
83        .iter()
84        .map(|p| running_squares(p, sr, scale))
85        .collect();
86    let samples = planes.first().map_or(0, |p| p.len());
87    let at = Blocks {
88        samples,
89        sr,
90        start_secs,
91        lift_db: 20.0 * scale.log10(),
92    };
93
94    let momentary = blocks(&squares, &at, MOMENTARY_SECS);
95    let short_term = blocks(&squares, &at, SHORT_TERM_SECS);
96
97    Loudness {
98        integrated_lufs: gated_mean(&momentary, INTEGRATED_GATE_LU),
99        range_lu: range(&short_term),
100        momentary_max_lufs: peak_of(&momentary),
101        short_term_max_lufs: peak_of(&short_term),
102        sample_peak,
103        sample_peak_dbfs: (sample_peak > 0.0).then(|| 20.0 * sample_peak.log10()),
104        peak_note: PEAK_NOTE,
105        momentary,
106        short_term,
107    }
108}
109
110/// A power of two over the peak: no square overflows, dividing is exact.
111fn headroom(peak: f64) -> f64 {
112    match peak > 1.0 {
113        true => 2f64.powi(peak.log2().ceil().min(f64::from(f64::MAX_EXP - 1)) as i32),
114        false => 1.0,
115    }
116}
117
118/// A prefix sum, so a block costs two reads rather than a pass of its own.
119fn running_squares(plane: &[f64], sr: f64, scale: f64) -> Vec<f64> {
120    let shelf = shelf(sr);
121    let cut = highpass(sr);
122    let mut first = State::default();
123    let mut second = State::default();
124    let mut out = Vec::with_capacity(plane.len() + 1);
125    let mut total = 0.0;
126    out.push(0.0);
127    for x in plane {
128        let k = second.step(&cut, first.step(&shelf, *x / scale));
129        total += k * k;
130        out.push(total);
131    }
132    out
133}
134
135struct Blocks {
136    samples: usize,
137    sr: f64,
138    start_secs: f64,
139    lift_db: f64,
140}
141
142fn blocks(squares: &[Vec<f64>], at: &Blocks, block_secs: f64) -> Vec<LoudnessFrame> {
143    let n = (block_secs * at.sr).round() as usize;
144    let step = ((STEP_SECS * at.sr).round() as usize).max(1);
145    if n == 0 || at.samples < n {
146        return Vec::new();
147    }
148    (0..=(at.samples - n))
149        .step_by(step)
150        .map(|k| LoudnessFrame {
151            t: at.start_secs + k as f64 / at.sr,
152            lufs: level(squares, k, n) + at.lift_db,
153        })
154        .collect()
155}
156
157fn level(squares: &[Vec<f64>], at: usize, n: usize) -> f64 {
158    let power: f64 = squares.iter().map(|s| (s[at + n] - s[at]) / n as f64).sum();
159    if power <= 0.0 {
160        return f64::NEG_INFINITY;
161    }
162    OFFSET_DB + 10.0 * power.log10()
163}
164
165/// The mean runs on POWER, not on the levels: averaging decibels is not averaging loudness.
166fn mean_above(frames: &[LoudnessFrame], threshold: f64) -> Option<f64> {
167    let kept: Vec<f64> = frames
168        .iter()
169        .map(|f| f.lufs)
170        .filter(|l| *l > threshold)
171        .collect();
172    if kept.is_empty() {
173        return None;
174    }
175    let from = kept.iter().copied().fold(OFFSET_DB, f64::max);
176    let power: f64 = kept.iter().map(|l| 10f64.powf((l - from) / 10.0)).sum();
177    Some(from + 10.0 * (power / kept.len() as f64).log10())
178}
179
180/// Both relative gates measure down from the ABSOLUTE-gated mean, never from each other.
181fn gated_mean(frames: &[LoudnessFrame], relative_lu: f64) -> Option<f64> {
182    mean_above(frames, mean_above(frames, ABSOLUTE_GATE)? + relative_lu)
183}
184
185/// EBU Tech 3342, whose gate is deliberately not the integrated measure's.
186fn range(short_term: &[LoudnessFrame]) -> Option<f64> {
187    let floor = mean_above(short_term, ABSOLUTE_GATE)? + RANGE_GATE_LU;
188    let mut kept: Vec<f64> = short_term
189        .iter()
190        .map(|f| f.lufs)
191        .filter(|l| *l > ABSOLUTE_GATE && *l > floor)
192        .collect();
193    if kept.is_empty() {
194        return None;
195    }
196    kept.sort_by(f64::total_cmp);
197    Some(percentile(&kept, 0.95) - percentile(&kept, 0.10))
198}
199
200fn percentile(sorted: &[f64], p: f64) -> f64 {
201    let at = ((sorted.len() as f64 - 1.0) * p).round() as usize;
202    sorted[at.min(sorted.len() - 1)]
203}
204
205fn peak_of(frames: &[LoudnessFrame]) -> Option<f64> {
206    frames
207        .iter()
208        .map(|f| f.lufs)
209        .filter(|l| l.is_finite())
210        .reduce(f64::max)
211}
sva_samples/measure/loudness.rs

sva_samples/measure/
loudness.rs
