Skip to content

Fix digital channel reading via high-level read() API #99

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
stevevanhooser merged 3 commits into from
Apr 10, 2026
Merged

Fix digital channel reading via high-level read() API #99

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5976556
Fix reading digital channels in neuropixelsGLX via high-level read() API
claude Apr 10, 2026
3b762cd
Expose neuropixelsGLX digital word as per-bit di channels
claude Apr 10, 2026
c17191a
Pre-compute digital line (col, bit) map and label in header
claude Apr 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 99 additions & 2 deletions +ndr/+format/+neuropixelsGLX/header.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,27 @@
% Number of saved channels
info.n_saved_chans = str2double(meta.nSavedChans);

% Parse snsApLfSy or snsMnMaXaDw to determine neural vs sync channels
% Parse snsApLfSy or snsMnMaXaDw to determine neural vs sync channels.
% Also compute, for digital lines:
% n_digital_word_cols : number of int16 columns in the .bin file
% that hold digital word data (stored last).
% n_digital_lines : number of single-bit digital lines exposed.
% digital_line_col : (n_digital_lines x 1) 0-based DW column
% offset (0 = first DW column).
% digital_line_bit : (n_digital_lines x 1) 0-based bit position
% within that column (0..15).
% digital_line_label : (n_digital_lines x 1) cellstr describing
% the underlying SpikeGLX line, e.g. 'XD0'
% for port-0 line 0, 'XD1.3' for port-1 line
% 3, or 'SY0.6' for sync col 0 bit 6.
%
% For NIDQ streams the count of active lines comes from niXDBytes1
% and niXDBytes2 (bytes captured per port). NI-DAQ hardware only
% enables digital input in whole-byte chunks, so every bit within a
% captured byte is electrically active even if the user only wired
% some of them; niXDChans1/2 is just informational and is not used
% to gate which lines are exposed. For IMEC streams there is no
% per-bit configuration; all 16 bits of each sync int16 are exposed.
if isfield(meta, 'snsApLfSy')
% imec stream: AP,LF,SY counts
counts = sscanf(meta.snsApLfSy, '%d,%d,%d');
Expand All @@ -82,21 +102,98 @@
info.stream_type = 'ap';
info.n_neural_chans = counts(1);
end
% IMEC sync word is an int16; each sync column provides 16 bits.
% In practice only bit 6 is the SMA sync input, but all 16 bits
% are exposed as independent digital lines so callers can pick
% whichever they need.
info.n_digital_word_cols = info.n_sync_chans;
n_lines = 16 * info.n_sync_chans;
info.n_digital_lines = n_lines;
info.digital_line_col = zeros(n_lines, 1);
info.digital_line_bit = zeros(n_lines, 1);
info.digital_line_label = cell(n_lines, 1);
idx = 0;
for c = 0:(info.n_sync_chans - 1)
for b = 0:15
idx = idx + 1;
info.digital_line_col(idx) = c;
info.digital_line_bit(idx) = b;
info.digital_line_label{idx} = sprintf('SY%d.%d', c, b);
end
end
elseif isfield(meta, 'snsMnMaXaDw')
% NI-DAQ stream: MN,MA,XA,DW
info.stream_type = 'nidq';
counts = sscanf(meta.snsMnMaXaDw, '%d,%d,%d,%d');
info.n_mn_chans = counts(1); % multiplexed neural
info.n_ma_chans = counts(2); % multiplexed analog
info.n_xa_chans = counts(3); % non-multiplexed analog
info.n_dw_chans = counts(4); % digital words
info.n_dw_chans = counts(4); % digital word int16 columns
info.n_neural_chans = counts(1) + counts(2) + counts(3);
info.n_sync_chans = counts(4);
info.n_digital_word_cols = counts(4);

% Bytes saved per port. NI hardware only enables digital input
% in whole-byte chunks, so each saved byte = 8 active lines
% regardless of how many of them the user actually wired up.
n_bytes_p0 = 0;
if isfield(meta, 'niXDBytes1')
n_bytes_p0 = str2double(meta.niXDBytes1);
end
n_bytes_p1 = 0;
if isfield(meta, 'niXDBytes2')
n_bytes_p1 = str2double(meta.niXDBytes2);
end

% If neither byte field is present, fall back to assuming every
% bit of every DW int16 column is in use (16 lines per column).
if n_bytes_p0 == 0 && n_bytes_p1 == 0
n_lines_p0 = 16 * info.n_dw_chans;
n_lines_p1 = 0;
else
n_lines_p0 = 8 * n_bytes_p0;
n_lines_p1 = 8 * n_bytes_p1;
end

