Aligning Concurrent Data Streams for Neuroscience

Many experiments now record from several systems at once: EEG or MEG, an eye tracker, a response box, a display.

Our guide Synchronization Strategies for Multimodal Neuroscience showed how to send event markers between systems to share event timing. Markers tell you when an event happened, but they cannot keep continuous recordings aligned between events.

This guide addresses continuous data stream alignment. The premise is simple; a single repeating synchronization pulse is sent to every recorder, including the system generating it. After the session, this pulse in each data stream allows you to align all signals to a single timeline.

The advantage of multiple aligned streams

VPixx hardware is often used to forward signals to other systems. Button presses, logged eye tracking, and triggers are passed to a receiver, be it an EEG amplifier, an MEG system, or some other third-party recording equipment.

This method collapses multiple data acquisition streams into a single master data stream, which is a reasonable approach. It also introduces a limitation: the number of signals the master data stream can handle is set by its own input channels, and signal resolution is set by the data stream's sampling rate. Even digital signals, which execute near instantaneously, get throttled to the sampling rate of the receiver (e.g., 600 Hz).

As of summer 2026, we have released a new utility called Digital Scope, which is an extension of our digital input log. Digital Scope allows researchers to use the DATAPixx3 as a recording stream that complements other acquisition systems rather than simply feeding them.

Records	Digital Input log	NEW Digital Scope log
Changes in first 16 digital inputs
Changes in last 8 digital inputs
Changes in all 24 digital outputs

With this method, you can leverage multiple recording systems working in parallel, as opposed to one single master data stream, without sacrificing temporal precision.

All hardware clocks differ slightly in their timing, due to differences in their crystals. Differences are typically expressed in ppm, or parts per million. Relative timing differences between clocks in the range of +/- 50 ppm is not uncommon.

In practice, that means that data streams on different clocks drift with respect to each other over time. Consider two events logged every two seconds on clock 1, which has a +50 ppm difference compared to clock 2:

Repeating synchronization pulses often allows for piecewise remapping of signals to frequent events. This eliminates drift that can add up when you use one-off event-alignment triggers that are spaced far apart.

While short trials and trigger-locked events like ERPs may not see a major difference compared to single alignment markers, longer blocks and time-series analyses like frequency tagging coherence analysis benefit a great deal from the correction.

General steps

Here is the general step-by-step process for generating a cross-system synchronization pulse:

Generate the pulse on the DATAPixx3 using the DOut schedule. A digital output, or Dout, schedule drives one output channel high for 30 ms, once per second. A triple beat marks the block onset, and a double beat logs the block end, automatically creating block epochs.
The exact pulse times and onset/offset signatures are entirely up to you. We recommend a pulse once per second (1 Hz) or less.
Connect that digital output to every external recorder. Connect a cable from the DATAPixx3 digital output channel to the receiver’s own TTL input.
Log events on the DATAPixx3. Use Digital Scope to record all 24 digital inputs and 24 digital outputs on the DATAPixx3 clock at a 1,000 kHz sampling rate. This data stream captures the synchronization pulse. It can also capture RESPONSEPixx button presses and other digital inputs like MRI triggers, and can record Pixel Mode video triggers and manual event triggers previously treated as output-only.

This same hardware records TRACKPixx3 eye tracking data, analog inputs, and microphone/audio inputs. Users can nominally synchronize analog and audio playback to video events using pixel sync, and save a digital output or Pixel Mode trigger on the same video frame to mark the playback time precisely
Fit the clocks after the session. Pair each pulse in the DATAPixx3 log with the same pulse on each recorder, fit the offset and drift, and map any event time from one clock to another.

Rinse and repeat. Add additional data streams via a split synchronization cable, and align all streams from the DATAPixx3 and other acquisition systems together in post-processing. Below is an expanded example that adds TRACKPixx3 gaze and blink flags on the DATAPixx3 stream, and acquisition on a second OPM-MEG system.

The animation below shows the stream alignment and how the TRACKPixx3 blink flags can be used to mask closed-eye instances in the neural trace.

The practical gain in signal strength

The above alignment method was used to align a rapid frequency-tagging reference signal timestamped on the DATAPixx3, with a neural recording on a dry-EEG headset (see our ECVP poster; paper in prep).

The two clocks differed by a highly consistent 49.3 ppm. Coherence analysis on the drift-corrected signal timing resulted in an 11-fold increase in statistical significance on a 2-minute trial.

