On device adaptation

.gitignore

@@ -84,3 +84,7 @@ __pycache__/
 *.pyc
 figures/
 .venv/
+
+# Generated CAP-window caches & trained model weights (large, reproducible)
+scripts/cap_pipeline/caches/
+*.pt

scripts/Arduino/ArduinoOnDeviceTrain/ArduinoOnDeviceTrain.ino

@@ -0,0 +1,275 @@
+/*
+  ArduinoOnDeviceTrain.ino
+
+  Detects CAPs live from the streamed nerve signal and TRAINS a small model
+  on-device (rather than running a frozen one). Task: pricking vs movement
+  (movement = dorsiflexion + plantarflexion).
+
+  For each detected CAP we extract 113 features, standardise them with running
+  stats, predict, and then do one SGD step using the label that came with the
+  stream. So the model learns while the data streams in.
+
+  Stream protocol (same as stream_raw.py, 119 bytes):
+    [0xAA][0xBB] sync
+    [label]      1 byte: 0=dorsi 1=plantar 2=pricking 254=eval-only 255=unknown
+    [ts]         4 bytes uint32 LE
+    [ch0..ch55]  56 x uint16 LE raw Intan ADC (uV = 0.195*(raw-32768))
+
+  Per CAP it prints:
+    CAP ts=N pred=P label=Y p=PROB acc=ACC
+      pred/label: 1=pricking 0=movement (label=-1 if unknown)
+      acc: running accuracy over labelled CAPs so far
+
+  LED: cyan idle, yellow armed, magenta computing, blue=pricking, green=movement.
+*/
+#include <Arduino.h>
+#include <math.h>
+#include <string.h>
+#include <stdint.h>
+#include "ondevice_model.h"
+
+// ── LED ──────────────────────────────────────────────────────────────────────
+#define PHASE_IDLE      0   // cyan
+#define PHASE_ARMED     1   // yellow
+#define PHASE_COMPUTE   2   // magenta
+#define PHASE_MOVEMENT  3   // green
+#define PHASE_PRICKING  4   // blue
+
+static uint32_t result_hold_until = 0;
+
+static void set_led(int phase) {
+  bool r = false, g = false, b = false;
+  switch (phase) {
+    case PHASE_IDLE:     g = true; b = true; break;   // cyan
+    case PHASE_ARMED:    r = true; g = true; break;   // yellow
+    case PHASE_COMPUTE:  r = true; b = true; break;   // magenta
+    case PHASE_MOVEMENT: g = true; break;             // green
+    case PHASE_PRICKING: b = true; break;             // blue
+  }
+  digitalWrite(LEDR, r ? LOW : HIGH);
+  digitalWrite(LEDG, g ? LOW : HIGH);
+  digitalWrite(LEDB, b ? LOW : HIGH);
+}
+
+// ── Signal processing (live CAP detection) ───────────────────────────────────
+#define N_CH      56
+#define N_RINGS    7
+#define CIRC_SIZE 150   // >= 100
+
+// 2 kHz notch, Q=30, fs=30 kHz: [b0,b1,b2,a1,a2]
+static const float NOTCH[5] = {
+   0.9930669720f, -1.8144236429f,  0.9930669720f,
+  -1.8144236429f,  0.9861339440f
+};
+// 1-3 kHz bandpass, 6th-order Butterworth, fs=30 kHz (6 SOS sections)
+static const float BP[6][5] = {
+  { 0.0000407482f,  0.0000814964f,  0.0000407482f, -1.4614934920f,  0.6136710879f },
+  { 1.0000000000f,  2.0000000000f,  1.0000000000f, -1.4091053901f,  0.6623962038f },
+  { 1.0000000000f,  2.0000000000f,  1.0000000000f, -1.6315199139f,  0.7102771058f },
+  { 1.0000000000f, -2.0000000000f,  1.0000000000f, -1.7840692078f,  0.8352552241f },
+  { 1.0000000000f, -2.0000000000f,  1.0000000000f, -1.5099810120f,  0.8562240048f },
+  { 1.0000000000f, -2.0000000000f,  1.0000000000f, -1.9018207709f,  0.9454426128f },
+};
+
+static float    circ_buf[N_CH][CIRC_SIZE];
+static int      circ_write = 0;
+static float    notch_s[N_RINGS][2];
+static float    bp_s[N_RINGS][6][2];
+static float    daa_hist[8][N_RINGS];
+static int      daa_write = 0;
+static float    ema_abs   = 5.0f;
+#define EMA_ALPHA 0.005f
+static int      refractory = 0;
+static int      countdown  = -1;
+static uint32_t peak_ts    = 0;
+static uint8_t  peak_label = 255;
+
+// ── Model + running stats ────────────────────────────────────────────────────
+static RunningStandardizer od_std;
+static LogRegSGD           od_model;
+static float               od_window[OD_N_CH * OD_WIN];   // 56*100, static (22 KB)
+static long                seen = 0;       // labelled CAPs
+static long                correct = 0;
+
+static inline float biquad(float x, const float* c, float* s) {
+  float w = x - c[3]*s[0] - c[4]*s[1];
+  float y = c[0]*w + c[1]*s[0] + c[2]*s[1];
+  s[1] = s[0];
+  s[0] = w;
+  return y;
+}
+
+static void classify_cap(uint32_t ts, uint8_t label);
+
+static void process_sample(const uint16_t* raw, uint32_t ts, uint8_t label) {
+  // 1. ADC -> uV
+  float uv[N_CH];
+  for (int c = 0; c < N_CH; c++) uv[c] = 0.195f * ((float)raw[c] - 32768.f);
+
+  // 2. Tripolar reference (mean of outer rings 0 and 6)
+  float ref = 0.f;
+  for (int c = 0;  c < 8;  c++) ref += uv[c];
+  for (int c = 48; c < 56; c++) ref += uv[c];
+  ref /= 16.f;
+  for (int c = 0; c < N_CH; c++) uv[c] -= ref;
+
+  // 3. Circular buffer of referenced signal (for window extraction)
+  for (int c = 0; c < N_CH; c++) circ_buf[c][circ_write] = uv[c];
+  circ_write = (circ_write + 1) % CIRC_SIZE;
+
+  // 4. Ring averages
+  float ring_avg[N_RINGS];
+  for (int r = 0; r < N_RINGS; r++) {
+    float s = 0.f;
+    for (int c = 0; c < 8; c++) s += uv[r*8 + c];
+    ring_avg[r] = s * 0.125f;
+  }
+
+  // 5. Causal notch + bandpass on each ring average
+  float filt[N_RINGS];
+  for (int r = 0; r < N_RINGS; r++) {
+    float x = biquad(ring_avg[r], NOTCH, notch_s[r]);
+    for (int s = 0; s < 6; s++) x = biquad(x, BP[s], bp_s[r][s]);
+    filt[r] = x;
+  }
+
+  // 6. Delay-and-add
+  for (int r = 0; r < N_RINGS; r++) daa_hist[daa_write][r] = filt[r];
+  float daa = 0.f;
+  for (int r = 0; r < N_RINGS; r++) {
+    int hidx = (daa_write - r + 8) % 8;
+    daa += daa_hist[hidx][r];
+  }
+  daa_write = (daa_write + 1) % 8;
+
+  // 7. EMA noise floor -> adaptive threshold (~4 sigma)
+  float abs_daa = fabsf(daa);
+  ema_abs = (1.f - EMA_ALPHA) * ema_abs + EMA_ALPHA * abs_daa;
+  float threshold = 5.01f * ema_abs;
+
+  // 8. Refractory
+  if (refractory > 0) refractory--;
+
+  // 9. Threshold crossing -> arm
+  if (countdown < 0 && refractory == 0 && abs_daa > threshold) {
+    countdown  = 50;
+    peak_ts    = ts;
+    peak_label = label;
+    refractory = 150;
+    set_led(PHASE_ARMED);
+  }
+
+  // 10. Countdown -> extract window and train/predict
+  if (countdown > 0) {
+    countdown--;
+  } else if (countdown == 0) {
+    countdown = -1;
+    classify_cap(peak_ts, peak_label);
+  }
+
+  // 11. Back to idle when the result colour expires
+  if (countdown < 0 && (uint32_t)millis() >= result_hold_until) set_led(PHASE_IDLE);
+}
+
+static void classify_cap(uint32_t ts, uint8_t label) {
+  set_led(PHASE_COMPUTE);
+
+  // Build the 56x100 window (row-major c*WIN+t, t=0 oldest .. 99 newest)
+  for (int t = 0; t < OD_WIN; t++) {
+    int pos = (circ_write - OD_WIN + t + 2 * CIRC_SIZE) % CIRC_SIZE;
+    for (int c = 0; c < N_CH; c++) od_window[c * OD_WIN + t] = circ_buf[c][pos];
+  }
+
+  // label 254 = evaluate only: predict + report, but do NOT train or touch the
+  // running stats. Lets the host measure the frozen model on a held-out set so an
+  // adaptation drop/recovery is visible (the running accuracy alone hides it).
+  bool eval = (label == 254);
+
+  // Features -> standardise -> predict
+  float  feat[OD_N_FEAT];
+  double z[OD_N_FEAT];
+  od_extract_features(od_window, feat);
+  if (!eval) od_std.update(feat);
+  od_std.transform(feat, z);
+  double p   = od_model.predict_proba(z);
+  int    pred = (p >= 0.5) ? 1 : 0;
+
+  set_led(pred == 1 ? PHASE_PRICKING : PHASE_MOVEMENT);
+  result_hold_until = (uint32_t)millis() + 300;
+
+  if (eval) {
+    Serial.print("EVAL ts="); Serial.print(ts);
+    Serial.print(" pred=");   Serial.print(pred);
+    Serial.print(" p=");      Serial.println(p, 3);
+    Serial.flush();
+    return;
+  }
+
+  // Label: pricking(2)=1, dorsi/plantar(0,1)=0, unknown(255)=-1
+  int y = (label == 255) ? -1 : (label == 2 ? 1 : 0);
+
+  // Train (only if we have a label) and track running accuracy
+  if (y >= 0) {
+    seen++;
+    if (pred == y) correct++;
+    od_model.update(z, y);
+  }
+
+  Serial.print("CAP ts=");  Serial.print(ts);
+  Serial.print(" pred=");   Serial.print(pred);
+  Serial.print(" label=");  Serial.print(y);
+  Serial.print(" p=");      Serial.print(p, 3);
+  Serial.print(" acc=");    Serial.println(seen ? (double)correct / (double)seen : 0.0, 3);
+  Serial.flush();
+}
+
+// ── Serial protocol (same as stream_raw.py) ──────────────────────────────────
+#define BAUD    115200
+#define PAYLOAD 117   // label(1) + ts(4) + 56*uint16(112)
+
+static uint8_t pkt_buf[PAYLOAD];
+static int     pkt_pos   = 0;
+static int     pkt_state = 0;   // 0=wait 0xAA  1=wait 0xBB  2=payload
+
+void setup() {
+  pinMode(LEDR, OUTPUT);
+  pinMode(LEDG, OUTPUT);
+  pinMode(LEDB, OUTPUT);
+  set_led(PHASE_IDLE);
+
+  Serial.begin(BAUD);
+  while (!Serial) {}
+
+  memset(circ_buf, 0, sizeof(circ_buf));
+  memset(notch_s,  0, sizeof(notch_s));
+  memset(bp_s,     0, sizeof(bp_s));
+  memset(daa_hist, 0, sizeof(daa_hist));
+
+  od_std.init();
+  od_model.init();   // lr=0.05, l2=1e-4
+
+  Serial.println("ArduinoOnDeviceTrain ready.");
+}
+
+void loop() {
+  while (Serial.available()) {
+    uint8_t b = (uint8_t)Serial.read();
+    if (pkt_state == 0) {
+      if (b == 0xAA) pkt_state = 1;
+    } else if (pkt_state == 1) {
+      if      (b == 0xBB) { pkt_state = 2; pkt_pos = 0; }
+      else if (b == 0xAA) pkt_state = 1;
+      else                pkt_state = 0;
+    } else {
+      pkt_buf[pkt_pos++] = b;
+      if (pkt_pos == PAYLOAD) {
+        uint8_t  label = pkt_buf[0];
+        uint32_t ts;     memcpy(&ts, pkt_buf + 1, 4);
+        uint16_t raw[N_CH]; memcpy(raw, pkt_buf + 5, N_CH * 2);
+        process_sample(raw, ts, label);
+        pkt_state = 0;
+        pkt_pos   = 0;
+      }
+    }
+  }
+}

