sc_neurocore_engine/analysis/
spectral.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 — Spectral analysis of spike trains
8
9use rustfft::{num_complex::Complex64, FftPlanner};
10
11/// Power spectral density of a binary spike train.
12///
13/// Returns `(psd, freqs_hz)` where `psd[k] = |FFT[k]|^2 / n` and
14/// `freqs_hz` are the corresponding one-sided frequency bins.
15pub fn power_spectrum(binary_train: &[i32], dt: f64) -> (Vec<f64>, Vec<f64>) {
16    let n = binary_train.len();
17    if n < 2 {
18        return (vec![], vec![]);
19    }
20    let mean: f64 = binary_train.iter().map(|&s| s as f64).sum::<f64>() / n as f64;
21    let mut buf: Vec<Complex64> = binary_train
22        .iter()
23        .map(|&s| Complex64::new(s as f64 - mean, 0.0))
24        .collect();
25
26    let mut planner = FftPlanner::new();
27    let fft = planner.plan_fft_forward(n);
28    fft.process(&mut buf);
29
30    let n_rfft = n / 2 + 1;
31    let nf = n as f64;
32    let psd: Vec<f64> = buf[..n_rfft]
33        .iter()
34        .map(|c| (c.re * c.re + c.im * c.im) / nf)
35        .collect();
36    let freqs: Vec<f64> = (0..n_rfft).map(|k| k as f64 / (nf * dt)).collect();
37    (psd, freqs)
38}
39
40#[cfg(test)]
41mod tests {
42    use super::*;
43
44    #[test]
45    fn test_power_spectrum_basic() {
46        let train = vec![1, 0, 1, 0, 1, 0, 1, 0, 1, 0];
47        let (psd, freqs) = power_spectrum(&train, 0.001);
48        assert_eq!(psd.len(), 6);
49        assert_eq!(freqs.len(), 6);
50        assert!((freqs[0]).abs() < 1e-12);
51        // Nyquist at 500 Hz
52        assert!((freqs[5] - 500.0).abs() < 1e-6);
53    }
54
55    #[test]
56    fn test_power_spectrum_empty() {
57        let (psd, freqs) = power_spectrum(&[1], 0.001);
58        assert!(psd.is_empty());
59        assert!(freqs.is_empty());
60    }
61
62    #[test]
63    fn test_power_spectrum_silence() {
64        let train = vec![0; 100];
65        let (psd, _) = power_spectrum(&train, 0.001);
66        // All zeros -> PSD all zero
67        assert!(psd.iter().all(|&v| v.abs() < 1e-12));
68    }
69
70    #[test]
71    fn test_power_spectrum_dc_removed() {
72        // Constant 1 -> after mean subtraction all zero
73        let train = vec![1; 64];
74        let (psd, _) = power_spectrum(&train, 0.001);
75        assert!(psd.iter().all(|&v| v.abs() < 1e-12));
76    }
77
78    #[test]
79    fn test_power_spectrum_peak_at_periodic() {
80        // Strong periodicity at bin 5 in 20-sample signal
81        let mut train = vec![0i32; 20];
82        for i in (0..20).step_by(4) {
83            train[i] = 1;
84        }
85        let (psd, _) = power_spectrum(&train, 0.001);
86        // Peak should be at the periodic frequency
87        let peak_idx = psd
88            .iter()
89            .enumerate()
90            .skip(1) // skip DC
91            .max_by(|a, b| a.1.partial_cmp(b.1).unwrap())
92            .unwrap()
93            .0;
94        // Period-4 pattern in 20 samples -> frequency at n/period = 20/4 = 5
95        // But FFT index = frequency * n * dt / (1/(n*dt))...
96        // Actually for n=20, period=4: peak at k = n/period = 5
97        // The actual peak may be at k=5 or a harmonic depending on DC removal
98        assert!(
99            peak_idx == 5 || peak_idx == 10,
100            "Expected peak at index 5 or 10, got {peak_idx}"
101        );
102    }
103
104    #[test]
105    fn test_power_spectrum_length() {
106        let train = vec![0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1];
107        let (psd, freqs) = power_spectrum(&train, 0.001);
108        assert_eq!(psd.len(), 13 / 2 + 1);
109        assert_eq!(psd.len(), freqs.len());
110    }
111}
sc_neurocore_engine/analysis/spectral.rs

sc_neurocore_engine/analysis/
spectral.rs
