sva_samples/measure/
alias.rs

1// Concern: scores a render against the same one oversampled, as ASR and per-band NMR | Non-concern: producing either render, the plain spectrum (spectrum.rs) | IO: (base, high, k, sr, start) -> Alias
2
3use crate::fft::fft;
4use crate::measure::spectrum::db;
5use crate::stft::hann_periodic as hann;
6
7/// Zwicker's critical bands; the last one runs to Nyquist.
8const BARK_EDGES: [f64; 25] = [
9    0.0, 100.0, 200.0, 300.0, 400.0, 510.0, 630.0, 770.0, 920.0, 1080.0, 1270.0, 1480.0, 1720.0,
10    2000.0, 2320.0, 2700.0, 3150.0, 3700.0, 4400.0, 5300.0, 6400.0, 7700.0, 9500.0, 12000.0,
11    15500.0,
12];
13
14/// The spreading triangle is steep down the Bark scale and shallow up it: a masker reaches
15/// more than twice as far above itself as below.
16const SPREAD_DOWN_DB_PER_BARK: f64 = 27.0;
17const SPREAD_UP_DB_PER_BARK: f64 = 12.0;
18
19/// What a full-scale sine reaches the ear at, so quiet has an energy.
20pub const PLAYBACK_DB_SPL: f64 = 90.0;
21
22/// 23 ms at 44.1 kHz, the length NMR is published at.
23pub const ALIAS_FRAME: usize = 1024;
24
25/// Under it a frame is fade or tail, its NMR two near-silences' ratio.
26const GATE_DB: f64 = -70.0;
27
28/// Above it the alias is audible (Brandenburg).
29pub const AUDIBLE_NMR_DB: f64 = -10.0;
30
31#[derive(Clone, Debug, PartialEq)]
32pub struct AliasBand {
33    pub lo_hz: f64,
34    pub hi_hz: f64,
35    pub signal_db: f64,
36    pub alias_db: f64,
37    pub nmr_db: f64,
38}
39
40#[derive(Clone, Debug, PartialEq)]
41pub struct Alias {
42    pub oversample: usize,
43    pub sample_rate: f64,
44    pub frame_size: usize,
45    pub frames: usize,
46    pub scored_frames: usize,
47    pub playback_db_spl: f64,
48    pub asr_db: f64,
49    pub nmr_db: f64,
50    pub nmr_peak_db: f64,
51    pub peak_at_secs: f64,
52    pub audible: bool,
53    /// Filled by whoever holds the graph: how many `instances` are not a function of `t` alone.
54    pub rate_dependent: usize,
55    pub instances: usize,
56    pub bands: Vec<AliasBand>,
57}
58
59pub fn worst(scored: impl IntoIterator<Item = Alias>) -> Option<Alias> {
60    scored
61        .into_iter()
62        .reduce(|held, next| match next.nmr_peak_db > held.nmr_peak_db {
63            true => next,
64            false => held,
65        })
66}
67
68/// The alias is what the base band holds and the `oversample`x render does not. One bin grid
69/// for both, so no decimation filter blurs the top octave this is measuring.
70pub fn measure_alias(
71    base: &[f64],
72    high: &[f64],
73    oversample: usize,
74    sample_rate: f64,
75    start_secs: f64,
76) -> Alias {
77    assert!(
78        oversample.is_power_of_two(),
79        "oversample must be a power of two so both transforms are radix-2, got {oversample}"
80    );
81    let frame = ALIAS_FRAME;
82    let wide = frame * oversample;
83    let hop = frame / 2;
84    let bins = frame / 2 + 1;
85    let bin_hz = sample_rate / frame as f64;
86    let edges = band_edges(sample_rate / 2.0);
87
88    let w_base = hann(frame);
89    let w_high = hann(wide);
90    let mut sig_bands = vec![0f64; edges.len() - 1];
91    let mut err_bands = vec![0f64; edges.len() - 1];
92    let mut sig_total = 0f64;
93    let mut err_total = 0f64;
94    let mut nmr_sum = 0f64;
95    let mut nmr_peak = f64::NEG_INFINITY;
96    let mut peak_at = start_secs;
97    let mut frames = 0usize;
98    let mut scored = 0usize;
99
100    let mut start = 0usize;
101    while start + frame <= base.len() && (start + frame) * oversample <= high.len() {
102        let (re_b, im_b) = transform(base, start, &w_base, 1);
103        let (re_h, im_h) = transform(high, start * oversample, &w_high, oversample);
104        frames += 1;
105
106        let mut sig = vec![0f64; bins];
107        let mut err = vec![0f64; bins];
108        for k in 0..bins {
109            sig[k] = re_h[k] * re_h[k] + im_h[k] * im_h[k];
110            let (dr, di) = (re_b[k] - re_h[k], im_b[k] - im_h[k]);
111            err[k] = dr * dr + di * di;
112        }
113        sig_total += sig.iter().sum::<f64>();
114        err_total += err.iter().sum::<f64>();
115
116        let s = group(&sig, bin_hz, &edges);
117        let e = group(&err, bin_hz, &edges);
118        for (acc, v) in sig_bands.iter_mut().zip(&s) {
119            *acc += v;
120        }
121        for (acc, v) in err_bands.iter_mut().zip(&e) {
122            *acc += v;
123        }
124
125        if db(s.iter().sum::<f64>().sqrt()) < GATE_DB {
126            start += hop;
127            continue;
128        }
129        let ratio = nmr_of(&s, &e, &edges);
130        nmr_sum += ratio;
131        scored += 1;
132        let frame_db = 10.0 * ratio.max(f64::MIN_POSITIVE).log10();
133        if frame_db > nmr_peak {
134            nmr_peak = frame_db;
135            peak_at = start_secs + start as f64 / sample_rate;
136        }
137        start += hop;
138    }
139
140    let mean = if scored > 0 {
141        nmr_sum / scored as f64
142    } else {
143        0.0
144    };
145    let nmr_db = 10.0 * mean.max(f64::MIN_POSITIVE).log10();
146    Alias {
147        oversample,
148        sample_rate,
149        frame_size: frame,
150        frames,
151        scored_frames: scored,
152        playback_db_spl: PLAYBACK_DB_SPL,
153        asr_db: 10.0
154            * (err_total / sig_total.max(f64::MIN_POSITIVE))
155                .max(f64::MIN_POSITIVE)
156                .log10(),
157        nmr_db,
158        nmr_peak_db: if scored > 0 { nmr_peak } else { nmr_db },
159        peak_at_secs: peak_at,
160        audible: scored > 0 && nmr_db > AUDIBLE_NMR_DB,
161        rate_dependent: 0,
162        instances: 0,
163        bands: bands(&sig_bands, &err_bands, &edges, frames.max(1)),
164    }
165}
166
167fn transform(x: &[f64], start: usize, window: &[f64], stride: usize) -> (Vec<f64>, Vec<f64>) {
168    let n = window.len();
169    let mut re = vec![0f64; n];
170    let mut im = vec![0f64; n];
171    for (i, w) in window.iter().enumerate() {
172        re[i] = x.get(start + i).copied().unwrap_or(0.0) * w;
173    }
174    fft(&mut re, &mut im);
175    let scale = 4.0 / n as f64;
176    let keep = n / (2 * stride) + 1;
177    re.truncate(keep);
178    im.truncate(keep);
179    for (r, i) in re.iter_mut().zip(im.iter_mut()) {
180        *r *= scale;
181        *i *= scale;
182    }
183    (re, im)
184}
185
186fn band_edges(nyquist: f64) -> Vec<f64> {
187    let mut edges: Vec<f64> = BARK_EDGES
188        .iter()
189        .copied()
190        .filter(|e| *e < nyquist)
191        .collect();
192    edges.push(nyquist);
193    edges
194}
195
196fn group(power: &[f64], bin_hz: f64, edges: &[f64]) -> Vec<f64> {
197    let mut out = vec![0f64; edges.len() - 1];
198    for (k, p) in power.iter().enumerate() {
199        let hz = k as f64 * bin_hz;
200        let b = edges
201            .windows(2)
202            .position(|w| hz >= w[0] && hz < w[1])
203            .unwrap_or(out.len() - 1);
204        out[b] += p;
205    }
206    out
207}
208
209/// Every NMR reads this one, so no band disagrees about what masks it. The offset it drops
210/// the spread by is our own, standing in for the tonality estimate this makes none of.
211fn thresholds(signal: &[f64], edges: &[f64]) -> Vec<f64> {
212    (0..signal.len())
213        .map(|j| {
214            let zj = j as f64 + 0.5;
215            let spread: f64 = signal
216                .iter()
217                .enumerate()
218                .map(|(i, s)| {
219                    let zi = i as f64 + 0.5;
220                    let slope = if zj >= zi {
221                        SPREAD_UP_DB_PER_BARK
222                    } else {
223                        SPREAD_DOWN_DB_PER_BARK
224                    };
225                    s * 10f64.powf(-slope * (zj - zi).abs() / 10.0)
226                })
227                .sum();
228            let offset = if zj <= 12.0 { 3.0 } else { 0.25 * zj };
229            let centre = (edges[j] + edges[j + 1]) / 2.0;
230            (spread * 10f64.powf(-offset / 10.0)).max(quiet_energy(centre))
231        })
232        .collect()
233}
234
235fn nmr_of(signal: &[f64], alias: &[f64], edges: &[f64]) -> f64 {
236    let masked = thresholds(signal, edges);
237    alias.iter().zip(&masked).map(|(a, m)| a / m).sum::<f64>() / signal.len() as f64
238}
239
240/// Terhardt's threshold in quiet, in a band energy's own units.
241fn quiet_energy(hz: f64) -> f64 {
242    let f = (hz / 1000.0).max(0.02);
243    let db_spl =
244        3.64 * f.powf(-0.8) - 6.5 * (-0.6 * (f - 3.3) * (f - 3.3)).exp() + 0.001 * f.powi(4);
245    10f64.powf((db_spl - PLAYBACK_DB_SPL) / 10.0)
246}
247
248fn bands(signal: &[f64], alias: &[f64], edges: &[f64], frames: usize) -> Vec<AliasBand> {
249    let mean: Vec<f64> = signal.iter().map(|s| s / frames as f64).collect();
250    let masked = thresholds(&mean, edges);
251    (0..mean.len())
252        .map(|j| {
253            let a = alias[j] / frames as f64;
254            AliasBand {
255                lo_hz: edges[j],
256                hi_hz: edges[j + 1],
257                signal_db: db(mean[j].sqrt()),
258                alias_db: db(a.sqrt()),
259                nmr_db: 10.0 * (a / masked[j]).max(f64::MIN_POSITIVE).log10(),
260            }
261        })
262        .collect()
263}
sva_samples/measure/alias.rs

sva_samples/measure/
alias.rs
