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/// The multiple of the rate an alias score's reference is read at.
29pub const ALIAS_OVERSAMPLE: usize = 4;
30
31/// Above it the alias is audible (Brandenburg).
32pub const AUDIBLE_NMR_DB: f64 = -10.0;
33
34#[derive(Clone, Debug, PartialEq)]
35pub struct AliasBand {
36    pub lo_hz: f64,
37    pub hi_hz: f64,
38    pub signal_db: f64,
39    pub alias_db: f64,
40    pub nmr_db: f64,
41}
42
43#[derive(Clone, Debug, PartialEq)]
44pub struct Alias {
45    pub oversample: usize,
46    pub sample_rate: f64,
47    pub frame_size: usize,
48    pub frames: usize,
49    pub scored_frames: usize,
50    pub playback_db_spl: f64,
51    pub asr_db: f64,
52    pub nmr_db: f64,
53    pub nmr_peak_db: f64,
54    pub peak_at_secs: f64,
55    pub audible: bool,
56    pub instances: usize,
57    pub bands: Vec<AliasBand>,
58}
59
60pub fn worst(scored: impl IntoIterator<Item = Alias>) -> Option<Alias> {
61    scored
62        .into_iter()
63        .reduce(|held, next| match next.nmr_peak_db > held.nmr_peak_db {
64            true => next,
65            false => held,
66        })
67}
68
69/// The alias is what the base band holds and the `oversample`x render does not. One bin grid
70/// for both, so no decimation filter blurs the top octave this is measuring.
71pub fn measure_alias(
72    base: &[f64],
73    high: &[f64],
74    oversample: usize,
75    sample_rate: f64,
76    start_secs: f64,
77) -> Alias {
78    assert!(
79        oversample.is_power_of_two(),
80        "oversample must be a power of two so both transforms are radix-2, got {oversample}"
81    );
82    let frame = ALIAS_FRAME;
83    let wide = frame * oversample;
84    let hop = frame / 2;
85    let bins = frame / 2 + 1;
86    let bin_hz = sample_rate / frame as f64;
87    let edges = band_edges(sample_rate / 2.0);
88
89    let w_base = hann(frame);
90    let w_high = hann(wide);
91    let mut sig_bands = vec![0f64; edges.len() - 1];
92    let mut err_bands = vec![0f64; edges.len() - 1];
93    let mut sig_total = 0f64;
94    let mut err_total = 0f64;
95    let mut nmr_sum = 0f64;
96    let mut nmr_peak = f64::NEG_INFINITY;
97    let mut peak_at = start_secs;
98    let mut frames = 0usize;
99    let mut scored = 0usize;
100
101    let mut start = 0usize;
102    while start + frame <= base.len() && (start + frame) * oversample <= high.len() {
103        let (re_b, im_b) = transform(base, start, &w_base, 1);
104        let (re_h, im_h) = transform(high, start * oversample, &w_high, oversample);
105        frames += 1;
106
107        let mut sig = vec![0f64; bins];
108        let mut err = vec![0f64; bins];
109        for k in 0..bins {
110            sig[k] = re_h[k] * re_h[k] + im_h[k] * im_h[k];
111            let (dr, di) = (re_b[k] - re_h[k], im_b[k] - im_h[k]);
112            err[k] = dr * dr + di * di;
113        }
114        sig_total += sig.iter().sum::<f64>();
115        err_total += err.iter().sum::<f64>();
116
117        let s = group(&sig, bin_hz, &edges);
118        let e = group(&err, bin_hz, &edges);
119        for (acc, v) in sig_bands.iter_mut().zip(&s) {
120            *acc += v;
121        }
122        for (acc, v) in err_bands.iter_mut().zip(&e) {
123            *acc += v;
124        }
125
126        if db(s.iter().sum::<f64>().sqrt()) < GATE_DB {
127            start += hop;
128            continue;
129        }
130        let ratio = nmr_of(&s, &e, &edges);
131        nmr_sum += ratio;
132        scored += 1;
133        let frame_db = 10.0 * ratio.max(f64::MIN_POSITIVE).log10();
134        if frame_db > nmr_peak {
135            nmr_peak = frame_db;
136            peak_at = start_secs + start as f64 / sample_rate;
137        }
138        start += hop;
139    }
140
141    let mean = if scored > 0 {
142        nmr_sum / scored as f64
143    } else {
144        0.0
145    };
146    let nmr_db = 10.0 * mean.max(f64::MIN_POSITIVE).log10();
147    Alias {
148        oversample,
149        sample_rate,
150        frame_size: frame,
151        frames,
152        scored_frames: scored,
153        playback_db_spl: PLAYBACK_DB_SPL,
154        asr_db: 10.0
155            * (err_total / sig_total.max(f64::MIN_POSITIVE))
156                .max(f64::MIN_POSITIVE)
157                .log10(),
158        nmr_db,
159        nmr_peak_db: if scored > 0 { nmr_peak } else { nmr_db },
160        peak_at_secs: peak_at,
161        audible: scored > 0 && nmr_db > AUDIBLE_NMR_DB,
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
