secs_engine/
pair_phase.rs

1//! Corrected phase and cancellation checks between main-channel pairs (ADR 0035).
2
3use crate::Note;
4use secs_dsp::Complex64;
5
6/// ADR 0035. Inspection range in hertz.
7pub const BAND_HZ: (f64, f64) = (40.0, 200.0);
8/// One-third-octave centre frequencies in the inspection range.
9pub const CENTERS_HZ: [f64; 6] = [50.0, 63.0, 80.0, 100.0, 125.0, 160.0];
10/// Cancellation threshold in decibels.
11pub const LOSS_DB: f64 = -6.0;
12/// Maximum number of reported pairs.
13pub const MAX_PAIRS: usize = 6;
14/// Lower bound for cancellation values in decibels.
15pub const FLOOR_DB: f64 = -30.0;
16
17/// One band's `(low Hz, high Hz, phase degrees, cancellation dB)` result.
18pub type Band = (f64, f64, f64, f64);
19
20/// Measure phase difference and cancellation in each one-third-octave band.
21///
22/// Phase is the angle of `sum(a * conj(b))` in degrees in `(-180, 180]`.
23/// Cancellation is
24/// `10 * log10(sum(|a + b|^2) / sum((|a| + |b|)^2))`, bounded below by
25/// [`FLOOR_DB`]. Bands without a bin in all three inputs are omitted.
26///
27/// Neighbouring bands share their edge, the geometric mean of the two
28/// centres, so a run of cancelling bands reads as one span after rounding.
29pub fn bands(f: &[f64], a: &[Complex64], b: &[Complex64]) -> Vec<Band> {
30    let edges = band_edges();
31    CENTERS_HZ
32        .iter()
33        .enumerate()
34        .filter_map(|(k, _)| {
35            let lo = edges[k];
36            let hi = edges[k + 1];
37            let bins = f
38                .iter()
39                .zip(a)
40                .zip(b)
41                .filter(|((frequency, _), _)| {
42                    let last = k + 1 == CENTERS_HZ.len();
43                    lo <= **frequency && (**frequency < hi || (last && **frequency <= hi))
44                })
45                .map(|((_, a), b)| (*a, *b))
46                .collect::<Vec<_>>();
47            if bins.is_empty() {
48                return None;
49            }
50
51            let cross = bins.iter().map(|(a, b)| *a * b.conj()).sum::<Complex64>();
52            let mut phase_deg = cross.arg().to_degrees();
53            if phase_deg <= -180.0 {
54                phase_deg = 180.0;
55            }
56
57            let combined = bins.iter().map(|(a, b)| (*a + *b).norm_sqr()).sum::<f64>();
58            let separate = bins
59                .iter()
60                .map(|(a, b)| (a.norm() + b.norm()).powi(2))
61                .sum::<f64>();
62            // A band whose energies do not compute (NaN, overflow in either
63            // sum) is not a cancellation: an overflowing denominator alone
64            // would read as a complete one. Minus infinity from a finite
65            // denominator and zero sum is a complete cancellation and clamps
66            // to the floor.
67            if !phase_deg.is_finite() || !combined.is_finite() || !separate.is_finite() {
68                return None;
69            }
70            let loss_db = if separate == 0.0 {
71                0.0
72            } else {
73                10.0 * (combined / separate).log10()
74            };
75            Some((lo, hi, phase_deg, loss_db.max(FLOOR_DB)))
76        })
77        .collect()
78}
79
80/// The band edges in hertz: `BAND_HZ.0`, the geometric means of neighbouring
81/// centres, then `BAND_HZ.1`, so the outer bands reach the ends of the
82/// inspection range.
83pub fn band_edges() -> Vec<f64> {
84    (0..=CENTERS_HZ.len())
85        .map(|k| {
86            if k == 0 {
87                BAND_HZ.0
88            } else if k == CENTERS_HZ.len() {
89                BAND_HZ.1
90            } else {
91                (CENTERS_HZ[k - 1] * CENTERS_HZ[k]).sqrt()
92            }
93        })
94        .collect()
95}
96
97/// Report main pairs with at least one band at or below [`LOSS_DB`].
98///
99/// `responses` is in design role order. Reports are sorted by each pair's
100/// worst cancellation, then by design order, and limited to [`MAX_PAIRS`].
101pub fn check(f: &[f64], responses: &[(String, Vec<Complex64>)]) -> Vec<Note> {
102    let mut candidates = Vec::new();
103    for left in 0..responses.len() {
104        for right in left + 1..responses.len() {
105            let affected = bands(f, &responses[left].1, &responses[right].1)
106                .into_iter()
107                .filter(|band| band.3 <= LOSS_DB)
108                .collect::<Vec<_>>();
109            if let Some(worst) = affected.iter().map(|band| band.3).min_by(f64::total_cmp) {
110                candidates.push((worst, left, right, affected));
111            }
112        }
113    }
114    candidates.sort_by(|a, b| {
115        a.0.total_cmp(&b.0)
116            .then_with(|| a.1.cmp(&b.1))
117            .then_with(|| a.2.cmp(&b.2))
118    });
119    candidates
120        .into_iter()
121        .take(MAX_PAIRS)
122        .map(|(_, left, right, bands)| Note::PairPhase {
123            a: responses[left].0.clone(),
124            b: responses[right].0.clone(),
125            bands,
126        })
127        .collect()
128}
129
130#[cfg(test)]
131mod tests {
132    use super::*;
133    use std::f64::consts::PI;
134
135    fn response(phase_deg: f64) -> Vec<Complex64> {
136        CENTERS_HZ
137            .iter()
138            .map(|_| Complex64::from_polar(1.0, -phase_deg * PI / 180.0))
139            .collect()
140    }
141
142    fn assert_all_close(values: &[Band], phase_deg: f64, loss_db: f64) {
143        assert_eq!(values.len(), CENTERS_HZ.len());
144        for band in values {
145            assert!((band.2 - phase_deg).abs() < 1e-10, "{band:?}");
146            assert!((band.3 - loss_db).abs() < 0.03, "{band:?}");
147        }
148    }
149
150    #[test]
151    fn pair_phase_bands_share_edges_between_neighbours() {
152        let edges = band_edges();
153        assert_eq!(edges.len(), CENTERS_HZ.len() + 1);
154        let same = response(0.0);
155        let bands = bands(&CENTERS_HZ, &same, &same);
156        assert_eq!(bands.len(), CENTERS_HZ.len());
157        for (k, band) in bands.iter().enumerate() {
158            assert!(band.0 < CENTERS_HZ[k] && CENTERS_HZ[k] < band.1, "{band:?}");
159            if k + 1 < bands.len() {
160                assert_eq!(band.1, bands[k + 1].0);
161            }
162        }
163        assert_eq!(edges[0], BAND_HZ.0);
164        assert_eq!(edges[CENTERS_HZ.len()], BAND_HZ.1);
165        assert!((edges[3] - 89.4427).abs() < 0.01);
166        assert!((edges[4] - 111.8034).abs() < 0.01);
167        assert!((edges[5] - 141.4214).abs() < 0.01);
168    }
169
170    #[test]
171    fn pair_phase_bands_skip_non_finite_and_include_the_top_edge() {
172        let one = vec![Complex64::new(1.0, 0.0)];
173        let minus = vec![Complex64::new(-1.0, 0.0)];
174        let top = bands(&[BAND_HZ.1], &one, &minus);
175        assert_eq!(top.len(), 1);
176        assert_eq!(top[0].3, FLOOR_DB);
177        assert!(bands(&[BAND_HZ.0 - 0.01], &one, &minus).is_empty());
178        assert!(bands(&[100.0], &[Complex64::new(f64::NAN, 0.0)], &one).is_empty());
179        let huge = vec![Complex64::new(1e200, 0.0)];
180        assert!(bands(&[100.0], &huge, &huge).is_empty());
181        // Only the denominator overflows: not a cancellation either.
182        let big = vec![Complex64::new(1.3e154, 0.0)];
183        let big_minus = vec![Complex64::new(-1e153, 0.0)];
184        assert!(bands(&[100.0], &big, &big_minus).is_empty());
185        assert!(check(&[100.0], &[("L".into(), huge.clone()), ("R".into(), huge)]).is_empty());
186    }
187
188    #[test]
189    fn pair_phase_bands_measure_phase_difference_and_loss() {
190        let same = response(0.0);
191        assert_all_close(&bands(&CENTERS_HZ, &same, &same), 0.0, 0.0);
192        assert_all_close(
193            &bands(&CENTERS_HZ, &same, &response(90.0)),
194            90.0,
195            -3.010_299_956_639_812,
196        );
197        assert_all_close(
198            &bands(&CENTERS_HZ, &same, &response(120.0)),
199            120.0,
200            -6.020_599_913_279_624,
201        );
202        assert_all_close(
203            &bands(&CENTERS_HZ, &same, &response(180.0)),
204            180.0,
205            FLOOR_DB,
206        );
207    }
208
209    #[test]
210    fn pair_phase_check_reports_only_cancelling_pairs_worst_first() {
211        let responses = vec![
212            ("A".into(), response(0.0)),
213            ("B".into(), response(121.0)),
214            ("C".into(), response(60.0)),
215        ];
216        let notes = check(&CENTERS_HZ, &responses);
217        assert!(matches!(
218            notes.as_slice(),
219            [Note::PairPhase { a, b, bands }]
220                if a == "A" && b == "B" && bands.len() == CENTERS_HZ.len()
221        ));
222
223        let responses = [0.0, 30.0, 60.0, 180.0, 210.0, 240.0]
224            .into_iter()
225            .enumerate()
226            .map(|(index, phase)| (format!("R{index}"), response(phase)))
227            .collect::<Vec<_>>();
228        let notes = check(&CENTERS_HZ, &responses);
229        assert_eq!(notes.len(), MAX_PAIRS);
230        let worst = notes
231            .iter()
232            .map(|note| match note {
233                Note::PairPhase { bands, .. } => bands
234                    .iter()
235                    .map(|band| band.3)
236                    .min_by(f64::total_cmp)
237                    .unwrap(),
238                _ => unreachable!(),
239            })
240            .collect::<Vec<_>>();
241        assert!(worst.windows(2).all(|pair| pair[0] <= pair[1]));
242        assert!(worst.iter().take(3).all(|value| *value == FLOOR_DB));
243    }
244}
secs_engine/pair_phase.rs

secs_engine/
pair_phase.rs
