Add Neuropixels SpikeGLX reader support to NDR

+ndr/+format/+neuropixelsGLX/NeuropixelsGLX_format.md

@@ -0,0 +1,151 @@
+# Neuropixels SpikeGLX Format — NDR Design Document
+
+## Overview
+
+This document describes the NDR support for Neuropixels data acquired with
+[SpikeGLX](https://billkarsh.github.io/SpikeGLX/), the open-source acquisition
+software for Neuropixels probes developed by Bill Karsh at the Allen Institute.
+
+## File Organization
+
+SpikeGLX saves data in a directory structure organized by run name, gate index,
+and trigger index. Each imec probe stream gets its own subdirectory:
+
+```
+runname_g0/                          # Gate 0
+  runname_g0_imec0/                  # Probe 0 subdirectory
+    runname_g0_t0.imec0.ap.bin       # AP-band binary data
+    runname_g0_t0.imec0.ap.meta      # AP-band metadata
+    runname_g0_t0.imec0.lf.bin       # LF-band binary data (Neuropixels 1.0)
+    runname_g0_t0.imec0.lf.meta      # LF-band metadata
+  runname_g0_imec1/                  # Probe 1 subdirectory (if multi-probe)
+    runname_g0_t0.imec1.ap.bin
+    runname_g0_t0.imec1.ap.meta
+    ...
+  runname_g0_t0.nidq.bin             # NI-DAQ auxiliary I/O (optional)
+  runname_g0_t0.nidq.meta
+```
+
+### Naming convention
+
+- **`_g<N>`** — Gate index (typically 0). Gates control when acquisition runs.
+- **`_t<N>`** — Trigger index (typically 0). Triggers segment data within a gate.
+- **`.imec<N>`** — Probe index (0-based). Each probe has independent clocks.
+- **`.ap`** — Action potential band (~30 kHz, high-pass filtered on-probe).
+- **`.lf`** — Local field potential band (~2.5 kHz, low-pass filtered). Only
+  available on Neuropixels 1.0 probes; Neuropixels 2.0 saves only AP data.
+
+## Binary File Format (.bin)
+
+The `.bin` files contain raw data as **interleaved int16** values with **no
+header**. Each time sample consists of one int16 value per saved channel,
+written consecutively:
+
+```
+[ch0_t0 ch1_t0 ch2_t0 ... chN_t0] [ch0_t1 ch1_t1 ... chN_t1] ...
+```
+
+- Data type: **int16** (signed 16-bit integer), little-endian
+- Byte order: **ieee-le** (Intel byte order)
+- No file header — the file begins immediately with sample data
+- Total samples = `fileSizeBytes / (nSavedChans * 2)`
+
+### Sync channel
+
+The last channel in each binary file is a **sync channel** (digital word).
+For a standard 384-channel Neuropixels 1.0 AP recording, the file contains
+385 channels: 384 neural + 1 sync.
+
+## Metadata File Format (.meta)
+
+Each `.bin` file has a companion `.meta` file containing acquisition
+parameters as `key=value` pairs, one per line.
+
+### Critical fields
+
+| Field               | Description                                        | Example          |
+|---------------------|----------------------------------------------------|------------------|
+| `imSampRate`        | Sampling rate in Hz                                | `30000`          |
+| `nSavedChans`       | Total channels per sample (neural + sync)          | `385`            |
+| `snsSaveChanSubset` | Which channels were saved (`all` or `0:383,768`)   | `all`            |
+| `snsApLfSy`         | Count of AP, LF, and sync channels                 | `384,0,1`        |
+| `fileSizeBytes`     | Size of the binary file in bytes                   | `231000000`      |
+| `fileTimeSecs`      | Duration of the recording in seconds               | `300.0`          |
+| `imAiRangeMax`      | Maximum voltage of ADC input range (V)             | `0.6`            |
+| `imAiRangeMin`      | Minimum voltage of ADC input range (V)             | `-0.6`           |
+| `imMaxInt`          | Maximum integer value for the ADC                  | `512`            |
+| `imDatPrb_type`     | Probe type (0=NP1.0, 21=NP2.0 single-shank, etc.) | `0`              |
+| `imDatPrb_sn`       | Probe serial number                                | `18005116102`    |
+| `imroTbl`           | Channel map with per-channel gains                 | `(0,384)(0 0...)` |
+| `typeThis`          | Stream type identifier                             | `imec`           |
+
+### Channel subsets
+
+When `snsSaveChanSubset` is not `all`, the user has selected a subset of
+channels to save. The format uses 0-based channel indices:
+
+- `0:383,768` — Save channels 0 through 383 and channel 768 (sync)
+- `0,2,4,6` — Save only even channels 0, 2, 4, 6
+- `0:95` — Save only the first 96 channels
+
+The NDR reader handles this by parsing the subset specification and mapping
+between file-order channel indices and original probe channel numbers.
+
+## Voltage Conversion
+
+Raw int16 values are converted to volts using the official SpikeGLX formula:
+
+```
+volts = int16_value * imAiRangeMax / imMaxInt / gain
+```
+
+where:
+- `imAiRangeMax` = maximum ADC voltage (typically 0.6 V)
+- `imMaxInt` = maximum ADC integer (512 for NP1.0, 8192 for NP2.0)
+- `gain` = per-channel gain from `imroTbl` (typically 500 for NP1.0 AP, 250 for LF; 80 for NP2.0)
+
+This yields ~2.34 uV/bit for NP1.0 AP and ~0.915 uV/bit for NP2.0 AP.
+
+## NDR Design Decisions
+
+### Scope of the reader
+
+The `ndr.reader.neuropixelsGLX` reader handles **one probe's AP-band stream**
+per instance. This means:
+
+- Each `.ap.bin`/`.ap.meta` file pair is one epoch stream
+- Multi-probe recordings use separate reader instances (one per `imec<N>`)
+- LF-band and NI-DAQ streams are not handled by this reader (future work)
+
+### Channel naming
+
+Neural channels are named `ai1` through `aiN` where N is the number of neural
+channels saved. The sync channel is exposed as `di1` (digital input). A time
+channel `t1` is always present.
+
+### Data type
+
+The reader returns raw int16 data from `readchannels_epochsamples` to preserve
+the native precision and enable efficient storage. Use
+`ndr.format.neuropixelsGLX.samples2volts` for voltage conversion.
+
+### Epoch structure
+
+Each `.ap.bin` file represents a single epoch. The `filenamefromepochfiles`
+method identifies the `.ap.meta` file from the epoch stream list, and the
+corresponding `.ap.bin` file is derived from it.
+
+## Format Functions
+
+| Function        | Purpose                                               |
+|-----------------|-------------------------------------------------------|
+| `readmeta`      | Parse a `.meta` file into a key-value structure       |
+| `header`        | Extract standardized recording parameters from `.meta`|
+| `read`          | Read binary data with time/channel subsetting         |
+| `samples2volts` | Convert raw int16 to voltage using gain parameters    |
+
+## References
+
+- [SpikeGLX Documentation](https://billkarsh.github.io/SpikeGLX/)
+- [SpikeGLX Data Format](https://billkarsh.github.io/SpikeGLX/Support/SpikeGLX_Datafile_Tools.html)
+- [Neuropixels](https://www.neuropixels.org/)

+ndr/+format/+neuropixelsGLX/header.m

@@ -0,0 +1,184 @@
+function info = header(metafilename)
+%HEADER Parse a SpikeGLX .meta file into a standardized header structure.
+%
+%   INFO = ndr.format.neuropixelsGLX.header(METAFILENAME)
+%
+%   Reads a SpikeGLX .meta file and extracts key recording parameters into
+%   a standardized structure with numeric fields suitable for data reading.
+%
+%   This function handles both AP-band and LF-band .meta files, and
+%   correctly parses the channel subset specification (snsSaveChanSubset)
+%   to determine exactly which channels were saved.
+%
+%   Inputs:
+%       METAFILENAME - Full path to the .meta file (char row vector).
+%
+%   Outputs:
+%       INFO - A scalar structure with the following fields:
+%           sample_rate      : Sampling rate in Hz (double).
+%           n_saved_chans    : Total number of saved channels including
+%                              the sync channel (double).
+%           n_neural_chans   : Number of neural (non-sync) channels saved
+%                              (double). For AP band this is the number of
+%                              AP channels; for LF band, the LF channels.
+%           n_sync_chans     : Number of sync channels (0 or 1) (double).
+%           saved_chan_list   : 1-based vector of saved channel indices
+%                              (double row vector). If all channels are
+%                              saved, this is 1:n_saved_chans.
+%           voltage_range    : Peak-to-peak voltage range [Vmin Vmax] in
+%                              volts (1x2 double).
+%           max_int          : Maximum integer value for the ADC (double).
+%           bits_per_sample  : Bits per sample value (always 16) (double).
+%           file_size_bytes  : Total file size in bytes (double).
+%           file_time_secs   : Recording duration in seconds (double).
+%           probe_type       : Probe type identifier (char).
+%           probe_sn         : Probe serial number (char).
+%           stream_type      : 'ap' or 'lf' (char), determined from the
+%                              meta file name.
+%           meta             : The raw meta structure from readmeta (struct).
+%
+%   Example:
+%       info = ndr.format.neuropixelsGLX.header('/data/run_g0_t0.imec0.ap.meta');
+%       fprintf('Sample rate: %g Hz\n', info.sample_rate);
+%       fprintf('Neural channels: %d\n', info.n_neural_chans);
+%       fprintf('Duration: %.2f s\n', info.file_time_secs);
+%
+%   See also: ndr.format.neuropixelsGLX.readmeta, ndr.format.neuropixelsGLX.read
+
+    arguments
+        metafilename (1,:) char {mustBeFile}
+    end
+
+    meta = ndr.format.neuropixelsGLX.readmeta(metafilename);
+
+    info = struct();
+    info.meta = meta;
+
+    % Sample rate
+    if isfield(meta, 'imSampRate')
+        info.sample_rate = str2double(meta.imSampRate);
+    elseif isfield(meta, 'niSampRate')
+        info.sample_rate = str2double(meta.niSampRate);
+    else
+        error('ndr:format:neuropixelsGLX:header:NoSampleRate', ...
+            'Could not find sample rate in meta file.');
+    end
+
+    % Number of saved channels
+    info.n_saved_chans = str2double(meta.nSavedChans);
+
+    % Parse snsApLfSy or snsMnMaXaDw to determine neural vs sync channels
+    if isfield(meta, 'snsApLfSy')
+        % imec stream: AP,LF,SY counts
+        counts = sscanf(meta.snsApLfSy, '%d,%d,%d');
+        % counts(1) = AP chans, counts(2) = LF chans, counts(3) = SY chans
+        info.n_sync_chans = counts(3);
+        % Determine if this is AP or LF from filename
+        [~, name, ~] = fileparts(metafilename);
+        if contains(name, '.lf')
+            info.stream_type = 'lf';
+            info.n_neural_chans = counts(2);
+        else
+            info.stream_type = 'ap';
+            info.n_neural_chans = counts(1);
+        end
+    elseif isfield(meta, 'snsMnMaXaDw')
+        % NI-DAQ stream
+        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;
+    else
+        % Fallback
+        info.stream_type = 'unknown';
+        info.n_neural_chans = info.n_saved_chans - 1;
+        info.n_sync_chans = 1;
+    end
+
+    % Parse saved channel subset
+    if isfield(meta, 'snsSaveChanSubset')
+        info.saved_chan_list = parse_channel_subset(meta.snsSaveChanSubset, info.n_saved_chans);
+    else
+        info.saved_chan_list = 1:info.n_saved_chans;
+    end
+
+    % Voltage range
+    if isfield(meta, 'imAiRangeMax')
+        vmax = str2double(meta.imAiRangeMax);
+        vmin = str2double(meta.imAiRangeMin);
+        info.voltage_range = [vmin vmax];
+    else
+        info.voltage_range = [-0.6 0.6]; % Neuropixels 1.0 default
+    end
+
+    % Max integer value
+    if isfield(meta, 'imMaxInt')
+        info.max_int = str2double(meta.imMaxInt);
+    else
+        info.max_int = 512; % Neuropixels 1.0 default
+    end
+
+    % Bits per sample
+    info.bits_per_sample = 16;
+
+    % File size and duration
+    if isfield(meta, 'fileSizeBytes')
+        info.file_size_bytes = str2double(meta.fileSizeBytes);
+    else
+        info.file_size_bytes = 0;
+    end
+
+    if isfield(meta, 'fileTimeSecs')
+        info.file_time_secs = str2double(meta.fileTimeSecs);
+    else
+        info.file_time_secs = 0;
+    end
+
+    % Probe information
+    if isfield(meta, 'imDatPrb_type')
+        info.probe_type = meta.imDatPrb_type;
+    else
+        info.probe_type = '';
+    end
+
+    if isfield(meta, 'imDatPrb_sn')
+        info.probe_sn = meta.imDatPrb_sn;
+    else
+        info.probe_sn = '';
+    end
+
+end
+
+
+function chan_list = parse_channel_subset(subset_str, n_saved_chans)
+%PARSE_CHANNEL_SUBSET Parse the snsSaveChanSubset field.
+%
+%   CHAN_LIST = PARSE_CHANNEL_SUBSET(SUBSET_STR, N_SAVED_CHANS)
+%
+%   The snsSaveChanSubset field can be:
+%     - 'all'         : All channels saved
+%     - '0:5,8,10:12' : Specific channels (0-based ranges/singles)
+%
+%   Returns a 1-based channel index vector.
+
+    if strcmpi(strtrim(subset_str), 'all')
+        chan_list = 1:n_saved_chans;
+        return;
+    end
+
+    chan_list = [];
+    parts = strsplit(strtrim(subset_str), ',');
+    for i = 1:numel(parts)
+        part = strtrim(parts{i});
+        if contains(part, ':')
+            range_vals = sscanf(part, '%d:%d');
+            chan_list = [chan_list, range_vals(1):range_vals(2)]; %#ok<AGROW>
+        else
+            chan_list = [chan_list, str2double(part)]; %#ok<AGROW>
+        end
+    end
+
+    % Convert from 0-based to 1-based
+    chan_list = chan_list + 1;
+
+end