sound_to_patch/audio/
stereo.rs

1//! The stereo screen: correlation, side over mid and the cross-spectrum
2//! phase, for the full band and three bands.
3//!
4//! The screen cannot see an anti-phase stereo chorus or a light reverb, and
5//! a mono stem can still get its width from a reverb tail. The confidence of
6//! the screen says so.
7
8use crate::cast;
9use crate::confidence::{Confidence, ConfidenceError, Level, NoSpread, Reason, ReasonCode, Signal};
10use crate::dsp::{ACTIVE_GATE_DB, ActiveFrames, Complex64, Stft, StftParams, frame_rms};
11use crate::units::{Correlation, Decibels, Degrees, Hz};
12
13use super::SampleRate;
14
15/// The upper edge of the sub band.
16pub const SUB_EDGE_HZ: Hz = Hz::from_const(120.0);
17/// The edge between the low-mid and the high band. The descriptors use it
18/// as a band edge too.
19pub const HIGH_EDGE_HZ: Hz = Hz::from_const(1000.0);
20/// The STFT frame length of the screen, for every f0.
21pub const SCREEN_FRAME_LEN: usize = 8192;
22/// The harmonics whose cross-spectrum phase the screen averages.
23pub const SCREEN_HARMONICS: usize = 20;
24/// The convention of the screen.
25pub const STEREO_CONVENTION: &str = "rho: Pearson correlation at zero lag; sm_db: \
2620 log10(rms(S) / rms(M)) with M = (L + R) / 2 and S = (L - R) / 2; per band from the \
278192-point Hann STFT over bins inside the band, DC excluded; ipd_deg: abs(angle(L conj R)) \
28at the largest cross-spectrum bin within 2 bins of h f0, h = 1..20 inside the band, \
29weighted by |L conj R|, then the mean over active frames, whose RMS over both channels, \
30sqrt((rms_L^2 + rms_R^2) / 2), lies within 40 dB of the loudest; bands: sub below 120 Hz, \
31low_mid 120 Hz to 1 kHz, high above 1 kHz";
32
33/// The bands of the screen.
34#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
35pub enum Band {
36    /// Every frequency.
37    Full,
38    /// Below 120 Hz.
39    Sub,
40    /// From 120 Hz to 1 kHz.
41    LowMid,
42    /// From 1 kHz up.
43    High,
44}
45
46impl Band {
47    /// Every band, full band first.
48    pub const ALL: [Self; 4] = [Self::Full, Self::Sub, Self::LowMid, Self::High];
49
50    /// Returns whether a frequency lies in the band: lower edge inclusive,
51    /// upper edge exclusive. 0 Hz lies in none.
52    #[must_use]
53    pub fn contains(self, frequency: f64) -> bool {
54        frequency > 0.0
55            && match self {
56                Self::Full => true,
57                Self::Sub => frequency < SUB_EDGE_HZ.get(),
58                Self::LowMid => (SUB_EDGE_HZ.get()..HIGH_EDGE_HZ.get()).contains(&frequency),
59                Self::High => frequency >= HIGH_EDGE_HZ.get(),
60            }
61    }
62}
63
64/// The statistics of one band. A value is `None` when it does not exist:
65/// a silent channel, a zero mid or side, or no f0.
66#[derive(Clone, Copy, Debug, Default, PartialEq)]
67pub struct StereoStats {
68    /// The correlation of the two channels.
69    pub rho: Option<Correlation>,
70    /// Side over mid in decibels.
71    pub sm_db: Option<Decibels>,
72    /// The mean phase difference at the harmonics of the f0.
73    pub ipd_deg: Option<Degrees>,
74}
75
76/// The stereo screen of a stem or of a sample range.
77#[derive(Clone, Debug, PartialEq)]
78pub struct StereoScreen {
79    bands: [StereoStats; 4],
80    confidence: Confidence<NoSpread>,
81}
82
83impl StereoScreen {
84    /// The screen of a one-channel file: every value is `None`.
85    ///
86    /// # Errors
87    ///
88    /// Never in practice; the reason is built from a fixed count.
89    pub fn one_channel() -> Result<Self, ConfidenceError> {
90        Ok(Self {
91            bands: [StereoStats::default(); 4],
92            confidence: Confidence::no_estimate(vec![Reason::new(
93                ReasonCode::OneChannel,
94                Signal::Count(1),
95            )?])?,
96        })
97    }
98
99    /// Returns the statistics of a band.
100    #[must_use]
101    pub fn band(&self, band: Band) -> StereoStats {
102        self.bands[match band {
103            Band::Full => 0,
104            Band::Sub => 1,
105            Band::LowMid => 2,
106            Band::High => 3,
107        }]
108    }
109
110    /// Returns the confidence of the screen.
111    #[must_use]
112    pub fn confidence(&self) -> &Confidence<NoSpread> {
113        &self.confidence
114    }
115}
116
117/// Runs the screen on two channels of equal length, or on a sample range
118/// of them. `f0` is the fundamental from the caller; without it every
119/// `ipd_deg` is `None`.
120///
121/// # Errors
122///
123/// Never in practice; the reasons are built from counts and flags.
124pub fn stereo_screen(
125    left: &[f64],
126    right: &[f64],
127    rate: SampleRate,
128    f0: Option<Hz>,
129) -> Result<StereoScreen, ConfidenceError> {
130    let len = left.len().min(right.len());
131    let (left, right) = (&left[..len], &right[..len]);
132    let params = StftParams::fixed(SCREEN_FRAME_LEN);
133    // Activity from the energy of both channels, so an anti-phase stem with
134    // no mid still counts as sounding.
135    let energy: Vec<f64> = frame_rms(left, params)
136        .iter()
137        .zip(frame_rms(right, params))
138        .map(|(l, r)| f64::midpoint(l * l, r * r).sqrt())
139        .collect();
140    let active = ActiveFrames::from_rms(&energy, ACTIVE_GATE_DB);
141    let bin_hz = params.bin_width(rate).get();
142    let nyquist = rate.nyquist().get();
143    let f0 = f0.map(Hz::get).filter(|f| *f > 0.0);
144
145    // Per band: sums of |L|^2, |R|^2 and Re(L conj R) over all frames.
146    let mut sums = [[0.0_f64; 3]; 4];
147    // Per band: the per-frame phase means of the active frames.
148    let mut phases: [Vec<f64>; 4] = Default::default();
149    let mut stft_left = Stft::new(params);
150    let mut stft_right = Stft::new(params);
151    for frame in 0..params.frame_count(len) {
152        let l = stft_left.frame(left, frame).to_vec();
153        let r = stft_right.frame(right, frame);
154        let cross: Vec<Complex64> = l.iter().zip(r).map(|(a, b)| a * b.conj()).collect();
155        for (k, ((a, c), b)) in l.iter().zip(&cross).zip(r).enumerate() {
156            let f = cast::usize_to_f64(k) * bin_hz;
157            for (index, band) in Band::ALL.iter().enumerate() {
158                if band.contains(f) {
159                    sums[index][0] += a.norm_sqr();
160                    sums[index][1] += b.norm_sqr();
161                    sums[index][2] += c.re;
162                }
163            }
164        }
165        if let (Some(f0), true) = (f0, active.is_active(frame)) {
166            for (index, band) in Band::ALL.iter().enumerate() {
167                if let Some(phase) = frame_phase(&cross, f0, bin_hz, nyquist, *band) {
168                    phases[index].push(phase);
169                }
170            }
171        }
172    }
173
174    let full = StereoStats {
175        rho: pearson(left, right),
176        sm_db: side_over_mid(left, right),
177        ipd_deg: mean_phase(&phases[0]),
178    };
179    let mut bands = [full; 4];
180    for index in 1..4 {
181        let [ll, rr, lr] = sums[index];
182        bands[index] = StereoStats {
183            rho: (ll > 0.0 && rr > 0.0).then(|| Correlation::clamped(lr / (ll * rr).sqrt())),
184            sm_db: band_side_over_mid(ll, rr, lr),
185            ipd_deg: mean_phase(&phases[index]),
186        };
187    }
188    let reasons = vec![
189        Reason::new(ReasonCode::StereoBlindSpot, Signal::Flag(true))?,
190        Reason::new(
191            ReasonCode::ActiveFrames,
192            Signal::Count(u64::try_from(active.count()).unwrap_or(u64::MAX)),
193        )?,
194    ];
195    let confidence = if active.count() == 0 {
196        Confidence::no_estimate(vec![Reason::new(
197            ReasonCode::SilentStem,
198            Signal::Flag(true),
199        )?])?
200    } else {
201        Confidence::estimate(Level::Medium, None, reasons)?
202    };
203    Ok(StereoScreen { bands, confidence })
204}
205
206/// The weighted mean of `abs(angle(L conj R))` at the harmonic peaks of one
207/// frame, or `None` when no harmonic lies in the band or every weight is 0.
208fn frame_phase(cross: &[Complex64], f0: f64, bin_hz: f64, nyquist: f64, band: Band) -> Option<f64> {
209    let (mut weighted, mut total) = (0.0, 0.0);
210    for h in 1..=SCREEN_HARMONICS {
211        let centre_hz = cast::usize_to_f64(h) * f0;
212        if centre_hz >= nyquist || !band.contains(centre_hz) {
213            continue;
214        }
215        let centre = cast::round_to_usize(centre_hz / bin_hz).unwrap_or(0);
216        let lo = centre.saturating_sub(2);
217        let hi = (centre + 2).min(cross.len().saturating_sub(1));
218        let Some(peak) = cross[lo..=hi]
219            .iter()
220            .max_by(|a, b| a.norm().total_cmp(&b.norm()))
221        else {
222            continue;
223        };
224        let weight = peak.norm();
225        weighted += weight * peak.arg().abs().to_degrees();
226        total += weight;
227    }
228    (total > 0.0).then(|| weighted / total)
229}
230
231fn mean_phase(values: &[f64]) -> Option<Degrees> {
232    (!values.is_empty())
233        .then(|| Degrees::clamped(values.iter().sum::<f64>() / cast::usize_to_f64(values.len())))
234}
235
236/// Pearson correlation at zero lag, or `None` when a channel has no
237/// variance.
238fn pearson(left: &[f64], right: &[f64]) -> Option<Correlation> {
239    if left.is_empty() {
240        return None;
241    }
242    let n = cast::usize_to_f64(left.len());
243    let mean_l = left.iter().sum::<f64>() / n;
244    let mean_r = right.iter().sum::<f64>() / n;
245    let (mut ll, mut rr, mut lr) = (0.0, 0.0, 0.0);
246    for (l, r) in left.iter().zip(right) {
247        let (dl, dr) = (l - mean_l, r - mean_r);
248        ll += dl * dl;
249        rr += dr * dr;
250        lr += dl * dr;
251    }
252    (ll > 0.0 && rr > 0.0).then(|| Correlation::clamped(lr / (ll * rr).sqrt()))
253}
254
255/// `20 log10(rms(S) / rms(M))`, or `None` when either is zero.
256fn side_over_mid(left: &[f64], right: &[f64]) -> Option<Decibels> {
257    let (mut mid, mut side) = (0.0, 0.0);
258    for (l, r) in left.iter().zip(right) {
259        mid += ((l + r) / 2.0).powi(2);
260        side += ((l - r) / 2.0).powi(2);
261    }
262    power_ratio_db(side, mid)
263}
264
265/// Side over mid from band sums of `|L|^2`, `|R|^2` and `Re(L conj R)`.
266fn band_side_over_mid(ll: f64, rr: f64, lr: f64) -> Option<Decibels> {
267    power_ratio_db((ll + rr - 2.0 * lr).max(0.0), (ll + rr + 2.0 * lr).max(0.0))
268}
269
270fn power_ratio_db(numerator: f64, denominator: f64) -> Option<Decibels> {
271    // Relative to the larger, a residue of rounding is no side signal.
272    let floor = 1e-24 * numerator.max(denominator);
273    (numerator > floor && denominator > floor)
274        .then(|| Decibels::from_power_ratio(numerator / denominator).ok())
275        .flatten()
276}
277
278#[cfg(test)]
279mod tests {
280    use super::*;
281
282    #[test]
283    fn the_constants_pass_the_unit_check() {
284        assert_eq!(SUB_EDGE_HZ, Hz::new(120.0).unwrap());
285        assert_eq!(HIGH_EDGE_HZ, Hz::new(1000.0).unwrap());
286    }
287
288    fn rate() -> SampleRate {
289        SampleRate::new(44_100).unwrap()
290    }
291
292    fn tone(len: usize) -> Vec<f64> {
293        (0..len)
294            .map(|n| {
295                (cast::usize_to_f64(n) * 0.031).sin() + 0.3 * (cast::usize_to_f64(n) * 0.17).sin()
296            })
297            .collect()
298    }
299
300    #[test]
301    fn equal_channels_are_fully_correlated_with_no_side_and_no_phase() {
302        let x = tone(20_000);
303        let f0 = Hz::new(0.031 * 44_100.0 / std::f64::consts::TAU).ok();
304        let screen = stereo_screen(&x, &x, rate(), f0).unwrap();
305        let full = screen.band(Band::Full);
306        assert!((full.rho.unwrap().get() - 1.0).abs() < 1e-12);
307        assert_eq!(full.sm_db, None);
308        assert!(full.ipd_deg.unwrap().get().abs() < 1e-9);
309    }
310
311    #[test]
312    fn opposite_channels_are_anti_correlated_with_no_mid() {
313        let x = tone(20_000);
314        let negated: Vec<f64> = x.iter().map(|v| -v).collect();
315        let screen = stereo_screen(&x, &negated, rate(), None).unwrap();
316        let full = screen.band(Band::Full);
317        assert!((full.rho.unwrap().get() + 1.0).abs() < 1e-12);
318        assert_eq!(full.sm_db, None);
319        assert_eq!(full.ipd_deg, None);
320        // No mid, but both channels sound: an estimate, and with an f0 the
321        // phase reads 180 degrees.
322        assert!(screen.confidence().is_estimate());
323        let f0 = Hz::new(0.031 * 44_100.0 / std::f64::consts::TAU).unwrap();
324        let with_f0 = stereo_screen(&x, &negated, rate(), Some(f0)).unwrap();
325        let phase = with_f0.band(Band::Full).ipd_deg.unwrap().get();
326        assert!((phase - 180.0).abs() < 1e-6, "{phase}");
327    }
328
329    #[test]
330    fn a_silent_channel_has_no_correlation_and_one_channel_no_values() {
331        let x = tone(20_000);
332        let screen = stereo_screen(&x, &vec![0.0; 20_000], rate(), None).unwrap();
333        assert_eq!(screen.band(Band::Full).rho, None);
334        let mono = StereoScreen::one_channel().unwrap();
335        assert_eq!(mono.band(Band::High), StereoStats::default());
336        assert!(!mono.confidence().is_estimate());
337    }
338
339    #[test]
340    fn bands_split_at_120_hz_and_1_khz() {
341        assert!(Band::Sub.contains(119.9) && !Band::Sub.contains(120.0));
342        assert!(Band::LowMid.contains(120.0) && !Band::LowMid.contains(1000.0));
343        assert!(Band::High.contains(1000.0));
344        assert!(!Band::Full.contains(0.0));
345    }
346}
sound_to_patch/audio/stereo.rs

sound_to_patch/audio/
stereo.rs
