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 * libm::exp(-0.6 * (f - 3.3) * (f - 3.3)) + 1e-3 * f * f * f * f
58}
59
60/// `bands`' Cam scale inverted, its power read through libm so it is one value native and in
61/// wasm.
62fn centre_hz(cam: f64) -> f64 {
63    (libm::pow(10.0, cam / 21.4) - 1.0) * 1000.0 / 4.37
64}
65
66fn intensity(db: f64) -> f64 {
67    libm::pow(10.0, db / 10.0)
68}
69
70struct Ear {
71    centres: Vec<f64>,
72    threshold: f64,
73    a: f64,
74    c: f64,
75}
76
77impl Ear {
78    fn new(sr: f64) -> Ear {
79        let centres: Vec<f64> = (0..BAND_COUNT)
80            .map(|k| centre_hz(FIRST_CAM + k as f64))
81            .filter(|hz| *hz < sr / 2.0)
82            .collect();
83        let threshold = intensity(quiet_db(CALIBRATION_HZ));
84        let mut ear = Ear {
85            centres,
86            threshold,
87            a: 2.0 * threshold,
88            c: 1.0,
89        };
90        let tone = [(CALIBRATION_HZ, intensity(CALIBRATION_DB))];
91        let total: f64 = ear
92            .specific(&ear.excitation(&tone, &[CALIBRATION_DB]))
93            .iter()
94            .sum();
95        ear.c = 1.0 / total;
96        ear
97    }
98
99    /// Each band's excitation off components (hz, intensity at the eardrum), each with the
100    /// level per ERB about it.
101    fn excitation(&self, parts: &[(f64, f64)], per_erb_db: &[f64]) -> Vec<f64> {
102        let p_1k = 4.0 * 1000.0 / erb_hz(1000.0);
103        self.centres
104            .iter()
105            .map(|&fc| {
106                let p = 4.0 * fc / erb_hz(fc);
107                parts
108                    .iter()
109                    .zip(per_erb_db)
110                    .map(|(&(hz, i), &x)| {
111                        let g = (hz - fc).abs() / fc;
112                        let slope = match hz < fc {
113                            true => (p - LOWER_SLOPE * (p / p_1k) * (x - LOWER_SLOPE_DB))
114                                .max(LEAST_LOWER_SHARE * p),
115                            false => p,
116                        };
117                        let pg = slope * g;
118                        i * (1.0 + pg) * libm::exp(-pg)
119                    })
120                    .sum()
121            })
122            .collect()
123    }
124
125    fn specific(&self, excitation: &[f64]) -> Vec<f64> {
126        excitation
127            .iter()
128            .map(|&e| {
129                let grown = libm::pow(e + self.a, ALPHA) - libm::pow(self.a, ALPHA);
130                let near = match e < self.threshold {
131                    true => libm::pow(2.0 * e / (e + self.threshold), 1.5),
132                    false => 1.0,
133                };
134                self.c * near * grown
135            })
136            .collect()
137    }
138
139    fn at_threshold(&self) -> f64 {
140        self.specific(&[self.threshold])[0]
141    }
142}
143
144/// Every frame centred on `start + k hop`, the target zero past its ends.
145pub fn measure(
146    planes: &[&[f64]],
147    sr: f64,
148    start: f64,
149    hop_secs: f64,
150    spl_db: f64,
151) -> LoudnessPattern {
152    let len = super::held(planes);
153    let ear = Ear::new(sr);
154    let n = super::scaled_pow2(WINDOW_SAMPLES_48K, sr);
155    let hop = ((hop_secs * sr).round() as usize).max(1);
156    let w: Vec<f64> = (0..n)
157        .map(|k| 0.5 - 0.5 * libm::cos(TAU * k as f64 / n as f64))
158        .collect();
159    let norm = 2.0 / (n as f64 * n as f64 * 0.375) / 0.5 * intensity(spl_db);
160    let bin_hz = sr / n as f64;
161    let weights: Vec<f64> = (0..=n / 2)
162        .map(|k| intensity(quiet_db(CALIBRATION_HZ) - quiet_db(k as f64 * bin_hz)))
163        .collect();
164    let frames = (0..=len.saturating_sub(1) / hop)
165        .map(|f| {
166            let centre = (f * hop) as i64;
167            let mut power = vec![0.0; n / 2 + 1];
168            for plane in planes {
169                let (mut re, mut im) = (vec![0.0; n], vec![0.0; n]);
170                for (i, wi) in w.iter().enumerate() {
171                    let at = centre + i as i64 - (n / 2) as i64;
172                    if at >= 0 && (at as usize) < len {
173                        re[i] = plane[at as usize] * wi;
174                    }
175                }
176                fft(&mut re, &mut im);
177                for k in 0..=n / 2 {
178                    power[k] += (re[k] * re[k] + im[k] * im[k]) / planes.len() as f64;
179                }
180            }
181            let parts: Vec<(f64, f64)> = (1..n / 2)
182                .map(|k| (k as f64 * bin_hz, power[k] * norm * weights[k]))
183                .collect();
184            let mut running = vec![0.0; parts.len() + 1];
185            for (i, p) in parts.iter().enumerate() {
186                running[i + 1] = running[i] + p.1;
187            }
188            let per_erb: Vec<f64> = parts
189                .iter()
190                .map(|&(hz, _)| {
191                    let half = erb_hz(hz) / 2.0;
192                    let first = ((hz - half) / bin_hz).ceil().max(1.0) as usize;
193                    let last = (((hz + half) / bin_hz).floor() as usize).min(parts.len());
194                    let held = running[last] - running[(first - 1).min(last)];
195                    10.0 * libm::log10(held + 1e-30)
196                })
197                .collect();
198            let specific = ear.specific(&ear.excitation(&parts, &per_erb));
199            PatternFrame {
200                t_secs: start + (f * hop) as f64 / sr,
201                total_sone: specific.iter().sum(),
202                specific,
203            }
204        })
205        .collect();
206    LoudnessPattern {
207        frame_secs: hop as f64 / sr,
208        window_secs: n as f64 / sr,
209        spl_db,
210        threshold_sone: ear.at_threshold(),
211        centres_hz: ear.centres,
212        frames,
213    }
214}
215
216/// DIN 45692's weighted first moment over Bark (Traunmüller's), in acum; `None` where silent.
217pub fn sharpness(centres_hz: &[f64], specific: &[f64]) -> Option<f64> {
218    let total: f64 = specific.iter().sum();
219    if total <= 0.0 {
220        return None;
221    }
222    let moment: f64 = centres_hz
223        .iter()
224        .zip(specific)
225        .map(|(&hz, &n)| {
226            let z = 26.81 * hz / (1960.0 + hz) - 0.53;
227            let g = match z > 15.8 {
228                true => 0.15 * libm::exp(0.42 * (z - 15.8)) + 0.85,
229                false => 1.0,
230            };
231            n * g * z
232        })
233        .sum();
234    Some(0.11 * moment / total)
235}
sva_samples/measure/loudness_pattern.rs

sva_samples/measure/
loudness_pattern.rs