% Compute the (col, bit) position of each active line.
% SpikeGLX storage layout: port0 lines occupy the first
% n_bytes_p0*8 bits of the concatenated digital bit stream,
% then port1 lines occupy the next n_bytes_p1*8 bits. The bit
% stream is laid out across n_dw_chans int16 columns (16 bits
% per column).
n_lines = n_lines_p0 + n_lines_p1;
info.n_digital_lines = n_lines;
info.digital_line_col = zeros(n_lines, 1);
info.digital_line_bit = zeros(n_lines, 1);
info.digital_line_label = cell(n_lines, 1);
idx = 0;
for k = 0:(n_lines_p0 - 1)
abs_bit = k;
idx = idx + 1;
info.digital_line_col(idx) = floor(abs_bit / 16);
info.digital_line_bit(idx) = mod(abs_bit, 16);
info.digital_line_label{idx} = sprintf('XD%d', k);
end
for k = 0:(n_lines_p1 - 1)
abs_bit = n_bytes_p0 * 8 + k;
idx = idx + 1;
info.digital_line_col(idx) = floor(abs_bit / 16);
info.digital_line_bit(idx) = mod(abs_bit, 16);
info.digital_line_label{idx} = sprintf('XD1.%d', k);
end
else
% Fallback
info.stream_type = 'unknown';
info.n_neural_chans = info.n_saved_chans - 1;
info.n_sync_chans = 1;
info.n_digital_word_cols = 1;
info.n_digital_lines = 16;
info.digital_line_col = zeros(16, 1);
info.digital_line_bit = (0:15)';
info.digital_line_label = cell(16, 1);
for b = 0:15
info.digital_line_label{b+1} = sprintf('bit%d', b);
end
end

% Parse saved channel subset
Expand Down
56 changes: 45 additions & 11 deletions +ndr/+reader/neuropixelsGLX.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,11 @@
%
% Channel mapping:
% - Neural channels are exposed as 'analog_in' (ai1..aiN)
% - The sync word is exposed as 'digital_in' (di1)
% - Digital lines are exposed as 'digital_in' (di1..diM), where each
% di channel is a single bit of the packed digital word(s). The
% number of lines is determined from metadata: for NIDQ streams
% it is 8 * (niXDBytes1 + niXDBytes2); for IMEC streams it is
% 16 * n_sync_chans (bit 6 is the SMA sync input in practice).
% - A single time channel 't1' is always present
%
% Data is returned as int16 to preserve native precision. Use
Expand Down Expand Up @@ -90,8 +94,9 @@
%
% Returns a structure array with fields 'name', 'type', and
% 'time_channel'. Neural channels are 'analog_in' (ai1..aiN),
% the sync channel is 'digital_in' (di1), and a time channel
% 't1' is always present.
% digital lines are 'digital_in' (di1..diM) with one entry
% per single-bit line in the packed digital word(s), and a
% time channel 't1' is always present.
%
% See also: ndr.format.neuropixelsGLX.header

Expand All @@ -109,10 +114,12 @@
'type', 'analog_in', 'time_channel', 1); %#ok<AGROW>
end

% Sync channel (digital_in)
if info.n_sync_chans > 0
channels(end+1) = struct('name', 'di1', ...
'type', 'digital_in', 'time_channel', 1);
% Digital lines (digital_in) — one per bit of the packed
% digital word(s). n_digital_lines comes from metadata
% (niXDBytes1/niXDBytes2 for NIDQ, 16*n_sync_chans for IMEC).
for i = 1:info.n_digital_lines
channels(end+1) = struct('name', ['di' int2str(i)], ...
'type', 'digital_in', 'time_channel', 1); %#ok<AGROW>
end
end

Expand Down Expand Up @@ -167,7 +174,9 @@
%
% For 'analog_in': returns int16 neural data.
% For 'time': returns double time stamps in seconds.
% For 'digital_in': returns int16 sync word values.
% For 'digital_in': returns int16 single-bit values (0 or 1)
% extracted from the packed digital word(s).
% CHANNEL gives the 1-based digital line(s).
%
% See also: ndr.format.neuropixelsGLX.read

Expand All @@ -189,9 +198,34 @@
data = read_samples(binfile, info, uint32(channel), s0, s1);

case {'digital_in', 'di'}
% Sync channel is the last channel in the file
sync_chan = info.n_saved_chans;
data = read_samples(binfile, info, uint32(sync_chan), s0, s1);
% CHANNEL is a vector of 1-based digital line indices
% (di1..di_n_digital_lines). The header pre-computes
% the (DW column, bit position) for each active line
% from niXDBytes1/niXDBytes2 (NIDQ) or n_sync_chans
% (IMEC), so the reader just looks up the mapping
% and extracts the requested bits with bitget.
line_idx = double(channel(:));
if any(line_idx < 1) || ...
any(line_idx > info.n_digital_lines)
error('ndr:reader:neuropixelsGLX:DigitalLineOutOfRange', ...
'Digital line out of range; valid lines are 1..%d.', ...
info.n_digital_lines);
end
first_dw_col = info.n_saved_chans - info.n_digital_word_cols + 1;
col_offsets = info.digital_line_col(line_idx);
bit_pos = info.digital_line_bit(line_idx);