There is a very real advantage to aligning streams as precisely as possible for these studies, and others where time-series analyses cross multiple recording systems.

Example code

Some DATAPixx3 units purchased before 2023 will need new firmware that has not been released yet. If you get an error trying to enable Digital Scope on your unit, please contact support@vpixx.com.

The VIEWPixx3 /OLED also supports Digital Scope but not digital output schedules; pulse trains must be driven manually or by Pixel Mode.

Below are examples using our APIs for Python (pypixxlib) and MATLAB/Psychtoolbox (Datapixx.mex), showing how to start the log, generate the start and stop markers, and play a continuous pulse train during a block.

MATLAB

% Requires Psychtoolbox, DATAPixx3 with firmware revision 29 or later, and VPixx software tools revision 1.12.0 or later.
SYNC_BIT   = 0;        % DOUT bit wired to the BNC (here DOUT 0)
LOG_ADDR   = 12e6;     % DATAPixx RAM address for the Digital Scope log
LOG_FRAMES = 250000;   % log capacity, in logged changes

cfg.SYNC_BIT   = SYNC_BIT;
cfg.RATE_HZ    = 1000;     % schedule sample rate: 1 ms resolution
cfg.PULSE_MS   = 30;       % pulse width
cfg.SCHED_ADDR = 8e6;      % DATAPixx RAM address for the pulse waveform

startBeat = pulseWave([0 0.1 0.2], 0.5, cfg);   % 3 pulses, 100 ms apart
train     = pulseWave(0, 1.0, cfg);             % one pulse, repeated every second
stopBeat  = pulseWave([0 0.3], 0.5, cfg);       % 2 pulses, 300 ms apart
Datapixx('Open');
Datapixx('StopAllSchedules');
Datapixx('SetDoutValues', 0, bitshift(1, SYNC_BIT));   % start low
Datapixx('RegWrRd');
Datapixx('SetDinLog', LOG_ADDR, LOG_FRAMES, 1);        % 1 = Full: activates the Digital Scope feature
Datapixx('StartDinLog');
Datapixx('RegWrRd');

playWave(startBeat, false, cfg);                       % 3 pulses 100 ms apart
WaitSecs(0.5 + 1.0);                                   % the beat, then 1 s of quiet
playWave(train, true, cfg);                            % 30 ms pulse every second, loops until stopped

WaitSecs(20);                                          % your trials go here

Datapixx('StopDoutSchedule');
Datapixx('SetDoutValues', 0, bitshift(1, SYNC_BIT));   % leave the bit low
Datapixx('RegWrRd');
WaitSecs(1);
playWave(stopBeat, false, cfg);                        % 2 pulses 300 ms apart
WaitSecs(0.5 + 1.0);

% Save the Digital Scope log: time, DIN word, DOUT word
Datapixx('StopDinLog');
Datapixx('RegWrRd');
st = Datapixx('GetDinStatus');
[logData, logTime] = Datapixx('ReadDinLog', st.newLogFrames);   % row 1 DIN, row 2 DOUT
logTable = array2table([logTime(:), double(logData(1, :))', double(logData(2, :))'], ...
                       'VariableNames', {'t_datapixx_s', 'din24', 'dout24'});
writetable(logTable, 'datapixx_log.csv');
Datapixx('Close');

function w = pulseWave(risingEdgesS, lengthS, cfg)
    w = zeros(1, round(lengthS*cfg.RATE_HZ));
    width = round(cfg.PULSE_MS*cfg.RATE_HZ/1000);
    for t = risingEdgesS
        i = round(t*cfg.RATE_HZ) + 1;
        w(i:min(end, i+width-1)) = bitshift(1, cfg.SYNC_BIT);
    end
end

function playWave(w, loop, cfg)
    Datapixx('StopDoutSchedule');
    Datapixx('RegWrRd');                           % stop is committed before the buffer is rewritten
    Datapixx('WriteDoutBuffer', w, cfg.SCHED_ADDR);
    if loop                                        % 0 frames = loop until stopped
        maxFrames = 0;
    else
        maxFrames = numel(w);
    end
    Datapixx('SetDoutSchedule', 0, cfg.RATE_HZ, maxFrames, cfg.SCHED_ADDR, numel(w));
    Datapixx('StartDoutSchedule');
    Datapixx('RegWrRd');
