sva_samples/measure/
loudness_pattern.rs

1// Concern: specific loudness per auditory band per frame, masking spread in | Non-concern: BS.1770 loudness (loudness.rs), comparing two (sva-core) | IO: (planes, rate, start, hop, spl) -> pattern
2
3use std::f64::consts::TAU;
4
5use super::bands::{BAND_COUNT, FIRST_CAM, erb_hz};
6use crate::fft::fft;
7
8pub const MODEL: &str = "specific loudness after Moore, Glasberg and Baer (1997), J. Audio Eng. \
9                         Soc. 45(4) 224-240, per frame: a periodic Hann of 2048 samples at 48 kHz \
10                         (its length scaled to the rate) centred on each frame time, every \
11                         channel's power per bin averaged, a full-scale sine at `spl_db`; the \
12                         outer and middle ear as Terhardt's (1979) threshold in quiet, 3.64 \
13                         f^-0.8 - 6.5 e^(-0.6 (f - 3.3)^2) + 0.001 f^4 dB (f in kHz), taken \
14                         against its value at 1 kHz; each band's excitation the power through a \
15                         rounded-exponential auditory filter (1 + p g) e^(-p g), g the \
16                         frequency's distance from the centre over the centre, p = 4 fc / \
17                         ERB(fc) above it and p - 0.35 (p / p(1 kHz)) (X - 51) below it, X the \
18                         level in dB SPL within one ERB of the component, its upward spread of \
19                         masking; N' = C ((E + A)^0.2 - A^0.2), A twice the excitation at \
20                         threshold, times (2 E / (E + E_T))^1.5 under threshold E_T, C so a 1 \
21                         kHz sine at 40 dB SPL is 1 sone; the total the sum over bands 1 Cam \
22                         apart. A cell at or under `threshold_sone` is inaudible: a comparison \
23                         counts it zero where both sides are. Measured, not ISO 532-2: no \
24                         binaural inhibition and no temporal integration, and the threshold in \
25                         quiet stands for the low-frequency internal noise.";
26
27const WINDOW_SAMPLES_48K: f64 = 2048.0;
28pub const DEFAULT_SPL_DB: f64 = 100.0;
29const ALPHA: f64 = 0.2;
30const LOWER_SLOPE_DB: f64 = 51.0;
31const LOWER_SLOPE: f64 = 0.35;
32const LEAST_LOWER_SHARE: f64 = 0.1;
33const CALIBRATION_HZ: f64 = 1000.0;
34const CALIBRATION_DB: f64 = 40.0;
35
36#[derive(Clone, Debug, PartialEq)]
37pub struct PatternFrame {
38    pub t_secs: f64,
39    pub total_sone: f64,
40    /// Sone per Cam, band by band.
41    pub specific: Vec<f64>,
42}
43
44#[derive(Clone, Debug, PartialEq)]
45pub struct LoudnessPattern {
46    pub frame_secs: f64,
47    pub window_secs: f64,
48    pub spl_db: f64,
49    pub centres_hz: Vec<f64>,
50    pub threshold_sone: f64,
51    pub frames: Vec<PatternFrame>,
52}
53
54/// Terhardt's threshold in quiet, dB SPL.
55fn quiet_db(hz: f64) -> f64 {
56    let f = hz.max(20.0) / 1000.0;
57    3.64 * libm::pow(f, -0.8) - 6.5 * sva_formula::exp(-0.6 * (f - 3.3) * (f - 3.3))
58        + 1e-3 * f * f * f * f
59}
60
61/// `bands`' Cam scale inverted, its power read through libm so it is one value native and in
62/// wasm.
63fn centre_hz(cam: f64) -> f64 {
64    (libm::pow(10.0, cam / 21.4) - 1.0) * 1000.0 / 4.37
65}
66
67fn intensity(db: f64) -> f64 {
68    libm::pow(10.0, db / 10.0)
69}
70
71struct Ear {
72    centres: Vec<f64>,
73    threshold: f64,
74    a: f64,
75    c: f64,
76}
77
78impl Ear {
79    fn new(sr: f64) -> Ear {
80        let centres: Vec<f64> = (0..BAND_COUNT)
81            .map(|k| centre_hz(FIRST_CAM + k as f64))
82            .filter(|hz| *hz < sr / 2.0)
83            .collect();
84        let threshold = intensity(quiet_db(CALIBRATION_HZ));
85        let mut ear = Ear {
86            centres,
87            threshold,
88            a: 2.0 * threshold,
89            c: 1.0,
90        };
91        let tone = [(CALIBRATION_HZ, intensity(CALIBRATION_DB))];
92        let total: f64 = ear
93            .specific(&ear.excitation(&tone, &[CALIBRATION_DB]))
94            .iter()
95            .sum();
96        ear.c = 1.0 / total;
97        ear
98    }
99
100    /// Each band's excitation off components (hz, intensity at the eardrum), each with the
101    /// level per ERB about it.
102    fn excitation(&self, parts: &[(f64, f64)], per_erb_db: &[f64]) -> Vec<f64> {
103        let p_1k = 4.0 * 1000.0 / erb_hz(1000.0);
104        self.centres
105            .iter()
106            .map(|&fc| {
107                let p = 4.0 * fc / erb_hz(fc);
108                parts
109                    .iter()
110                    .zip(per_erb_db)
111                    .map(|(&(hz, i), &x)| {
112                        let g = (hz - fc).abs() / fc;
113                        let slope = match hz < fc {
114                            true => (p - LOWER_SLOPE * (p / p_1k) * (x - LOWER_SLOPE_DB))
115                                .max(LEAST_LOWER_SHARE * p),
116                            false => p,
117                        };
118                        let pg = slope * g;
119                        i * (1.0 + pg) * sva_formula::exp(-pg)
120                    })
121                    .sum()
122            })
123            .collect()
124    }
125
126    fn specific(&self, excitation: &[f64]) -> Vec<f64> {
127        excitation
128            .iter()
129            .map(|&e| {
130                let grown = libm::pow(e + self.a, ALPHA) - libm::pow(self.a, ALPHA);
131                let near = match e < self.threshold {
132                    true => libm::pow(2.0 * e / (e + self.threshold), 1.5),
133                    false => 1.0,
134                };
135                self.c * near * grown
136            })
137            .collect()
138    }
139
140    fn at_threshold(&self) -> f64 {
141        self.specific(&[self.threshold])[0]
142    }
143}
144
145/// Every frame centred on `start + k hop`, the target zero past its ends.
146pub fn measure(
147    planes: &[&[f64]],
148    sr: f64,
149    start: f64,
150    hop_secs: f64,
151    spl_db: f64,
152) -> LoudnessPattern {
153    let len = super::held(planes);
154    let ear = Ear::new(sr);
155    let n = super::scaled_pow2(WINDOW_SAMPLES_48K, sr);
156    let hop = ((hop_secs * sr).round() as usize).max(1);
157    let w: Vec<f64> = (0..n)
158        .map(|k| 0.5 - 0.5 * libm::cos(TAU * k as f64 / n as f64))
159        .collect();
160    let norm = 2.0 / (n as f64 * n as f64 * 0.375) / 0.5 * intensity(spl_db);
161    let bin_hz = sr / n as f64;
162    let weights: Vec<f64> = (0..=n / 2)
163        .map(|k| intensity(quiet_db(CALIBRATION_HZ) - quiet_db(k as f64 * bin_hz)))
164        .collect();
165    let frames = (0..=len.saturating_sub(1) / hop)
166        .map(|f| {
167            let centre = (f * hop) as i64;
168            let mut power = vec![0.0; n / 2 + 1];
169            for plane in planes {
170                let (mut re, mut im) = (vec![0.0; n], vec![0.0; n]);
171                for (i, wi) in w.iter().enumerate() {
172                    let at = centre + i as i64 - (n / 2) as i64;
173                    if at >= 0 && (at as usize) < len {
174                        re[i] = plane[at as usize] * wi;
175                    }
176                }
177                fft(&mut re, &mut im);
178                for k in 0..=n / 2 {
179                    power[k] += (re[k] * re[k] + im[k] * im[k]) / planes.len() as f64;
180                }
181            }
182            let parts: Vec<(f64, f64)> = (1..n / 2)
183                .map(|k| (k as f64 * bin_hz, power[k] * norm * weights[k]))
184                .collect();
185            let mut running = vec![0.0; parts.len() + 1];
186            for (i, p) in parts.iter().enumerate() {
187                running[i + 1] = running[i] + p.1;
188            }
189            let per_erb: Vec<f64> = parts
190                .iter()
191                .map(|&(hz, _)| {
192                    let half = erb_hz(hz) / 2.0;
193                    let first = ((hz - half) / bin_hz).ceil().max(1.0) as usize;
194                    let last = (((hz + half) / bin_hz).floor() as usize).min(parts.len());
195                    let held = running[last] - running[(first - 1).min(last)];
196                    10.0 * libm::log10(held + 1e-30)
197                })
198                .collect();
199            let specific = ear.specific(&ear.excitation(&parts, &per_erb));
200            PatternFrame {
201                t_secs: start + (f * hop) as f64 / sr,
202                total_sone: specific.iter().sum(),
203                specific,
204            }
205        })
206        .collect();
207    LoudnessPattern {
208        frame_secs: hop as f64 / sr,
209        window_secs: n as f64 / sr,
210        spl_db,
211        threshold_sone: ear.at_threshold(),
212        centres_hz: ear.centres,
213        frames,
214    }
215}
sva_samples/measure/loudness_pattern.rs

sva_samples/measure/
loudness_pattern.rs
