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    pub instances: usize,
54    pub bands: Vec<AliasBand>,
55}
56
57pub fn worst(scored: impl IntoIterator<Item = Alias>) -> Option<Alias> {
58    scored
59        .into_iter()
60        .reduce(|held, next| match next.nmr_peak_db > held.nmr_peak_db {
61            true => next,
62            false => held,
63        })
64}
65
66/// The alias is what the base band holds and the `oversample`x render does not. One bin grid
67/// for both, so no decimation filter blurs the top octave this is measuring.
68pub fn measure_alias(
69    base: &[f64],
70    high: &[f64],
71    oversample: usize,
72    sample_rate: f64,
73    start_secs: f64,
74) -> Alias {
75    assert!(
76        oversample.is_power_of_two(),
77        "oversample must be a power of two so both transforms are radix-2, got {oversample}"
78    );
79    let frame = ALIAS_FRAME;
80    let wide = frame * oversample;
81    let hop = frame / 2;
82    let bins = frame / 2 + 1;
83    let bin_hz = sample_rate / frame as f64;
84    let edges = band_edges(sample_rate / 2.0);
85
86    let w_base = hann(frame);
87    let w_high = hann(wide);
88    let mut sig_bands = vec![0f64; edges.len() - 1];
89    let mut err_bands = vec![0f64; edges.len() - 1];
90    let mut sig_total = 0f64;
91    let mut err_total = 0f64;
92    let mut nmr_sum = 0f64;
93    let mut nmr_peak = f64::NEG_INFINITY;
94    let mut peak_at = start_secs;
95    let mut frames = 0usize;
96    let mut scored = 0usize;
97
98    let mut start = 0usize;
99    while start + frame <= base.len() && (start + frame) * oversample <= high.len() {
100        let (re_b, im_b) = transform(base, start, &w_base, 1);
101        let (re_h, im_h) = transform(high, start * oversample, &w_high, oversample);
102        frames += 1;
103
104        let mut sig = vec![0f64; bins];
105        let mut err = vec![0f64; bins];
106        for k in 0..bins {
107            sig[k] = re_h[k] * re_h[k] + im_h[k] * im_h[k];
108            let (dr, di) = (re_b[k] - re_h[k], im_b[k] - im_h[k]);
109            err[k] = dr * dr + di * di;
110        }
111        sig_total += sig.iter().sum::<f64>();
112        err_total += err.iter().sum::<f64>();
113
114        let s = group(&sig, bin_hz, &edges);
115        let e = group(&err, bin_hz, &edges);
116        for (acc, v) in sig_bands.iter_mut().zip(&s) {
117            *acc += v;
118        }
119        for (acc, v) in err_bands.iter_mut().zip(&e) {
120            *acc += v;
121        }
122
123        if db(s.iter().sum::<f64>().sqrt()) < GATE_DB {
124            start += hop;
125            continue;
126        }
127        let ratio = nmr_of(&s, &e, &edges);
128        nmr_sum += ratio;
129        scored += 1;
130        let frame_db = 10.0 * ratio.max(f64::MIN_POSITIVE).log10();
131        if frame_db > nmr_peak {
132            nmr_peak = frame_db;
133            peak_at = start_secs + start as f64 / sample_rate;
134        }
135        start += hop;
136    }
137
138    let mean = if scored > 0 {
139        nmr_sum / scored as f64
140    } else {
141        0.0
142    };
143    let nmr_db = 10.0 * mean.max(f64::MIN_POSITIVE).log10();
144    Alias {
145        oversample,
146        sample_rate,
147        frame_size: frame,
148        frames,
149        scored_frames: scored,
150        playback_db_spl: PLAYBACK_DB_SPL,
151        asr_db: 10.0
152            * (err_total / sig_total.max(f64::MIN_POSITIVE))
153                .max(f64::MIN_POSITIVE)
154                .log10(),
155        nmr_db,
156        nmr_peak_db: if scored > 0 { nmr_peak } else { nmr_db },
157        peak_at_secs: peak_at,
158        audible: scored > 0 && nmr_db > AUDIBLE_NMR_DB,
159        instances: 0,
160        bands: bands(&sig_bands, &err_bands, &edges, frames.max(1)),
161    }
162}
163
164fn transform(x: &[f64], start: usize, window: &[f64], stride: usize) -> (Vec<f64>, Vec<f64>) {
165    let n = window.len();
166    let mut re = vec![0f64; n];
167    let mut im = vec![0f64; n];
168    for (i, w) in window.iter().enumerate() {
169        re[i] = x.get(start + i).copied().unwrap_or(0.0) * w;
170    }
171    fft(&mut re, &mut im);
172    let scale = 4.0 / n as f64;
173    let keep = n / (2 * stride) + 1;
174    re.truncate(keep);
175    im.truncate(keep);
176    for (r, i) in re.iter_mut().zip(im.iter_mut()) {
177        *r *= scale;
178        *i *= scale;
179    }
180    (re, im)
181}
182
183fn band_edges(nyquist: f64) -> Vec<f64> {
184    let mut edges: Vec<f64> = BARK_EDGES
185        .iter()
186        .copied()
187        .filter(|e| *e < nyquist)
188        .collect();
189    edges.push(nyquist);
190    edges
191}
192
193fn group(power: &[f64], bin_hz: f64, edges: &[f64]) -> Vec<f64> {
194    let mut out = vec![0f64; edges.len() - 1];
195    for (k, p) in power.iter().enumerate() {
196        let hz = k as f64 * bin_hz;
197        let b = edges
198            .windows(2)
199            .position(|w| hz >= w[0] && hz < w[1])
200            .unwrap_or(out.len() - 1);
201        out[b] += p;
202    }
203    out
204}
205
206/// Every NMR reads this one, so no band disagrees about what masks it. The offset it drops
207/// the spread by is our own, standing in for the tonality estimate this makes none of.
208fn thresholds(signal: &[f64], edges: &[f64]) -> Vec<f64> {
209    (0..signal.len())
210        .map(|j| {
211            let zj = j as f64 + 0.5;
212            let spread: f64 = signal
213                .iter()
214                .enumerate()
215                .map(|(i, s)| {
216                    let zi = i as f64 + 0.5;
217                    let slope = if zj >= zi {
218                        SPREAD_UP_DB_PER_BARK
219                    } else {
220                        SPREAD_DOWN_DB_PER_BARK
221                    };
222                    s * 10f64.powf(-slope * (zj - zi).abs() / 10.0)
223                })
224                .sum();
225            let offset = if zj <= 12.0 { 3.0 } else { 0.25 * zj };
226            let centre = (edges[j] + edges[j + 1]) / 2.0;
227            (spread * 10f64.powf(-offset / 10.0)).max(quiet_energy(centre))
228        })
229        .collect()
230}
231
232fn nmr_of(signal: &[f64], alias: &[f64], edges: &[f64]) -> f64 {
233    let masked = thresholds(signal, edges);
234    alias.iter().zip(&masked).map(|(a, m)| a / m).sum::<f64>() / signal.len() as f64
235}
236
237/// Terhardt's threshold in quiet, in a band energy's own units.
238fn quiet_energy(hz: f64) -> f64 {
239    let f = (hz / 1000.0).max(0.02);
240    let db_spl =
241        3.64 * f.powf(-0.8) - 6.5 * (-0.6 * (f - 3.3) * (f - 3.3)).exp() + 0.001 * f.powi(4);
242    10f64.powf((db_spl - PLAYBACK_DB_SPL) / 10.0)
243}
244
245fn bands(signal: &[f64], alias: &[f64], edges: &[f64], frames: usize) -> Vec<AliasBand> {
246    let mean: Vec<f64> = signal.iter().map(|s| s / frames as f64).collect();
247    let masked = thresholds(&mean, edges);
248    (0..mean.len())
249        .map(|j| {
250            let a = alias[j] / frames as f64;
251            AliasBand {
252                lo_hz: edges[j],
253                hi_hz: edges[j + 1],
254                signal_db: db(mean[j].sqrt()),
255                alias_db: db(a.sqrt()),
256                nmr_db: 10.0 * (a / masked[j]).max(f64::MIN_POSITIVE).log10(),
257            }
258        })
259        .collect()
260}
sva_samples/measure/alias.rs

sva_samples/measure/
alias.rs
