autoeq 0.4.36

//! Microphone phase calibration loader — GD-Opt v2 Phase GD-1f.
//!
//! Below ~50 Hz a USB measurement mic's own phase can drift ±30°, so
//! any bass-phase extraction that doesn't subtract the mic's response
//! attributes mic artefacts to the room. See
//! `docs/gd_opt_v2_plan.md` §2.6 and §2.8 (`"mic_phase_uncalibrated"`
//! advisory).
//!
//! This module provides:
//! 1. [`MicPhaseCalibration`] — the loaded 4-column calibration
//!    `(freq, mag_db, phase_deg, coherence)`.
//! 2. [`load_mic_phase_calibration`] — CSV loader with the same
//!    header-driven column discovery as
//!    [`crate::read::load_driver_measurement`], so callers can drop a
//!    calibration file authored by any tool that produces those four
//!    named columns.
//! 3. [`MicPhaseCalibration::apply_to_curve`] — in-place correction
//!    that subtracts the mic's magnitude and phase from a measured
//!    `Curve` and attenuates the curve's own coherence by the mic
//!    coherence.
//!
//! The struct, loader, and application are all read-only with
//! respect to the filesystem — nothing is mutated outside of the
//! `&mut Curve` the caller passes in.

use std::error::Error;
use std::fs::File;
use std::io::{BufRead, BufReader};
use std::path::Path;

use ndarray::Array1;

use crate::Curve;

/// Microphone calibration carrying both magnitude **and** phase
/// deviations at each frequency, plus a per-bin coherence that
/// bounds how much we trust each correction.
///
/// A "flat" mic has `mag_db[k] == 0.0` and `phase_deg[k] == 0.0`
/// across all bins; any deviation reports the mic's bias relative
/// to an ideal reference.
///
/// Frequencies must be strictly increasing so that
/// [`MicPhaseCalibration::apply_to_curve`] can do monotonic
/// interpolation. The loader enforces this.
#[derive(Debug, Clone, PartialEq)]
pub struct MicPhaseCalibration {
    /// Frequency points in Hz, strictly increasing.
    pub freq: Array1<f64>,
    /// Magnitude deviation in dB (positive = mic is louder than flat).
    pub mag_db: Array1<f64>,
    /// Phase deviation in degrees (positive = mic leads the reference).
    pub phase_deg: Array1<f64>,
    /// Per-bin coherence γ² from the calibration capture, in `[0, 1]`.
    /// Callers use this to down-weight corrections where the cal
    /// itself was noisy.
    pub coherence: Array1<f64>,
}

impl MicPhaseCalibration {
    /// Build an identity calibration (no deviation) on a given grid.
    /// Only useful for tests and as a no-op when a calibration file is
    /// missing — the caller still has to **choose** to apply this
    /// rather than treat the mic as uncalibrated.
    pub fn identity(freq: Array1<f64>) -> Self {
        let n = freq.len();
        Self {
            freq,
            mag_db: Array1::zeros(n),
            phase_deg: Array1::zeros(n),
            coherence: Array1::ones(n),
        }
    }

    /// Interpolate the calibration at a single frequency. Returns
    /// `(mag_db, phase_deg, coherence)` at that frequency, or `None`
    /// if the cal is empty. Uses piecewise linear interpolation in
    /// linear frequency with flat extrapolation at both ends (same
    /// policy as the CSV reader's `read_curve_from_csv`).
    pub fn sample_at(&self, freq_hz: f64) -> Option<(f64, f64, f64)> {
        let n = self.freq.len();
        if n == 0 {
            return None;
        }
        if !freq_hz.is_finite() {
            return None;
        }
        if freq_hz <= self.freq[0] {
            return Some((self.mag_db[0], self.phase_deg[0], self.coherence[0]));
        }
        if freq_hz >= self.freq[n - 1] {
            return Some((
                self.mag_db[n - 1],
                self.phase_deg[n - 1],
                self.coherence[n - 1],
            ));
        }
        // Binary search for the bracket.
        let idx = self
            .freq
            .as_slice()
            .expect("contiguous")
            .partition_point(|&f| f < freq_hz);
        let x0 = self.freq[idx - 1];
        let x1 = self.freq[idx];
        let dx = x1 - x0;
        if dx.abs() < f64::EPSILON {
            return Some((self.mag_db[idx], self.phase_deg[idx], self.coherence[idx]));
        }
        let t = (freq_hz - x0) / dx;
        Some((
            self.mag_db[idx - 1] * (1.0 - t) + self.mag_db[idx] * t,
            self.phase_deg[idx - 1] * (1.0 - t) + self.phase_deg[idx] * t,
            self.coherence[idx - 1] * (1.0 - t) + self.coherence[idx] * t,
        ))
    }

