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

sva_samples/measure/
alias.rs
