multi-probe handling

calc_CellMetrics/loadOpenEphysDigital.m

@@ -1,57 +1,85 @@
-function openephysDig = loadOpenEphysDigital(session)
+function openephysDig = loadOpenEphysDigital(session, varargin)
 % function to load digital inputs from Open Ephys and align the data properly
 %
 % INPUT
 % session: session struct
+% Optional inputs:
+%   probeLetter: 'A' or 'B' for specific probe TTL data. If empty, tries ProbeA first, then legacy format
 %
 % OUTPUTS
 % intanDig.on: on-state changes for each channel [in seconds]
 % intanDig.off: off-state changes for each channel [in seconds]
 %
-% By Peter Petersen
+% By Peter Petersen, Mingze Dou
 % [email protected]
 
-
-% Load TTL data
-% Each path must contain either:
-% Legacy format:
-%   1. timestamps.npy
-%   2. channel_states.npy
-%   3. channels.npy
-%   4. full_words.npy
-% New format:
-%   1. timestamps.npy
-%   2. states.npy
-%   3. sample_numbers.npy
-%   4. full_words.npy
+p = inputParser;
+addParameter(p,'probeLetter','',@ischar);
+parse(p,varargin{:});
+parameters = p.Results;
 
 TTL_paths = {};
 epochs_startTime = [];
 ephys_t0 = [];
 
 % Determine file format and set paths
 for i = 1:numel(session.epochs)