end

Python

import csv
import time

from pypixxlib import _libdpx as dp

SYNC_BIT = 0                                    # DOUT bit wired to the BNC (here DOUT 0)
RATE_HZ = 1000                                  # schedule sample rate: 1 ms resolution
PULSE_MS = 30                                   # pulse width
SCHED_ADDR = 8_000_000                          # DATAPixx RAM address for the pulse waveform
LOG_ADDR = 12_000_000                           # DATAPixx RAM address for the Digital Scope log
LOG_FRAMES = 250_000                            # log capacity, in logged changes


def pulse_wave(rising_edges_s, length_s):
    """One DOUT value per millisecond, SYNC_BIT high for PULSE_MS at each rising edge."""
    wave = [0] * round(length_s * RATE_HZ)
    width = round(PULSE_MS * RATE_HZ / 1000)
    for t in rising_edges_s:
        start = round(t * RATE_HZ)
        for i in range(start, start + width):
            wave[i] = 1 << SYNC_BIT
    return wave


START_BEAT = pulse_wave([0.0, 0.1, 0.2], 0.5)  # 3 pulses, 100 ms apart
TRAIN = pulse_wave([0.0], 1.0)                 # one pulse, repeated every second
STOP_BEAT = pulse_wave([0.0, 0.3], 0.5)        # 2 pulses, 300 ms apart

dp.DPxOpen()
dp.DPxStopAllScheds()
dp.DPxSetDoutValue(0, 1 << SYNC_BIT)            # start low
dp.DPxUpdateRegCache()
dp.DPxEnableDinDio48()                          # activates the Digital Scope feature: log DOUT alongside DIN
log = dp.DPxSetDinLog(buffBaseAddr=LOG_ADDR, numBuffFrames=LOG_FRAMES, mode="Full")
dp.DPxStartDinLog()
dp.DPxUpdateRegCache()


def play(wave, loop=False):
    dp.DPxStopDoutSched()
    dp.DPxUpdateRegCache()
    dp.DPxWriteDoutBuffer(wave, SCHED_ADDR)
    frames = 0 if loop else len(wave)           # 0 = loop until stopped
    dp.DPxSetDoutSchedule(0.0, RATE_HZ, frames, SCHED_ADDR, len(wave))
    dp.DPxStartDoutSched()
    dp.DPxUpdateRegCache()


play(START_BEAT)                                # 3 pulses 100 ms apart
time.sleep(0.5 + 1.0)                           # the beat, then 1 s of quiet
play(TRAIN, loop=True)                          # 30 ms pulse every second

time.sleep(20)                                  # your trials go here

dp.DPxStopDoutSched()
dp.DPxSetDoutValue(0, 1 << SYNC_BIT)            # leave the bit low
dp.DPxUpdateRegCache()
time.sleep(1.0)
play(STOP_BEAT)                                 # 2 pulses 300 ms apart
time.sleep(0.5 + 1.0)

