exg 0.0.5

EXG (EEG/ECG/EMG) preprocessing — native Rust DSP + FIF reader, numerical parity with MNE-Python
Documentation
//! Parity tests: verify exg filter design matches MNE-Python exactly.
//!
//! Ground-truth vectors generated by scripts/verify_firwin_design.py using
//! MNE-Python 1.11.0 + scipy.signal.firwin.
//!
//! Tests verify:
//! - Filter lengths match MNE
//! - Filter sums match MNE
//! - Individual coefficients match MNE to f32 precision

use exg::filter::design::*;

// ── MNE ground-truth coefficients ───────────────────────────────────────────

// From firwin_design_coeffs.json:
// lowpass_75_256: N=47, h[0]=0, h[1]=0.00115616...
const LP_75_256_FIRST3: [f64; 3] = [0.0, 0.0011561643652309693, -0.0006078593677176805];

// highpass_0p5_256: N=1691, h[0]=2.6964792e-05
const HP_05_256_FIRST3: [f64; 3] = [2.6964792235451508e-05, 2.7081967561582585e-05, 2.720056823081876e-05];

// bandpass_0p1_75_256: N=8449
const BP_01_75_FIRST3: [f64; 3] = [5.3976340822146335e-06, 5.402292600631122e-06, 5.406962396556493e-06];

// bandpass_1_40_256: N=845
const BP_1_40_FIRST3: [f64; 3] = [5.417010928228019e-05, 5.4638797362817064e-05, 5.5119126933146115e-05];

// notch_60_256: N=1691 (from create_filter, NOT _firwin_design standalone)
const NOTCH_60_FIRST3: [f64; 3] = [-5.21386110830481306179e-05, -2.09709410768981113283e-05, 4.86110882540915340520e-05];

// ── LOWPASS ─────────────────────────────────────────────────────────────────

#[test]
fn lowpass_75hz_256hz_length() {
    assert_eq!(design_lowpass(75.0, 256.0).len(), 47);
}

#[test]
fn lowpass_75hz_256hz_sum() {
    let h = design_lowpass(75.0, 256.0);
    let s: f64 = h.iter().map(|&v| v as f64).sum();
    approx::assert_abs_diff_eq!(s, 1.0, epsilon = 1e-6);
}

#[test]
fn lowpass_75hz_256hz_coefficients() {
    let h = design_lowpass(75.0, 256.0);
    for (i, &expected) in LP_75_256_FIRST3.iter().enumerate() {
        approx::assert_abs_diff_eq!(h[i] as f64, expected, epsilon = 1e-10);
    }
}

#[test]
fn lowpass_40hz_256hz_length() {
    assert_eq!(design_lowpass(40.0, 256.0).len(), 85);
}

#[test]
fn lowpass_100hz_512hz_length() {
    assert_eq!(design_lowpass(100.0, 512.0).len(), 69);
}

#[test]
fn lowpass_trans_bandwidth() {
    approx::assert_abs_diff_eq!(auto_trans_bandwidth_lowpass(75.0, 256.0), 18.75, epsilon = 1e-6);
    approx::assert_abs_diff_eq!(auto_trans_bandwidth_lowpass(40.0, 256.0), 10.0, epsilon = 1e-6);
    approx::assert_abs_diff_eq!(auto_trans_bandwidth_lowpass(100.0, 512.0), 25.0, epsilon = 1e-6);
    approx::assert_abs_diff_eq!(auto_trans_bandwidth_lowpass(120.0, 256.0), 8.0, epsilon = 1e-6);
}

// ── HIGHPASS ────────────────────────────────────────────────────────────────

#[test]
fn highpass_0p5_256_length() {
    assert_eq!(design_highpass(0.5, 256.0).len(), 1691);
}

#[test]
fn highpass_0p5_256_coefficients() {
    let h = design_highpass(0.5, 256.0);
    for (i, &expected) in HP_05_256_FIRST3.iter().enumerate() {
        approx::assert_abs_diff_eq!(h[i] as f64, expected, epsilon = 1e-10);
    }
}

