Timestamp compensation: Physio16 realignment, filter recalibration, system delay, mock fidelity

.gitignore

@@ -20,3 +20,13 @@ settings.local.json
 
 # macOS
 .DS_Store
+
+# misc
+.codegraph/
+delay_capture*/
+
+# Audio latency test outputs (large XDF recordings, plots, CLI logs)
+scripts/audio_latency_test/recordings/
+
+# Python cache
+__pycache__/

CMakeLists.txt

@@ -27,7 +27,7 @@ option(EGIAMP_BUILD_MOCK "Build the mock EGI Amp device for testing" ON)
 # By default, liblsl is fetched automatically from GitHub.
 # To use a pre-installed liblsl, set LSL_INSTALL_ROOT.
 set(LSL_INSTALL_ROOT "" CACHE PATH "Path to installed liblsl (optional)")
-set(LSL_FETCH_REF "dev" CACHE STRING "liblsl version to fetch from GitHub")
+set(LSL_FETCH_REF "v1.18.0.b2" CACHE STRING "liblsl version to fetch from GitHub")
 
 if(LSL_INSTALL_ROOT)
     # Use pre-installed liblsl

README.md

@@ -35,7 +35,6 @@ The CLI provides a lightweight alternative to the GUI:
 - `--amp-id <id>` - Amplifier ID (default: 0)
 - `--sample-rate <hz>` - Sample rate in Hz (default: 1000). Forces amplifier to this rate if already running at a different rate. Valid rates: 250, 500, 1000 (decimated) or 500, 1000, 2000, 4000, 8000 (native).
 - `--fast-recovery` - Use native rate mode (no FPGA anti-alias filter) for lower latency. See [Sample Rate Modes](#sample-rate-modes).
-- `--align-timestamps` - Adjust timestamps to compensate for anti-alias filter delay. See [Timestamp Alignment](#timestamp-alignment).
 - `--impedance` - Enable impedance testing mode
 - `--native-format` - Transmit raw int32 ADC counts instead of float microvolts
 - `--shutdown` - Shutdown the Amp Server (terminates all connections)
@@ -52,11 +51,11 @@ The CLI provides a lightweight alternative to the GUI:
 # With native format (int32 ADC counts)
 ./EGIAmpServerCLI --address 10.10.10.51 --native-format
 
-# Low-latency mode (fast recovery)
+# Low-latency mode (fast recovery / native)
 ./EGIAmpServerCLI --address 10.10.10.51 --sample-rate 1000 --fast-recovery
 
-# With timestamp alignment for filter delay
-./EGIAmpServerCLI --address 10.10.10.51 --sample-rate 1000 --align-timestamps
+# Decimated mode — timestamp compensation is applied automatically
+./EGIAmpServerCLI --address 10.10.10.51 --sample-rate 1000
 ```
 
 ## Connecting to an Already-Running Amplifier
@@ -78,9 +77,8 @@ This is safe to use alongside Net Station — it will not stop or reconfigure th
 
 The following flags cause the CLI to stop, reconfigure, and restart the amplifier — even if it was started by another application:
 
-- `--sample-rate <hz>` — reinitializes if the detected rate differs from the requested rate
+- `--sample-rate <hz>` — reinitializes if the detected rate differs from the requested rate. At 500/1000 Hz, where native and decimated are indistinguishable in the data stream, requesting a decimated rate also forces reinitialization to guarantee decimated mode (so the automatic [timestamp compensation](#timestamp-compensation) is correct)
 - `--fast-recovery` — reinitializes to ensure native (unfiltered) mode
-- `--align-timestamps` — reinitializes at 500/1000 Hz to ensure decimated (filtered) mode, since the operating mode cannot be distinguished from the data stream alone
 
 **Warning**: Reinitialization will interrupt any active Net Station recording. If you need to coexist with Net Station, omit these flags and let the CLI match the existing configuration.
 
@@ -100,14 +98,14 @@ The NA400/NA410 amplifiers support two operating modes that affect anti-aliasing
 
 Uses the FPGA's digital anti-aliasing filter to downsample from the ADC's native rate. This provides:
 - Better frequency response (~400 Hz bandwidth at 1000 Hz sample rate)
-- Higher latency due to filter delay (36-112 samples depending on rate)
+- Higher latency due to the filter group delay (36-111 ms depending on rate), which the app compensates for automatically — see [Timestamp Compensation](#timestamp-compensation)
 
 Available decimated rates: 250, 500, 1000 Hz
 
 ### Native Mode (Fast Recovery)
 
 Bypasses the FPGA filter and samples directly at the requested rate. This provides:
-- Lower latency (~3 samples)
+- Lower latency (no filter group delay)
 - Reduced bandwidth (~1/4 of sample rate, e.g., 250 Hz at 1000 Hz sample rate)
 - Optimized for EEG-TMS and real-time BCI applications
 
@@ -117,32 +115,52 @@ Use `--fast-recovery` to enable native mode for rates that support both modes (5
 
 ### Filter Delay by Sample Rate
 
-| Mode | Sample Rate | Filter Delay (samples) | Filter Delay (ms) |
-|------|-------------|------------------------|-------------------|
-| Decimated | 250 Hz | 112 | 448 ms |
-| Decimated | 500 Hz | 66 | 132 ms |
-| Decimated | 1000 Hz | 36 | 36 ms |
-| Native | 500-8000 Hz | ~3 | ~3 ms |
+DIN→EEG group delay of the FPGA anti-alias filter (decimated mode only), measured with `scripts/delay_capture_sweep.py`:
 
-## Timestamp Alignment
+| Mode | Sample Rate | Filter delay (ms) | (samples) |
+|------|-------------|-------------------|-----------|
+| Decimated | 250 Hz | 111 | ~28 |
+| Decimated | 500 Hz | 61 | ~30 |
+| Decimated | 1000 Hz | 36 | 36 |
+| Native | 500-8000 Hz | 0 | 0 |
 
-When using decimated mode, the FPGA anti-aliasing filter introduces a delay between when brain activity occurs and when it appears in the data stream. The `--align-timestamps` option compensates for this by adjusting LSL timestamps backward by the filter delay amount.
+## Timestamp Compensation
 
-### When to Use
+The application adjusts the LSL timestamp it assigns to each sample so that the timestamp reflects **when the signal actually occurred**, not when the bytes happened to arrive at this client. This makes EEG, Physio16, and DIN events line up with each other and with external event markers. Three corrections are applied; together they are what we call timestamp compensation.
 
-- **ERP analysis**: Enable `--align-timestamps` to align EEG data with event markers
-- **Real-time BCI**: Use `--fast-recovery` instead (no filter delay to compensate)
-- **Raw recording**: Disable alignment if you prefer unmodified timestamps
+There is **no flag to toggle this** — it follows the configured mode automatically. (The old `--align-timestamps` flag was removed.) The app trusts the mode you configure: a decimated rate means the filter corrections apply; a native rate (`--fast-recovery`, or any rate above 1000 Hz) means they do not. At 500/1000 Hz, where native and decimated are indistinguishable in the data stream, the app assumes the mode you asked for (and, when forcing a rate, reinitializes to guarantee it).
 
-### Limitations
+### 1. System (pipeline) delay — always applied
 
-**Important**: Timestamp alignment only works correctly when this application initializes the amplifier. If Net Station or another application previously initialized the amplifier, the current operating mode (decimated vs native) cannot be queried from AmpServer. In this case:
+Every pushed timestamp (EEG, Physio16, **and** DIN) is moved earlier by `systemdelayms` to account for the time between digitization and the sample becoming available here: device firmware + network transmission + our read path. Measured at ~5 ms in native mode using `scripts/audio_latency_test`.
 
-1. The application will reinitialize the amplifier to ensure the correct mode
-2. This will interrupt any existing Net Station recording
-3. To avoid this, start EGIAmpServer before Net Station, or restart the amplifier
+The default is **4.5 ms** — deliberately a little *under* the measured value. Under-compensating is the safe direction: it guarantees a sample is never back-dated to *before* the event that produced it (a response must not precede its stimulus). It is user-configurable; raise it toward the measured latency if you prefer tighter alignment and can accept that risk.
 
-If you need to join an existing Net Station session without reinitialization, do not use `--align-timestamps` unless you are certain of the current mode.
+### 2. FPGA filter offset — decimated mode only
+
+In decimated mode the anti-alias filter delays the EEG relative to the (unfiltered) DIN by the group delay in the table above (111/61/36 ms at 250/500/1000 Hz). The app subtracts this from the EEG/Physio timestamps so filtered EEG lines up with DIN. Native mode has no filter, so nothing is subtracted. DIN itself is never filter-shifted (only the system delay applies to it).
+
+### 3. Physio16 realignment — decimated mode only
+
+The Physio16 (PNS/PIB) acquisition path runs **ahead** of the FPGA-filtered EEG by a fixed ~33 ms in decimated mode, so the same event lands at different sample indices in the EEG vs physio channels of the combined stream. To keep them sample-aligned, the physio channels are buffered (delayed) by `physioaligndelayms` (default **33 ms**, applied as `floor(33 ms × rate)` samples). The physio channels therefore read **zeros for the first ~33 ms** after streaming starts, then real data. Native mode applies no physio delay.
+
+### Configuration
+
+Both values live in the config file (and persist via File → Save Configuration):
+
+```xml
+<settings>
+  ...
+  <physioaligndelayms>33</physioaligndelayms>   <!-- Physio16 realignment, decimated only -->
+  <systemdelayms>4.5</systemdelayms>            <!-- pipeline-latency compensation, always -->
+</settings>
+```
+
+The filter-offset values (111/61/36 ms) are built in. All three were measured on real hardware; see `scripts/audio_latency_test/README.md` (system delay) and `notebooks/delay_inspection.ipynb` (filter and physio delays).
+
+### Coexisting with Net Station
+
+The filter corrections are only correct when the operating mode is known. If another application (e.g. Net Station) already started the amplifier at an ambiguous rate (500/1000 Hz), the app applies the corrections for the mode **you configured**, assuming it matches. If you are unsure of the running mode, attach with a native configuration (no filter corrections) or let the app reinitialize to a known mode (which interrupts the existing session).
 
 ## Impedance Testing
 

ampserver_config.cfg

@@ -9,4 +9,16 @@
   <amplifierid>0</amplifierid>
   <samplingrate>1000</samplingrate>
   <channelcount>280</channelcount>
+  <!-- Physio16 (PNS/PIB) delay vs FPGA-filtered EEG in decimated mode, in ms.
+       Fixed 33 ms physical skew (300 s DIN+EEG+Physio16 sweep, corroborated across
+       device versions; scripts/delay_capture_sweep.py + analyze_delay_sweep.py).
+       Applied as floor(33 ms x rate) samples -> 32 ms @ 250/500 Hz, 33 ms @ 1000 Hz.
+       Decimated rates only; native rates have ~0 skew. -->
+  <physioaligndelayms>33</physioaligndelayms>
+  <!-- System/pipeline latency (ms) subtracted from every pushed timestamp
+       (EEG, physio, DIN) to compensate digitization->client delay (firmware +
+       network + read path). Measured ~5 ms via scripts/audio_latency_test;
+       default 4.5 deliberately under-compensates so a sample is never back-dated
+       before its event. User-tunable. -->
+  <systemdelayms>4.5</systemdelayms>
 </settings>

mock/include/MockAmplifier.h

@@ -28,6 +28,14 @@ struct AmplifierState {
     // Sampling configuration
     int nativeRate = 1000;
     int decimatedRate = 1000;
+    bool decimated = true;  // true after cmd_SetDecimatedRate (FPGA filter on),
+                            // false after cmd_SetNativeRate (no filter, no skew)
+
+    // Active streaming rate = whichever mode was last selected.
+    int activeRate() const {
+        int r = decimated ? decimatedRate : nativeRate;
+        return r > 0 ? r : 1000;
+    }
 
     // Calibration/test signal state
     bool all10KOhms = false;
@@ -165,11 +173,6 @@ class MockAmplifier {
     // For synthetic data generation
     double phase_ = 0.0;
 
-    // DIN counter - cycles 0x0000 to 0xFFFF at 1kHz base rate
-    // At higher sample rates, values are duplicated:
-    // 2kHz = 2x duplicates, 4kHz = 4x, 8kHz = 8x
-    uint32_t dinCounter_ = 0;
-
     int64_t getCurrentTimeMicros() const;
 };
 

mock/src/DataStreamGenerator.cpp

@@ -119,16 +119,31 @@ void DataStreamGenerator::stopListening(int64_t ampId) {
 }
 
 void DataStreamGenerator::streamingThread() {
-    // We send packets at ~200 Hz (every 5ms) and adjust samples per packet
-    // based on sample rate to achieve the correct rate
-    const auto packetInterval = std::chrono::milliseconds(5);
+    // Pace to a steady deadline so the effective sample rate is accurate. Each
+    // packet carries `samplesPerPacket` samples (~5 ms worth) and must therefore
+    // span exactly samplesPerPacket/rate seconds. Computing the interval that way
+    // (rather than a fixed 5 ms) keeps non-1000 rates exact too: e.g. 250 Hz ->
+    // 1 sample / 4 ms, 500 Hz -> 2 samples / 4 ms. sleep_until compensates for the
+    // time spent generating/sending, which sleep_for did not.
+    auto nextDeadline = std::chrono::steady_clock::now();
 
     while (running_) {
         if (amplifier_->isStreaming() && listening_) {
             sendDataToClients();
         }
 
-        std::this_thread::sleep_for(packetInterval);
+        const int rate = amplifier_->state().activeRate();
+        int samplesPerPacket = (rate * 5) / 1000;
+        if (samplesPerPacket < 1) samplesPerPacket = 1;
+        const auto interval = std::chrono::nanoseconds(
+            static_cast<int64_t>(1'000'000'000LL * samplesPerPacket / rate));
+
+        nextDeadline += interval;
+        const auto now = std::chrono::steady_clock::now();
+        if (nextDeadline < now) {
+            nextDeadline = now;  // fell behind (e.g. stall) — resync, don't burst
+        }
+        std::this_thread::sleep_until(nextDeadline);
     }
 }
 
@@ -158,7 +173,7 @@ void DataStreamGenerator::generateSyntheticData(std::vector<uint8_t>& buffer) {
     // Calculate samples per packet based on sample rate
     // We send packets every 5ms, so samples = rate * 0.005
     // 1000 Hz = 5 samples, 2000 Hz = 10 samples, 4000 Hz = 20 samples, 8000 Hz = 40 samples
-    int sampleRate = state.decimatedRate > 0 ? state.decimatedRate : 1000;
+    int sampleRate = state.activeRate();
     int samplesPerPacket = (sampleRate * 5) / 1000;
     if (samplesPerPacket < 1) samplesPerPacket = 1;
 

mock/src/MockAmplifier.cpp

@@ -4,6 +4,40 @@
 
 namespace mock {
 
+namespace {
+// DIN->EEG anti-alias filter group delay (ms) in decimated mode. Mirrors
+// getFilterDelaySeconds() in the app (measured values). Native mode = no filter.
+int filterDelayMs(int rate) {
+    switch (rate) {
+        case 250:  return 111;
+        case 500:  return 61;
+        case 1000: return 36;
+        default:   return 0;
+    }
+}
+
+// The PIB/physio path leads the FPGA-filtered EEG by this fixed physical delay in
+// decimated mode (see notebooks/delay_inspection.ipynb). Quantized to whole
+// samples at the active rate, matching the device (and the app's PhysioDelayLine).
+constexpr int PHYSIO_LEAD_MS = 33;
+
+// Common reference event injected into DIN pin-1 / EEG ch1 / physio ch1 (emulates
+// the test jig). One pulse per second, 250 ms wide.
+constexpr double REF_PERIOD_S = 1.0;
+constexpr double REF_PULSE_S = 0.25;
+constexpr double REF_EEG_PULSE_UV = 500.0;     // EEG ch1 deflection during the pulse
+constexpr double REF_PHYSIO_PULSE_UV = 1000.0; // physio ch1 deflection during the pulse
+
+// Is the reference pulse active at sample index k (handles negative k from delays)?
+bool refActiveAtSample(long long k, int sampleRate) {
+    const long long period = static_cast<long long>(REF_PERIOD_S * sampleRate);
+    const long long width = static_cast<long long>(REF_PULSE_S * sampleRate);
+    if (period <= 0) return false;
+    const long long m = ((k % period) + period) % period;
+    return m < width;
+}
+}  // namespace
+
 MockAmplifier::MockAmplifier(int64_t id) {
     state_.ampId = id;
     state_.startTime = getCurrentTimeMicros();
@@ -92,6 +126,7 @@ bool MockAmplifier::setNativeRate(int rate) {
         }
     }
     state_.nativeRate = rate;
+    state_.decimated = false;  // native mode: no FPGA anti-alias filter, no skew
     return true;
 }
 
@@ -100,6 +135,7 @@ bool MockAmplifier::setDecimatedRate(int rate) {
     // Real amplifier accepts any rate value without validation
     // Valid rates: 250, 500, 1000
     state_.decimatedRate = rate;
+    state_.decimated = true;  // decimated mode: FPGA filter introduces EEG/physio skew
     return true;
 }
 
@@ -369,16 +405,31 @@ void MockAmplifier::generatePacketFormat2(PacketFormat2_SamplePacket& packet) {
     std::lock_guard<std::mutex> lock(mutex_);
     std::memset(&packet, 0, sizeof(packet));
 
-    // Digital inputs - cycle 0x0000 to 0xFFFF at 1kHz base rate
-    // At higher sample rates, duplicate the values:
-    // 1kHz: each value once, 2kHz: twice, 4kHz: 4x, 8kHz: 8x
-    int sampleRate = state_.decimatedRate > 0 ? state_.decimatedRate : 1000;
-    int duplicateFactor = sampleRate / 1000;
-    if (duplicateFactor < 1) duplicateFactor = 1;
-
-    // Calculate which 1kHz "tick" we're on
-    uint32_t din1kHzTick = static_cast<uint32_t>(state_.packetCounter / duplicateFactor);
-    packet.digitalInputs = static_cast<uint16_t>(din1kHzTick & 0xFFFF);
+    // Active streaming rate (decimated or native mode).
+    int sampleRate = state_.activeRate();
+
+    // Common reference event (emulates the test jig pulsing DIN pin-1 while
+    // injecting the same pulse into EEG ch1 and physio ch1). DIN reflects the
+    // event immediately; EEG and physio see it delayed in decimated mode so the
+    // device's EEG-vs-DIN filter delay and physio-leads-EEG-by-33ms skews appear.
+    const long long n = static_cast<long long>(state_.packetCounter);
+    // EEG filter delay rounds to the nearest sample (the app applies it as a
+    // continuous timestamp offset, so nearest-sample is the closest the mock's
+    // sample grid can get). The physio lead floors, matching the device's own
+    // whole-sample quantization (32 ms @ 250/500, 33 ms @ 1000).
+    const int eegDelaySamples = state_.decimated
+        ? (filterDelayMs(sampleRate) * sampleRate + 500) / 1000 : 0;
+    const int physioLeadSamples = state_.decimated
+        ? PHYSIO_LEAD_MS * sampleRate / 1000 : 0;
+    int physioDelaySamples = eegDelaySamples - physioLeadSamples;
+    if (physioDelaySamples < 0) physioDelaySamples = 0;
+
+    // DIN is active-low: 0xFFFF idle, pin-1 (bit 0) cleared while the event is on.
+    uint16_t din = 0xFFFF;
+    if (refActiveAtSample(n, sampleRate)) {
+        din &= static_cast<uint16_t>(~0x0001);
+    }
+    packet.digitalInputs = din;
 
     // TR byte (255 = no GTEN activity)
     packet.tr = state_.gtenTrainRunning ? 251 : 255;
@@ -455,6 +506,12 @@ void MockAmplifier::generatePacketFormat2(PacketFormat2_SamplePacket& packet) {
         packet.eegData[ch] = static_cast<int32_t>(value / scaleFactor);
     }
 
+    // Inject the (filter-delayed) reference pulse into EEG ch1 so it lags the DIN
+    // pin-1 edge by the anti-alias filter delay (skip in impedance mode).
+    if (!impedanceMode && refActiveAtSample(n - eegDelaySamples, sampleRate)) {
+        packet.eegData[0] += static_cast<int32_t>(REF_EEG_PULSE_UV / scaleFactor);
+    }
+
     // Aux data
     for (int i = 0; i < 3; ++i) {
         packet.auxData[i] = 0;
@@ -494,6 +551,14 @@ void MockAmplifier::generatePacketFormat2(PacketFormat2_SamplePacket& packet) {
         }
     }
 
+    // Inject the reference pulse into physio ch1 (port 1). It uses a smaller delay
+    // than EEG (leads it by PHYSIO_LEAD_MS), so the same event lands at an earlier
+    // sample index than in EEG — exactly the skew PhysioDelayLine corrects.
+    if (!impedanceMode && (state_.physioConnectionStatus & 0x01) &&
+        refActiveAtSample(n - physioDelaySamples, sampleRate)) {
+        packet.pib1_Data[0] += static_cast<int32_t>(REF_PHYSIO_PULSE_UV / PHYSIO_SCALE_1_8);
+    }
+
     ++state_.packetCounter;
 }