diff --git a/+ndr/+format/+neuropixelsGLX/header.m b/+ndr/+format/+neuropixelsGLX/header.m index 7b95c9e..6582a37 100644 --- a/+ndr/+format/+neuropixelsGLX/header.m +++ b/+ndr/+format/+neuropixelsGLX/header.m @@ -83,11 +83,15 @@ info.n_neural_chans = counts(1); end elseif isfield(meta, 'snsMnMaXaDw') - % NI-DAQ stream + % NI-DAQ stream: MN,MA,XA,DW info.stream_type = 'nidq'; counts = sscanf(meta.snsMnMaXaDw, '%d,%d,%d,%d'); - info.n_neural_chans = counts(1); % MN channels - info.n_sync_chans = 0; + 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_neural_chans = counts(1) + counts(2) + counts(3); + info.n_sync_chans = counts(4); else % Fallback info.stream_type = 'unknown'; @@ -107,6 +111,10 @@ vmax = str2double(meta.imAiRangeMax); vmin = str2double(meta.imAiRangeMin); info.voltage_range = [vmin vmax]; + elseif isfield(meta, 'niAiRangeMax') + vmax = str2double(meta.niAiRangeMax); + vmin = str2double(meta.niAiRangeMin); + info.voltage_range = [vmin vmax]; else info.voltage_range = [-0.6 0.6]; % Neuropixels 1.0 default end @@ -114,10 +122,20 @@ % Max integer value if isfield(meta, 'imMaxInt') info.max_int = str2double(meta.imMaxInt); + elseif isfield(meta, 'niMaxInt') + info.max_int = str2double(meta.niMaxInt); else info.max_int = 512; % Neuropixels 1.0 default end + % NI-DAQ gains + if isfield(meta, 'niMNGain') + info.ni_mn_gain = str2double(meta.niMNGain); + end + if isfield(meta, 'niMAGain') + info.ni_ma_gain = str2double(meta.niMAGain); + end + % Bits per sample info.bits_per_sample = 16; diff --git a/+ndr/+format/+neuropixelsGLX/read.m b/+ndr/+format/+neuropixelsGLX/read.m index cadd301..dea1a5e 100644 --- a/+ndr/+format/+neuropixelsGLX/read.m +++ b/+ndr/+format/+neuropixelsGLX/read.m @@ -23,10 +23,13 @@ % SR - Sampling rate in Hz (positive double, default: from .meta). % channels - 1-based vector of channel indices to read (default: [] % meaning all channels). Must be within [1, numChans]. +% scale - If true (default), convert raw int16 data to volts +% using ndr.format.neuropixelsGLX.samples2volts. Requires +% a companion .meta file. If false, return raw int16. % % Outputs: -% DATA - N x C int16 matrix, where N is the number of time samples -% and C is the number of channels read. +% DATA - N x C matrix of samples. Double volts if scale is true, +% int16 raw values if scale is false. % T - N x 1 double vector of time points in seconds. % T0_T1 - 1x2 double vector [startTime endTime] for the full file. % @@ -47,21 +50,27 @@ options.numChans (1,1) {mustBeInteger, mustBeNonnegative} = 0 options.SR (1,1) {mustBeNumeric, mustBeNonnegative} = 0 options.channels (1,:) {mustBeNumeric, mustBeInteger, mustBePositive} = [] + options.scale (1,1) logical = true end - % If numChans or SR not provided, read from companion .meta file - if options.numChans == 0 || options.SR == 0 - metafilename = [binfilename(1:end-3) 'meta']; + % Read companion .meta file if needed for numChans/SR or scaling + metafilename = [binfilename(1:end-3) 'meta']; + info = []; + if options.numChans == 0 || options.SR == 0 || options.scale if ~isfile(metafilename) - error('ndr:format:neuropixelsGLX:read:NoMetaFile', ... - 'No .meta file found at %s and numChans/SR not specified.', metafilename); - end - info = ndr.format.neuropixelsGLX.header(metafilename); - if options.numChans == 0 - options.numChans = info.n_saved_chans; - end - if options.SR == 0 - options.SR = info.sample_rate; + if options.numChans == 0 || options.SR == 0 + error('ndr:format:neuropixelsGLX:read:NoMetaFile', ... + 'No .meta file found at %s and numChans/SR not specified.', metafilename); + end + % scale requested but no meta file; will skip scaling below + else + info = ndr.format.neuropixelsGLX.header(metafilename); + if options.numChans == 0 + options.numChans = info.n_saved_chans; + end + if options.SR == 0 + options.SR = info.sample_rate; + end end end @@ -124,12 +133,21 @@ 'byteOrder', 'ieee-le', ... 'headerSkip', uint64(0)); + % Scale to volts if requested and header info is available + if options.scale && ~isempty(info) && ~isempty(data) + data = ndr.format.neuropixelsGLX.samples2volts(data, info, double(channelsToRead)); + end + % Generate time vector if ~isempty(data) t = ndr.time.fun.samples2times((s0_actual:s1_actual)', t0_t1, SR); else t = []; - data = int16([]); + if ~options.scale + data = int16([]); + else + data = double([]); + end end end diff --git a/+ndr/+format/+neuropixelsGLX/samples2volts.m b/+ndr/+format/+neuropixelsGLX/samples2volts.m index 793ecce..b4bdcfd 100644 --- a/+ndr/+format/+neuropixelsGLX/samples2volts.m +++ b/+ndr/+format/+neuropixelsGLX/samples2volts.m @@ -1,7 +1,8 @@ -function volts = samples2volts(data, info) +function volts = samples2volts(data, info, channels) %SAMPLES2VOLTS Convert raw int16 samples to voltage in volts. % % VOLTS = ndr.format.neuropixelsGLX.samples2volts(DATA, INFO) +% VOLTS = ndr.format.neuropixelsGLX.samples2volts(DATA, INFO, CHANNELS) % % Converts raw int16 Neuropixels data to voltage using the gain and range % parameters from the meta file header. @@ -16,16 +17,22 @@ % Per-channel gains are automatically parsed from the imroTbl field. % If imroTbl is absent, default gains are used (500 for AP, 250 for LF). % +% For NI-DAQ streams, gains are determined from niMNGain and niMAGain +% fields in the meta file. +% % Inputs: -% DATA - N x C int16 matrix of raw samples. -% INFO - Header structure from ndr.format.neuropixelsGLX.header. +% DATA - N x C int16 matrix of raw samples. +% INFO - Header structure from ndr.format.neuropixelsGLX.header. +% CHANNELS - Optional 1-based channel indices corresponding to the +% columns of DATA. If omitted, columns are assumed to be +% channels 1:C in order. % % Outputs: % VOLTS - N x C double matrix of voltages in volts. % % Example: % info = ndr.format.neuropixelsGLX.header('run_g0_t0.imec0.ap.meta'); -% [data, ~] = ndr.format.neuropixelsGLX.read('run_g0_t0.imec0.ap.bin', 0, 1); +% [data, ~] = ndr.format.neuropixelsGLX.read('run_g0_t0.imec0.ap.bin', 0, 1, 'scale', false); % volts = ndr.format.neuropixelsGLX.samples2volts(data(:,1:info.n_neural_chans), info); % % See also: ndr.format.neuropixelsGLX.header, ndr.format.neuropixelsGLX.read @@ -33,23 +40,40 @@ arguments data {mustBeNumeric} info (1,1) struct + channels (1,:) {mustBeNumeric, mustBeInteger, mustBePositive} = [] end - vmax = info.voltage_range(2); % imAiRangeMax + vmax = info.voltage_range(2); + n_chans = size(data, 2); - % Parse per-channel gains from imroTbl if available - if isfield(info.meta, 'imroTbl') - gains = parse_imro_gains(info.meta.imroTbl, info.stream_type); - n_chans = size(data, 2); - if numel(gains) >= n_chans - gains = gains(1:n_chans); + if strcmpi(info.stream_type, 'nidq') + % NI-DAQ stream: build gain vector for all channels, then select + all_gains = build_nidq_gains(info, info.n_saved_chans); + if ~isempty(channels) + gains = all_gains(channels); + else + gains = all_gains(1:n_chans); + end + volts = double(data) .* (vmax ./ (info.max_int .* gains)); + elseif isfield(info.meta, 'imroTbl') + % Imec stream with per-channel gains from imroTbl + all_gains = parse_imro_gains(info.meta.imroTbl, info.stream_type); + if ~isempty(channels) + % Pad if needed, then index + if numel(all_gains) < max(channels) + all_gains = [all_gains, repmat(all_gains(end), 1, max(channels) - numel(all_gains))]; + end + gains = all_gains(channels); else - gains = [gains, repmat(gains(end), 1, n_chans - numel(gains))]; + if numel(all_gains) >= n_chans + gains = all_gains(1:n_chans); + else + gains = [all_gains, repmat(all_gains(end), 1, n_chans - numel(all_gains))]; + end end - % Official formula: V = int16 * imAiRangeMax / imMaxInt / gain volts = double(data) .* (vmax ./ (info.max_int .* gains)); else - % Default gain for Neuropixels 1.0 AP band + % Default gain for Neuropixels 1.0 if strcmpi(info.stream_type, 'ap') default_gain = 500; else @@ -61,6 +85,41 @@ end +function gains = build_nidq_gains(info, n_chans) +%BUILD_NIDQ_GAINS Build per-channel gain vector for NI-DAQ streams. +% +% NI-DAQ channels are ordered: MN (neural), MA (auxiliary analog), +% XA (non-multiplexed analog), DW (digital words). +% MN channels use niMNGain, MA channels use niMAGain, XA channels +% have gain=1 (already in volts), DW channels have gain=1. + + mn_gain = 1; + ma_gain = 1; + if isfield(info, 'ni_mn_gain') + mn_gain = info.ni_mn_gain; + end + if isfield(info, 'ni_ma_gain') + ma_gain = info.ni_ma_gain; + end + + n_mn = 0; n_ma = 0; n_xa = 0; + if isfield(info, 'n_mn_chans'), n_mn = info.n_mn_chans; end + if isfield(info, 'n_ma_chans'), n_ma = info.n_ma_chans; end + if isfield(info, 'n_xa_chans'), n_xa = info.n_xa_chans; end + + gains = [repmat(mn_gain, 1, n_mn), ... + repmat(ma_gain, 1, n_ma), ... + ones(1, n_xa)]; + + % Pad or trim to match the number of data columns + if numel(gains) >= n_chans + gains = gains(1:n_chans); + else + gains = [gains, ones(1, n_chans - numel(gains))]; + end +end + + function gains = parse_imro_gains(imroTbl_str, stream_type) %PARSE_IMRO_GAINS Extract per-channel gains from imroTbl string. % diff --git a/.github/badges/tests.svg b/.github/badges/tests.svg index 7dfcaa8..8102a6a 100644 --- a/.github/badges/tests.svg +++ b/.github/badges/tests.svg @@ -1 +1 @@ -teststests108 passed108 passed \ No newline at end of file +teststests121 passed121 passed \ No newline at end of file diff --git a/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m index aac1339..4bef82f 100644 --- a/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m +++ b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX.m @@ -351,14 +351,50 @@ function testFormatReadFunction(testCase) testCase.BinFilename, 0, 0.001, ... 'numChans', testCase.NumTotalChans, ... 'SR', testCase.SR, ... - 'channels', 1:2); + 'channels', 1:2, ... + 'scale', false); - testCase.verifyClass(data, 'int16', 'Format read should return int16.'); + testCase.verifyClass(data, 'int16', 'Format read should return int16 when scale is false.'); testCase.verifySize(data, [size(data,1), 2], 'Should have 2 columns.'); testCase.verifyGreaterThan(numel(t), 0, 'Time vector should not be empty.'); testCase.verifyEqual(t0_t1_range(1), 0, 'AbsTol', 1e-9, 'File should start at t=0.'); end + function testFormatReadScaled(testCase) + %TESTFORMATREADSCALED Verify read returns volts when scale is true. + [data_scaled, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.001, ... + 'numChans', testCase.NumTotalChans, ... + 'SR', testCase.SR, ... + 'channels', 1:2, ... + 'scale', true); + + testCase.verifyClass(data_scaled, 'double', 'Scaled data should be double.'); + + % Compare against manual conversion + [data_raw, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.001, ... + 'numChans', testCase.NumTotalChans, ... + 'SR', testCase.SR, ... + 'channels', 1:2, ... + 'scale', false); + + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + expected = ndr.format.neuropixelsGLX.samples2volts(data_raw, info, [1 2]); + testCase.verifyEqual(data_scaled, expected, 'AbsTol', 1e-15, ... + 'Scaled read should match manual samples2volts conversion.'); + end + + function testFormatReadScaleDefault(testCase) + %TESTFORMATREADSCALEDEFAULT Verify read defaults to scaled output. + [data, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.001, ... + 'channels', 1:2); + + testCase.verifyClass(data, 'double', ... + 'Default read should return double (scaled).'); + end + function testChannelSubsetParsing(testCase) %TESTCHANNELSUBSETPARSING Verify channel subset field in header. info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); diff --git a/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX_nidq.m b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX_nidq.m new file mode 100644 index 0000000..2bda6be --- /dev/null +++ b/tools/tests/+ndr/+unittest/+reader/TestNeuropixelsGLX_nidq.m @@ -0,0 +1,181 @@ +classdef TestNeuropixelsGLX_nidq < matlab.unittest.TestCase +%TESTNEUROPIXELSGLX_NIDQ Unit tests for NI-DAQ stream support in neuropixelsGLX. +% +% Tests header parsing, voltage scaling, and read for NIDQ-format files +% with niMNGain, niMAGain, and snsMnMaXaDw fields. +% +% Example: +% results = runtests('ndr.unittest.reader.TestNeuropixelsGLX_nidq'); + + properties (Constant) + SR = 25000; + NumSamples = 500; + NumMN = 8; % multiplexed neural + NumMA = 2; % multiplexed analog + NumXA = 1; % non-multiplexed analog + NumDW = 1; % digital words + MNGain = 200; + MAGain = 100; + VMax = 5.0; + MaxInt = 32768; + end + + properties (SetAccess=protected) + TempDir char = '' + MetaFilename char = '' + BinFilename char = '' + NumTotalChans double = NaN + end + + methods (TestClassSetup) + function setupOnce(testCase) + testCase.TempDir = fullfile(tempdir, ['ndr_nidq_test_' char(java.util.UUID.randomUUID)]); + mkdir(testCase.TempDir); + + nMN = testCase.NumMN; + nMA = testCase.NumMA; + nXA = testCase.NumXA; + nDW = testCase.NumDW; + testCase.NumTotalChans = nMN + nMA + nXA + nDW; + + testCase.MetaFilename = fullfile(testCase.TempDir, 'test_nidq.nidq.meta'); + testCase.BinFilename = fullfile(testCase.TempDir, 'test_nidq.nidq.bin'); + + % Generate test data + nSamples = testCase.NumSamples; + nTotal = testCase.NumTotalChans; + data = zeros(nSamples, nTotal, 'int16'); + for c = 1:nTotal + data(:, c) = int16((c-1)*10 + (1:nSamples))'; + end + + % Write binary file + fid = fopen(testCase.BinFilename, 'w', 'ieee-le'); + interleaved = reshape(data', 1, []); + fwrite(fid, interleaved, 'int16'); + fclose(fid); + + % Write meta file + fileSizeBytes = nSamples * nTotal * 2; + fileTimeSecs = nSamples / testCase.SR; + fid = fopen(testCase.MetaFilename, 'w'); + fprintf(fid, 'niSampRate=%g\n', testCase.SR); + fprintf(fid, 'nSavedChans=%d\n', nTotal); + fprintf(fid, 'snsMnMaXaDw=%d,%d,%d,%d\n', nMN, nMA, nXA, nDW); + fprintf(fid, 'snsSaveChanSubset=all\n'); + fprintf(fid, 'fileSizeBytes=%d\n', fileSizeBytes); + fprintf(fid, 'fileTimeSecs=%.6f\n', fileTimeSecs); + fprintf(fid, 'niAiRangeMax=%.1f\n', testCase.VMax); + fprintf(fid, 'niAiRangeMin=-%.1f\n', testCase.VMax); + fprintf(fid, 'niMaxInt=%d\n', testCase.MaxInt); + fprintf(fid, 'niMNGain=%d\n', testCase.MNGain); + fprintf(fid, 'niMAGain=%d\n', testCase.MAGain); + fclose(fid); + end + end + + methods (TestClassTeardown) + function teardownOnce(testCase) + if ~isempty(testCase.TempDir) && isfolder(testCase.TempDir) + try + rmdir(testCase.TempDir, 's'); + catch + end + end + end + end + + methods (Test) + + function testHeaderNidq(testCase) + %TESTHEADERNIDQ Verify header parses NIDQ fields correctly. + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + + testCase.verifyEqual(info.stream_type, 'nidq'); + testCase.verifyEqual(info.sample_rate, testCase.SR); + testCase.verifyEqual(info.n_saved_chans, testCase.NumTotalChans); + testCase.verifyEqual(info.n_mn_chans, testCase.NumMN); + testCase.verifyEqual(info.n_ma_chans, testCase.NumMA); + testCase.verifyEqual(info.n_xa_chans, testCase.NumXA); + testCase.verifyEqual(info.n_dw_chans, testCase.NumDW); + testCase.verifyEqual(info.voltage_range, [-testCase.VMax testCase.VMax]); + testCase.verifyEqual(info.max_int, testCase.MaxInt); + testCase.verifyEqual(info.ni_mn_gain, testCase.MNGain); + testCase.verifyEqual(info.ni_ma_gain, testCase.MAGain); + end + + function testSamples2VoltsNidqMN(testCase) + %TESTSAMPLES2VOLTSNIDQMN Verify scaling for MN channels. + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + raw = int16([1000; -1000; 0]); + volts = ndr.format.neuropixelsGLX.samples2volts(raw, info, 1); + + expected = double(raw) * testCase.VMax / (testCase.MaxInt * testCase.MNGain); + testCase.verifyEqual(volts, expected, 'AbsTol', 1e-15, ... + 'MN channel scaling incorrect.'); + end + + function testSamples2VoltsNidqMA(testCase) + %TESTSAMPLES2VOLTSNIDQMA Verify scaling for MA channels. + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + raw = int16([500; -500]); + ma_chan = testCase.NumMN + 1; % first MA channel + volts = ndr.format.neuropixelsGLX.samples2volts(raw, info, ma_chan); + + expected = double(raw) * testCase.VMax / (testCase.MaxInt * testCase.MAGain); + testCase.verifyEqual(volts, expected, 'AbsTol', 1e-15, ... + 'MA channel scaling incorrect.'); + end + + function testSamples2VoltsNidqMixed(testCase) + %TESTSAMPLES2VOLTSNIDQMIXED Verify mixed MN+MA channel scaling. + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + raw = int16([1000 500; -1000 -500]); + channels = [1, testCase.NumMN + 1]; % one MN, one MA + volts = ndr.format.neuropixelsGLX.samples2volts(raw, info, channels); + + scale_mn = testCase.VMax / (testCase.MaxInt * testCase.MNGain); + scale_ma = testCase.VMax / (testCase.MaxInt * testCase.MAGain); + expected = [double(raw(:,1)) * scale_mn, double(raw(:,2)) * scale_ma]; + testCase.verifyEqual(volts, expected, 'AbsTol', 1e-15, ... + 'Mixed MN+MA scaling incorrect.'); + end + + function testReadNidqScaled(testCase) + %TESTREADNIDQSCALED Verify read with default scaling for NIDQ. + [data, t, t0_t1] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.001); + + testCase.verifyClass(data, 'double', 'Scaled NIDQ data should be double.'); + testCase.verifyGreaterThan(numel(t), 0); + testCase.verifyEqual(t0_t1(1), 0, 'AbsTol', 1e-9); + end + + function testReadNidqUnscaled(testCase) + %TESTREADNIDQUNSCALED Verify read with scale=false for NIDQ. + [data, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.001, 'scale', false); + + testCase.verifyClass(data, 'int16', 'Unscaled NIDQ data should be int16.'); + end + + function testReadNidqScaledMatchesManual(testCase) + %TESTREADNIDQSCALEDMATCHESMANUAL Verify scaled read matches manual conversion. + channels = [1, testCase.NumMN + 1]; % MN and MA channel + + [data_scaled, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.01, 'channels', channels); + + [data_raw, ~, ~] = ndr.format.neuropixelsGLX.read(... + testCase.BinFilename, 0, 0.01, 'channels', channels, 'scale', false); + + info = ndr.format.neuropixelsGLX.header(testCase.MetaFilename); + expected = ndr.format.neuropixelsGLX.samples2volts(data_raw, info, channels); + + testCase.verifyEqual(data_scaled, expected, 'AbsTol', 1e-15, ... + 'Scaled read should match manual samples2volts for NIDQ.'); + end + + end + +end