sc_neurocore_engine/analysis/
lfp.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later
2// Commercial license available
3// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
4// © Code 2020–2026 Miroslav Šotek. All rights reserved.
5// ORCID: 0009-0009-3560-0851
6// Contact: www.anulum.li | protoscience@anulum.li
7// SC-NeuroCore — Spike-LFP coupling: phase locking, coherence, phase
8
9use rustfft::{num_complex::Complex64, FftPlanner};
10use std::f64::consts::PI;
11
12/// Compute instantaneous phase via analytic signal (one-sided FFT Hilbert).
13fn instantaneous_phase(signal: &[f64]) -> Vec<f64> {
14    let n = signal.len();
15    if n == 0 {
16        return vec![];
17    }
18    let mut buf: Vec<Complex64> = signal.iter().map(|&s| Complex64::new(s, 0.0)).collect();
19    let mut planner = FftPlanner::new();
20    let fft_fwd = planner.plan_fft_forward(n);
21    fft_fwd.process(&mut buf);
22
23    // Zero negative frequencies, double positive (analytic signal)
24    buf[0] = Complex64::new(0.0, 0.0); // DC = 0 after centering isn't required here
25    for k in 1..n {
26        if k <= n / 2 {
27            buf[k] *= 2.0;
28        } else {
29            buf[k] = Complex64::new(0.0, 0.0);
30        }
31    }
32    if n > 0 {
33        buf[0] = Complex64::new(signal.iter().sum::<f64>(), 0.0); // preserve DC? No — match Python
34    }
35    // Actually match the Python: fft * 2 * (arange(n) > 0)
36    // So: bin 0 is zeroed, bins 1..n are doubled
37    // Let me redo:
38    let mut buf2: Vec<Complex64> = signal.iter().map(|&s| Complex64::new(s, 0.0)).collect();
39    fft_fwd.process(&mut buf2);
40    buf2[0] = Complex64::new(0.0, 0.0);
41    for k in 1..n {
42        buf2[k] *= 2.0;
43    }
44    let fft_inv = planner.plan_fft_inverse(n);
45    fft_inv.process(&mut buf2);
46    let scale = 1.0 / n as f64;
47    buf2.iter().map(|c| (c * scale).arg()).collect()
48}
49
50/// Phase locking value (PLV) between spikes and LFP phase.
51///
52/// PLV = |mean(exp(j * phase_at_spikes))|.
53pub fn phase_locking_value(binary_train: &[i32], lfp_signal: &[f64]) -> f64 {
54    let n = binary_train.len().min(lfp_signal.len());
55    if n == 0 {
56        return 0.0;
57    }
58    let phase = instantaneous_phase(&lfp_signal[..n]);
59    let spike_idx: Vec<usize> = (0..n).filter(|&i| binary_train[i] > 0).collect();
60    if spike_idx.is_empty() {
61        return 0.0;
62    }
63    let mut sum_re = 0.0;
64    let mut sum_im = 0.0;
65    for &i in &spike_idx {
66        sum_re += phase[i].cos();
67        sum_im += phase[i].sin();
68    }
69    let count = spike_idx.len() as f64;
70    ((sum_re / count).powi(2) + (sum_im / count).powi(2)).sqrt()
71}
72
73/// Spike-field coherence (SFC).
74///
75/// Returns `(coherence, freqs_hz)`. SFC = |S_xy|^2 / (S_xx * S_yy).
76pub fn spike_field_coherence(
77    binary_train: &[i32],
78    lfp_signal: &[f64],
79    dt: f64,
80) -> (Vec<f64>, Vec<f64>) {
81    let n = binary_train.len().min(lfp_signal.len());
82    if n < 2 {
83        return (vec![], vec![]);
84    }
85    let mean_a: f64 = binary_train[..n].iter().map(|&s| s as f64).sum::<f64>() / n as f64;
86    let mean_b: f64 = lfp_signal[..n].iter().sum::<f64>() / n as f64;
87
88    let mut buf_a: Vec<Complex64> = binary_train[..n]
89        .iter()
90        .map(|&s| Complex64::new(s as f64 - mean_a, 0.0))
91        .collect();
92    let mut buf_b: Vec<Complex64> = lfp_signal[..n]
93        .iter()
94        .map(|&s| Complex64::new(s - mean_b, 0.0))
95        .collect();
96
97    let mut planner = FftPlanner::new();
98    let fft = planner.plan_fft_forward(n);
99    fft.process(&mut buf_a);
100    fft.process(&mut buf_b);
101
102    let n_rfft = n / 2 + 1;
103    let mut sfc = vec![0.0f64; n_rfft];
104    for k in 0..n_rfft {
105        let sab = buf_a[k] * buf_b[k].conj();
106        let saa = buf_a[k].norm_sqr();
107        let sbb = buf_b[k].norm_sqr();
108        let denom = saa * sbb;
109        sfc[k] = if denom > 1e-30 {
110            sab.norm_sqr() / denom
111        } else {
112            0.0
113        };
114    }
115    let freqs: Vec<f64> = (0..n_rfft).map(|k| k as f64 / (n as f64 * dt)).collect();
116    (sfc, freqs)
117}
118
119/// Histogram of LFP phase at spike times.
120///
121/// Returns `(counts, bin_centres_rad)` with bins spanning [-pi, pi].
122pub fn spike_phase_histogram(
123    binary_train: &[i32],
124    lfp_signal: &[f64],
125    n_bins: usize,
126) -> (Vec<i64>, Vec<f64>) {
127    let n = binary_train.len().min(lfp_signal.len());
128    let phase = instantaneous_phase(&lfp_signal[..n]);
129
130    let edges: Vec<f64> = (0..=n_bins)
131        .map(|k| -PI + k as f64 * 2.0 * PI / n_bins as f64)
132        .collect();
133    let centres: Vec<f64> = (0..n_bins)
134        .map(|k| (edges[k] + edges[k + 1]) / 2.0)
135        .collect();
136
137    let mut hist = vec![0i64; n_bins];
138    for i in 0..n {
139        if binary_train[i] > 0 {
140            let p = phase[i];
141            let mut k = ((p + PI) / (2.0 * PI) * n_bins as f64).floor() as usize;
142            if k >= n_bins {
143                k = n_bins - 1;
144            }
145            hist[k] += 1;
146        }
147    }
148    (hist, centres)
149}
150
151#[cfg(test)]
152mod tests {
153    use super::*;
154
155    fn make_lfp(n: usize, freq_hz: f64, dt: f64) -> Vec<f64> {
156        (0..n)
157            .map(|i| (2.0 * PI * freq_hz * i as f64 * dt).sin())
158            .collect()
159    }
160
161    #[test]
162    fn test_instantaneous_phase_length() {
163        let lfp = make_lfp(100, 10.0, 0.001);
164        let phase = instantaneous_phase(&lfp);
165        assert_eq!(phase.len(), 100);
166    }
167
168    #[test]
169    fn test_phase_locking_value_basic() {
170        let n = 1000;
171        let lfp = make_lfp(n, 10.0, 0.001);
172        // Spikes locked to peaks (phase ~ 0)
173        let phase = instantaneous_phase(&lfp);
174        let mut train = vec![0i32; n];
175        for i in 0..n {
176            if phase[i].abs() < 0.3 {
177                train[i] = 1;
178            }
179        }
180        let plv = phase_locking_value(&train, &lfp);
181        assert!(
182            plv > 0.5,
183            "PLV={plv} should be high for phase-locked spikes"
184        );
185    }
186
187    #[test]
188    fn test_phase_locking_value_no_spikes() {
189        let lfp = make_lfp(100, 10.0, 0.001);
190        let train = vec![0i32; 100];
191        assert_eq!(phase_locking_value(&train, &lfp), 0.0);
192    }
193
194    #[test]
195    fn test_spike_field_coherence_shape() {
196        let n = 200;
197        let lfp = make_lfp(n, 20.0, 0.001);
198        let mut train = vec![0i32; n];
199        for i in (0..n).step_by(10) {
200            train[i] = 1;
201        }
202        let (sfc, freqs) = spike_field_coherence(&train, &lfp, 0.001);
203        assert_eq!(sfc.len(), n / 2 + 1);
204        assert_eq!(freqs.len(), sfc.len());
205        assert!(sfc.iter().all(|&v| (0.0..=1.0 + 1e-10).contains(&v)));
206    }
207
208    #[test]
209    fn test_spike_field_coherence_empty() {
210        let (sfc, freqs) = spike_field_coherence(&[1], &[0.5], 0.001);
211        assert!(sfc.is_empty());
212        assert!(freqs.is_empty());
213    }
214
215    #[test]
216    fn test_spike_phase_histogram_shape() {
217        let n = 500;
218        let lfp = make_lfp(n, 10.0, 0.001);
219        let mut train = vec![0i32; n];
220        for i in (0..n).step_by(5) {
221            train[i] = 1;
222        }
223        let (hist, centres) = spike_phase_histogram(&train, &lfp, 36);
224        assert_eq!(hist.len(), 36);
225        assert_eq!(centres.len(), 36);
226        let total: i64 = hist.iter().sum();
227        let expected: i64 = train.iter().map(|&v| v as i64).sum();
228        assert_eq!(total, expected);
229    }
230
231    #[test]
232    fn test_spike_phase_histogram_no_spikes() {
233        let lfp = make_lfp(100, 10.0, 0.001);
234        let train = vec![0i32; 100];
235        let (hist, _) = spike_phase_histogram(&train, &lfp, 12);
236        assert!(hist.iter().all(|&v| v == 0));
237    }
238
239    #[test]
240    fn test_phase_locking_random_spikes_low() {
241        // Random spikes -> PLV should be low
242        let n = 2000;
243        let lfp = make_lfp(n, 15.0, 0.001);
244        let mut train = vec![0i32; n];
245        let mut rng = 42u64;
246        for i in 0..n {
247            rng = rng.wrapping_mul(6364136223846793005).wrapping_add(1);
248            if rng.is_multiple_of(20) {
249                train[i] = 1;
250            }
251        }
252        let plv = phase_locking_value(&train, &lfp);
253        assert!(plv < 0.3, "PLV={plv} should be low for random spikes");
254    }
255}
sc_neurocore_engine/analysis/lfp.rs

sc_neurocore_engine/analysis/
lfp.rs
