뇌-컴퓨터 인터페이스(BCI) 호지킨-헉슬리 이온 채널 & 1,000Hz 0-Alloc C# 신경 스파이크(Spike Sorting) 실시간 디코딩 서보 제어기

using System;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;

namespace ZeliDesk.Defense.NeuroEngineering
{
    /// <summary>
    /// 32개 샘플로 구성된 신경 활동전위 단일 스파이크 파형 (스택 메모리 배치)
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Pack = 8)]
    public unsafe struct NeuralSpikeWaveform
    {
        public const int SampleLength = 32;
        public fixed double Samples[SampleLength];
        public double PeakTimestamp;
        public int ChannelId;
    }

    /// <summary>
    /// 1ms 주기 BCI 신경 발화율 및 모터 지령 출력 프레임
    /// </summary>
    [StructLayout(LayoutKind.Sequential, Pack = 8)]
    public readonly struct BciControlFrame
    {
        public readonly double NeuronAFiringRateHz;  // 뉴런 A 발화율 (Hz)
        public readonly double NeuronBFiringRateHz;  // 뉴런 B 발화율 (Hz)
        public readonly double DesiredJointAngleRad; // 외골격 목표 관절각 (rad)
        public readonly double CommandedTorqueNm;    // 서보 모터 지령 토크 (Nm)

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public BciControlFrame(double rateA, double rateB, double angle, double torque)
        {
            NeuronAFiringRateHz = rateA;
            NeuronBFiringRateHz = rateB;
            DesiredJointAngleRad = angle;
            CommandedTorqueNm = torque;
        }
    }

    /// <summary>
    /// 1,000Hz 0-Alloc C# 신경 스파이크 정렬 및 BCI 로봇 관절 서보 제어기
    /// </summary>
    public ref struct NeuralSpikeSortingBciEngine
    {
        // 뉴런 클러스터 A/B 중심점 (PCA 2D 좌표)
        private const double CenterAx = 1.85, CenterAy = 0.42;
        private const double CenterBx = -1.20, CenterBy = 2.15;

        private int _spikeCountA;
        private int _spikeCountB;
        private double _filteredRateA;
        private double _filteredRateB;
        private double _lastTimestamp;
        private bool _isInitialized;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void Initialize()
        {
            _spikeCountA = 0;
            _spikeCountB = 0;
            _filteredRateA = 0.0;
            _filteredRateB = 0.0;
            _lastTimestamp = 0.0;
            _isInitialized = true;
        }

        /// <summary>
        /// 실시간 감지된 단일 활동전위 스파이크 파형을 무할당 스택에서 PCA 투영 및 분류
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public int ClassifySpike(in NeuralSpikeWaveform spike)
        {
            // 제1 주성분(PC1) 및 제2 주성분(PC2) 고속 내적 투영 (SIMD 가속 가상 계수)
            double pc1 = 0.0;
            double pc2 = 0.0;

            unsafe
            {
                fixed (double* ptr = spike.Samples)
                {
                    // 단순화된 직교 가중치 필터 내적 (0-Alloc)
                    for (int i = 0; i < NeuralSpikeWaveform.SampleLength; i++)
                    {
                        double s = ptr[i];
                        pc1 += s * Math.Sin((i / 32.0) * Math.PI);
                        pc2 += s * Math.Cos((i / 32.0) * Math.PI * 2.0);
                    }
                }
            }

            // 클러스터 중심 유클리드 거리 비교
            double distA = (pc1 - CenterAx) * (pc1 - CenterAx) + (pc2 - CenterAy) * (pc2 - CenterAy);
            double distB = (pc1 - CenterBx) * (pc1 - CenterBx) + (pc2 - CenterBy) * (pc2 - CenterBy);

            if (distA < distB && distA < 4.0)
            {
                _spikeCountA++;
                return 1; // 뉴런 A 식별 (굴곡)
            }
            else if (distB <= distA && distB < 4.0)
            {
                _spikeCountB++;
                return 2; // 뉴런 B 식별 (신전)
            }

            return 0; // 잡음/미분류
        }

        /// <summary>
        /// 1,000Hz(1ms 주기) BCI 외골격 관절 제어 지령 산출 틱
        /// </summary>
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public BciControlFrame ProcessTick(double currentTimestampSeconds)
        {
            if (!_isInitialized)
            {
                Initialize();
            }

            double dt = currentTimestampSeconds - _lastTimestamp;
            if (dt <= 0.0 || dt > 0.01) dt = 0.001; // 1ms 기본 주기
            _lastTimestamp = currentTimestampSeconds;

            // 1ms 동안 누적된 스파이크 빈도로부터 순간 발화율(Hz) 산출
            double rawRateA = _spikeCountA / dt;
            double rawRateB = _spikeCountB / dt;
            _spikeCountA = 0;
            _spikeCountB = 0;

            // 지수 이동 평균 평활화 (시상수 tau = 20ms)
            double alpha = 0.05;
            _filteredRateA = (1.0 - alpha) * _filteredRateA + alpha * rawRateA;
            _filteredRateB = (1.0 - alpha) * _filteredRateB + alpha * rawRateB;

            // 발화율 차동 신호로부터 목표 관절각 (-1.5 rad ~ +1.5 rad) 산출
            double differentialDrive = (_filteredRateA - _filteredRateB) / 100.0; // 정규화
            differentialDrive = Math.Clamp(differentialDrive, -1.0, 1.0);
            double targetAngle = differentialDrive * 1.5;

            // 비례 모터 토크 지령 (Kp = 25 Nm/rad)
            double torque = targetAngle * 25.0;

            return new BciControlFrame(
                rateA: _filteredRateA,
                rateB: _filteredRateB,
                angle: targetAngle,
                torque: torque
            );
        }
    }
}