scripts/Arduino/ArduinoOnDeviceTrain/README.md

@@ -0,0 +1,34 @@
+# ArduinoOnDeviceTrain
+
+Detects CAPs live from the streamed nerve data and **trains** a small model on
+the board (pricking vs movement). It runs the full detection front-end and then
+does one SGD step per detected CAP, using the label that comes with the stream.
+
+## Files
+- `ArduinoOnDeviceTrain.ino` — the sketch (detection + on-device training)
+- `ondevice_model.h` / `ondevice_model.cpp` — the model (features, standardizer,
+  logistic regression + SGD). Same code the host test in `../ondevice_train/` checks.
+
+## Run it
+1. Open the sketch in the Arduino IDE and upload to the Nano 33 BLE Sense.
+2. Stream the data to it from the laptop:
+   ```
+   python scripts/cap_pipeline/ondevice/stream_raw.py
+   ```
+3. Watch the serial monitor. Per detected CAP it prints:
+   ```
+   CAP ts=.. pred=.. label=.. p=.. acc=..
+   ```
+   `pred`/`label`: 1=pricking, 0=movement. `acc` is the running accuracy, which
+   should climb as it learns.
+
+LED: cyan = idle, yellow = CAP detected, magenta = computing, blue = pricking,
+green = movement.
+
+## Notes
+- Stream slower than real time if the board can't keep up (30 kHz x 56 ch is a lot
+  over serial) — it's a replay, so the rate doesn't have to be realtime.
+- `stream_raw.py` drives the adaptation demo: it streams a sequence of `PHASES`
+  (e.g. RAT4 then RAT9), parses the `pred=`/`label=` reply, prints a sliding-window
+  accuracy so the shift and recovery are visible, and logs the learning curve to
+  `caches/learning_curve.csv`. Edit `PHASES`/`SERIAL_PORT` at the top of the file.