dspk_synapse/features/
conduction.rs

1//! Muscle Fiber Conduction Velocity (MFCV) estimation for planar HD-EMG grids (`conduction.rs`).
2//!
3//! Estimates action potential propagation speed $v$ (physiologically 2–7 m/s) along longitudinal
4//! muscle fiber columns of a 2D HD-EMG array using sub-sample parabolic cross-correlation between
5//! adjacent single-differential channels.
6
7use serde::{Deserialize, Serialize};
8use dspk_ops::math::{cross_correlation, peak_lag};
9
10/// Largest delay searched between adjacent channels, as a fraction of the snippet length.
11const MAX_LAG_FRACTION: usize = 3;
12/// Channel pairs with less energy product than this are skipped (flat signals).
13const MIN_PAIR_ENERGY: f32 = 1e-8;
14/// Normalized correlation a channel pair needs to count.
15const MIN_PAIR_CORRELATION: f32 = 0.2;
16/// Delays below this (samples) carry no direction and are skipped.
17const MIN_DELAY_SAMPLES: f32 = 1e-3;
18
19impl ConductionVelocityEstimate {
20    /// No estimate (too few channels or samples, or no correlated pair): every field NaN.
21    pub const UNDEFINED: Self = Self { velocity_m_per_s: f32::NAN, mean_delay_ms: f32::NAN, correlation_r: f32::NAN };
22}
23
24/// Estimated Muscle Fiber Conduction Velocity (MFCV) along a longitudinal HD-EMG electrode column.
25#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
26pub struct ConductionVelocityEstimate {
27    /// Propagation speed in meters per second ($\text{m/s}$, physiological human range $3.0 - 6.0\,\text{m/s}$).
28    pub velocity_m_per_s: f32,
29    /// Mean inter-electrode delay in milliseconds ($\text{ms}$).
30    pub mean_delay_ms: f32,
31    /// Normalized peak cross-correlation coefficient $r \in [0, 1]$ across longitudinal channel pairs.
32    pub correlation_r: f32,
33}
34
35/// Estimates Muscle Fiber Conduction Velocity (MFCV) from a multi-channel HD-EMG template or snippet
36/// (`[channels, samples]`, organized as `rows x cols` where `cols` runs along the muscle fiber axis)
37/// with inter-electrode spacing `pitch_um` ($\mu\text{m}$) and `sample_rate_hz`.
38pub fn estimate_hdemg_conduction_velocity(
39    waveform: &[f32],
40    rows: usize,
41    cols: usize,
42    samples: usize,
43    sample_rate_hz: f64,
44    pitch_um: f32,
45) -> ConductionVelocityEstimate {
46    assert_eq!(waveform.len(), rows * cols * samples);
47    if rows == 0 || cols < 3 || samples < 5 {
48        return ConductionVelocityEstimate::UNDEFINED;
49    }
50
51    // 1. Find the row with maximal RMS energy
52    let mut best_row = 0usize;
53    let mut best_energy = -1.0f32;
54    for r in 0..rows {
55        let mut e = 0.0f32;
56        for c in 0..cols {
57            let ch = r * cols + c;
58            let slice = &waveform[ch * samples..(ch + 1) * samples];
59            e += slice.iter().map(|v| v * v).sum::<f32>();
60        }
61        if e > best_energy {
62            best_energy = e;
63            best_row = r;
64        }
65    }
66
67    // 2. Compute longitudinal single-differential signals along `best_row` (cols - 1 pairs)
68    let num_sd = cols - 1;
69    let mut sd = vec![0.0f32; num_sd * samples];
70    for c in 0..num_sd {
71        let ch0 = best_row * cols + c;
72        let ch1 = best_row * cols + c + 1;
73        for t in 0..samples {
74            sd[c * samples + t] = waveform[ch1 * samples + t] - waveform[ch0 * samples + t];
75        }
76    }
77
78    // 3. Cross-correlate adjacent single-differential channels with sub-sample parabolic refinement
79    let max_lag = (samples / MAX_LAG_FRACTION).max(2);
80    let mut weighted_delay_samples = 0.0f32;
81    let mut weight_sum = 0.0f32;
82    let mut corr_sum = 0.0f32;
83    let mut pair_count = 0usize;
84
85    for c in 0..(num_sd - 1) {
86        let x = &sd[c * samples..(c + 1) * samples];
87        let y = &sd[(c + 1) * samples..(c + 2) * samples];
88        let ex: f32 = x.iter().map(|v| v * v).sum();
89        let ey: f32 = y.iter().map(|v| v * v).sum();
90        let denom = (ex * ey).sqrt();
91        if denom < MIN_PAIR_ENERGY {
92            continue;
93        }
94
95        let corrs: Vec<f32> = cross_correlation(x, y, max_lag).into_iter().map(|c| c / denom).collect();
96        let Some(peak) = peak_lag(&corrs, max_lag) else { continue };
97        let best_r = peak.value;
98        let lag_samples = peak.fractional_lag().abs();
99        if best_r > MIN_PAIR_CORRELATION && lag_samples > MIN_DELAY_SAMPLES {
100            let w = best_r * denom;
101            weighted_delay_samples += w * lag_samples;
102            weight_sum += w;
103            corr_sum += best_r;
104            pair_count += 1;
105        }
106    }
107
108    if weight_sum <= MIN_PAIR_ENERGY || pair_count == 0 {
109        return ConductionVelocityEstimate::UNDEFINED;
110    }
111
112    let mean_lag_samples = weighted_delay_samples / weight_sum;
113    let delay_sec = (mean_lag_samples as f64) / sample_rate_hz.max(1.0);
114    let pitch_m = (pitch_um as f64) * 1e-6;
115    // Delays are at least MIN_DELAY_SAMPLES, so the velocity is finite
116    let velocity_m_per_s = (pitch_m / delay_sec) as f32;
117
118    ConductionVelocityEstimate {
119        velocity_m_per_s,
120        mean_delay_ms: (delay_sec * 1e3) as f32,
121        correlation_r: corr_sum / (pair_count as f32),
122    }
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128
129    #[test]
130    fn test_hdemg_conduction_velocity_estimation() {
131        // 4x8 grid, 8000 um (8 mm) inter-electrode pitch, 4000 Hz sample rate
132        // True velocity = 4.0 m/s -> delay per 8 mm pitch = 2.0 ms = 8.0 samples
133        let (rows, cols, samples) = (4, 8, 96);
134        let pitch_um = 8000.0f32;
135        let fs = 4000.0f64;
136        let mut wave = vec![0.0f32; rows * cols * samples];
137
138        for c in 0..cols {
139            let ch = 1 * cols + c;
140            let center = 20.0 + (c as f32) * 8.0;
141            for t in 0..samples {
142                let u = (t as f32 - center) / 3.5;
143                // Biphasic propagating action potential
144                wave[ch * samples + t] = -100.0 * u * (-0.5 * u * u).exp() * (c as f32 * 0.05 + 1.0);
145            }
146        }
147
148        let est = estimate_hdemg_conduction_velocity(&wave, rows, cols, samples, fs, pitch_um);
149        assert!(
150            (est.velocity_m_per_s - 4.0).abs() < 0.3,
151            "est velocity = {} m/s, expected 4.0 m/s",
152            est.velocity_m_per_s
153        );
154        assert!(est.correlation_r > 0.9);
155    }
156}
dspk_synapse/features/conduction.rs

dspk_synapse/features/
conduction.rs
