1use super::resample::{bessel_i0, kaiser, sinc};
2use crate::media::{AudioBuffer, Error, Result};
3use rayon::prelude::*;
4use std::f64::consts::PI;
5
6const HIST_MIN: f64 = -70.0;
7const HIST_MAX: f64 = 20.0;
8const HIST_STEP: f64 = 0.01;
9const HIST_BINS: usize = ((HIST_MAX - HIST_MIN) / HIST_STEP) as usize;
10const SHORT: usize = 30;
11const MOMENTARY: usize = 4;
12pub(crate) const TAPS: usize = 32;
13const CUTOFF: f64 = 0.97;
14const BETA: f64 = 9.0;
15const CHUNK: usize = 1024;
16const PARALLEL: usize = 8192;
17
18pub(crate) fn lufs(e: f64) -> f64 {
19 if e <= 0.0 {
20 f64::NEG_INFINITY
21 } else {
22 -0.691 + 10.0 * e.log10()
23 }
24}
25
26#[derive(Debug, Clone, Copy)]
27struct Coeffs {
28 b0: f64,
29 b1: f64,
30 b2: f64,
31 a1: f64,
32 a2: f64,
33}
34
35fn prefilter(rate: f64) -> Coeffs {
36 let f0 = 1_681.974_450_955_533;
37 let gain_db = 3.999_843_853_973_347;
38 let q = 0.707_175_236_955_419_6;
39 let k = (PI * f0 / rate).tan();
40 let vh = 10f64.powf(gain_db / 20.0);
41 let vb = vh.powf(0.499_666_774_154_541_6);
42 let a0 = 1.0 + k / q + k * k;
43 Coeffs {
44 b0: (vh + vb * k / q + k * k) / a0,
45 b1: 2.0 * (k * k - vh) / a0,
46 b2: (vh - vb * k / q + k * k) / a0,
47 a1: 2.0 * (k * k - 1.0) / a0,
48 a2: (1.0 - k / q + k * k) / a0,
49 }
50}
51
52fn rlb(rate: f64) -> Coeffs {
53 let f0 = 38.135_470_876_024_44;
54 let q = 0.500_327_037_323_877_3;
55 let k = (PI * f0 / rate).tan();
56 let a0 = 1.0 + k / q + k * k;
57 Coeffs {
58 b0: 1.0,
59 b1: -2.0,
60 b2: 1.0,
61 a1: 2.0 * (k * k - 1.0) / a0,
62 a2: (1.0 - k / q + k * k) / a0,
63 }
64}
65
66fn weights(channels: usize) -> Vec<f64> {
67 match channels {
68 5 => vec![1.0, 1.0, 1.0, 1.41, 1.41],
69 6 => vec![1.0, 1.0, 1.0, 0.0, 1.41, 1.41],
70 8 => vec![1.0, 1.0, 1.0, 0.0, 1.41, 1.41, 1.41, 1.41],
71 n => vec![1.0; n],
72 }
73}
74
75pub(crate) fn oversampling(rate: f64) -> usize {
76 if rate < 88_200.0 {
77 4
78 } else if rate < 176_400.0 {
79 2
80 } else {
81 1
82 }
83}
84
85pub(crate) fn polyphase(factor: usize) -> Vec<[f32; TAPS]> {
86 if factor < 2 {
87 return Vec::new();
88 }
89 let half = TAPS as f64 / 2.0;
90 let i0 = bessel_i0(BETA);
91 (0..factor)
92 .map(|p| {
93 let mut w = [0.0f64; TAPS];
94 for (j, v) in w.iter_mut().enumerate() {
95 let t = j as f64 - (half - 1.0) - p as f64 / factor as f64;
96 *v = sinc(CUTOFF * t) * kaiser(t / half, BETA, i0);
97 }
98 let sum: f64 = w.iter().sum();
99 w.map(|v| (v / sum) as f32)
100 })
101 .collect()
102}
103
104#[derive(Debug, Clone)]
105struct Histogram {
106 count: Vec<u64>,
107 energy: Vec<f64>,
108}
109
110impl Histogram {
111 fn new() -> Self {
112 Histogram {
113 count: vec![0; HIST_BINS],
114 energy: vec![0.0; HIST_BINS],
115 }
116 }
117
118 fn bin(l: f64) -> usize {
119 (((l - HIST_MIN) / HIST_STEP).floor().max(0.0) as usize).min(HIST_BINS - 1)
120 }
121
122 fn centre(i: usize) -> f64 {
123 HIST_MIN + (i as f64 + 0.5) * HIST_STEP
124 }
125
126 fn add(&mut self, energy: f64) {
127 let l = lufs(energy);
128 if l >= HIST_MIN {
129 let b = Self::bin(l);
130 self.count[b] += 1;
131 self.energy[b] += energy;
132 }
133 }
134
135 fn mean_from(&self, from: usize) -> Option<f64> {
136 let e: f64 = self.energy[from..].iter().sum();
137 let n: u64 = self.count[from..].iter().sum();
138 (n > 0).then(|| e / n as f64)
139 }
140
141 fn gate(&self, relative: f64) -> Option<usize> {
142 let t = lufs(self.mean_from(0)?) + relative;
143 Some(if t < HIST_MIN { 0 } else { Self::bin(t) })
144 }
145}
146
147#[derive(Debug, Clone, Copy, PartialEq)]
148pub struct Loudness {
149 pub integrated: f64,
150 pub range: f64,
151 pub true_peak: f64,
152 pub sample_peak: f64,
153 pub momentary_max: f64,
154 pub short_term_max: f64,
155}
156
157#[derive(Debug, Clone)]
158struct Channel {
159 s: [f64; 4],
160 carry: f64,
161 sums: Vec<f64>,
162 buffer: Vec<f32>,
163 out: Vec<f32>,
164 sample_peak: f32,
165 true_peak: f32,
166}
167
168impl Channel {
169 fn new() -> Self {
170 Channel {
171 s: [0.0; 4],
172 carry: 0.0,
173 sums: Vec::new(),
174 buffer: vec![0.0; TAPS - 1],
175 out: vec![0.0; CHUNK],
176 sample_peak: 0.0,
177 true_peak: 0.0,
178 }
179 }
180
181 fn peaks(&mut self, bank: &[[f32; TAPS]]) {
182 let peak = peak(&self.buffer[TAPS - 1..]);
183 self.sample_peak = self.sample_peak.max(peak);
184 let tp = true_peak(&self.buffer, bank, &mut self.out);
185 self.true_peak = self.true_peak.max(tp).max(peak);
186 let n = self.buffer.len();
187 self.buffer.copy_within(n - (TAPS - 1).., 0);
188 self.buffer.truncate(TAPS - 1);
189 }
190}
191
192fn run(chans: &mut [Channel], k: &Filters, segments: &[usize], bank: &[[f32; TAPS]]) {
193 if chans.len() == 2 {
194 weigh::<2>(chans, k, segments);
195 } else {
196 weigh::<1>(chans, k, segments);
197 }
198 chans.iter_mut().for_each(|c| c.peaks(bank));
199}
200
201fn weigh<const N: usize>(chans: &mut [Channel], k: &Filters, segments: &[usize]) {
202 let mut s: [[f64; 4]; N] = std::array::from_fn(|i| chans[i].s);
203 let mut sums: [Vec<f64>; N] = std::array::from_fn(|i| std::mem::take(&mut chans[i].sums));
204 sums.iter_mut().for_each(Vec::clear);
205 let data: [&[f32]; N] = std::array::from_fn(|i| &chans[i].buffer[TAPS - 1..]);
206 let (pre, rlb) = (k.pre, k.rlb);
207 let mut at = 0;
208 for &len in segments {
209 let mut acc = [0.0f64; N];
210 let part: [&[f32]; N] = std::array::from_fn(|i| &data[i][at..at + len]);
211 (0..len).for_each(|t| {
212 for ((st, a), p) in s.iter_mut().zip(acc.iter_mut()).zip(&part) {
213 let x = p[t] as f64;
214 let [p1, p2, r1, r2] = st;
215 let y = pre.b0 * x + *p1;
216 *p1 = pre.b1 * x - pre.a1 * y + *p2;
217 *p2 = pre.b2 * x - pre.a2 * y;
218 let z = y + *r1;
219 *r1 = -2.0 * y - rlb.a1 * z + *r2;
220 *r2 = y - rlb.a2 * z;
221 *a += z * z;
222 }
223 });
224 for i in 0..N {
225 sums[i].push(acc[i]);
226 s[i] = s[i].map(|v| if v.abs() < 1e-30 { 0.0 } else { v });
227 }
228 at += len;
229 }
230 for (c, (st, sum)) in chans.iter_mut().zip(s.into_iter().zip(sums)) {
231 c.s = st;
232 c.sums = sum;
233 }
234}
235
236pub(crate) fn peak(x: &[f32]) -> f32 {
237 let mut acc = [0.0f32; 16];
238 let mut chunks = x.chunks_exact(16);
239 for c in &mut chunks {
240 for k in 0..16 {
241 acc[k] = acc[k].max(c[k].abs());
242 }
243 }
244 let tail = chunks
245 .remainder()
246 .iter()
247 .fold(0.0f32, |m, v| m.max(v.abs()));
248 acc.iter().fold(tail, |m, &v| m.max(v))
249}
250
251pub(crate) fn true_peak(buffer: &[f32], bank: &[[f32; TAPS]], scratch: &mut [f32]) -> f32 {
252 if bank.is_empty() {
253 return 0.0;
254 }
255 let n = buffer.len().saturating_sub(TAPS - 1);
256 let mut peak = 0.0f32;
257 let mut start = 0;
258 while start < n {
259 let len = (n - start).min(scratch.len());
260 let out = &mut scratch[..len];
261 interpolated(&buffer[start..start + len + TAPS - 1], bank, out, peak);
262 peak = peak.max(self::peak(out));
263 start += len;
264 }
265 peak
266}
267
268pub(crate) fn interpolated(buffer: &[f32], bank: &[[f32; TAPS]], out: &mut [f32], floor: f32) {
269 #[cfg(target_arch = "x86_64")]
270 if std::arch::is_x86_feature_detected!("avx2") && std::arch::is_x86_feature_detected!("fma") {
271 unsafe { interpolated_fma(buffer, bank, out, floor) };
272 return;
273 }
274 interpolated_plain(buffer, bank, out, floor)
275}
276
277#[cfg(target_arch = "x86_64")]
278#[target_feature(enable = "avx2,fma")]
279unsafe fn interpolated_fma(buffer: &[f32], bank: &[[f32; TAPS]], out: &mut [f32], floor: f32) {
280 interpolated_plain(buffer, bank, out, floor)
281}
282
283#[inline(always)]
284fn interpolated_plain(buffer: &[f32], bank: &[[f32; TAPS]], out: &mut [f32], floor: f32) {
285 const B: usize = 64;
286 let n = out.len();
287 assert!(buffer.len() >= n + TAPS - 1);
288 let gain = bank
289 .iter()
290 .map(|c| c.iter().map(|v| v.abs()).sum::<f32>())
291 .fold(0.0f32, f32::max);
292 let mut o = 0;
293 while o + B <= n {
294 if peak(&buffer[o..o + B + TAPS - 1]) * gain <= floor {
295 out[o..o + B].fill(0.0);
296 o += B;
297 continue;
298 }
299 let mut m = [0.0f32; B];
300 for c in bank {
301 let mut acc = [0.0f32; B];
302 for (j, &w) in c.iter().enumerate() {
303 let x: &[f32; B] = buffer[o + TAPS - 1 - j..o + TAPS - 1 - j + B]
304 .try_into()
305 .unwrap();
306 for k in 0..B {
307 acc[k] += w * x[k];
308 }
309 }
310 for k in 0..B {
311 m[k] = m[k].max(acc[k].abs());
312 }
313 }
314 out[o..o + B].copy_from_slice(&m);
315 o += B;
316 }
317 for (i, v) in out.iter_mut().enumerate().skip(o) {
318 *v = bank
319 .iter()
320 .map(|c| {
321 c.iter()
322 .enumerate()
323 .map(|(j, &w)| w * buffer[i + TAPS - 1 - j])
324 .sum::<f32>()
325 .abs()
326 })
327 .fold(0.0, f32::max);
328 }
329}
330
331#[derive(Debug, Clone, Copy)]
332struct Filters {
333 pre: Coeffs,
334 rlb: Coeffs,
335}
336
337#[derive(Debug, Clone)]
338pub struct LoudnessMeter {
339 rate: u32,
340 weights: Vec<f64>,
341 filters: Filters,
342 bank: Vec<[f32; TAPS]>,
343 channels: Vec<Channel>,
344 segments: Vec<usize>,
345 block: usize,
346 position: usize,
347 ring: [f64; SHORT],
348 ring_at: usize,
349 seen: u64,
350 momentary: f64,
351 short_term: f64,
352 momentary_max: f64,
353 short_term_max: f64,
354 integrated: Histogram,
355 range: Histogram,
356}
357
358impl LoudnessMeter {
359 pub fn new(rate: u32, channels: u16) -> Result<Self> {
360 if rate == 0 || channels == 0 {
361 return Err(Error::Invalid(format!(
362 "loudness meter at {rate} Hz with {channels} channels"
363 )));
364 }
365 let r = rate as f64;
366 Ok(LoudnessMeter {
367 rate,
368 weights: weights(channels as usize),
369 filters: Filters {
370 pre: prefilter(r),
371 rlb: rlb(r),
372 },
373 bank: polyphase(oversampling(r)),
374 channels: vec![Channel::new(); channels as usize],
375 segments: Vec::new(),
376 block: ((r * 0.1).round() as usize).max(1),
377 position: 0,
378 ring: [0.0; SHORT],
379 ring_at: 0,
380 seen: 0,
381 momentary: f64::NEG_INFINITY,
382 short_term: f64::NEG_INFINITY,
383 momentary_max: f64::NEG_INFINITY,
384 short_term_max: f64::NEG_INFINITY,
385 integrated: Histogram::new(),
386 range: Histogram::new(),
387 })
388 }
389
390 pub fn set_weights(&mut self, weights: &[f64]) {
391 for (i, w) in self.weights.iter_mut().enumerate() {
392 *w = weights.get(i).copied().unwrap_or(1.0);
393 }
394 }
395
396 pub fn process(&mut self, input: &AudioBuffer) -> Result<()> {
397 let nch = self.channels.len();
398 if input.rate != self.rate || input.channels as usize != nch {
399 return Err(Error::Invalid(format!(
400 "meter expects {} Hz {} channels, got {} Hz {} channels",
401 self.rate, nch, input.rate, input.channels
402 )));
403 }
404 let frames = input.frames();
405 if frames == 0 {
406 return Ok(());
407 }
408 self.segments.clear();
409 let mut left = frames;
410 let mut room = self.block - self.position;
411 while left > 0 {
412 let take = left.min(room);
413 self.segments.push(take);
414 left -= take;
415 room = self.block;
416 }
417 for (c, ch) in self.channels.iter_mut().enumerate() {
418 let start = ch.buffer.len();
419 ch.buffer.resize(start + frames, 0.0);
420 for (o, f) in ch.buffer[start..]
421 .iter_mut()
422 .zip(input.samples.chunks_exact(nch))
423 {
424 *o = f[c];
425 }
426 }
427 let (filters, segments, bank) = (self.filters, &self.segments, &self.bank);
428 if frames * nch >= PARALLEL && nch > 1 {
429 self.channels
430 .par_chunks_mut(2)
431 .for_each(|pair| run(pair, &filters, segments, bank));
432 } else {
433 self.channels
434 .chunks_mut(2)
435 .for_each(|pair| run(pair, &filters, segments, bank));
436 }
437 for k in 0..self.segments.len() {
438 let len = self.segments[k];
439 for ch in self.channels.iter_mut() {
440 ch.carry += ch.sums[k];
441 }
442 self.position += len;
443 if self.position == self.block {
444 let e = self
445 .channels
446 .iter_mut()
447 .zip(&self.weights)
448 .map(|(ch, w)| w * std::mem::take(&mut ch.carry))
449 .sum();
450 self.finish(e);
451 }
452 }
453 Ok(())
454 }
455
456 fn finish(&mut self, e: f64) {
457 self.ring[self.ring_at] = e;
458 self.ring_at = (self.ring_at + 1) % SHORT;
459 self.position = 0;
460 self.seen += 1;
461 let window = |k: usize| {
462 let s: f64 = (0..k)
463 .map(|i| self.ring[(self.ring_at + SHORT - 1 - i) % SHORT])
464 .sum();
465 s / (k * self.block) as f64
466 };
467 let em = window(MOMENTARY);
468 let es = window(SHORT);
469 self.momentary = lufs(em);
470 self.short_term = lufs(es);
471 if self.seen >= MOMENTARY as u64 {
472 self.integrated.add(em);
473 self.momentary_max = self.momentary_max.max(self.momentary);
474 }
475 if self.seen >= SHORT as u64 {
476 self.range.add(es);
477 self.short_term_max = self.short_term_max.max(self.short_term);
478 }
479 }
480
481 pub fn momentary(&self) -> f64 {
482 self.momentary
483 }
484
485 pub fn short_term(&self) -> f64 {
486 self.short_term
487 }
488
489 pub fn integrated(&self) -> f64 {
490 let h = &self.integrated;
491 h.gate(-10.0)
492 .and_then(|b| h.mean_from(b))
493 .map_or(f64::NEG_INFINITY, lufs)
494 }
495
496 pub fn range(&self) -> f64 {
497 let h = &self.range;
498 let Some(start) = h.gate(-20.0) else {
499 return 0.0;
500 };
501 let n: u64 = h.count[start..].iter().sum();
502 if n == 0 {
503 return 0.0;
504 }
505 let percentile = |p: f64| {
506 let rank = ((p * (n - 1) as f64).round() as u64).min(n - 1);
507 let mut c = 0u64;
508 for i in start..HIST_BINS {
509 c += h.count[i];
510 if c > rank {
511 return Histogram::centre(i);
512 }
513 }
514 Histogram::centre(HIST_BINS - 1)
515 };
516 (percentile(0.95) - percentile(0.10)).max(0.0)
517 }
518
519 pub fn true_peak(&self) -> f64 {
520 let p = self.channels.iter().fold(0.0f32, |m, c| m.max(c.true_peak));
521 20.0 * (p as f64).log10()
522 }
523
524 pub fn sample_peak(&self) -> f64 {
525 let p = self
526 .channels
527 .iter()
528 .fold(0.0f32, |m, c| m.max(c.sample_peak));
529 20.0 * (p as f64).log10()
530 }
531
532 pub fn summary(&self) -> Loudness {
533 Loudness {
534 integrated: self.integrated(),
535 range: self.range(),
536 true_peak: self.true_peak(),
537 sample_peak: self.sample_peak(),
538 momentary_max: self.momentary_max,
539 short_term_max: self.short_term_max,
540 }
541 }
542}