    /// Subtract the mic's magnitude and phase contributions from
    /// `curve`, in place.
    ///
    /// - `curve.spl[k] -= mag_db(curve.freq[k])`
    /// - `curve.phase[k] -= phase_deg(curve.freq[k])` (when phase is present)
    /// - `curve.coherence[k] *= cal_coherence(curve.freq[k])`
    ///   (when coherence is present) — noisy cal bins drag the
    ///   curve's own coherence down, which is exactly what the
    ///   confidence gate (GD-1g) needs.
    ///
    /// Silently no-op when `curve.freq.len() != curve.spl.len()` —
    /// malformed curves are rejected by the preceding CSV reader and
    /// shouldn't reach this path, but we fail closed rather than
    /// panic.
    pub fn apply_to_curve(&self, curve: &mut Curve) {
        let n = curve.freq.len();
        if n == 0 || curve.spl.len() != n {
            return;
        }
        for i in 0..n {
            let (mag, phase, coh) = match self.sample_at(curve.freq[i]) {
                Some(tuple) => tuple,
                None => continue,
            };
            curve.spl[i] -= mag;
            if let Some(ref mut phase_arr) = curve.phase
                && phase_arr.len() == n
            {
                phase_arr[i] -= phase;
            }
            if let Some(ref mut coh_arr) = curve.coherence
                && coh_arr.len() == n
            {
                coh_arr[i] *= coh;
            }
        }
    }
}

