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
14/// At or under it the shelf stands at or past Nyquist.
15pub const LEAST_RATE_HZ: f64 = 2.0 * SHELF_HZ;
16
17fn shelf(sr: f64) -> Coeffs {
18    let k = libm::tan(std::f64::consts::PI * SHELF_HZ / sr);
19    let high = libm::pow(10f64, SHELF_GAIN_DB / 20.0);
20    let mid = libm::pow(high, SHELF_MID);
21    let a0 = 1.0 + k / SHELF_Q + k * k;
22    Coeffs {
23        b0: (high + mid * k / SHELF_Q + k * k) / a0,
24        b1: 2.0 * (k * k - high) / a0,
25        b2: (high - mid * k / SHELF_Q + k * k) / a0,
26        a1: 2.0 * (k * k - 1.0) / a0,
27        a2: (1.0 - k / SHELF_Q + k * k) / a0,
28    }
29}
30
31fn highpass(sr: f64) -> Coeffs {
32    let k = libm::tan(std::f64::consts::PI * HIGHPASS_HZ / sr);
33    let a0 = 1.0 + k / HIGHPASS_Q + k * k;
34    Coeffs {
35        b0: 1.0,
36        b1: -2.0,
37        b2: 1.0,
38        a1: 2.0 * (k * k - 1.0) / a0,
39        a2: (1.0 - k / HIGHPASS_Q + k * k) / a0,
40    }
41}
42
43const OFFSET_DB: f64 = -0.691;
44const ABSOLUTE_GATE: f64 = -70.0;
45const INTEGRATED_GATE_LU: f64 = -10.0;
46
47const RANGE_GATE_LU: f64 = -20.0;
48
49pub const MOMENTARY_SECS: f64 = 0.4;
50const SHORT_TERM_SECS: f64 = 3.0;
51
52const STEP_SECS: f64 = 0.1;
53
54pub const PEAK_NOTE: &str = "sample peak, not true peak: no oversampling exists here, so an inter-sample peak above \
55     this figure is not measured";
56
57#[derive(Clone, Copy, Debug, PartialEq)]
58pub struct LoudnessFrame {
59    /// The block's START, the global seconds every other representation reports.
60    pub t: f64,
61    pub lufs: f64,
62}
63
64#[derive(Clone, Debug, PartialEq)]
65pub struct Loudness {
66    pub integrated_lufs: Option<f64>,
67    pub range_lu: Option<f64>,
68    pub momentary_max_lufs: Option<f64>,
69    pub short_term_max_lufs: Option<f64>,
70    pub sample_peak: f64,
71    pub sample_peak_dbfs: Option<f64>,
72    pub peak_note: &'static str,
73    pub momentary_secs: f64,
74    pub hop_secs: f64,
75    pub momentary: Vec<LoudnessFrame>,
76    pub short_term: Vec<LoudnessFrame>,
77}
78
79/// A rate with no K-weighting.
80#[derive(Clone, Copy, Debug, PartialEq)]
81pub struct RateUnderWeighting {
82    pub rate_hz: f64,
83    pub least_hz: f64,
84}
85
86/// So a quarter block is a whole step.
87pub const LEAST_BLOCK_SAMPLES: usize = 4;
88
89#[derive(Clone, Copy, Debug, PartialEq)]
90pub enum Unweighted {
91    Rate(RateUnderWeighting),
92    Block { window_secs: f64, samples: usize },
93}
94
95/// Every component weighs 1.0: BS.1770's 1.41 surround lift needs a channel identity an
96/// indexed component does not carry.
97pub fn analyze(
98    planes: &[&[f64]],
99    sr: f64,
100    start_secs: f64,
101    window_secs: f64,
102) -> Result<Loudness, Unweighted> {
103    if sr <= LEAST_RATE_HZ {
104        return Err(Unweighted::Rate(RateUnderWeighting {
105            rate_hz: sr,
106            least_hz: LEAST_RATE_HZ,
107        }));
108    }
109    let samples = (window_secs * sr).round();
110    if samples.is_nan() || samples < LEAST_BLOCK_SAMPLES as f64 {
111        return Err(Unweighted::Block {
112            window_secs,
113            samples: samples.max(0.0) as usize,
114        });
115    }
116    let (block, step) = (samples as usize, (samples / 4.0).round() as usize);
117    let sample_peak = planes
118        .iter()
119        .flat_map(|p| p.iter())
120        .fold(0.0f64, |acc, v| acc.max(v.abs()));
121    let scale = headroom(sample_peak);
122    let squares = weighted_squares(planes, sr, scale);
123    let at = Blocks {
124        sr,
125        start_secs,
126        lift_db: 20.0 * libm::log10(scale),
127    };
128
129    let seconds = |secs: f64| (secs * sr).round() as usize;
130    let momentary = blocks(&squares, &at, (block, step));
131    let short_term = blocks(
132        &squares,
133        &at,
134        (seconds(SHORT_TERM_SECS), seconds(STEP_SECS)),
135    );
136
137    Ok(Loudness {
138        integrated_lufs: gated_mean(&momentary, INTEGRATED_GATE_LU),
139        range_lu: range(&short_term),
140        momentary_max_lufs: peak_of(&momentary),
141        short_term_max_lufs: peak_of(&short_term),
142        sample_peak,
143        sample_peak_dbfs: (sample_peak > 0.0).then(|| 20.0 * libm::log10(sample_peak)),
144        peak_note: PEAK_NOTE,
145        momentary_secs: block as f64 / sr,
146        hop_secs: step as f64 / sr,
147        momentary,
148        short_term,
149    })
150}
151
152/// A power of two at the peak: dividing is exact; only 3000 dB under it does a square underflow.
153fn headroom(peak: f64) -> f64 {
154    match peak > 0.0 && peak.is_finite() {
155        true => {
156            let exponent = libm::log2(peak)
157                .ceil()
158                .clamp(f64::from(f64::MIN_EXP - 1), f64::from(f64::MAX_EXP - 1));
159            2f64.powi(exponent as i32)
160        }
161        false => 1.0,
162    }
163}
164
165fn weighted_squares(planes: &[&[f64]], sr: f64, scale: f64) -> Vec<f64> {
166    let samples = planes.first().map_or(0, |p| p.len());
167    let mut out = vec![0.0; samples];
168    let (shelf, cut) = (shelf(sr), highpass(sr));
169    for plane in planes {
170        let (mut first, mut second) = (State::default(), State::default());
171        for (sum, x) in out.iter_mut().zip(plane.iter()) {
172            let k = second.step(&cut, first.step(&shelf, *x / scale));
173            *sum += k * k;
174        }
175    }
176    out
177}
178
179struct Blocks {
180    sr: f64,
181    start_secs: f64,
182    lift_db: f64,
183}
184
185/// A block sums its own samples, piece by piece between edges, whatever loud stretch precedes it.
186fn blocks(squares: &[f64], at: &Blocks, (n, step): (usize, usize)) -> Vec<LoudnessFrame> {
187    let step = step.max(1);
188    if n == 0 || squares.len() < n {
189        return Vec::new();
190    }
191    let starts: Vec<usize> = (0..=(squares.len() - n)).step_by(step).collect();
192    let mut edges: Vec<usize> = starts.iter().flat_map(|k| [*k, k + n]).collect();
193    edges.sort_unstable();
194    edges.dedup();
195    let pieces: Vec<f64> = edges
196        .windows(2)
197        .map(|w| squares[w[0]..w[1]].iter().sum())
198        .collect();
199    let piece = |sample: usize| edges.binary_search(&sample).expect("a block edge");
200    starts
201        .iter()
202        .map(|&k| {
203            let power = pieces[piece(k)..piece(k + n)].iter().sum::<f64>() / n as f64;
204            LoudnessFrame {
205                t: at.start_secs + k as f64 / at.sr,
206                lufs: level(power) + at.lift_db,
207            }
208        })
209        .collect()
210}
211
212fn level(power: f64) -> f64 {
213    if power <= 0.0 {
214        return f64::NEG_INFINITY;
215    }
216    OFFSET_DB + 10.0 * libm::log10(power)
217}
218
219/// The mean runs on POWER, not on the levels: averaging decibels is not averaging loudness.
220fn mean_above(frames: &[LoudnessFrame], threshold: f64) -> Option<f64> {
221    let kept: Vec<f64> = frames
222        .iter()
223        .map(|f| f.lufs)
224        .filter(|l| *l > threshold)
225        .collect();
226    if kept.is_empty() {
227        return None;
228    }
229    let from = kept.iter().copied().fold(OFFSET_DB, f64::max);
230    let power: f64 = kept
231        .iter()
232        .map(|l| libm::pow(10f64, (l - from) / 10.0))
233        .sum();
234    Some(from + 10.0 * libm::log10(power / kept.len() as f64))
235}
236
237/// Both relative gates measure down from the ABSOLUTE-gated mean, never from each other.
238fn gated_mean(frames: &[LoudnessFrame], relative_lu: f64) -> Option<f64> {
239    mean_above(frames, mean_above(frames, ABSOLUTE_GATE)? + relative_lu)
240}
241
242/// EBU Tech 3342, whose gate is deliberately not the integrated measure's.
243fn range(short_term: &[LoudnessFrame]) -> Option<f64> {
244    let floor = mean_above(short_term, ABSOLUTE_GATE)? + RANGE_GATE_LU;
245    let mut kept: Vec<f64> = short_term
246        .iter()
247        .map(|f| f.lufs)
248        .filter(|l| *l > ABSOLUTE_GATE && *l > floor)
249        .collect();
250    if kept.is_empty() {
251        return None;
252    }
253    kept.sort_by(f64::total_cmp);
254    Some(percentile(&kept, 0.95) - percentile(&kept, 0.10))
255}
256
257fn percentile(sorted: &[f64], p: f64) -> f64 {
258    let at = ((sorted.len() as f64 - 1.0) * p).round() as usize;
259    sorted[at.min(sorted.len() - 1)]
260}
261
262fn peak_of(frames: &[LoudnessFrame]) -> Option<f64> {
263    frames
264        .iter()
265        .map(|f| f.lufs)
266        .filter(|l| l.is_finite())
267        .reduce(f64::max)
268}
sva_samples/measure/loudness.rs

sva_samples/measure/
loudness.rs