-    % Try new format path first
-    newFormatPath = fullfile(session.epochs{i}.name,'events','Neuropix-PXI-100.ProbeA','TTL');
-    legacyFormatPath = fullfile(session.epochs{i}.name,'events','Neuropix-PXI-100.0','TTL_1');
-
-    % Check which path exists and set accordingly
-    if exist(fullfile(session.general.basePath, newFormatPath), 'dir')
-        TTL_paths{i} = newFormatPath;
-        timestampPath = fullfile(session.general.basePath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.ProbeA','timestamps.npy');
-        isNewFormat = true;  %format flag
+
+    if ~isempty(parameters.probeLetter)
+        % Specific probe requested
+        TTL_path = fullfile(session.epochs{i}.name,'events',['Neuropix-PXI-100.Probe' parameters.probeLetter],'TTL');
+        disp(['Checking TTL path for epoch ' num2str(i) ': ' fullfile(session.general.basePath, TTL_path)])
+        if exist(fullfile(session.general.basePath, TTL_path), 'dir')
+            TTL_paths{i} = TTL_path;
+            isNewFormat = true;
+            disp('Found TTL directory')
+        else
+            error(['TTL path not found for Probe' parameters.probeLetter ' in epoch ' num2str(i)]);
+        end
     else
-        TTL_paths{i} = legacyFormatPath;
-        timestampPath = fullfile(session.general.basePath,session.epochs{i}.name,'continuous','Neuropix-PXI-100.1','timestamps.npy');
-        isNewFormat = false;  
+        % Try ProbeA format first, then legacy format
+        probeA_path = fullfile(session.epochs{i}.name,'events','Neuropix-PXI-100.ProbeA','TTL');
+        legacy_path = fullfile(session.epochs{i}.name,'events','Neuropix-PXI-100.0','TTL_1');
+
+        if exist(fullfile(session.general.basePath, probeA_path), 'dir')
+            TTL_paths{i} = probeA_path;
+            isNewFormat = true;
+        elseif exist(fullfile(session.general.basePath, legacy_path), 'dir')
+            TTL_paths{i} = legacy_path;
+            isNewFormat = false;
+        else
+            error(['No valid TTL path found for epoch ' num2str(i)]);
+        end
     end
-
+
+    % Set timestamp path based on format
+    if isNewFormat
+        probeStr = parameters.probeLetter;
+        if isempty(probeStr)
+            probeStr = 'A';
+        end
+        timestampPath = fullfile(session.general.basePath, session.epochs{i}.name, 'continuous', ...
+            ['Neuropix-PXI-100.Probe' probeStr], 'timestamps.npy');
+    else
+        timestampPath = fullfile(session.general.basePath, session.epochs{i}.name, 'continuous', ...
+            'Neuropix-PXI-100.1', 'timestamps.npy');
+    end
+
+    disp(['Checking timestamp path for epoch ' num2str(i) ': ' timestampPath])
+    if exist(timestampPath, 'file')
+        temp = readNPY(timestampPath);
+        disp(['Found ' num2str(length(temp)) ' timestamps, first value: ' num2str(temp(1))])
+    else
+        disp('Timestamp file not found')
+    end
+
     epochs_startTime(i) = session.epochs{i}.startTime;
-    temp = readNPY(timestampPath);
     if isNewFormat
-        ephys_t0(i) = double(temp(1));  % timestamps in second according to new format of open ephys
+        ephys_t0(i) = double(temp(1));
+        disp(['ephys_t0 for epoch ' num2str(i) ': ' num2str(ephys_t0(i))])
     else
-        ephys_t0(i) = double(temp(1))/session.extracellular.sr;  
+        ephys_t0(i) = double(temp(1))/session.extracellular.sr;
+        disp(['ephys_t0 for epoch ' num2str(i) ': ' num2str(ephys_t0(i))])
     end
 end
 
@@ -60,23 +88,31 @@
 
 % Load and process first epoch data
 basePath = fullfile(session.general.basePath, TTL_paths{1});
-timestamps = readNPY(fullfile(basePath,'timestamps.npy'));
+ttlPath = fullfile(basePath,'timestamps.npy');
+disp(['Checking TTL timestamps path: ' ttlPath])
+if exist(ttlPath, 'file')
+    timestamps = readNPY(ttlPath);
+    disp(['Found ' num2str(length(timestamps)) ' TTL timestamps'])
+else
+    disp('TTL timestamp file not found')
+end
 if isNewFormat
-    openephysDig.timestamps = epochs_startTime(1) + double(timestamps) - ephys_t0(1);  % timestamps in second
+    openephysDig.timestamps = epochs_startTime(1) + double(timestamps) - ephys_t0(1); % timestamps in second
 else
-    openephysDig.timestamps = epochs_startTime(1) + double(timestamps)/session.extracellular.sr - ephys_t0(1);  % Legacy format
+    openephysDig.timestamps = epochs_startTime(1) + double(timestamps)/session.extracellular.sr - ephys_t0(1); % legacy format
 end
 
 % Check which format to use for first epoch
 if exist(fullfile(basePath,'states.npy'), 'file')
+    % New format files
     openephysDig.channel_states = readNPY(fullfile(basePath,'states.npy'));
     openephysDig.channels = readNPY(fullfile(basePath,'sample_numbers.npy'));
-    openephysDig.full_words = readNPY(fullfile(basePath,'full_words.npy'));
 else
+    % Legacy format files
     openephysDig.channel_states = readNPY(fullfile(basePath,'channel_states.npy'));
     openephysDig.channels = readNPY(fullfile(basePath,'channels.npy'));
-    openephysDig.full_words = readNPY(fullfile(basePath,'full_words.npy'));
 end
+openephysDig.full_words = readNPY(fullfile(basePath,'full_words.npy'));
 
 openephysDig.on{1} = double(openephysDig.timestamps(openephysDig.channel_states == 1));
 openephysDig.off{1} = double(openephysDig.timestamps(openephysDig.channel_states == -1));
@@ -92,18 +128,20 @@
         timestamps = readNPY(fullfile(basePath,'timestamps.npy'));
 
         if isNewFormat
-            timestamps = epochs_startTime(i) + double(timestamps) - ephys_t0(i);  % timestamps in second
+            timestamps = epochs_startTime(i) + double(timestamps) - ephys_t0(i); % timestamps in second
         else
-            timestamps = epochs_startTime(i) + double(timestamps)/session.extracellular.sr - ephys_t0(i);  % Legacy format
+            timestamps = epochs_startTime(i) + double(timestamps)/session.extracellular.sr - ephys_t0(i); % legacy format
         end
         openephysDig.timestamps = [openephysDig.timestamps; timestamps];
 
         % Load states and channels based on format
         if exist(fullfile(basePath,'states.npy'), 'file')
+            % New format files
             channel_states = readNPY(fullfile(basePath,'states.npy'));
             openephysDig.channel_states = [openephysDig.channel_states; channel_states];
             openephysDig.channels = [openephysDig.channels; readNPY(fullfile(basePath,'sample_numbers.npy'))];
         else
+            % Legacy format files
             channel_states = readNPY(fullfile(basePath,'channel_states.npy'));
             openephysDig.channel_states = [openephysDig.channel_states; channel_states];
             openephysDig.channels = [openephysDig.channels; readNPY(fullfile(basePath,'channels.npy'))];
@@ -134,4 +172,4 @@
 
 % Saving data
 saveStruct(openephysDig,'digitalseries','session',session);
-end
+end

calc_CellMetrics/loadOpenEphysDigitalNidaq.m

@@ -0,0 +1,189 @@
+function openephysDig = loadOpenEphysDigitalNidaq(session, varargin)
+% loadOpenEphysDigitalNidaq - Load and process digital TTL signals from OpenEphys NI-DAQ Board
+%
+% This function processes TTL signals recorded by the NI-DAQ board in OpenEphys, handling
+% multiple recording epochs and aligning timestamps to the probe's timeline.
+%
+% INPUTS:
+%   session        - Session struct containing recording metadata and paths
+%
+% Optional Name-Value Pairs:
+%   'channelNum'   - Integer. Channel number on NI-DAQ board to process
+%                    0: Clock signal (1 Hz)
+%                    1: Camera trigger (120 Hz)
+%                    Default: [] (processes all channels)
+%   'probeLetter'  - Character. Specifies which probe's timeline to use for alignment
+%                    'A': Use ProbeA timestamps (default)
+%                    'B': Use ProbeB timestamps
+%
+% OUTPUTS:
+%   openephysDig   - Structure containing processed TTL data with fields:
+%     .timestamps  - [Nx1] Vector of all event timestamps (in seconds)
+%     .states      - [Nx1] Vector of binary states (0/1) for each timestamp
+%     .epochNum    - [Nx1] Vector indicating the epoch number for each timestamp
+%     .on          - {1x1} Cell containing timestamps of rising edges (0→1)
+%     .off         - {1x1} Cell containing timestamps of falling edges (1→0)
+%     .diagnostics - Structure containing timing analysis and data quality info
+%                   including timestamp duplicates and state conflicts
+%
+% NOTES:
+%   - Rising and falling edges are paired based on channel-specific transitions
+%   - Timestamps are aligned to probe recording start time
+%   - Function warns if rising edges lack corresponding falling edges
+%   - Empty epochs or missing files are handled gracefully with warnings
+%   - Duplicate timestamps are analyzed with their states for conflicts
+%
+% EXAMPLE:
+%   % Load camera trigger signals from channel 1
+%   ttl = loadOpenEphysDigitalNidaq(session, 'channelNum', 1)
+%
+% See also: readNPY, saveStruct
+%
+% By Mingze Dou
+
+p = inputParser;
+addParameter(p,'channelNum', [], @isnumeric);
+addParameter(p,'probeLetter', 'A', @ischar);
+parse(p,varargin{:});
+parameters = p.Results;
+
+% Check for existing digitalseries file
+existingFile = fullfile(session.general.basePath, [session.general.name '.openephysDig.digitalseries.mat']);
+if exist(existingFile, 'file')
+    load(existingFile, 'openephysDig');
+    return
+end
+
+% Initialize
+ttlPaths = {};
+epochsStartTime = [];
+ephysT0 = [];
+validEpochs = [];
+openephysDig.on = {[]};
+openephysDig.off = {[]};
+openephysDig.timestamps = [];
+openephysDig.states = [];
+openephysDig.epochNum = [];
+openephysDig.diagnostics = struct();
+
+% Find valid epochs
+for i = 1:numel(session.epochs)
+    ttlPath = fullfile(session.epochs{i}.name, 'events', 'NI-DAQmx-116.PXIe-6341', 'TTL');
+    fullWordsPath = fullfile(session.general.basePath, ttlPath, 'full_words.npy');
+
+    if exist(fullWordsPath, 'file') && (~isempty(parameters.channelNum) && ...
+            any(bitand(readNPY(fullWordsPath), 2^parameters.channelNum)) || isempty(parameters.channelNum))
+        validEpochs = [validEpochs i];
+        ttlPaths{end+1} = ttlPath;
+
+        % Get probe timestamps
+        probeStr = parameters.probeLetter;
+        if ~isempty(parameters.probeLetter), probeStr = parameters.probeLetter; end
+
+        timestampPath = fullfile(session.general.basePath, session.epochs{i}.name, 'continuous', ...
+            ['Neuropix-PXI-103.Probe' probeStr], 'timestamps.npy');
+
+        if ~exist(timestampPath, 'file')
+            error(['Timestamp file not found for Probe' probeStr ' in epoch ' num2str(i)]);
+        end
+
+        timestampData = readNPY(timestampPath);
+        epochsStartTime(end+1) = session.epochs{i}.startTime;
+        ephysT0(end+1) = double(timestampData(1));
+    end
+end
+
+if isempty(validEpochs)
+    error('No valid epochs found with the specified channel');
+end
+
+% Process epochs
+for i = 1:length(validEpochs)
+    currentEpoch = validEpochs(i);
+    basePath = fullfile(session.general.basePath, ttlPaths{i});
+    timestamps = readNPY(fullfile(basePath,'timestamps.npy'));
+    fullWords = readNPY(fullfile(basePath,'full_words.npy'));
+    if isempty(timestamps) || isempty(fullWords)
+        warning('Empty data found in epoch %d', currentEpoch);
+        continue;
+    end
+
+    % Align timestamps and get states
+    alignedTimestamps = epochsStartTime(i) + double(timestamps) - ephysT0(i);
+    channel_states = bitget(fullWords, parameters.channelNum + 1);
+
+    % Find unique timestamps and duplicates
+    [uniqueTimestamps, uniqueIdx] = unique(alignedTimestamps);
+    uniqueStates = channel_states(uniqueIdx);
+
+    % Analyze duplicates and their states
+    [~, ~, dupGroups] = unique(alignedTimestamps);  % Get grouping for all timestamps
+    duplicateInfo = struct('timestamps', [], 'states', [], 'hasConflict', []);
+
+    % Find timestamps that appear more than once
+    [dupCounts, dupValues] = hist(alignedTimestamps, unique(alignedTimestamps));
+    duplicateTimestamps = dupValues(dupCounts > 1);
+
+    if ~isempty(duplicateTimestamps)
+        for j = 1:length(duplicateTimestamps)
+            dupTime = duplicateTimestamps(j);
+            dupIndices = find(alignedTimestamps == dupTime);
+            dupStates = channel_states(dupIndices);
+
+            % Store duplicate information
+            duplicateInfo.timestamps(j) = dupTime;
+            duplicateInfo.states{j} = dupStates;
+            duplicateInfo.hasConflict(j) = length(unique(dupStates)) > 1;
+        end
+
+        % Store in diagnostics
+        openephysDig.diagnostics.epoch(i).duplicates = duplicateInfo;
+    end
+
+    % Initialize arrays for this epoch
+    rising_edges = [];
+    falling_edges = [];
+
+    % Process the state sequence
+    state_idx = 1;
+    while state_idx < length(uniqueStates)
+        % Find next rising edge
+        while state_idx < length(uniqueStates) && uniqueStates(state_idx) == 0
+            state_idx = state_idx + 1;
+        end
+
+        if state_idx < length(uniqueStates)
+            % Found a rising edge
+            rising_edge_idx = state_idx;
+
+            % Look for the next falling edge for this specific channel
+            state_idx = state_idx + 1;
+            while state_idx < length(uniqueStates) && uniqueStates(state_idx) == 1
+                state_idx = state_idx + 1;
+            end
+
+            if state_idx < length(uniqueStates)
+                % Found the corresponding falling edge
+                falling_edge_idx = state_idx;
+
+                % Store the edges
+                rising_edges = [rising_edges; rising_edge_idx];
+                falling_edges = [falling_edges; falling_edge_idx];
+            else
+                warning('Found rising edge without corresponding falling edge at timestamp %f in epoch %d', ...
+                    uniqueTimestamps(rising_edge_idx), currentEpoch);
+            end
+        end
+    end
+
+    % Store results using the correct timestamps
+    openephysDig.on{1} = [openephysDig.on{1}; uniqueTimestamps(rising_edges)];
+    openephysDig.off{1} = [openephysDig.off{1}; uniqueTimestamps(falling_edges)];
+    openephysDig.timestamps = [openephysDig.timestamps; uniqueTimestamps];
+    openephysDig.states = [openephysDig.states; uniqueStates];
+    openephysDig.epochNum = [openephysDig.epochNum; repmat(currentEpoch, length(uniqueTimestamps), 1)];
+end
+
+% Save
+saveStruct(openephysDig,'digitalseries','session',session);
+end

calc_CellMetrics/loadOpenEphysSettingsFile.m

@@ -5,13 +5,13 @@
 % https://open-ephys.github.io/gui-docs/User-Manual/Recording-data/Binary-format.html
 
 p = inputParser;
-addParameter(p,'probe_letter','A',@ischar);
+addParameter(p,'probeLetter','A',@ischar);
 
 % Parsing inputs
 parse(p,varargin{:})
 parameters = p.Results;
 
-switch parameters.probe_letter
+switch parameters.probeLetter
     case 'A'
         stream_id = 1;
     case 'B'