#[test]
fn highpass_0p1_256_length() {
    assert_eq!(design_highpass(0.1, 256.0).len(), 8449);
}

#[test]
fn highpass_1p0_256_length() {
    assert_eq!(design_highpass(1.0, 256.0).len(), 845);
}

// ── BANDPASS ────────────────────────────────────────────────────────────────

#[test]
fn bandpass_0p1_75_256_length() {
    assert_eq!(design_bandpass(0.1, 75.0, 256.0).len(), 8449);
}

#[test]
fn bandpass_0p1_75_256_sum() {
    let h = design_bandpass(0.1, 75.0, 256.0);
    let s: f64 = h.iter().map(|&v| v as f64).sum();
    assert!(s.abs() < 1e-6, "bandpass sum should be ~0, got {s}");
}

#[test]
fn bandpass_0p1_75_256_coefficients() {
    let h = design_bandpass(0.1, 75.0, 256.0);
    for (i, &expected) in BP_01_75_FIRST3.iter().enumerate() {
        approx::assert_abs_diff_eq!(h[i] as f64, expected, epsilon = 1e-10);
    }
}

#[test]
fn bandpass_1_40_256_length() {
    assert_eq!(design_bandpass(1.0, 40.0, 256.0).len(), 845);
}

#[test]
fn bandpass_1_40_256_sum() {
    let h = design_bandpass(1.0, 40.0, 256.0);
    let s: f64 = h.iter().map(|&v| v as f64).sum();
    assert!(s.abs() < 1e-6, "bandpass sum = {s}");
}

#[test]
fn bandpass_1_40_256_coefficients() {
    let h = design_bandpass(1.0, 40.0, 256.0);
    for (i, &expected) in BP_1_40_FIRST3.iter().enumerate() {
        approx::assert_abs_diff_eq!(h[i] as f64, expected, epsilon = 1e-10);
    }
}

#[test]
fn bandpass_0p5_100_512_length() {
    assert_eq!(design_bandpass(0.5, 100.0, 512.0).len(), 3381);
}

// ── NOTCH ───────────────────────────────────────────────────────────────────

#[test]
fn notch_60hz_256hz_length() {
    assert_eq!(design_notch(60.0, 256.0, None, None).len(), 1691);
}

#[test]
fn notch_60hz_256hz_sum() {
    let h = design_notch(60.0, 256.0, None, None);
    let s: f64 = h.iter().map(|&v| v as f64).sum();
    approx::assert_abs_diff_eq!(s, 1.0, epsilon = 1e-3);
}

#[test]
fn notch_60hz_256hz_coefficients() {
    let h = design_notch(60.0, 256.0, None, None);
    for (i, &expected) in NOTCH_60_FIRST3.iter().enumerate() {
        approx::assert_abs_diff_eq!(h[i] as f64, expected, epsilon = 1e-10);
    }
}

#[test]
fn notch_50hz_256hz_length() {
    assert_eq!(design_notch(50.0, 256.0, None, None).len(), 1691);
}

#[test]
fn notch_60hz_512hz_length() {
    assert_eq!(design_notch(60.0, 512.0, None, None).len(), 3381);
}

// ── EDF READER ──────────────────────────────────────────────────────────────
// EDF reader tests use a test file if available

#[test]
fn edf_parity_test_file() {
    let path = "scripts/parity_vectors/test.edf";
    if !std::path::Path::new(path).exists() {
        eprintln!("Skipping EDF parity test (no test file)");
        return;
    }
    let raw = exg::edf::open_raw_edf(path).unwrap();
    // Should match MNE: 3 channels, 2560 samples, 256 Hz
    assert_eq!(raw.num_channels(), 3);
    assert_eq!(raw.num_samples(), 2560);
    approx::assert_abs_diff_eq!(raw.header.sample_rate, 256.0, epsilon = 0.1);
}