diff --git a/+ndr/+format/+neuropixelsGLX/header.m b/+ndr/+format/+neuropixelsGLX/header.m index 6582a37..f2b3299 100644 --- a/+ndr/+format/+neuropixelsGLX/header.m +++ b/+ndr/+format/+neuropixelsGLX/header.m @@ -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'); @@ -82,6 +102,25 @@ 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'; @@ -89,14 +128,72 @@ 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 diff --git a/+ndr/+reader/neuropixelsGLX.m b/+ndr/+reader/neuropixelsGLX.m index dfef9c5..c6a01df 100644 --- a/+ndr/+reader/neuropixelsGLX.m +++ b/+ndr/+reader/neuropixelsGLX.m @@ -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 @@ -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 @@ -109,10 +114,12 @@ 'type', 'analog_in', 'time_channel', 1); %#ok 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 end end @@ -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 @@ -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', ... diff --git a/+ndr/+test/+reader/+neuropixelsGLX/test.m b/+ndr/+test/+reader/+neuropixelsGLX/test.m index eb7bdf4..eda3e26 100644 --- a/+ndr/+test/+reader/+neuropixelsGLX/test.m +++ b/+ndr/+test/+reader/+neuropixelsGLX/test.m @@ -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'); @@ -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 diff --git a/+ndr/reader.m b/+ndr/reader.m index e86e694..5d38382 100644 --- a/+ndr/reader.m +++ b/+ndr/reader.m @@ -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'}, if ~useSamples, % must compute the samples to be read s0 = round(1+t0*channelstruct(1).samplerate); s1 = round(1+t1*channelstruct(1).samplerate); diff --git a/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m index 4bef82f..dddec10 100644 --- a/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m +++ b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m @@ -160,8 +160,10 @@ function testGetChannelsEpoch(testCase) %TESTGETCHANNELSEPOCH Verify channel listing. channels = testCase.Reader.getchannelsepoch({testCase.MetaFilename}, 1); - % Should have: 1 time + N neural + 1 sync = N+2 - expectedTotal = testCase.NumNeuralChansActual + 2; + % IMEC sync word is a single int16 column = 16 single-bit + % digital lines (di1..di16). Total: 1 time + N neural + 16 di. + nDigitalLines = 16; + expectedTotal = testCase.NumNeuralChansActual + 1 + nDigitalLines; testCase.verifyNumElements(channels, expectedTotal, 'Wrong number of channels.'); % First channel should be time @@ -176,9 +178,15 @@ function testGetChannelsEpoch(testCase) ['Wrong type for neural channel ' int2str(i)]); end - % Sync channel (last) - testCase.verifyEqual(channels(end).name, 'di1', 'Last channel should be di1.'); - testCase.verifyEqual(channels(end).type, 'digital_in', 'Last channel type should be digital_in.'); + % Digital lines: di1..di16 follow the neural channels. + di_start = 1 + testCase.NumNeuralChansActual + 1; % after time + neural + for i = 1:nDigitalLines + idx = di_start + i - 1; + testCase.verifyEqual(channels(idx).name, ['di' int2str(i)], ... + ['Wrong name for digital line ' int2str(i)]); + testCase.verifyEqual(channels(idx).type, 'digital_in', ... + ['Wrong type for digital line ' int2str(i)]); + end end function testSampleRate(testCase)