/// Load a 4-column microphone phase calibration CSV.
///
/// # CSV format
/// Header row required. Column discovery is header-name driven
/// (case-insensitive) so authors can reorder columns or name them in
/// any of the canonical SOTF spellings:
///
/// | Column | Recognised header names |
/// |---|---|
/// | freq | `frequency_hz`, `frequency`, `freq`, `hz` |
/// | mag_db | `mag_db`, `magnitude_db`, `magnitude`, `spl`, `spl_db`, `db` |
/// | phase_deg | `phase_deg`, `phase` |
/// | coherence | `coherence` |
///
/// Missing columns trigger an error — this is a strict 4-column
/// loader. For 2-column magnitude-only calibrations use the
/// pre-existing [`math_audio_dsp::analysis::MicrophoneCompensation`].
///
/// Frequencies must be strictly increasing. Rows with non-finite
/// values are dropped.
pub fn load_mic_phase_calibration(path: &Path) -> Result<MicPhaseCalibration, Box<dyn Error>> {
    let file = File::open(path)?;
    let reader = BufReader::new(file);

    let mut freq_col: Option<usize> = None;
    let mut mag_col: Option<usize> = None;
    let mut phase_col: Option<usize> = None;
    let mut coh_col: Option<usize> = None;
    let mut header_parsed = false;

    let mut freqs = Vec::new();
    let mut mags = Vec::new();
    let mut phases = Vec::new();
    let mut cohs = Vec::new();

    for (line_num, line) in reader.lines().enumerate() {
        let line = line?;
        let trimmed = line.trim();
        if trimmed.is_empty() || trimmed.starts_with('#') || trimmed.starts_with("//") {
            continue;
        }

        let parts: Vec<&str> = if trimmed.contains(',') {
            trimmed.split(',').map(|s| s.trim()).collect()
        } else {
            trimmed.split_whitespace().collect()
        };

        if !header_parsed {
            // A valid 4-column calibration MUST start with a named
            // header — without it we can't tell `phase_deg` from
            // `coherence`.
            for (idx, col_name) in parts.iter().enumerate() {
                let lower = col_name.to_lowercase();
                if coh_col.is_none() && lower == "coherence" {
                    coh_col = Some(idx);
                } else if phase_col.is_none()
                    && (lower.contains("phase") || lower == "phase_deg")
                {
                    phase_col = Some(idx);
                } else if freq_col.is_none()
                    && (lower.contains("freq") || lower == "hz" || lower == "frequency_hz")
                {
                    freq_col = Some(idx);
                } else if mag_col.is_none()
                    && (lower.contains("mag_db")
                        || lower.contains("magnitude")
                        || lower.contains("spl")
                        || lower == "db"
                        || lower == "spl_db")
                {
                    mag_col = Some(idx);
                }
            }
            if freq_col.is_none()
                || mag_col.is_none()
                || phase_col.is_none()
                || coh_col.is_none()
            {
                return Err(format!(
                    "mic phase calibration at {path:?} must have named columns for \
                     frequency / magnitude_db / phase_deg / coherence; got header {parts:?}"
                )
                .into());
            }
            header_parsed = true;
            continue;
        }

        let freq_idx = freq_col.unwrap();
        let mag_idx = mag_col.unwrap();
        let phase_idx = phase_col.unwrap();
        let coh_idx = coh_col.unwrap();
        let max_idx = freq_idx.max(mag_idx).max(phase_idx).max(coh_idx);
        if parts.len() <= max_idx {
            continue; // short row; skip
        }

        let (Ok(f), Ok(m), Ok(p), Ok(c)) = (
            parts[freq_idx].parse::<f64>(),
            parts[mag_idx].parse::<f64>(),
            parts[phase_idx].parse::<f64>(),
            parts[coh_idx].parse::<f64>(),
        ) else {
            log::debug!(
                "[mic_phase_calibration] Skipping malformed row {} in {path:?}: {:?}",
                line_num + 1,
                parts
            );
            continue;
        };
        if !(f.is_finite() && m.is_finite() && p.is_finite() && c.is_finite()) {
            continue;
        }
        freqs.push(f);
        mags.push(m);
        phases.push(p);
        cohs.push(c);
    }

    if freqs.is_empty() {
        return Err(format!(
            "mic phase calibration at {path:?} contained no valid data rows"
        )
        .into());
    }

    for pair in freqs.windows(2) {
        if pair[0] >= pair[1] {
            return Err(format!(
                "mic phase calibration at {path:?} frequencies are not strictly increasing \
                 (found {} before {})",
                pair[0], pair[1]
            )
            .into());
        }
    }

    Ok(MicPhaseCalibration {
        freq: Array1::from_vec(freqs),
        mag_db: Array1::from_vec(mags),
        phase_deg: Array1::from_vec(phases),
        coherence: Array1::from_vec(cohs),
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::io::Write;
    use tempfile::NamedTempFile;

    fn write_cal_csv(csv: &str) -> NamedTempFile {
        let mut f = NamedTempFile::new().unwrap();
        f.write_all(csv.as_bytes()).unwrap();
        f.flush().unwrap();
        f
    }

    #[test]
    fn loads_canonical_four_column_csv() {
        let csv = "\
frequency_hz,mag_db,phase_deg,coherence
20.0,2.0,-10.0,0.95
50.0,1.0,-5.0,0.98
200.0,0.0,0.0,0.99
2000.0,-0.5,2.0,0.99
20000.0,-3.0,15.0,0.90
";
        let f = write_cal_csv(csv);
        let cal = load_mic_phase_calibration(f.path()).unwrap();
        assert_eq!(cal.freq.len(), 5);
        assert!((cal.freq[0] - 20.0).abs() < 1e-9);
        assert!((cal.mag_db[0] - 2.0).abs() < 1e-9);
        assert!((cal.phase_deg[1] + 5.0).abs() < 1e-9);
        assert!((cal.coherence[4] - 0.90).abs() < 1e-9);
    }

    #[test]
    fn column_order_is_header_driven() {
        // Shuffled columns: the loader must key off names, not position.
        let csv = "\
phase_deg,coherence,frequency_hz,mag_db
-10.0,0.95,20.0,2.0
-5.0,0.98,50.0,1.0
";
        let f = write_cal_csv(csv);
        let cal = load_mic_phase_calibration(f.path()).unwrap();
        assert!((cal.freq[0] - 20.0).abs() < 1e-9);
        assert!((cal.phase_deg[0] + 10.0).abs() < 1e-9);
        assert!((cal.mag_db[0] - 2.0).abs() < 1e-9);
    }

    #[test]
    fn missing_column_rejects_file() {
        // No coherence column → the 4-col contract is violated.
        let csv = "\
frequency_hz,mag_db,phase_deg
20.0,2.0,-10.0
";
        let f = write_cal_csv(csv);
        assert!(load_mic_phase_calibration(f.path()).is_err());
    }

    #[test]
    fn non_monotonic_frequencies_reject_file() {
        let csv = "\
frequency_hz,mag_db,phase_deg,coherence
200.0,0.0,0.0,0.99
20.0,2.0,-10.0,0.95
";
        let f = write_cal_csv(csv);
        let err = load_mic_phase_calibration(f.path())
            .unwrap_err()
            .to_string();
        assert!(err.contains("not strictly increasing"), "got: {err}");
    }

    #[test]
    fn malformed_row_is_dropped_but_loader_still_succeeds() {
        let csv = "\
frequency_hz,mag_db,phase_deg,coherence
20.0,2.0,-10.0,0.95
50.0,not-a-number,-5.0,0.98
200.0,0.0,0.0,0.99
";
        let f = write_cal_csv(csv);
        let cal = load_mic_phase_calibration(f.path()).unwrap();
        assert_eq!(cal.freq.len(), 2);
        assert!((cal.freq[1] - 200.0).abs() < 1e-9);
    }

    #[test]
    fn sample_at_exact_node_matches_stored_value() {
        let cal = MicPhaseCalibration {
            freq: Array1::from_vec(vec![20.0, 200.0, 2000.0]),
            mag_db: Array1::from_vec(vec![2.0, 0.0, -3.0]),
            phase_deg: Array1::from_vec(vec![-10.0, 0.0, 5.0]),
            coherence: Array1::from_vec(vec![0.95, 0.99, 0.90]),
        };
        let (m, p, c) = cal.sample_at(200.0).unwrap();
        assert!((m - 0.0).abs() < 1e-9);
        assert!((p - 0.0).abs() < 1e-9);
        assert!((c - 0.99).abs() < 1e-9);
    }

    #[test]
    fn sample_at_midpoint_interpolates() {
        let cal = MicPhaseCalibration {
            freq: Array1::from_vec(vec![100.0, 200.0]),
            mag_db: Array1::from_vec(vec![0.0, 4.0]),
            phase_deg: Array1::from_vec(vec![0.0, 20.0]),
            coherence: Array1::from_vec(vec![1.0, 0.5]),
        };
        // At 150 Hz we're halfway between 100 and 200 → all three
        // values should interpolate linearly.
        let (m, p, c) = cal.sample_at(150.0).unwrap();
        assert!((m - 2.0).abs() < 1e-9, "mag @ 150 Hz should be 2 dB, got {m}");
        assert!(
            (p - 10.0).abs() < 1e-9,
            "phase @ 150 Hz should be 10°, got {p}"
        );
        assert!(
            (c - 0.75).abs() < 1e-9,
            "coherence @ 150 Hz should be 0.75, got {c}"
        );
    }

    #[test]
    fn sample_at_below_min_returns_first() {
        let cal = MicPhaseCalibration {
            freq: Array1::from_vec(vec![50.0, 500.0]),
            mag_db: Array1::from_vec(vec![1.0, 2.0]),
            phase_deg: Array1::from_vec(vec![-10.0, 5.0]),
            coherence: Array1::from_vec(vec![0.9, 0.95]),
        };
        // Below the cal's min frequency: flat extrapolation.
        let (m, _, _) = cal.sample_at(10.0).unwrap();
        assert!((m - 1.0).abs() < 1e-9);
    }

    #[test]
    fn sample_at_above_max_returns_last() {
        let cal = MicPhaseCalibration {
            freq: Array1::from_vec(vec![50.0, 500.0]),
            mag_db: Array1::from_vec(vec![1.0, 2.0]),
            phase_deg: Array1::from_vec(vec![-10.0, 5.0]),
            coherence: Array1::from_vec(vec![0.9, 0.95]),
        };
        let (m, _, _) = cal.sample_at(50_000.0).unwrap();
        assert!((m - 2.0).abs() < 1e-9);
    }

    #[test]
    fn identity_is_transparent() {
        let freq = Array1::from_vec(vec![20.0, 200.0, 2000.0]);
        let cal = MicPhaseCalibration::identity(freq.clone());
        let mut curve = Curve {
            freq: freq.clone(),
            spl: Array1::from_vec(vec![60.0, 80.0, 90.0]),
            phase: Some(Array1::from_vec(vec![-30.0, 0.0, 45.0])),
            coherence: Some(Array1::from_vec(vec![0.98, 0.99, 0.97])),
            ..Default::default()
        };
        let before = curve.clone();
        cal.apply_to_curve(&mut curve);
        assert_eq!(before.spl, curve.spl);
        assert_eq!(before.phase, curve.phase);
        assert_eq!(before.coherence, curve.coherence);
    }

    #[test]
    fn apply_subtracts_mag_and_phase() {
        // Curve freq grid = cal freq grid → exact subtraction.
        let freq = Array1::from_vec(vec![20.0, 200.0, 2000.0]);
        let cal = MicPhaseCalibration {
            freq: freq.clone(),
            mag_db: Array1::from_vec(vec![2.0, 0.0, -3.0]),
            phase_deg: Array1::from_vec(vec![-10.0, 0.0, 5.0]),
            coherence: Array1::from_vec(vec![0.5, 1.0, 0.8]),
        };
        let mut curve = Curve {
            freq: freq.clone(),
            spl: Array1::from_vec(vec![60.0, 80.0, 90.0]),
            phase: Some(Array1::from_vec(vec![0.0, 0.0, 0.0])),
            coherence: Some(Array1::from_vec(vec![1.0, 1.0, 1.0])),
            ..Default::default()
        };
        cal.apply_to_curve(&mut curve);
        // spl: subtract mic's magnitude
        assert!((curve.spl[0] - 58.0).abs() < 1e-9); // 60 - 2
        assert!((curve.spl[1] - 80.0).abs() < 1e-9); // 80 - 0
        assert!((curve.spl[2] - 93.0).abs() < 1e-9); // 90 - (-3)
        // phase: subtract mic's phase
        let phase = curve.phase.as_ref().unwrap();
        assert!((phase[0] - 10.0).abs() < 1e-9); // 0 - (-10)
        assert!((phase[1] - 0.0).abs() < 1e-9);
        assert!((phase[2] + 5.0).abs() < 1e-9); // 0 - 5
        // coherence: multiply by mic's coherence
        let coh = curve.coherence.as_ref().unwrap();
        assert!((coh[0] - 0.5).abs() < 1e-9);
        assert!((coh[1] - 1.0).abs() < 1e-9);
        assert!((coh[2] - 0.8).abs() < 1e-9);
    }

    #[test]
    fn apply_skips_when_phase_or_coherence_absent() {
        // Without phase/coherence on the Curve, the cal should only
        // modify `spl` and not panic.
        let freq = Array1::from_vec(vec![20.0, 200.0]);
        let cal = MicPhaseCalibration {
            freq: freq.clone(),
            mag_db: Array1::from_vec(vec![2.0, 0.0]),
            phase_deg: Array1::from_vec(vec![-10.0, 0.0]),
            coherence: Array1::from_vec(vec![0.5, 1.0]),
        };
        let mut curve = Curve {
            freq,
            spl: Array1::from_vec(vec![60.0, 80.0]),
            phase: None,
            coherence: None,
            ..Default::default()
        };
        cal.apply_to_curve(&mut curve);
        assert_eq!(curve.phase, None);
        assert_eq!(curve.coherence, None);
        assert!((curve.spl[0] - 58.0).abs() < 1e-9);
    }
}