# Save the Digital Scope log: time, DIN word, DOUT word
dp.DPxStopDinLog()
dp.DPxUpdateRegCache()
dp.DPxGetDinStatus(log)
rows = dp.DPxReadDinLog(log, nReadFrames=0)
with open("datapixx_log.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["t_datapixx_s", "din24", "dout24"])
    for t, din, dout in rows:
        writer.writerow([t, int(din) & 0xFFFFFF, int(dout) & 0xFFFFFF])
dp.DPxClose()

Aligning the data in post-processing

Below is a simple example of how to align data in post-processing. You need two files: the Digital Scope log and your EEG trigger channel, exported as one value per sample.

Toy versions of both data files can be downloaded [here] as .csv, so you can run the code and see results before recording anything. Download the file and unzip it into the same folder as your scripts.
The alignment method follows six discrete steps:

Find the pulses. In each recording, find the moment each synchronization pulse switches on. On the DATAPixx3 these times come from the Digital Scope log. On the EEG they come from the trigger channel.
Find the beats. Pulses less than half a second apart belong to one beat. Three pulses 100 ms apart mark the start of a block, and two pulses 300 ms apart mark the end.
Line up the start. Matching the first start beat on the EEG with a start beat in the DATAPixx3 log gives a first estimate of the offset between the two clocks.
Match the pulses outward. From that starting point, predict where each DATAPixx3 pulse should appear on the EEG and take the EEG pulse within 10 ms of that spot. The code refits a straight line, offset plus drift, each time the window widens, until the whole recording is covered.
Check the result. Rather than return a wrong answer, the code stops with an error if the EEG matches two blocks equally well, if fewer than 10 pulses match, or if the drift is implausibly large. Otherwise, it reports how many pulses matched, the drift in parts per million and the pulse-fit residual.
Move events between clocks. To place a button press on the EEG clock, the code interpolates between the matched pulses on either side. Times outside the recording, or inside a gap of more than 2.5 s with no pulses (for example, between blocks), return NaN rather than a guess.

With the toy data, the EEG clock reads 7.3 s ahead of the DATAPixx3 clock, runs 40 ppm fast, and missed one pulse. The output should read:

matched 159 of 160 DATAPixx pulses in 2 part(s); drift +39.73 ppm; pulse-fit residual rms 0.289 ms, max 0.518 ms
button at   12.250 s on the DATAPixx3 =   19.551 s on the EEG
button at   30.875 s on the DATAPixx3 =   38.177 s on the EEG
button at   47.500 s on the DATAPixx3 =   54.802 s on the EEG
button at   95.125 s on the DATAPixx3 =  102.429 s on the EEG

MATLAB

% Align an EEG recording to the DATAPixx3 clock using one synchronization
% pulse, then move button presses onto the EEG clock. Needs MATLAB R2019a
% or later, no toolboxes.

% ---- Your settings ------------------------------------------------------
SYNC_BIT        = 0;      % DOUT bit that carried the pulse (same as the setup)
EEG_FS          = 1000;   % EEG sampling rate, Hz
EEG_TRIGGER_BIT = 0;      % bit of the EEG trigger channel the BNC is wired to

logTable   = readmatrix('datapixx_log.csv');   % written by the setup example
logT       = logTable(:, 1);
logDout    = logTable(:, 3);
eegTrigger = readmatrix('eeg_trigger.csv');    % EEG trigger channel exported from your EEG software, one value per sample

% Step 1: the pulse times on each clock
dpx  = risingEdges(logT, logDout, SYNC_BIT);
eegT = (0:numel(eegTrigger)-1)' / EEG_FS;
eeg  = risingEdges(eegT, eegTrigger, EEG_TRIGGER_BIT);

% Steps 2 to 5: match the pulses and fit the EEG clock to the DATAPixx3 clock
fit = fitClock(dpx, eeg);

% Step 6: move button presses from the DATAPixx3 clock to the EEG clock
buttonTDpx = [12.250; 30.875; 47.500; 95.125];   % fake presses, DATAPixx3 clock (s)
buttonTEeg = mapToRecorder(fit, buttonTDpx);
for k = 1:numel(buttonTDpx)
    fprintf('button at %8.3f s on the DATAPixx3 = %8.3f s on the EEG\n', buttonTDpx(k), buttonTEeg(k));
end


% ---- Step 1: find the pulses --------------------------------------------
function t = risingEdges(tSamples, word, bit)
% Times where one bit of word switches from 0 to 1.
    level = bitand(bitshift(uint32(word(:)), -bit), 1);
    rise = find(diff(double(level)) == 1) + 1;
    t = tSamples(rise);
    t = t(:);
end


% ---- Step 2: find the beats ---------------------------------------------
function beats = findBeats(edges)
% Groups of close pulses: 'start' (3 pulses, 100 ms), 'stop' (2 pulses,
% 300 ms), 'start?' (a start beat that lost one pulse) or 'unknown'.
    edges = sort(edges(:));
    beats = struct('kind', {}, 't', {});
    if isempty(edges)
        return;
    end
    breaks = [0; find(diff(edges) > 0.5); numel(edges)];
    tol = 0.04;
    for k = 1:numel(breaks)-1
        g = edges(breaks(k)+1:breaks(k+1));
        if numel(g) < 2
            continue;
        end
        gaps = diff(g);
        if numel(g) == 3 && all(abs(gaps - 0.1) < tol)
            kind = 'start';
        elseif numel(g) == 2 && abs(gaps - 0.3) < tol
            kind = 'stop';
        elseif numel(g) == 2 && (abs(gaps - 0.1) < tol || abs(gaps - 0.2) < tol)
            kind = 'start?';
        else
            kind = 'unknown';
        end
        beats(end+1) = struct('kind', kind, 't', g(1)); %#ok<AGROW>
    end
end


% ---- Steps 3 to 5: match the pulses and fit the clocks ------------------
function fit = fitClock(dpxEdges, recEdges, offset)
% Match every DATAPixx pulse with the recorder pulse it produced, then fit
% recorder = slope * datapixx + intercept. Pass offset (recorder minus
% DATAPixx time, s) only if the start beats are ambiguous.
    pairTol   = 0.010;   % a recorder pulse must sit within 10 ms of where it is expected
    maxGap    = 2.5;     % no mapping across gaps with no pulses longer than this
    maxAbsPpm = 1000;    % larger drift means the pairing went wrong
    minPairs  = 10;      % fewer matched pulses than this is refused
    dpx = checkEdges(dpxEdges, 'DATAPixx pulses');
    rec = checkEdges(recEdges, 'recorder pulses');
    if numel(dpx) < 2 || numel(rec) < 2
        error('fitClock:tooFew', 'Need at least two pulses on each side.');
    end
    if nargin < 3
        [offsets, anchors] = candidates(dpx, rec);
    else
        offsets = offset;
        [~, k] = min(abs(dpx + offset - rec(1)));
        anchors = dpx(k);
    end
    res = cell(numel(offsets), 1);
    n = zeros(numel(offsets), 1);
    for c = 1:numel(offsets)
        [i, j, p] = growFit(dpx, rec, offsets(c), anchors(c), pairTol);
        res{c} = {i, j, p};
        n(c) = numel(i);
    end
    [~, order] = sort(-n);
    best = res{order(1)};
    for c = order(2:end)'
        if n(c) >= n(order(1)) && abs(res{c}{3}(2) - best{3}(2)) > pairTol
            error('fitClock:ambiguous', ['Ambiguous: the recording fits the blocks starting at ' ...
                  '%.3f s and %.3f s equally well. Pass an offset, or give blocks different ' ...
                  'lengths.'], anchors(order(1)), anchors(c));
        end
    end
    i = best{1};
    j = best{2};
    p = best{3};
    if numel(i) < max(2, minPairs)
        error('fitClock:tooFewPairs', ['Only %d pulses matched (need %d); check the bit, the ' ...
              'cable and that the two recordings overlap.'], numel(i), max(2, minPairs));
    end
    driftPpm = (p(1) - 1) * 1e6;
    if abs(driftPpm) > maxAbsPpm
        error('fitClock:implausible', 'Fitted drift %+.0f ppm is implausible; the matching is wrong.', ...
              driftPpm);
    end
    fit.slope     = p(1);
    fit.intercept = p(2);
    fit.dpx       = dpx(i);
    fit.rec       = rec(j);
    fit.maxGap    = maxGap;
    r = fit.rec - (p(1)*fit.dpx + p(2));
    fprintf(['matched %d of %d DATAPixx pulses in %d part(s); drift %+.2f ppm; ' ...
             'pulse-fit residual rms %.3f ms, max %.3f ms\n'], numel(i), numel(dpx), ...
            1 + sum(diff(fit.dpx) > maxGap), driftPpm, 1e3*sqrt(mean(r.^2)), 1e3*max(abs(r)));
end


function x = checkEdges(x, name)
    if ~isvector(x) || ~isnumeric(x)
        error('fitClock:badInput', '%s: expected a vector of pulse times.', name);
    end
    x = double(x(:));
    if any(~isfinite(x))
        error('fitClock:badInput', '%s: contains NaN or infinite times.', name);
    end
    k = find(diff(x) <= 0, 1);
    if ~isempty(k)
        error('fitClock:notIncreasing', ['%s: times must increase, but pulse %d (%.6f s) is not ' ...
              'after pulse %d. Split a clock reset or joined sessions and align each part.'], ...
              name, k+1, x(k+1), k);
    end
end


function [offsets, anchors] = candidates(dpx, rec)
% Every (offset, anchor) the start beats allow: the first start beat on the
% recorder against each start beat on the DATAPixx3.
    a = findBeats(dpx);
    a = [a(strcmp({a.kind}, 'start')).t];
    b = findBeats(rec);
    b = b(strncmp({b.kind}, 'start', 5));
    if isempty(a) || isempty(b)
        error('fitClock:noStart', 'No start beat found on one side; check the cable and the bit.');
    end
    shifts = 0;
    if strcmp(b(1).kind, 'start?')
        g = rec(rec >= b(1).t & rec < b(1).t + 0.5);
        if abs(g(2) - g(1) - 0.1) < 0.04
            shifts = [0 0.1];                    % the first or the last pulse was lost
        end
    end
    [A, S] = meshgrid(a, shifts);
    offsets = b(1).t - S(:) - A(:);
    anchors = A(:);
end


function [i, j] = pairEdges(dpx, rec, p, tol)
% For each DATAPixx pulse, the nearest recorder pulse within tol of the
% line. Each recorder pulse is used once, by its closest match.
    pred = p(1)*dpx + p(2);
    j = interp1(rec, (1:numel(rec))', pred, 'nearest', 'extrap');
    err = abs(rec(j) - pred);
    cand = find(err < tol);
    [~, s] = sort(err(cand));
    cand = cand(s);
    [~, first] = unique(j(cand), 'first');
    i = sort(cand(first));
    j = j(i);
end


function [i, j, p] = growFit(dpx, rec, offset, anchor, tol)
% Match pulses outward from the anchor beat, refitting the line each time
% the window doubles, so drift never pushes a true match out of reach.
    p = [1 offset];
    half = 10;
    reach = max(dpx(end) - anchor, anchor - dpx(1));
    while true
        win = find(abs(dpx - anchor) <= half);
        [i, j] = pairEdges(dpx(win), rec, p, tol);
        i = win(i);
        if numel(i) >= 2
            p = polyfit(dpx(i), rec(j), 1);
        end
        if half >= reach
            break;
        end
        half = 2*half;
    end
    [i, j] = pairEdges(dpx, rec, p, tol);
    if numel(i) >= 2
        p = polyfit(dpx(i), rec(j), 1);
    end
end


% ---- Step 6: move events between clocks ---------------------------------
function out = mapToRecorder(fit, t)
% DATAPixx3 times onto the recorder clock, interpolating between matched
% pulses. NaN outside them, or across a gap longer than fit.maxGap.
    out = interp1(fit.dpx, fit.rec, t, 'linear', NaN);
    x = fit.dpx(:);
    for k = 1:numel(t)
        i = find(x <= t(k), 1, 'last');
        ok = ~isempty(i) && (x(i) == t(k) || (i < numel(x) && x(i+1) - x(i) <= fit.maxGap));
        if ~ok
            out(k) = NaN;
        end
    end
end

Python

# Align an EEG recording to the DATAPixx3 clock using one synchronization
# pulse, then move button presses onto the EEG clock. Needs Python 3 and NumPy.
import numpy as np

# ---- Your settings --------------------------------------------------------
SYNC_BIT = 0                # DOUT bit that carried the pulse (same as the setup)
EEG_FS = 1000               # EEG sampling rate, Hz
EEG_TRIGGER_BIT = 0         # bit of the EEG trigger channel the BNC is wired to

PAIR_TOL_S = 0.010          # a recorder pulse must sit within 10 ms of where it is expected
MAX_GAP_S = 2.5             # no mapping across gaps with no pulses longer than this
MAX_ABS_PPM = 1000.0        # larger drift means the pairing went wrong
MIN_PAIRS = 10              # fewer matched pulses than this is refused

CLUSTER_GAP_S = 0.5         # pulses closer than this belong to one beat
START_SPACING_S = 0.100     # start beat: 3 pulses 100 ms apart
STOP_SPACING_S = 0.300      # stop beat: 2 pulses 300 ms apart
SPACING_TOL_S = 0.040


class AlignmentError(ValueError):
    """The pulses do not give one clear, trustworthy match."""


# ---- Step 1: find the pulses ----------------------------------------------
def rising_edges(t, word, bit):
    """Times where one bit of `word` switches from 0 to 1."""
    level = (np.asarray(word, np.int64) >> bit) & 1
    rises = np.flatnonzero(np.diff(level) == 1) + 1
    return np.asarray(t, float)[rises]


# ---- Step 2: find the beats -----------------------------------------------
def find_beats(edges):
    """List of (kind, time of first pulse) for each group of close pulses.
    "start?" is a start beat that lost one pulse; it still never looks like
    the 300 ms stop beat."""
    edges = np.asarray(edges, float)
    if edges.size == 0:
        return []
    groups = np.split(edges, np.flatnonzero(np.diff(edges) > CLUSTER_GAP_S) + 1)
    beats = []
    for g in groups:
        if g.size < 2:
            continue
        gaps = np.diff(g)
        if g.size == 3 and np.all(np.abs(gaps - START_SPACING_S) < SPACING_TOL_S):
            beats.append(("start", g[0]))
        elif g.size == 2 and abs(gaps[0] - STOP_SPACING_S) < SPACING_TOL_S:
            beats.append(("stop", g[0]))
        elif g.size == 2 and (abs(gaps[0] - START_SPACING_S) < SPACING_TOL_S
                              or abs(gaps[0] - 2 * START_SPACING_S) < SPACING_TOL_S):
            beats.append(("start?", g[0]))
        else:
            beats.append(("unknown", g[0]))
    return beats


# ---- Steps 3 to 5: match the pulses and fit the clocks --------------------
def check_edges(x, name):
    x = np.asarray(x, float)
    if x.ndim != 1:
        raise AlignmentError(f"{name}: expected a 1-D list of pulse times")
    if not np.all(np.isfinite(x)):
        raise AlignmentError(f"{name}: contains NaN or infinite times")
    bad = np.flatnonzero(np.diff(x) <= 0)
    if bad.size:
        k = int(bad[0])
        raise AlignmentError(f"{name}: times must increase, but pulse {k + 1} ({x[k + 1]:.6f} s) "
                             f"is not after pulse {k}. Split a clock reset or joined sessions "
                             "and align each part.")
    return x


def pair(dpx, rec, slope, intercept):
    """For each DATAPixx pulse, the nearest recorder pulse within PAIR_TOL_S
    of the line. Each recorder pulse is used once, by its closest match."""
    pred = slope * dpx + intercept
    j = np.clip(np.searchsorted(rec, pred), 1, len(rec) - 1)
    j = np.where(np.abs(rec[j - 1] - pred) < np.abs(rec[j] - pred), j - 1, j)
    err = np.abs(rec[j] - pred)
    cand = np.flatnonzero(err < PAIR_TOL_S)
    cand = cand[np.argsort(err[cand], kind="stable")]
    _, first = np.unique(j[cand], return_index=True)
    i = np.sort(cand[first])
    return i, j[i]


def candidates(dpx, rec):
    """Every (offset, anchor) the start beats allow: the first start beat on
    the recorder against each start beat on the DATAPixx3."""
    a = [t for k, t in find_beats(dpx) if k == "start"]
    b = [(k, t) for k, t in find_beats(rec) if k.startswith("start")]
    if not a or not b:
        raise AlignmentError("no start beat found on one side; check the cable and the bit")
    kind, t_rec = b[0]
    shifts = [0.0]
    if kind == "start?":
        g = rec[(rec >= t_rec) & (rec < t_rec + CLUSTER_GAP_S)]
        if abs(g[1] - g[0] - START_SPACING_S) < SPACING_TOL_S:
            shifts = [0.0, START_SPACING_S]       # the first or the last pulse was lost
    return [(t_rec - s - t_dpx, t_dpx) for t_dpx in a for s in shifts]


def grow_fit(dpx, rec, offset, anchor):
    """Match pulses outward from the anchor beat, refitting the line each
    time the window doubles, so drift never pushes a true match out of reach."""
    slope = 1.0
    intercept = offset
    half = 10.0
    reach = max(dpx[-1] - anchor, anchor - dpx[0])
    while True:
        win = np.flatnonzero(np.abs(dpx - anchor) <= half)
        i, j = pair(dpx[win], rec, slope, intercept)
        i = win[i]
        if len(i) >= 2:
            slope, intercept = np.polyfit(dpx[i], rec[j], 1)
        if half >= reach:
            break
        half *= 2
    i, j = pair(dpx, rec, slope, intercept)
    if len(i) >= 2:
        slope, intercept = np.polyfit(dpx[i], rec[j], 1)
    return i, j, slope, intercept


class ClockFit:
    """Matched pulses and the fitted line recorder = slope * datapixx + intercept."""

    def __init__(self, dpx, rec, slope, intercept, n_dpx):
        self.dpx = dpx
        self.rec = rec
        self.slope = slope
        self.intercept = intercept
        self.n_dpx = n_dpx

    def supported(self, t_dpx):
        """True on a matched pulse, or between two matched pulses at most MAX_GAP_S apart."""
        t = np.asarray(t_dpx, float)
        i = np.searchsorted(self.dpx, t, side="right")
        inside = (i > 0) & (i < len(self.dpx))
        k = np.clip(i, 1, len(self.dpx) - 1)
        short = (self.dpx[k] - self.dpx[k - 1]) <= MAX_GAP_S
        return (inside & short) | np.isin(t, self.dpx)

    def to_recorder(self, t_dpx):
        """DATAPixx3 times onto the recorder clock, interpolating between matched pulses."""
        t = np.asarray(t_dpx, float)
        out = np.interp(t, self.dpx, self.rec)
        return np.where(self.supported(t), out, np.nan)

    def report(self):
        drift = (self.slope - 1.0) * 1e6
        r = self.rec - (self.slope * self.dpx + self.intercept)
        n_parts = 1 + int(np.sum(np.diff(self.dpx) > MAX_GAP_S))
        return (f"matched {len(self.dpx)} of {self.n_dpx} DATAPixx pulses in {n_parts} part(s); "
                f"drift {drift:+.2f} ppm; pulse-fit residual rms {np.sqrt(np.mean(r ** 2)) * 1e3:.3f} ms, "
                f"max {np.max(np.abs(r)) * 1e3:.3f} ms")


def fit_clock(dpx_edges, rec_edges, offset=None):
    """Match every DATAPixx pulse with the recorder pulse it produced, then fit.
    Pass `offset` (recorder minus DATAPixx time, s) only if the start beats
    are ambiguous."""
    dpx = check_edges(dpx_edges, "DATAPixx pulses")
    rec = check_edges(rec_edges, "recorder pulses")
    if len(dpx) < 2 or len(rec) < 2:
        raise AlignmentError("need at least two pulses on each side")
    if offset is None:
        cands = candidates(dpx, rec)
    else:
        cands = [(offset, dpx[np.argmin(np.abs(dpx + offset - rec[0]))])]
    fits = [grow_fit(dpx, rec, o, a) for o, a in cands]
    order = sorted(range(len(fits)), key=lambda k: -len(fits[k][0]))
    i, j, slope, intercept = fits[order[0]]
    for k in order[1:]:
        if len(fits[k][0]) >= len(i) and abs(fits[k][3] - intercept) > PAIR_TOL_S:
            raise AlignmentError(f"ambiguous: the recording fits the blocks starting at "
                                 f"{cands[order[0]][1]:.3f} s and {cands[k][1]:.3f} s equally well. "
                                 "Pass offset=, or give blocks different lengths.")
    if len(i) < max(2, MIN_PAIRS):
        raise AlignmentError(f"only {len(i)} pulses matched (need {max(2, MIN_PAIRS)}); check the bit, "
                             "the cable and that the two recordings overlap")
    drift = (slope - 1.0) * 1e6
    if abs(drift) > MAX_ABS_PPM:
        raise AlignmentError(f"fitted drift {drift:+.0f} ppm is implausible; the matching is wrong")
    return ClockFit(dpx[i], rec[j], float(slope), float(intercept), len(dpx))


# ---- Run it ---------------------------------------------------------------
log = np.genfromtxt("datapixx_log.csv", delimiter=",", names=True)   # written by the setup example
log_t = log["t_datapixx_s"]
log_dout = log["dout24"].astype(int)
eeg_trigger = np.loadtxt("eeg_trigger.csv").astype(int)   # EEG trigger channel exported from your EEG software, one value per sample

# Step 1: the pulse times on each clock
dpx = rising_edges(log_t, log_dout, SYNC_BIT)
eeg_t = np.arange(len(eeg_trigger)) / EEG_FS
eeg = rising_edges(eeg_t, eeg_trigger, EEG_TRIGGER_BIT)

# Steps 2 to 5: match the pulses and fit the EEG clock to the DATAPixx3 clock
fit = fit_clock(dpx, eeg)
print(fit.report())

# Step 6: move button presses from the DATAPixx3 clock to the EEG clock
button_t_dpx = np.array([12.250, 30.875, 47.500, 95.125])   # fake presses, DATAPixx3 clock (s)
button_t_eeg = fit.to_recorder(button_t_dpx)
for t_dpx, t_eeg in zip(button_t_dpx, button_t_eeg):
    print(f"button at {t_dpx:8.3f} s on the DATAPixx3 = {t_eeg:8.3f} s on the EEG")