n_samples = double(s1) - double(s0) + 1;
data = zeros(n_samples, numel(channel), 'int16');
unique_cols = unique(col_offsets);
for u = 1:numel(unique_cols)
file_col = first_dw_col + unique_cols(u);
raw = read_samples(binfile, info, uint32(file_col), s0, s1);
idx = find(col_offsets == unique_cols(u));
for k = 1:numel(idx)
data(:, idx(k)) = int16(bitget(raw, bit_pos(idx(k)) + 1));
end
end

otherwise
error('ndr:reader:neuropixelsGLX:UnknownChannelType', ...
Expand Down
128 changes: 128 additions & 0 deletions +ndr/+test/+reader/+neuropixelsGLX/test.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,12 @@ function test(varargin)
for c = 1:nNeuralChans
data(:, c) = int16(round(500 * sin(2 * pi * c * t_vec)));
end
% Sync channel is a 16-bit packed digital word. Each bit is exposed as
% a separate digital line (di1..di16). Use a counter pattern that
% exercises bits 0..14 (mod 2^15 keeps values non-negative so the int16
% representation is unambiguous).
sync_pattern = int16(mod((0:nSamples-1), 2^15));
data(:, nTotalChans) = sync_pattern(:);

% Write binary
fid = fopen(binfile, 'w', 'ieee-le');
Expand Down Expand Up @@ -102,6 +108,128 @@ function test(varargin)
disp(['Max error vs expected: ' num2str(max_error)]);
assert(max_error == 0, 'Data mismatch detected!');

% Verify that the IMEC sync word is exposed as 16 single-bit digital
% lines (di1..di16), one per bit of the int16 sync column.
di_idx = find(strcmp({channels.type}, 'digital_in'));
assert(numel(di_idx) == 16, ...
sprintf('Expected 16 IMEC digital lines, got %d.', numel(di_idx)));

% Read individual digital lines through the high-level read() API.
% This is the code path that previously returned [] because
% ndr.reader.read() routed 'digital_in' to readevents_epochsamples
% (which is abstract for this format).
t1_end = (nSamples-1) / SR;
[d_di1, t_di] = r.read({metafile}, 'di1', 't0', 0, 't1', t1_end);
disp(['Read ' int2str(size(d_di1, 1)) ' samples from channel di1 via r.read().']);
assert(~isempty(d_di1), 'Digital di1 read returned empty data.');
assert(~isempty(t_di), 'Digital di1 read returned empty time.');
assert(size(d_di1, 1) == nSamples, ...
sprintf('Digital sample count mismatch: got %d, expected %d.', ...
size(d_di1, 1), nSamples));
% di1 is bit 0 of the sync word (alternates 0,1,0,1,...).
expected_di1 = int16(bitget(sync_pattern(:), 1));
assert(isequal(d_di1, expected_di1), 'di1 (bit 0) extraction mismatch.');

% di8 is bit 7 of the sync word.
[d_di8, ~] = r.read({metafile}, 'di8', 't0', 0, 't1', t1_end);
expected_di8 = int16(bitget(sync_pattern(:), 8));
assert(isequal(d_di8, expected_di8), 'di8 (bit 7) extraction mismatch.');

% di12 is bit 11 of the sync word (2^11 = 2048; the pattern reaches
% 2999 so this bit is non-trivially exercised).
[d_di12, ~] = r.read({metafile}, 'di12', 't0', 0, 't1', t1_end);
expected_di12 = int16(bitget(sync_pattern(:), 12));
assert(any(expected_di12 ~= 0), ...
'Test pattern should exercise bit 11; check sync_pattern.');
assert(isequal(d_di12, expected_di12), 'di12 (bit 11) extraction mismatch.');

% === NIDQ format test ===
% Mirrors the user-reported configuration: snsMnMaXaDw=0,0,8,1 with
% niXDBytes1=1, so 8 analog inputs plus 8 digital lines packed into
% the low byte of a single int16 DW column.
disp('--- NIDQ format test ---');

nXA = 8;
nDW = 1;
nNidqChans = nXA + nDW;
nNidqSamples = 500;
SR_ni = 10593.220339;

nidq_subdir = fullfile(tempdir_path, 'nidq_g0');
mkdir(nidq_subdir);
nidq_metafile = fullfile(nidq_subdir, 'nidq_g0_t0.nidq.meta');
nidq_binfile = fullfile(nidq_subdir, 'nidq_g0_t0.nidq.bin');

% Build data: XA0..XA7 sine waves + DW column with an 8-bit pattern.
ni_data = zeros(nNidqSamples, nNidqChans, 'int16');
t_vec_ni = (0:nNidqSamples-1)' / SR_ni;
for c = 1:nXA
ni_data(:, c) = int16(round(1000 * sin(2 * pi * c * t_vec_ni)));
end
ni_sync = int16(mod((0:nNidqSamples-1), 2^8)); % low byte: 0..255
ni_data(:, end) = ni_sync(:);

fid = fopen(nidq_binfile, 'w', 'ieee-le');
fwrite(fid, reshape(ni_data', 1, []), 'int16');
fclose(fid);

fid = fopen(nidq_metafile, 'w');
fprintf(fid, 'niSampRate=%.6f\n', SR_ni);
fprintf(fid, 'nSavedChans=%d\n', nNidqChans);
fprintf(fid, 'snsMnMaXaDw=0,0,%d,%d\n', nXA, nDW);
fprintf(fid, 'snsSaveChanSubset=all\n');
fprintf(fid, 'fileSizeBytes=%d\n', nNidqSamples * nNidqChans * 2);
fprintf(fid, 'fileTimeSecs=%.6f\n', nNidqSamples / SR_ni);
fprintf(fid, 'niAiRangeMax=5\n');
fprintf(fid, 'niAiRangeMin=-5\n');
fprintf(fid, 'niMaxInt=32768\n');
fprintf(fid, 'niXDBytes1=1\n');
fprintf(fid, 'niXDChans1=0:7\n');
fprintf(fid, 'niXAChans1=0:7\n');
fprintf(fid, 'typeThis=nidq\n');
fclose(fid);

r_ni = ndr.reader('neuropixelsGLX');
ni_channels = r_ni.getchannelsepoch({nidq_metafile}, 1);
ni_di_count = sum(strcmp({ni_channels.type}, 'digital_in'));
assert(ni_di_count == 8, ...
sprintf('Expected 8 NIDQ digital lines (niXDBytes1=1), got %d.', ...
ni_di_count));
ni_ai_count = sum(strcmp({ni_channels.type}, 'analog_in'));
assert(ni_ai_count == nXA, ...
sprintf('Expected %d NIDQ analog lines, got %d.', nXA, ni_ai_count));

% Read di1 via r.read() — the exact call the user reported failing.
t1_ni = (nNidqSamples-1) / SR_ni;
[d_ni_di1, t_ni_di1] = r_ni.read({nidq_metafile}, 'di1', 't0', 0, 't1', t1_ni);
assert(~isempty(d_ni_di1), 'NIDQ di1 returned empty data (user-reported bug).');
assert(~isempty(t_ni_di1), 'NIDQ di1 returned empty time.');
assert(size(d_ni_di1, 1) == nNidqSamples, ...
sprintf('NIDQ di1 sample count mismatch: got %d, expected %d.', ...
size(d_ni_di1, 1), nNidqSamples));
expected_ni_di1 = int16(bitget(ni_sync(:), 1));
assert(isequal(d_ni_di1, expected_ni_di1), 'NIDQ di1 (bit 0) mismatch.');

% Read di8 (bit 7, the high bit of the byte).
[d_ni_di8, ~] = r_ni.read({nidq_metafile}, 'di8', 't0', 0, 't1', t1_ni);
expected_ni_di8 = int16(bitget(ni_sync(:), 8));
assert(isequal(d_ni_di8, expected_ni_di8), 'NIDQ di8 (bit 7) mismatch.');

% di9 must be out of range for an 8-line NIDQ configuration. Call
% readchannels_epochsamples directly because the high-level read()
% would fail earlier at channel-name lookup (di9 isn't listed).
threw = false;
try
r_ni.readchannels_epochsamples('digital_in', 9, ...
{nidq_metafile}, 1, 1, 10);
catch ME
threw = strcmp(ME.identifier, ...
'ndr:reader:neuropixelsGLX:DigitalLineOutOfRange');
end
assert(threw, 'Expected DigitalLineOutOfRange error for line 9 in 8-line NIDQ.');

disp('NIDQ digital line read: OK.');

disp('All checks passed.');

if plotit
Expand Down
3 changes: 2 additions & 1 deletion +ndr/reader.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,8 @@
if b,
switch (channelstruct(1).ndr_type),
% readchannels_epochsamples
case {'analog_input','analog_output','analog_in','analog_out','ai','ao'},
case {'analog_input','analog_output','analog_in','analog_out','ai','ao', ...
'digital_input','digital_output','digital_in','digital_out','di','do'},

Check notice

Code scanning / Code Analyzer

Extra comma is unnecessary. Note

Extra comma is unnecessary.
if ~useSamples, % must compute the samples to be read
s0 = round(1+t0*channelstruct(1).samplerate);
s1 = round(1+t1*channelstruct(1).samplerate);
Expand Down
Loading
Loading