diff --git a/.gitignore b/.gitignore index 0c0c833..2aef8d1 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ *.asv *.pyc +*.mexa64 *.mexw64 *.mexmaci64 ._* @@ -7,8 +8,4 @@ behavior/dlc/Run_startdlc.bat preProcessing/FilterM/FilterX.exp preProcessing/FilterM/FilterX.lib preProcessing/FilterM/FilterX.mexw64 -preProcessing/FilterM/FilterX.exp -preProcessing/FilterM/FilterX.mexw64 -preProcessing/FilterM/FilterX.lib -preProcessing/FilterM/FilterX.exp spikeSorting/KiloSort1/CUDA/CUDA.xml diff --git a/behavior/general_behavior_file.m b/behavior/general_behavior_file.m index 14e61dc..013fbbc 100644 --- a/behavior/general_behavior_file.m +++ b/behavior/general_behavior_file.m @@ -20,6 +20,12 @@ function general_behavior_file(varargin) % optitrack .csv file % *_Behavior*.mat file from crcns_pfc-2_data % +% OPTIONAL PARAMETERS: +% 'video_overlay' - logical, enable video overlay in manual_trackerjumps (default: false) +% 'clean_tracker_jumps' - logical, enable manual tracker jump cleaning GUI (default: true) +% 'convert_xy_to_cm' - logical, convert coordinates from pixels to cm (default: true) +% 'maze_sizes' - vector, maze dimensions in cm for coordinate conversion (required if convert_xy_to_cm=true) +% % TODO: % -make so you can choose (w/ varargin) which method to use (some sessions meet several) % This will require making each method into a sub/local function @@ -41,6 +47,7 @@ function general_behavior_file(varargin) addParameter(p, 'clean_tracker_jumps', true, @islogical); % option to manually clean tracker jumps addParameter(p, 'convert_xy_to_cm', true, @islogical); % option to convert xy to cm (best if used with clean_tracker_jumps) addParameter(p, 'maze_sizes', [], @isvector); % list of maze sizes (x-dim, in cm) per non-sleep epoch (if same maze & cam pos over epochs use single number) +addParameter(p, 'video_overlay', false, @islogical); % option to overlay tracking on video frame in manual_trackerjumps parse(p, varargin{:}); basepaths = p.Results.basepath; @@ -54,12 +61,12 @@ function general_behavior_file(varargin) clean_tracker_jumps = p.Results.clean_tracker_jumps; convert_xy_to_cm = p.Results.convert_xy_to_cm; maze_sizes = p.Results.maze_sizes; +video_overlay = p.Results.video_overlay; if convert_xy_to_cm && isempty(maze_sizes) error('Cannot run with "convert_xy_to_cm" option without providing maze_sizes.'); end - if ~iscell(basepaths) basepaths = {basepaths}; end @@ -152,6 +159,11 @@ function general_behavior_file(varargin) epoch_df = epoch_df(epoch_df.environment ~= "sleep", :); start = epoch_df.startTime; stop = epoch_df.stopTime; + % make full paths using epoch_df.basepath and epoch_df.name + session_paths = cellfun(@(x,y) fullfile(x,y), ... + epoch_df.basepath, ... + epoch_df.name, ... + 'UniformOutput', false); good_idx = manual_trackerjumps(behavior.timestamps, ... behavior.position.x, ... @@ -159,7 +171,9 @@ function general_behavior_file(varargin) start, ... stop, ... basepath, ... - 'darkmode', true); + session_paths, ... + 'darkmode', true, ... + 'video_overlay', p.Results.video_overlay); behavior.position.x(~good_idx) = NaN; behavior.position.y(~good_idx) = NaN; @@ -598,9 +612,9 @@ function general_behavior_file(varargin) source = 'basename.posTrials.mat'; % posTrials is sometimes moved -elseif exist([basepath, filesep, ['oldfiles\posTrials.mat']], 'file') +elseif exist([basepath, filesep, ['oldfiles', filesep, 'posTrials.mat']], 'file') disp('detected posTrials.mat') - load([basepath, filesep, ['oldfiles\posTrials.mat']], 'posTrials'); + load([basepath, filesep, ['oldfiles', filesep, 'posTrials.mat']], 'posTrials'); positions = [posTrials{1}; posTrials{2}]; [~, idx] = sort(positions(:, 1)); positions = positions(idx, :); diff --git a/behavior/manual_trackerjumps.m b/behavior/manual_trackerjumps.m index 66da971..0794b31 100644 --- a/behavior/manual_trackerjumps.m +++ b/behavior/manual_trackerjumps.m @@ -1,4 +1,4 @@ -function good_idx = manual_trackerjumps(ts,x,y,StartofRec,EndofRec,basepath,varargin) +function good_idx = manual_trackerjumps(ts,x,y,StartofRec,EndofRec,basepath,session_paths, varargin) % Manually cut out xy coordinates that are outside the bounds of your maze. % @@ -11,6 +11,13 @@ % y % StartofRec: ts indicating the start points of your event % EndofRec: ts indicating the end points of your event +% basepath: path to session folder +% session_paths: array of full paths for each subsession (for video overlay) +% +% Optional Parameters: +% 'video_overlay': logical, overlay tracking on video frame (default: false) +% 'darkmode': logical, use dark theme for plots (default: true) +% 'save_boundary_file': logical, save boundary selection (default: false) % % Copyright (C) 2018 Ryan E Harvey @@ -22,6 +29,7 @@ addParameter(p,'alpha',.2,@isnumeric); addParameter(p,'add_scatter',true,@islogical); addParameter(p,'save_boundary_file',false,@islogical); +addParameter(p,'video_overlay',false,@islogical); parse(p,varargin{:}); darkmode = p.Results.darkmode; @@ -30,6 +38,8 @@ alpha = p.Results.alpha; add_scatter = p.Results.add_scatter; save_boundary_file = p.Results.save_boundary_file; +video_overlay = p.Results.video_overlay; +video_frame_time = []; savets=[]; for i=1:length(StartofRec) @@ -37,9 +47,73 @@ xtemp=x(ts>=StartofRec(i) & ts<=EndofRec(i)); ytemp=y(ts>=StartofRec(i) & ts<=EndofRec(i)); tstemp=ts(ts>=StartofRec(i) & ts<=EndofRec(i)); + + % Check if we have valid data + if isempty(xtemp) || isempty(ytemp) + warning('Event %d has no data points', i); + continue; + end + + % Remove NaN values for visualization + valid_idx = ~isnan(xtemp) & ~isnan(ytemp); + if sum(valid_idx) == 0 + warning('Event %d has only NaN values', i); + continue; + end + + fprintf('Event %d: %d total points, %d valid points\n', i, length(xtemp), sum(valid_idx)); + + % Auto-detect video path if video_overlay is enabled but path is empty + current_video_path = ''; + basename = basenameFromBasepath(session_paths{i}); + if p.Results.video_overlay && isempty(current_video_path) + fprintf('Video overlay enabled, searching for video file...\n'); + + % Look for video files in MergePoints folders first + if exist(fullfile(session_paths{i}, [basename, '.MergePoints.events.mat']), 'file') + load(fullfile(session_paths{i}, [basename, '.MergePoints.events.mat']), 'MergePoints'); + + % Check subsession folders for MP4/AVI files + for k = 1:length(MergePoints.foldernames) + mp4_files = dir(fullfile(session_paths{i}, MergePoints.foldernames{k}, '*.mp4')); + if ~isempty(mp4_files) + current_video_path = fullfile(session_paths{i}, MergePoints.foldernames{k}, mp4_files(1).name); + fprintf('Found video: %s\n', current_video_path); + break; + end + + avi_files = dir(fullfile(session_paths{i}, MergePoints.foldernames{k}, '*.avi')); + if ~isempty(avi_files) + current_video_path = fullfile(session_paths{i}, MergePoints.foldernames{k}, avi_files(1).name); + fprintf('Found video: %s\n', current_video_path); + break; + end + end + end + + % If still not found, check main session folder + if isempty(current_video_path) + mp4_files = dir(fullfile(session_paths{i}, '*.mp4')); + if ~isempty(mp4_files) + current_video_path = fullfile(session_paths{i}, mp4_files(1).name); + fprintf('Found video: %s\n', current_video_path); + else + avi_files = dir(fullfile(session_paths{i}, '*.avi')); + if ~isempty(avi_files) + current_video_path = fullfile(session_paths{i}, avi_files(1).name); + fprintf('Found video: %s\n', current_video_path); + end + end + end + + if isempty(current_video_path) + fprintf('Warning: No video file found for overlay. Proceeding without video.\n'); + end + end + % use the gui to cut out points [~,~,in]=restrictMovement(xtemp,ytemp,restic_dir,darkmode,axis_equal,... - alpha,add_scatter); + alpha,add_scatter,video_overlay,current_video_path,video_frame_time,tstemp); % save the ts where the tracker error exists savets=[savets,tstemp(in)]; end @@ -55,7 +129,7 @@ end end -function [x,y,in]=restrictMovement(x,y,direction,darkmode,axis_equal,alpha,add_scatter) +function [x,y,in]=restrictMovement(x,y,direction,darkmode,axis_equal,alpha,add_scatter,video_overlay,video_path,video_frame_time,timestamps) % restrictMovement allows you to draw a line around xy coordinates in order % to eliminate certain points you don't want...ie when the rat jumps out of % maze or tracker errors. @@ -67,6 +141,10 @@ % % Input x,y: coordinates % direction: 1 (default) to remove outside points; 0 to remove inside points +% video_overlay: logical to overlay video frame as background +% video_path: path to video file +% video_frame_time: specific time to extract frame, or empty for middle frame +% timestamps: timestamps for current event % % % Output x,y: retained coordinates @@ -79,32 +157,227 @@ if nargin<3 direction=1; end +if nargin<8 + video_overlay=false; +end +if nargin<9 + video_path=''; +end +if nargin<10 + video_frame_time=[]; +end +if nargin<11 + timestamps=[]; +end -% set up figure -if darkmode - fig=figure;plot(x,y,'Color',[1,1,1,alpha]);hold on - if add_scatter - scatter(x,y,3,'w','filled');hold on +% Validate input data +if isempty(x) || isempty(y) + error('Empty x or y coordinates'); +end + +% Remove NaN values for visualization +valid_idx = ~isnan(x) & ~isnan(y); +if sum(valid_idx) == 0 + error('All coordinates are NaN'); +end + +x_plot = x(valid_idx); +y_plot = y(valid_idx); + +fprintf('Plotting %d valid points (out of %d total)\n', sum(valid_idx), length(x)); +fprintf('X range: [%.2f, %.2f], Y range: [%.2f, %.2f]\n', ... + min(x_plot), max(x_plot), min(y_plot), max(y_plot)); + +% Load video frame if requested +video_frame = []; +if video_overlay && ~isempty(video_path) && exist(video_path, 'file') + try + fprintf('Loading video frame from: %s\n', video_path); + + % Determine which frame to extract + if isempty(video_frame_time) && ~isempty(timestamps) + % Use middle timestamp of current event + video_frame_time = median(timestamps) - timestamps(1); + fprintf('Using median timestamp: %.3f s\n', video_frame_time); + elseif isempty(video_frame_time) + % Use middle of video + video_frame_time = []; + fprintf('Using middle frame of video\n'); + else + fprintf('Using specified timestamp: %.3f s\n', video_frame_time); + end + + % Try to read video frame + try + v = VideoReader(video_path); + + % Check if video_frame_time exceeds video duration + if ~isempty(video_frame_time) && video_frame_time > v.Duration + fprintf('Warning: Requested time (%.3f s) exceeds video duration (%.3f s). Skipping video overlay.\n', ... + video_frame_time, v.Duration); + % use previous frame if available + if exist('video_frame', 'var') && ~isempty(video_frame) + fprintf('Using previously loaded frame\n'); + else + video_frame = []; + end + elseif ~isempty(video_frame_time) + v.CurrentTime = video_frame_time; + if hasFrame(v) + video_frame = readFrame(v); + fprintf('Successfully loaded video frame\n'); + else + fprintf('Warning: No frame available at specified time\n'); + % use previous frame if available + if exist('video_frame', 'var') && ~isempty(video_frame) + fprintf('Using previously loaded frame\n'); + else + video_frame = []; + end + end + else + v.CurrentTime = v.Duration / 2; % Middle of video + if hasFrame(v) + video_frame = readFrame(v); + fprintf('Successfully loaded video frame\n'); + else + fprintf('Warning: No frame available at specified time\n'); + % use previous frame if available + if exist('video_frame', 'var') && ~isempty(video_frame) + fprintf('Using previously loaded frame\n'); + else + video_frame = []; + end + end + end + catch video_err + fprintf('Warning: VideoReader failed: %s\n', video_err.message); + fprintf('Attempting ffmpeg method...\n'); + + % Get video duration using ffprobe for validation + if ~isempty(video_frame_time) + duration_cmd = sprintf('ffprobe -v error -show_entries format=duration -of default=noprint_wrappers=1:nokey=1 "%s" 2>/dev/null', video_path); + [status_dur, duration_str] = system(duration_cmd); + if status_dur == 0 + video_duration = str2double(strtrim(duration_str)); + if video_frame_time > video_duration + fprintf('Warning: Requested time (%.3f s) exceeds video duration (%.3f s). Skipping video overlay.\n', ... + video_frame_time, video_duration); + video_frame = []; + % Skip ffmpeg extraction + continue_to_ffmpeg = false; + else + continue_to_ffmpeg = true; + end + else + % Could not get duration, try ffmpeg anyway + continue_to_ffmpeg = true; + end + else + continue_to_ffmpeg = true; + end + + % Fallback: try to extract frame using system ffmpeg + if continue_to_ffmpeg + temp_frame_path = tempname; + temp_frame_path = [temp_frame_path, '.jpg']; + + if ~isempty(video_frame_time) + cmd = sprintf('ffmpeg -i "%s" -ss %.3f -vframes 1 -y "%s" 2>/dev/null', ... + video_path, video_frame_time, temp_frame_path); + else + cmd = sprintf('ffmpeg -i "%s" -vframes 1 -y "%s" 2>/dev/null', ... + video_path, temp_frame_path); + end + + [status, ~] = system(cmd); + if status == 0 && exist(temp_frame_path, 'file') + video_frame = imread(temp_frame_path); + delete(temp_frame_path); + fprintf('Successfully extracted frame using ffmpeg\n'); + else + fprintf('Warning: ffmpeg extraction also failed\n'); + if exist(temp_frame_path, 'file') + delete(temp_frame_path); + end + end + end + end + + catch err + fprintf('Warning: Could not load video frame: %s\n', err.message); + video_frame = []; end - title('Click around the points you want to keep') - xlabel('X') - ylabel('Y') - axis tight - if axis_equal - axis equal +elseif video_overlay && isempty(video_path) + fprintf('Warning: video_overlay requested but no video_path provided\n'); +elseif video_overlay && ~exist(video_path, 'file') + fprintf('Warning: video file not found: %s\n', video_path); +end + +% set up figure +figure; + +% Display video frame as background if available +if ~isempty(video_frame) + % Display video frame + imshow(video_frame); + hold on; + + % Get image dimensions for coordinate scaling + [img_height, img_width, ~] = size(video_frame); + + % Scale tracking coordinates to match image dimensions if needed + % Assume tracking coordinates are already in pixel coordinates + x_scaled = x_plot; + y_scaled = y_plot; + + % Plot tracking data over video + if darkmode + plot(x_scaled, y_scaled, 'Color', [1,1,1,alpha], 'LineWidth', 1.5); + if add_scatter + scatter(x_scaled, y_scaled, 8, 'w', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.5); + end + else + plot(x_scaled, y_scaled, 'Color', [1,0,0,alpha], 'LineWidth', 1.5); + if add_scatter + scatter(x_scaled, y_scaled, 8, 'r', 'filled', 'MarkerEdgeColor', 'k', 'LineWidth', 0.5); + end end - darkBackground(gcf,[0.1 0.1 0.1],[0.7 0.7 0.7]) + + title('Click around the points you want to keep (Video overlay)') + + % Set axis limits to match image + xlim([0, img_width]); + ylim([0, img_height]); + axis equal; + else - fig=figure;plot(x,y,'Color',[0,0,0,alpha]);hold on - if add_scatter - scatter(x,y,3,'k','filled');hold on - end - title('Click around the points you want to keep') - xlabel('X') - ylabel('Y') - axis tight - if axis_equal - axis equal + % Original plotting without video overlay + if darkmode + plot(x_plot,y_plot,'Color',[1,1,1,alpha]);hold on + if add_scatter + scatter(x_plot,y_plot,3,'w','filled');hold on + end + title('Click around the points you want to keep') + xlabel('X') + ylabel('Y') + axis tight + if axis_equal + axis equal + end + darkBackground(gcf,[0.1 0.1 0.1],[0.7 0.7 0.7]) + else + plot(x_plot,y_plot,'Color',[0,0,0,alpha]);hold on + if add_scatter + scatter(x_plot,y_plot,3,'k','filled');hold on + end + title('Click around the points you want to keep') + xlabel('X') + ylabel('Y') + axis tight + if axis_equal + axis equal + end end end disp('PRESS "ENTER" TO EXIT') diff --git a/behavior/process_and_sync_dlc.m b/behavior/process_and_sync_dlc.m index 2568604..1c190ea 100644 --- a/behavior/process_and_sync_dlc.m +++ b/behavior/process_and_sync_dlc.m @@ -51,12 +51,39 @@ end % read video to get framerate - obj = VideoReader(fullfile(video_file.folder, video_file.name)); - fs = obj.FrameRate; + % Try VideoReader first, if it fails use ffprobe + try + obj = VideoReader(fullfile(video_file.folder, video_file.name)); + fs = obj.FrameRate; + if isempty(fs) || fs == 0 + error('Frame rate is empty or zero'); + end + catch ME + warning(ME.identifier, 'VideoReader failed: %s. Attempting to use ffprobe...', ME.message); + % Use ffprobe to get frame rate + video_path = fullfile(video_file.folder, video_file.name); + [status, result] = system(sprintf('ffprobe -v error -select_streams v:0 -show_entries stream=r_frame_rate -of default=noprint_wrappers=1:nokey=1 "%s"', video_path)); + if status == 0 + % Parse the frame rate (usually in format "30/1" or "29.97") + rate_parts = strsplit(strtrim(result), '/'); + if length(rate_parts) == 2 + fs = str2double(rate_parts{1}) / str2double(rate_parts{2}); + else + fs = str2double(result); + end + fprintf('Frame rate extracted using ffprobe: %.2f fps\n', fs); + else + error('Could not extract frame rate from video using VideoReader or ffprobe. Please install video codecs or provide frame rate manually.'); + end + end % load csv with proper header df = load_dlc_csv(fullfile(file(1).folder, file(1).name)); + fprintf('DLC CSV loaded: %d rows (tracking frames)\n', size(df, 1)); + fprintf('Video frame rate: %.2f fps\n', fs); + fprintf('Expected duration: %.2f seconds\n', size(df, 1) / fs); + % locate columns with [x,y,likelihood] field_names = df.Properties.VariableNames; @@ -76,6 +103,9 @@ x = df{:, x_col}; y = df{:, y_col}; + fprintf('Before sync - X: %d rows, Y: %d rows, ts: %d rows\n', ... + size(x, 1), size(y, 1), length(ts)); + tempTracking{count} = sync_ttl(file.folder, x, y, ts, fs, pulses_delta_range); trackFolder(count) = ii; count = count + 1; @@ -185,9 +215,43 @@ end load(fullfile(folder, 'digitalIn.events.mat'), 'digitalIn') +% Find the correct TTL channel by matching expected frame rate Len = cellfun(@length, digitalIn.timestampsOn, 'UniformOutput', false); -[~, idx] = max(cell2mat(Len)); -bazlerTtl = digitalIn.timestampsOn{idx}; +fprintf('Digital input channels found with pulse counts:\n'); +for ch = 1:length(Len) + if Len{ch} > 0 + pulse_freq = 1 / median(diff(digitalIn.timestampsOn{ch})); + fprintf(' Channel %d: %d pulses, median frequency: %.2f Hz\n', ch, Len{ch}, pulse_freq); + end +end + +% Find channel whose frequency is closest to expected video frame rate (fs) +best_idx = []; +min_freq_diff = inf; +for ch = 1:length(digitalIn.timestampsOn) + if ~isempty(digitalIn.timestampsOn{ch}) && length(digitalIn.timestampsOn{ch}) > 10 + pulse_freq = 1 / median(diff(digitalIn.timestampsOn{ch})); + freq_diff = abs(pulse_freq - fs); + if freq_diff < min_freq_diff && abs(pulse_freq - fs) / fs < 0.2 + min_freq_diff = freq_diff; + best_idx = ch; + end + end +end + +if isempty(best_idx) + warning('No digital input channel matches expected frame rate (%.2f Hz). Using channel with most pulses.', fs); + [~, best_idx] = max(cell2mat(Len)); +else + fprintf('Selected channel %d for camera TTLs (frequency matches expected %.2f Hz)\n', best_idx, fs); +end + +bazlerTtl = digitalIn.timestampsOn{best_idx}; + +% Debug: Check TTL timestamps +fprintf('=== DEBUG: TTL timestamp analysis ===\n'); +fprintf('bazlerTtl range: %.6f to %.6f seconds\n', min(bazlerTtl), max(bazlerTtl)); +fprintf('bazlerTtl sample values: [%.6f, %.6f, %.6f, ...]\n', bazlerTtl(1), bazlerTtl(2), bazlerTtl(3)); %check for extra pulses of much shorter distance than they should extra_pulses = diff(bazlerTtl) < ((1 / fs) - (1 / fs) * pulses_delta_range); @@ -203,7 +267,6 @@ warning('MyComponent:MyWarningID', 'First ttl %f may indicate video was started before intan', bazlerTtl(1)); end - if isempty(bazlerTtl) warning('no real ttls, something went wrong') disp('fix the issue and type dbcont to continue') @@ -212,6 +275,16 @@ basler_intan_diff = length(bazlerTtl) - size(x, 1); +fprintf('\n=== TTL vs DLC Frame Mismatch Analysis ===\n'); +fprintf('Number of TTL pulses: %d\n', length(bazlerTtl)); +fprintf('Number of DLC frames: %d\n', size(x, 1)); +if basler_intan_diff > 0 + fprintf('Difference: %d more TTL pulses than DLC frames\n', basler_intan_diff); +else + fprintf('Difference: %d more DLC frames than TTL pulses\n', abs(basler_intan_diff)); +end +fprintf('==========================================\n\n'); + [x, y, ts, bazlerTtl, fs, notes] = match_basler_frames_to_ttl(bazlerTtl, basler_intan_diff, x, y, ts, fs); [~, folder_name] = fileparts(folder); diff --git a/preProcessing/acqID.m b/preProcessing/acqID.m index 3ff8d89..1b9c716 100644 --- a/preProcessing/acqID.m +++ b/preProcessing/acqID.m @@ -1,6 +1,6 @@ -function [datpaths, recordingnames] = acqID(basepath, sortFiles, altSort, ignoreFolders) +function [datpaths, recordingnames] = acqID(basepath, sortFiles, altSort, ignoreFolders, behaviorOnly) -%% [datpaths, recordingnames] = acqID(basepath, sortFiles, altSort) +%% [datpaths, recordingnames] = acqID(basepath, sortFiles, altSort, ignoreFolders, behaviorOnly) % Function to acquire the ID of dat files (Intan or openEphys) within the % given session folder (basepath) and sort according to the indicated @@ -29,6 +29,10 @@ % folder containing original copies of the data. % Example input may look like: ["backup", % "ignore"]. +% behaviorOnly Logical flag to only search for digitalin.dat +% files (behavior-only mode without electrophysiology). +% Default is false (searches for continuous.dat, +% amplifier.dat, and digitalin.dat). %% OUTPUTS % datpaths Returns cell array of ordered datpaths in "PATH " @@ -36,6 +40,10 @@ % recordingnames Returns cell array of ordered subsession folder % names. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if nargin < 5 + behaviorOnly = false; % Default to standard mode +end + if sortFiles && (~isempty(altSort)) error('sortFiles cannot be empty while altSort provides an order. Please choose to either sort by time (sortFiles=true) or designate a manual order of concatenation (altSort)'); end @@ -51,7 +59,18 @@ listing = dir("**"); allFolders = struct2table(listing); -useIDX = find(contains(allFolders.name, 'continuous.dat') | contains(allFolders.name, 'amplifier.dat')); + +% Choose which files to search for based on behaviorOnly mode +if behaviorOnly + % Behavior-only mode: only search for digitalin.dat + useIDX = find(contains(allFolders.name, 'digitalin.dat')); + fprintf('Running in behavior-only mode: searching only for digitalin.dat files\n'); +else + % Standard mode: search for all recording files + useIDX = find(contains(allFolders.name, 'continuous.dat') | ... + contains(allFolders.name, 'amplifier.dat')); +end + folderNames = allFolders(useIDX, :).folder; %this assumes that your recording folders with the date are within your @@ -82,11 +101,22 @@ expNum(i) = '1'; recNum(i) = '1'; else - recordingnames{i} = relParts{1}; - expIdx = find(startsWith(relParts, 'experiment'), 1); - recIdx = find(startsWith(relParts, 'recording'), 1); - expNum(i) = ternary(~isempty(expIdx), extractAfter(relParts{expIdx}, 'experiment'), '1'); - recNum(i) = ternary(~isempty(recIdx), extractAfter(relParts{recIdx}, 'recording'), '1'); + % OpenEphys format: multiple folders (experiment/recording structure) + recordingnames{i} = parts{1}; + expIdx = find(startsWith(parts, 'experiment'), 1); + recIdx = find(startsWith(parts, 'recording'), 1); + + if ~isempty(expIdx) + expNum(i) = extractAfter(parts{expIdx}, 'experiment'); + else + expNum(i) = '1'; + end + + if ~isempty(recIdx) + recNum(i) = extractAfter(parts{recIdx}, 'recording'); + else + recNum(i) = '1'; + end end else fprintf('.dat file found nested in a folder labeled %s . Skipping: \n', ignoreFolders(f)); diff --git a/preProcessing/concatenateDats.m b/preProcessing/concatenateDats.m index 9da9198..5f65c5b 100755 --- a/preProcessing/concatenateDats.m +++ b/preProcessing/concatenateDats.m @@ -28,6 +28,10 @@ function concatenateDats(varargin) % folder containing original copies of the data. % Example input may look like: ["backup", % "ignore"]. +% behaviorOnly Set to true for behavior-only sessions (no +% amplifier.dat). Reads metadata from info.rhd or +% settings.xml. Skips LFP, Kilosort, and SleepScore. +% Default false p = inputParser; addParameter(p, 'basepath', pwd, @isfolder); % by default, current folder @@ -35,6 +39,7 @@ function concatenateDats(varargin) addParameter(p, 'sortFiles', true, @islogical); addParameter(p, 'altSort', [], @isnumeric); addParameter(p, 'ignoreFolders', "", @isstring); +addParameter(p, 'behaviorOnly', false, @islogical); parse(p, varargin{:}); @@ -43,13 +48,14 @@ function concatenateDats(varargin) sortFiles = p.Results.sortFiles; altSort = p.Results.altSort; ignoreFolders = p.Results.ignoreFolders; +behaviorOnly = p.Results.behaviorOnly; if sortFiles && (~isempty(altSort)) error('sortFiles cannot be empty while altSort provides an order. Please choose to either sort by time (sortFiles=true) or designate a manual order of concatenation (altSort)'); end basename = basenameFromBasepath(basepath); -[datpaths, recordingnames] = acqID(basepath, sortFiles, altSort, ignoreFolders); +[datpaths, recordingnames] = acqID(basepath, sortFiles, altSort, ignoreFolders, behaviorOnly); if isempty(datpaths) disp('no subsessions detected, exiting concatenation'); return @@ -60,7 +66,14 @@ function concatenateDats(varargin) end end -otherdattypes = {'analogin'; 'digitalin'; 'auxiliary'; 'time'; 'supply'}; +if behaviorOnly + % In behavior-only mode, digitalin is the primary file (replaces amplifier) + % so we don't include it in otherdattypes + otherdattypes = {'analogin'; 'auxiliary'; 'time'; 'supply'}; +else + % Standard mode: amplifier is primary, digitalin is in otherdattypes + otherdattypes = {'analogin'; 'digitalin'; 'auxiliary'; 'time'; 'supply'}; +end toFill = zeros(size(otherdattypes)); fileBase = cell(size(datpaths)); @@ -79,7 +92,7 @@ function concatenateDats(varargin) end for ii = 1:length(toFill) if toFill(ii) == 1 - fillMissingDats('basepath', basepath, 'fileType', otherdattypes{ii}); + fillMissingDats('basepath', basepath, 'fileType', otherdattypes{ii}, 'behaviorOnly', behaviorOnly); end end else @@ -99,7 +112,7 @@ function concatenateDats(varargin) if datCount(j) < length(datpaths) toFill(j) = false; disp([otherdattypes{j}, ' present in ', num2str(datCount(j)), ... - '/', num2str(length(datpaths)), ' subfolders, will not concatenate.']); + filesep, num2str(length(datpaths)), ' subfolders, will not concatenate.']); %If receiving this message, can make 'fillMissingDatFiles' true %to fill the subfolders missing the files with zero'd data else @@ -232,12 +245,61 @@ function concatenateDats(varargin) end %% Create MergePoints file based on concatenation -sessionInfo = LoadXml(fullfile(basepath, [basename, '.xml'])); -nSamp = []; -for didx = 1:length(datpaths) - t = dir(datpaths{didx}(2:end - 2)); - dataTypeNBytes = numel(typecast(cast(0, 'int16'), 'uint8')); - nSamp(didx) = t.bytes / (sessionInfo.nChannels * dataTypeNBytes); +% For behavior-only sessions, read metadata from info.rhd or session file +% instead of requiring XML file +if behaviorOnly + % Try to load session file to get sampling rate + sessionFile = fullfile(basepath, [basename, '.session.mat']); + if exist(sessionFile, 'file') + load(sessionFile, 'session'); + sessionInfo.SampleRate = session.extracellular.sr; + sessionInfo.nChannels = 1; % For digitalin.dat, we use bytes directly + else + % Try to read from info.rhd + rhdFile = dir(fullfile(basepath, '**', '*.rhd')); + if ~isempty(rhdFile) + try + [~, ~, ~, ~, ~, ~, frequency_parameters, ~] = ... + read_Intan_RHD2000_file_bz('basepath', rhdFile(1).folder); + sessionInfo.SampleRate = frequency_parameters.board_dig_in_sample_rate; + sessionInfo.nChannels = 1; % For digitalin.dat + catch + warning('Could not read sampling rate from info.rhd. Using default 20000 Hz'); + sessionInfo.SampleRate = 20000; + sessionInfo.nChannels = 1; + end + else + warning('No session file or info.rhd found. Using default sampling rate 20000 Hz'); + sessionInfo.SampleRate = 20000; + sessionInfo.nChannels = 1; + end + end + + % For behavior-only, calculate samples from digitalin.dat (uint16, 2 bytes per sample) + nSamp = []; + for didx = 1:length(datpaths) + % Find digitalin.dat in this subsession folder + subsessionPath = datpaths{didx}(2:end - 2); % Remove quotes and spaces + [parentPath, ~, ~] = fileparts(subsessionPath); + digitalinPath = fullfile(parentPath, 'digitalin.dat'); + + if exist(digitalinPath, 'file') + t = dir(digitalinPath); + % digitalin.dat: uint16, 2 bytes per sample + nSamp(didx) = t.bytes / 2; + else + error(['digitalin.dat not found in ', parentPath]); + end + end +else + % Standard mode: use XML file + sessionInfo = LoadXml(fullfile(basepath, [basename, '.xml'])); + nSamp = []; + for didx = 1:length(datpaths) + t = dir(datpaths{didx}(2:end - 2)); + dataTypeNBytes = numel(typecast(cast(0, 'int16'), 'uint8')); + nSamp(didx) = t.bytes / (sessionInfo.nChannels * dataTypeNBytes); + end end cumsum_nSamp = cumsum(nSamp); diff --git a/preProcessing/fillMissingDats.m b/preProcessing/fillMissingDats.m index c47b0ea..c04177a 100644 --- a/preProcessing/fillMissingDats.m +++ b/preProcessing/fillMissingDats.m @@ -33,10 +33,12 @@ function fillMissingDats(varargin) isFileType = @(x) sum(strcmp(x, otherdattypes)) == 1; addParameter(p, 'basepath', pwd, @isfolder) addParameter(p, 'fileType', [], isFileType) +addParameter(p, 'behaviorOnly', false, @islogical); parse(p, varargin{:}) basepath = p.Results.basepath; fileType = p.Results.fileType; +behaviorOnly = p.Results.behaviorOnly; % get file types and data types files_table = table(); @@ -49,31 +51,83 @@ function fillMissingDats(varargin) %% check location of each file and what to concat typeFiles = dir([basepath, filesep, '*', filesep, fileType, '.dat']); -ampFiles = dir([basepath, filesep, '*', filesep, 'amplifier.dat']); -contFiles = dir([basepath, filesep, '**', filesep, 'continuous.dat']); %check for openEphys -ampFiles = cat(1, ampFiles, contFiles); + +% For behavior-only sessions, use digitalin.dat as reference instead of amplifier.dat +if behaviorOnly || ampNch == 0 + disp('Behavior-only mode: using digitalin.dat as reference'); + ampFiles = dir([basepath, filesep, '*', filesep, 'digitalin.dat']); + if isempty(ampFiles) + warning('No digitalin.dat files found in subsessions. Cannot fill missing files.'); + return + end +else + % Standard mode: use amplifier.dat or continuous.dat + ampFiles = dir([basepath, filesep, '*', filesep, 'amplifier.dat']); + contFiles = dir([basepath, filesep, '**', filesep, 'continuous.dat']); %check for openEphys + ampFiles = cat(1, ampFiles, contFiles); +end + typeFolders = {typeFiles.folder}; ampFolders = {ampFiles.folder}; fillInds = find(~ismember(ampFolders, typeFolders)); %index of folders that need fill -% calculate number of channels to for fill file +if isempty(fillInds) + disp(['All subsessions already have ', fileType, '.dat files. No filling needed.']); + return +end + +% calculate number of channels for fill file refTypeInd = find(ismember(typeFolders, ampFolders), 1); + +if isempty(refTypeInd) + warning(['Cannot find a subsession with both ', fileType, '.dat and reference file. Cannot determine channel count.']); + return +end + refTypeSize = typeFiles(refTypeInd).bytes; -refAmpInd = strcmp(ampFolders, typeFiles(refTypeInd).folder); +refAmpInd = find(strcmp(ampFolders, typeFiles(refTypeInd).folder), 1); + +if isempty(refAmpInd) + warning('Cannot find matching reference file. Cannot fill missing files.'); + return +end + refAmpSize = ampFiles(refAmpInd).bytes; -typeNch = refTypeSize * ampNch / refAmpSize; + +% Calculate number of channels based on file type +if behaviorOnly || ampNch == 0 + % For behavior-only, use digitalin as reference (16 channels, uint16) + digitalNch = 16; % Standard Intan digital input channels + refNch = digitalNch; + refBytesPerSample = 2; % uint16 + typeNch = refTypeSize * refNch / refAmpSize; +else + % Standard mode with amplifier channels + typeNch = refTypeSize * ampNch / refAmpSize; +end %% Fill in missing files for ii = 1:length(fillInds) fillIdx = fillInds(ii); localAmpSize = ampFiles(fillIdx).bytes; - nPoints = localAmpSize / (ampNch * 2); + + % Calculate number of samples based on reference file + if behaviorOnly || ampNch == 0 + % For behavior-only, reference is digitalin.dat (16 channels, uint16, 2 bytes) + nPoints = localAmpSize / (refNch * refBytesPerSample); + else + % Standard mode with amplifier + nPoints = localAmpSize / (ampNch * 2); + end zeroData = zeros(typeNch, nPoints); filepath = [ampFiles(fillIdx).folder, filesep, fileType, '.dat']; + disp(['Creating missing file: ', filepath]); fid = fopen(filepath, 'w'); fwrite(fid, zeroData, files_table.data_type{contains(files_table.files, fileType)}); fclose(fid); end + +disp(['Successfully filled ', num2str(length(fillInds)), ' missing ', fileType, '.dat files']); diff --git a/preProcessing/preprocessSession.m b/preProcessing/preprocessSession.m index 3457d5c..346fc36 100644 --- a/preProcessing/preprocessSession.m +++ b/preProcessing/preprocessSession.m @@ -15,6 +15,10 @@ function preprocessSession(varargin) % ------------------------------------------------------------------------- % basepath Basepath for experiment. It contains all session % folders. If not provided takes pwd +% behaviorOnly Set to true for behavior-only sessions (no amplifier.dat, +% no spike sorting). Reads metadata from info.rhd or +% settings.xml. Skips LFP, Kilosort, and SleepScore. +% Default false % analogChannels List of analog channels with pulses to be detected (it % supports Intan Buzsaki Edition) % digitalChannels List of digital channels with pulses to be detected (it @@ -94,6 +98,7 @@ function preprocessSession(varargin) p = inputParser; addParameter(p, 'basepath', pwd, @isfolder); % by default, current folder +addParameter(p, 'behaviorOnly', false, @islogical); % behavior-only mode (no amplifier.dat) addParameter(p, 'fillMissingDatFiles', false, @islogical); addParameter(p, 'analogInputs', false, @islogical); addParameter(p, 'analogChannels', [], @isnumeric); @@ -130,6 +135,7 @@ function preprocessSession(varargin) parse(p, varargin{:}); basepath = p.Results.basepath; +behaviorOnly = p.Results.behaviorOnly; fillMissingDatFiles = p.Results.fillMissingDatFiles; analogInputs = p.Results.analogInputs; analogChannels = p.Results.analogChannels; @@ -180,12 +186,16 @@ function preprocessSession(varargin) [~, basename] = fileparts(basepath); % Get xml file in order -if ~exist([basepath, filesep, basename, '.xml'], 'file') - xmlFile = checkFile('fileType', '.xml', 'searchSubdirs', true); - xmlFile = xmlFile(1); - if ~(strcmp(xmlFile.folder, basepath) && strcmp(xmlFile.name(1:end-4), basename)) - copyfile([xmlFile.folder, filesep, xmlFile.name], [basepath, filesep, basename, '.xml']) +if ~behaviorOnly + if ~exist([basepath, filesep, basename, '.xml'], 'file') + xmlFile = checkFile('fileType', '.xml', 'searchSubdirs', true); + xmlFile = xmlFile(1); + if ~(strcmp(xmlFile.folder, basepath) && strcmp(xmlFile.name(1:end-4), basename)) + copyfile([xmlFile.folder, filesep, xmlFile.name], [basepath, filesep, basename, '.xml']) + end end +else + disp('Behavior-only mode: skipping XML file requirement'); end % Check info.rhd @@ -202,16 +212,25 @@ function preprocessSession(varargin) %% Make SessionInfo % Manually ID bad channels at this point. automating it would be good -session = sessionTemplate(basepath, 'showGUI', false); +if behaviorOnly + disp('Creating behavior-only session metadata...'); + session = sessionTemplate_behaviorOnly(basepath, 'showGUI', false); +else + session = sessionTemplate(basepath, 'showGUI', false); +end save(fullfile(basepath, [basename, '.session.mat']), 'session'); %% Concatenate sessions disp('Concatenate session folders...'); concatenateDats('basepath', basepath, 'fillMissingDatFiles', fillMissingDatFiles, ... - 'sortFiles', sortFiles, 'altSort', altSort, 'ignoreFolders', ignoreFolders); + 'sortFiles', sortFiles, 'altSort', altSort, 'ignoreFolders', ignoreFolders, 'behaviorOnly', behaviorOnly); %% run again to add epochs from basename.MergePoints.mat -session = sessionTemplate(basepath, 'showGUI', false); +if behaviorOnly + session = sessionTemplate_behaviorOnly(basepath, 'showGUI', false); +else + session = sessionTemplate(basepath, 'showGUI', false); +end save(fullfile(basepath, [basename, '.session.mat']), 'session'); %% Process additional inputs - CHECK FOR OUR LAB @@ -251,6 +270,28 @@ function preprocessSession(varargin) cleanPulses(pulses.ints{1}(:)); end +% Skip electrophysiology processing in behavior-only mode +if behaviorOnly + disp('Behavior-only mode: Skipping LFP, spike sorting, and brain state detection'); + + % Still process tracking if requested + if getPos + % check for pre existing deeplab cut + if exist(path_to_dlc_bat_file, 'file') + system(path_to_dlc_bat_file, '-echo') + end + % put tracking into standard format + general_behavior_file('basepath', basepath) + end + + % Log and exit + results = p.Results; + save(fullfile(basepath, 'preprocessSession_params.mat'), 'results') + targetFile = fullfile(basepath, 'preprocessSession.log'); + copyfile(which('preprocessSession.m'), targetFile); + return +end + if LFPbeforeKilo %Make LFP runLFP(basepath, basename, session); diff --git a/preProcessing/read_Intan_RHD2000_file_bz.m b/preProcessing/read_Intan_RHD2000_file_bz.m index c73e653..f80f6f0 100644 --- a/preProcessing/read_Intan_RHD2000_file_bz.m +++ b/preProcessing/read_Intan_RHD2000_file_bz.m @@ -54,25 +54,39 @@ %______________________________________________________________________ p = inputParser; -p.addParameter('basepath',pwd,@isfolder); +p.addParameter('basepath',pwd,@isfolder); p.parse(varargin{:}); path = p.Results.basepath; -% path = [pwd,filesep]; -file = [ls(fullfile(path,'*.rhd'))]; - -if (file == 0) +% Cross-platform file discovery using dir +d = dir(fullfile(path,'*.rhd')); +if isempty(d) + warning('No .rhd files found in path: %s', path); return; -elseif size(file,1) > 1 +end + +% Reconstruct full absolute paths consistently across platforms +files = fullfile({d.folder}, {d.name}); +files = string(files); + +if length(files) == 1 + file = char(files(1)); +else + % Multiple .rhd files found, select the one matching basename basename = basenameFromBasepath(path); - i=1; - while i<=size(file,1) - tempFile = file(i,:); - if contains(tempFile, basename) - file = []; - file = tempFile; + file_found = false; + for i = 1:length(files) + [~, fname, ~] = fileparts(files(i)); + if contains(fname, basename) + file = char(files(i)); + file_found = true; + break; end - i=i+1; + end + if ~file_found + % If no basename match, use the first file + file = char(files(1)); + warning('Multiple .rhd files found, but none matched basename "%s". Using: %s', basename, file); end end @@ -82,7 +96,9 @@ % file = d(end).name; tic; -filename = fullfile(path,file); +file = strip(file); +% file is already a full path, no need to use fullfile again +filename = file; fid = fopen(filename, 'r'); s = dir(filename); @@ -139,19 +155,19 @@ 'note1', fread_QString(fid), ... 'note2', fread_QString(fid), ... 'note3', fread_QString(fid) ); - + % If data file is from GUI v1.1 or later, see if temperature sensor data % was saved. num_temp_sensor_channels = 0; if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 1) ... - || (data_file_main_version_number > 1)) + || (data_file_main_version_number > 1)) num_temp_sensor_channels = fread(fid, 1, 'int16'); end % If data file is from GUI v1.3 or later, load eval board mode. eval_board_mode = 0; if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 3) ... - || (data_file_main_version_number > 1)) + || (data_file_main_version_number > 1)) eval_board_mode = fread(fid, 1, 'int16'); end @@ -250,7 +266,7 @@ new_trigger_channel(1).digital_edge_polarity = fread(fid, 1, 'int16'); new_channel(1).electrode_impedance_magnitude = fread(fid, 1, 'single'); new_channel(1).electrode_impedance_phase = fread(fid, 1, 'single'); - + if (channel_enabled) switch (signal_type) case 0 @@ -276,7 +292,7 @@ error('Unknown channel type'); end end - + end end end @@ -326,7 +342,7 @@ end % Temp sensor is sampled once per data block if (num_temp_sensor_channels > 0) - bytes_per_block = bytes_per_block + 1 * 2 * num_temp_sensor_channels; + bytes_per_block = bytes_per_block + 1 * 2 * num_temp_sensor_channels; end % How many data blocks remain in this file? @@ -358,7 +374,7 @@ end if (data_present) - + % Pre-allocate memory for data. fprintf(1, 'Allocating memory for data...\n'); @@ -391,7 +407,7 @@ % integeters to signed integers to accomidate negative (adjusted) % timestamps for pretrigger data. if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 2) ... - || (data_file_main_version_number > 1)) + || (data_file_main_version_number > 1)) t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'int32'); else t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint32'); @@ -444,28 +460,28 @@ fclose(fid); if (data_present) - + fprintf(1, 'Parsing data...\n'); % Extract digital input channels to separate variables. for i=1:num_board_dig_in_channels - mask = 2^(board_dig_in_channels(i).native_order) * ones(size(board_dig_in_raw)); - board_dig_in_data(i, :) = (bitand(board_dig_in_raw, mask) > 0); + mask = 2^(board_dig_in_channels(i).native_order) * ones(size(board_dig_in_raw)); + board_dig_in_data(i, :) = (bitand(board_dig_in_raw, mask) > 0); end for i=1:num_board_dig_out_channels - mask = 2^(board_dig_out_channels(i).native_order) * ones(size(board_dig_out_raw)); - board_dig_out_data(i, :) = (bitand(board_dig_out_raw, mask) > 0); + mask = 2^(board_dig_out_channels(i).native_order) * ones(size(board_dig_out_raw)); + board_dig_out_data(i, :) = (bitand(board_dig_out_raw, mask) > 0); end % Scale voltage levels appropriately. amplifier_data = amplifier_data - 32768; -% amplifier_data = 0.195 * (amplifier_data - 32768); % converting units into microvolts. Raly: commented this out because want the raw units as this is what is recorded when .dat file outputs are used directly + % amplifier_data = 0.195 * (amplifier_data - 32768); % converting units into microvolts. Raly: commented this out because want the raw units as this is what is recorded when .dat file outputs are used directly aux_input_data = 37.4e-6 * aux_input_data; % units = volts supply_voltage_data = 74.8e-6 * supply_voltage_data; % units = volts if (eval_board_mode == 1) board_adc_data = 152.59e-6 * (board_adc_data - 32768); % units = volts elseif (eval_board_mode == 13) % Intan Recording Controller - board_adc_data = 312.5e-6 * (board_adc_data - 32768); % units = volts + board_adc_data = 312.5e-6 * (board_adc_data - 32768); % units = volts else board_adc_data = 50.354e-6 * board_adc_data; % units = volts end @@ -587,99 +603,99 @@ function a = fread_QString(fid) -% a = read_QString(fid) -% -% Read Qt style QString. The first 32-bit unsigned number indicates -% the length of the string (in bytes). If this number equals 0xFFFFFFFF, -% the string is null. + % a = read_QString(fid) + % + % Read Qt style QString. The first 32-bit unsigned number indicates + % the length of the string (in bytes). If this number equals 0xFFFFFFFF, + % the string is null. -a = ''; -length = fread(fid, 1, 'uint32'); -if length == hex2num('ffffffff') - return; -end -% convert length from bytes to 16-bit Unicode words -length = length / 2; + a = ''; + length = fread(fid, 1, 'uint32'); + if length == hex2num('ffffffff') + return; + end + % convert length from bytes to 16-bit Unicode words + length = length / 2; -for i=1:length - a(i) = fread(fid, 1, 'uint16'); -end + for i=1:length + a(i) = fread(fid, 1, 'uint16'); + end -return + return function s = plural(n) -% s = plural(n) -% -% Utility function to optionally plurailze words based on the value -% of n. + % s = plural(n) + % + % Utility function to optionally plurailze words based on the value + % of n. -if (n == 1) - s = ''; -else - s = 's'; -end + if (n == 1) + s = ''; + else + s = 's'; + end -return + return function out = notch_filter(in, fSample, fNotch, Bandwidth) -% out = notch_filter(in, fSample, fNotch, Bandwidth) -% -% Implements a notch filter (e.g., for 50 or 60 Hz) on vector 'in'. -% fSample = sample rate of data (in Hz or Samples/sec) -% fNotch = filter notch frequency (in Hz) -% Bandwidth = notch 3-dB bandwidth (in Hz). A bandwidth of 10 Hz is -% recommended for 50 or 60 Hz notch filters; narrower bandwidths lead to -% poor time-domain properties with an extended ringing response to -% transient disturbances. -% -% Example: If neural data was sampled at 30 kSamples/sec -% and you wish to implement a 60 Hz notch filter: -% -% out = notch_filter(in, 30000, 60, 10); - -tstep = 1/fSample; -Fc = fNotch*tstep; - -L = length(in); - -% Calculate IIR filter parameters -d = exp(-2*pi*(Bandwidth/2)*tstep); -b = (1 + d*d)*cos(2*pi*Fc); -a0 = 1; -a1 = -b; -a2 = d*d; -a = (1 + d*d)/2; -b0 = 1; -b1 = -2*cos(2*pi*Fc); -b2 = 1; - -out = zeros(size(in)); -out(1) = in(1); -out(2) = in(2); -% (If filtering a continuous data stream, change out(1) and out(2) to the -% previous final two values of out.) - -% Run filter -for i=3:L - out(i) = (a*b2*in(i-2) + a*b1*in(i-1) + a*b0*in(i) - a2*out(i-2) - a1*out(i-1))/a0; -end + % out = notch_filter(in, fSample, fNotch, Bandwidth) + % + % Implements a notch filter (e.g., for 50 or 60 Hz) on vector 'in'. + % fSample = sample rate of data (in Hz or Samples/sec) + % fNotch = filter notch frequency (in Hz) + % Bandwidth = notch 3-dB bandwidth (in Hz). A bandwidth of 10 Hz is + % recommended for 50 or 60 Hz notch filters; narrower bandwidths lead to + % poor time-domain properties with an extended ringing response to + % transient disturbances. + % + % Example: If neural data was sampled at 30 kSamples/sec + % and you wish to implement a 60 Hz notch filter: + % + % out = notch_filter(in, 30000, 60, 10); + + tstep = 1/fSample; + Fc = fNotch*tstep; + + L = length(in); + + % Calculate IIR filter parameters + d = exp(-2*pi*(Bandwidth/2)*tstep); + b = (1 + d*d)*cos(2*pi*Fc); + a0 = 1; + a1 = -b; + a2 = d*d; + a = (1 + d*d)/2; + b0 = 1; + b1 = -2*cos(2*pi*Fc); + b2 = 1; + + out = zeros(size(in)); + out(1) = in(1); + out(2) = in(2); + % (If filtering a continuous data stream, change out(1) and out(2) to the + % previous final two values of out.) + + % Run filter + for i=3:L + out(i) = (a*b2*in(i-2) + a*b1*in(i-1) + a*b0*in(i) - a2*out(i-2) - a1*out(i-1))/a0; + end -return + return function move_to_base_workspace(variable) -% move_to_base_workspace(variable) -% -% Move variable from function workspace to base MATLAB workspace so -% user will have access to it after the program ends. + % move_to_base_workspace(variable) + % + % Move variable from function workspace to base MATLAB workspace so + % user will have access to it after the program ends. -variable_name = inputname(1); -assignin('base', variable_name, variable); + variable_name = inputname(1); + assignin('base', variable_name, variable); -return; + return; diff --git a/preProcessing/sessionTemplate_behaviorOnly.m b/preProcessing/sessionTemplate_behaviorOnly.m new file mode 100644 index 0000000..1fb8f9b --- /dev/null +++ b/preProcessing/sessionTemplate_behaviorOnly.m @@ -0,0 +1,1093 @@ +function session = sessionTemplate_behaviorOnly(input1, varargin) + % sessionTemplate_behaviorOnly - Create session metadata for behavior-only recordings + % + % This function creates a minimal session struct for recordings that only contain + % behavioral data (digital/analog inputs, tracking) without electrophysiology data + % (no amplifier.dat, no spike sorting needed). + % + % It reads metadata from: + % - info.rhd (Intan) or settings.xml (OpenEphys) + % - time.dat (to calculate duration) + % - digitalin.dat, analogin.dat, digitalout.dat (for channel info) + % + % USAGE: + % session = sessionTemplate_behaviorOnly(basepath) + % session = sessionTemplate_behaviorOnly(basepath, 'showGUI', true) + % session = sessionTemplate_behaviorOnly(session) % update existing session + % + % INPUTS: + % input1 - Either basepath (string) or existing session struct + % showGUI - Show session GUI editor (default: false) + % basename - Override basename detection + % addBehaviorFiles - Auto-detect tracking/video/DLC files (default: false) + % nDigitalChannels - Number of digital input channels (default: 16 for Intan) + % nAnalogChannels - Number of analog input channels (default: 8 for Intan) + % + % OUTPUTS: + % session - Session metadata structure + % + % EXAMPLE: + % session = sessionTemplate_behaviorOnly('/path/to/behavior/session'); + % + % % For a system with 32 digital channels: + % session = sessionTemplate_behaviorOnly(basepath, 'nDigitalChannels', 32); + % + % See also: sessionTemplate, preprocessSession + + % Parse inputs + p = inputParser; + addRequired(p, 'input1', @(X) (ischar(X) && exist(X, 'dir')) || isstruct(X)); + addParameter(p, 'basename', [], @ischar); + addParameter(p, 'showGUI', false, @islogical); + addParameter(p, 'addBehaviorFiles', false, @islogical); + addParameter(p, 'nDigitalChannels', 16, @isnumeric); % Intan default + addParameter(p, 'nAnalogChannels', 8, @isnumeric); % Intan default + + parse(p, input1, varargin{:}) + basename = p.Results.basename; + showGUI = p.Results.showGUI; + addBehaviorFiles = p.Results.addBehaviorFiles; + nDigitalChannels = p.Results.nDigitalChannels; + nAnalogChannels = p.Results.nAnalogChannels; + + % Initialize session struct and basepath + if ischar(input1) + basepath = input1; + cd(basepath) + elseif isstruct(input1) + session = input1; + + if isfield(session.general, 'basePath') && exist(session.general.basePath, 'dir') + basepath = session.general.basePath; + cd(basepath) + else + basepath = pwd; + end + + end + + % Get basename + if isempty(basename) + basename = basenameFromBasepath(basepath); + end + + % Load existing session file if it exists + if ~exist('session', 'var') && exist(fullfile(basepath, [basename, '.session.mat']), 'file') + disp('Loading existing basename.session.mat file') + session = loadSession(basepath, basename); + elseif ~exist('session', 'var') + session = []; + end + + pathPieces = regexp(basepath, filesep, 'split'); + + %% General metadata + session.general.version = 5; + session.general.name = basename; + session.general.basePath = basepath; + session.general.sessionType = 'Behavior'; + session.general.experimenters = ''; % Add experimenter names if known + + %% Animal metadata + if ~isfield(session, 'animal') + session.animal.name = pathPieces{end - 1}; + session.animal.sex = 'Unknown'; + session.animal.species = 'Unknown'; + session.animal.strain = 'Unknown'; + session.animal.geneticLine = ''; + end + + %% Detect recording system and read metadata + rhdFile = dir(fullfile(basepath, '**', '*.rhd')); + settingsFile = dir(fullfile(basepath, '**', 'settings.xml')); + + if ~isempty(rhdFile) + % Intan system detected + disp('Intan system detected (info.rhd found)') + session = readIntanBehaviorMetadata(session, basepath, rhdFile(1)); + + elseif ~isempty(settingsFile) + % OpenEphys system detected + disp('OpenEphys system detected (settings.xml found)') + session = readOpenEphysBehaviorMetadata(session, basepath, settingsFile(1)); + + else + % No metadata files found - use defaults and warn user + warning('No info.rhd or settings.xml found. Using default parameters.'); + session = setDefaultBehaviorMetadata(session, basepath); + end + + %% Calculate duration from time.dat or data files + session = calculateSessionDuration(session, basepath, basename); + + %% Create epochs from MergePoints or default + session = createBehaviorEpochs(session, basepath, basename); + + %% Extract date/time from folder or files + session = extractSessionDateTime(session, basepath, basename); + + %% Add behavioral data paths (if requested) + if addBehaviorFiles + session = addBehavioralPaths(session, basepath, basename); + else + % Remove behavioral file fields if they exist (cleanup from previous runs) + fieldsToRemove = {'videos', 'deepLabCut', 'optitrack', 'behavioralTracking'}; + + for i = 1:length(fieldsToRemove) + + if isfield(session, fieldsToRemove{i}) + session = rmfield(session, fieldsToRemove{i}); + end + + end + + end + + %% Show GUI if requested + if showGUI + session = gui_session(session); + end + +end + +%% Helper Functions + +function session = readIntanBehaviorMetadata(session, basepath, rhdFile) + % readIntanBehaviorMetadata - Extract metadata from Intan RHD2000 files + % + % Reads recording parameters from an Intan info.rhd file including: + % - Sampling rates (main and LFP) + % - Digital input channels (board_dig_in) + % - Analog input channels (board_adc) + % - Auxiliary channels (accelerometer, etc.) + % - Recording notes + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % rhdFile - Dir struct entry for the .rhd file + % + % OUTPUTS: + % session - Session struct with Intan metadata populated + % + % NOTES: + % - Uses read_Intan_RHD2000_file_bz() to parse binary RHD format + % - Falls back to setDefaultBehaviorMetadata() on error + % - Sets extracellular.nChannels = 0 (behavior-only, no amplifier data) + + try + % Use existing Intan reader + rhdPath = fullfile(rhdFile.folder, rhdFile.name); + + % Read RHD file + [amplifier_channels, notes, aux_input_channels, ~, ... + board_dig_in_channels, supply_voltage_channels, frequency_parameters, board_adc_channels] = ... + read_Intan_RHD2000_file_bz('basepath', rhdFile.folder); + + % Sampling rate + session.extracellular.sr = frequency_parameters.amplifier_sample_rate; + session.extracellular.srLfp = 1250; % Standard LFP rate (not used in behavior-only) + + % Digital channels + if ~isempty(board_dig_in_channels) + nDigitalChannels = length(board_dig_in_channels); + session.inputs.nDigitalChannels = nDigitalChannels; + session.inputs.digitalChannels = 1:nDigitalChannels; + + for i = 1:nDigitalChannels + session.inputs.digitalChannelNames{i} = board_dig_in_channels(i).native_channel_name; + end + + else + session.inputs.nDigitalChannels = 0; + end + + % Analog channels (board ADC) + if ~isempty(board_adc_channels) + nAnalogChannels = length(board_adc_channels); + session.inputs.nAnalogChannels = nAnalogChannels; + session.inputs.analogChannels = 1:nAnalogChannels; + + for i = 1:nAnalogChannels + session.inputs.analogChannelNames{i} = board_adc_channels(i).native_channel_name; + end + + else + session.inputs.nAnalogChannels = 0; + end + + % Auxiliary channels (accelerometer, etc.) + if ~isempty(aux_input_channels) + nAuxChannels = length(aux_input_channels); + session.inputs.nAuxChannels = nAuxChannels; + + for i = 1:nAuxChannels + session.inputs.auxChannelNames{i} = aux_input_channels(i).native_channel_name; + end + + else + session.inputs.nAuxChannels = 0; + end + + % Notes + if ~isempty(notes) && isstruct(notes) + % Handle different possible field names for notes + if isfield(notes, 'note1') + session.general.notes = notes.note1; + elseif isfield(notes, 'notes') + session.general.notes = notes.notes; + else + % Use first available field + noteFields = fieldnames(notes); + + if ~isempty(noteFields) + session.general.notes = notes.(noteFields{1}); + end + + end + + elseif ischar(notes) || isstring(notes) + % Handle case where notes is directly a string + session.general.notes = char(notes); + end + + % Set minimal extracellular fields (even though no spikes) + session.extracellular.nChannels = 0; % No amplifier channels in behavior-only + session.extracellular.nElectrodeGroups = 0; + session.extracellular.nSpikeGroups = 0; + session.extracellular.precision = 'int16'; % Standard precision for raw data + session.extracellular.leastSignificantBit = 0.195; % Intan default (in µV) + + disp([' Sampling rate: ', num2str(session.extracellular.sr), ' Hz']) + disp([' Digital channels: ', num2str(session.inputs.nDigitalChannels)]) + disp([' Analog channels: ', num2str(session.inputs.nAnalogChannels)]) + + catch ME + warning(['Failed to read Intan metadata: ', ME.message]); + session = setDefaultBehaviorMetadata(session, basepath); + end + +end + +function session = readOpenEphysBehaviorMetadata(session, basepath, settingsFile) + % readOpenEphysBehaviorMetadata - Extract metadata from OpenEphys settings.xml + % + % Parses OpenEphys XML configuration to extract recording parameters. + % Attempts to find sampling rate from SIGNALCHAIN->PROCESSOR->EDITOR hierarchy. + % Falls back to detectChannelsFromFiles() for channel detection since XML + % may not contain complete channel information. + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % settingsFile - Dir struct entry for settings.xml + % + % OUTPUTS: + % session - Session struct with OpenEphys metadata populated + % + % NOTES: + % - Default sampling rate is 30000 Hz if not found in XML + % - Uses xml2struct() for XML parsing + % - Calls detectChannelsFromFiles() to count actual data channels + % - Falls back to setDefaultBehaviorMetadata() on error + + try + % Parse XML file + settingsPath = fullfile(settingsFile.folder, settingsFile.name); + xmlStruct = xml2struct(settingsPath); + + % Extract sampling rate (typically in SIGNALCHAIN > PROCESSOR > EDITOR) + % OpenEphys default is 30000 Hz + session.extracellular.sr = 30000; % Default, will try to parse from XML + + % Try to extract actual sampling rate from XML + try + + if isfield(xmlStruct, 'SETTINGS') && isfield(xmlStruct.SETTINGS, 'SIGNALCHAIN') + signalChain = xmlStruct.SETTINGS.SIGNALCHAIN; + + if isfield(signalChain, 'PROCESSOR') + processors = signalChain.PROCESSOR; + + if ~iscell(processors) + processors = {processors}; + end + + for i = 1:length(processors) + + if isfield(processors{i}, 'Attributes') && isfield(processors{i}.Attributes, 'name') + + if contains(processors{i}.Attributes.name, 'Sources', 'IgnoreCase', true) + + if isfield(processors{i}, 'EDITOR') && isfield(processors{i}.EDITOR, 'Attributes') + + if isfield(processors{i}.EDITOR.Attributes, 'SampleRate') + session.extracellular.sr = str2double(processors{i}.EDITOR.Attributes.SampleRate); + end + + end + + end + + end + + end + + end + + end + + catch + warning('Could not parse sampling rate from settings.xml, using default 30000 Hz'); + end + + session.extracellular.srLfp = 1250; + + % Set minimal extracellular fields + session.extracellular.nChannels = 0; + session.extracellular.nElectrodeGroups = 0; + session.extracellular.nSpikeGroups = 0; + session.extracellular.precision = 'int16'; % Standard precision for raw data + session.extracellular.leastSignificantBit = 0.195; % Default (Intan = 0.195, Amplipex = 0.3815) + + % Detect digital/analog channels from actual files + session = detectChannelsFromFiles(session, basepath, nDigitalChannels, nAnalogChannels); + + disp([' Sampling rate: ', num2str(session.extracellular.sr), ' Hz']) + disp([' Digital channels detected: ', num2str(session.inputs.nDigitalChannels)]) + disp([' Analog channels detected: ', num2str(session.inputs.nAnalogChannels)]) + + catch ME + warning(['Failed to read OpenEphys metadata: ', ME.message]); + session = setDefaultBehaviorMetadata(session, basepath); + end + +end + +function session = setDefaultBehaviorMetadata(session, basepath) + % setDefaultBehaviorMetadata - Apply fallback metadata when no info files exist + % + % Used when neither info.rhd (Intan) nor settings.xml (OpenEphys) are found. + % Sets conservative default values and relies on detectChannelsFromFiles() + % to determine actual channel counts from data files. + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % + % OUTPUTS: + % session - Session struct with default metadata + % + % DEFAULTS: + % - Sampling rate: 20000 Hz (common Intan default) + % - LFP rate: 1250 Hz + % - No electrode groups or spike groups (behavior-only) + % - Channel counts determined from actual files + % + % NOTES: + % - Displays warning that user should verify parameters + % - Should only be used as last resort + + warning('Using default behavior metadata. Please verify sampling rate and channel counts!'); + + session.extracellular.sr = 20000; % Intan default + session.extracellular.srLfp = 1250; + session.extracellular.nChannels = 0; + session.extracellular.nElectrodeGroups = 0; + session.extracellular.nSpikeGroups = 0; + + % Detect channels from files + session = detectChannelsFromFiles(session, basepath, nDigitalChannels, nAnalogChannels); + +end + +function session = detectChannelsFromFiles(session, basepath, defaultDigitalChannels, defaultAnalogChannels) + % detectChannelsFromFiles - Determine channel counts from actual data files + % + % Searches for digitalin.dat, analogin.dat, and digitalout.dat files to + % determine what input channels are present. Uses configurable defaults + % for channel counts since files don't self-describe their channel count. + % + % INPUTS: + % session - Session struct being built + % basepath - Path to session directory + % defaultDigitalChannels - Number of digital input channels (default: 16 for Intan) + % defaultAnalogChannels - Number of analog input channels (default: 8 for Intan) + % + % OUTPUTS: + % session - Session struct with populated channel information: + % - inputs.nDigitalChannels, digitalChannels, digitalChannelNames + % - inputs.nAnalogChannels, analogChannels, analogChannelNames + % - inputs.hasDigitalOut + % + % FILE FORMATS: + % digitalin.dat - uint16, multiplexed digital channels + % analogin.dat - uint16, multiplexed analog channels + % digitalout.dat - uint16, digital output channels + % + % NOTES: + % - Searches recursively with dir(fullfile(basepath, '**', filename)) + % - Channel names are auto-generated (e.g., 'Digital-0', 'Analog-0') + % - Defaults can be overridden via function parameters + + % Check for digitalin.dat + digitalFiles = dir(fullfile(basepath, '**', 'digitalin.dat')); + + if ~isempty(digitalFiles) + digitalFile = digitalFiles(1); + fileInfo = dir(fullfile(digitalFile.folder, digitalFile.name)); + % Use the configurable default + session.inputs.nDigitalChannels = defaultDigitalChannels; + session.inputs.digitalChannels = 1:defaultDigitalChannels; + + % Generate default channel names + for i = 1:defaultDigitalChannels + session.inputs.digitalChannelNames{i} = sprintf('Digital-%d', i-1); + end + else + session.inputs.nDigitalChannels = 0; + end + + % Check for analogin.dat + analogFiles = dir(fullfile(basepath, '**', 'analogin.dat')); + + if ~isempty(analogFiles) + analogFile = analogFiles(1); + fileInfo = dir(fullfile(analogFile.folder, analogFile.name)); + % Use the configurable default + session.inputs.nAnalogChannels = defaultAnalogChannels; + session.inputs.analogChannels = 1:defaultAnalogChannels; + + % Generate default channel names + for i = 1:defaultAnalogChannels + session.inputs.analogChannelNames{i} = sprintf('Analog-%d', i-1); + end + else + session.inputs.nAnalogChannels = 0; + end + + % Check for digitalout.dat + digitalOutFiles = dir(fullfile(basepath, '**', 'digitalout.dat')); + + if ~isempty(digitalOutFiles) + session.inputs.hasDigitalOut = true; + else + session.inputs.hasDigitalOut = false; + end + +end + +function session = calculateSessionDuration(session, basepath, basename) + % calculateSessionDuration - Determine recording duration from data files + % + % Attempts to calculate total session duration using available files, + % prioritizing time.dat for accuracy, then falling back to data files. + % + % INPUTS: + % session - Session struct being built + % basepath - Path to session directory + % basename - Base name for session files (unused but kept for compatibility) + % + % OUTPUTS: + % session - Session struct with populated: + % - extracellular.nSamples (total samples) + % - general.duration (seconds) + % + % FILE PRIORITY: + % 1. time.dat - Most accurate (contains actual timestamps) + % 2. digitalin.dat - Good fallback (file size / 2 bytes = samples) + % 3. analogin.dat - Last resort (file size / (nChannels * 2 bytes) = samples) + % + % NOTES: + % - time.dat contains int32 timestamps in sample units + % - Duration = nSamples / samplingRate + % - Searches recursively for files in subdirectories + % - Displays warning if no usable files found + + % Try to use time.dat first + timeFiles = dir(fullfile(basepath, '**', 'time.dat')); + + if ~isempty(timeFiles) + timeFile = fullfile(timeFiles(1).folder, timeFiles(1).name); + + try + fid = fopen(timeFile, 'r'); + timeStamps = fread(fid, inf, 'int32'); + fclose(fid); + + if ~isempty(timeStamps) + % time.dat contains timestamps in samples + session.extracellular.nSamples = double(max(timeStamps)); + session.general.duration = session.extracellular.nSamples / session.extracellular.sr; + disp([' Duration (from time.dat): ', num2str(session.general.duration), ' seconds']) + return + end + + catch + warning('Could not read time.dat'); + end + + end + + % Fall back to digitalin.dat or analogin.dat + digitalFiles = dir(fullfile(basepath, '**', 'digitalin.dat')); + + if ~isempty(digitalFiles) + digitalFile = fullfile(digitalFiles(1).folder, digitalFiles(1).name); + fileInfo = dir(digitalFile); + % digitalin.dat: uint16, 2 bytes per sample + nSamples = fileInfo.bytes / 2; + session.extracellular.nSamples = nSamples; + session.general.duration = nSamples / session.extracellular.sr; + disp([' Duration (from digitalin.dat): ', num2str(session.general.duration), ' seconds']) + return + end + + % Try analogin.dat + analogFiles = dir(fullfile(basepath, '**', 'analogin.dat')); + + if ~isempty(analogFiles) + analogFile = fullfile(analogFiles(1).folder, analogFiles(1).name); + fileInfo = dir(analogFile); + % analogin.dat: uint16, 8 channels, 2 bytes per sample per channel + nSamples = fileInfo.bytes / (session.inputs.nAnalogChannels * 2); + session.extracellular.nSamples = nSamples; + session.general.duration = nSamples / session.extracellular.sr; + disp([' Duration (from analogin.dat): ', num2str(session.general.duration), ' seconds']) + return + end + + % If nothing found + warning('Could not calculate session duration. No time.dat, digitalin.dat, or analogin.dat found.'); + session.extracellular.nSamples = 0; + session.general.duration = 0; + +end + +function s = xml2struct(xmlFile) + % xml2struct - Convert XML file to MATLAB struct + % + % Recursively parses XML document into nested struct with element names + % as field names and attributes stored in .Attributes subfields. + % + % INPUTS: + % xmlFile - Path to XML file + % + % OUTPUTS: + % s - Struct representation of XML hierarchy + % + % EXAMPLE: + % For XML: ... + % Output: s.SETTINGS.PROCESSOR.Attributes.name = 'test' + % s.SETTINGS.PROCESSOR.Attributes.rate = '30000' + % + % NOTES: + % - Uses MATLAB's built-in xmlread() for parsing + % - Calls parseChildNodes() recursively for nested elements + % - Replaces '-' and ':' in element names with '_' + % - Element attributes stored in .Attributes substruct + % - Basic implementation; consider external tools for complex XML + + try + tree = xmlread(xmlFile); + s = parseChildNodes(tree); + catch ME + error(['Failed to parse XML file: ', ME.message]); + end + +end + +function children = parseChildNodes(theNode) + % parseChildNodes - Recursively parse XML DOM nodes into struct + % + % Helper function for xml2struct(). Traverses XML document tree and + % converts nodes and their attributes into nested MATLAB struct. + % + % INPUTS: + % theNode - XML DOM node object (from xmlread) + % + % OUTPUTS: + % children - Struct containing parsed child elements and attributes + % + % PARSING RULES: + % - Only ELEMENT_NODE types are processed (ignores text, comments, etc.) + % - Node names with '-' or ':' are converted to '_' for valid fieldnames + % - Node attributes stored in .Attributes substruct + % - Recursively processes nested child nodes + % + % NOTES: + % - Called internally by xml2struct() + % - Uses Java DOM API (getChildNodes, getAttributes, etc.) + + children = struct; + + if theNode.hasChildNodes + childNodes = theNode.getChildNodes; + numChildNodes = childNodes.getLength; + + for count = 1:numChildNodes + theChild = childNodes.item(count - 1); + + if theChild.getNodeType == theChild.ELEMENT_NODE + childName = char(theChild.getNodeName); + childName = strrep(childName, '-', '_'); + childName = strrep(childName, ':', '_'); + + % Get attributes + if theChild.hasAttributes + attributes = theChild.getAttributes; + numAttributes = attributes.getLength; + attributeStruct = struct; + + for i = 1:numAttributes + attr = attributes.item(i - 1); + attrName = char(attr.getName); + attrValue = char(attr.getValue); + attributeStruct.(attrName) = attrValue; + end + + children.(childName).Attributes = attributeStruct; + end + + % Recurse on child nodes + if theChild.hasChildNodes + childData = parseChildNodes(theChild); + + if ~isempty(fieldnames(childData)) + children.(childName) = childData; + end + + end + + end + + end + + end + +end + +function session = createBehaviorEpochs(session, basepath, basename) + % createBehaviorEpochs - Define temporal epochs for the session + % + % Creates epoch definitions from available information, prioritizing + % MergePoints.events.mat (from concatenateDats), then subsession folders, + % finally falling back to single default epoch. + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % basename - Base name for session files + % + % OUTPUTS: + % session - Session struct with populated epochs cell array + % Each epoch contains: name, startTime, stopTime, + % behavioralParadigm, environment, manipulation + % + % EPOCH SOURCES (in priority order): + % 1. MergePoints.events.mat - Most accurate (from concatenateDats) + % 2. Subsession folders (session_01, session_02, etc.) - Approximate + % 3. Single default epoch - Last resort + % + % BEHAVIORAL INFERENCE: + % - Attempts to infer behavioralParadigm from folder names + % - Keywords: 'sleep'/'rest' → Sleep, 'maze'/'task' → Maze, 'open'/'field' → Open field + % + % NOTES: + % - Skips update if existing epochs already match MergePoints + % - Warns user if MergePoints missing with multiple subsessions + % - Subsession epoch times are approximate until concatenateDats runs + + % Check for MergePoints file (created by concatenateDats) + mergePointsFile = fullfile(basepath, [basename, '.MergePoints.events.mat']); + + if exist(mergePointsFile, 'file') + + try + load(mergePointsFile, 'MergePoints'); + + % Clear existing epochs if they don't match MergePoints + shouldUpdate = true; + + if isfield(session, 'epochs') && ~isempty(session.epochs) + % Check if epochs already match MergePoints + if length(session.epochs) == size(MergePoints.foldernames, 2) + % Check if timestamps match + allMatch = true; + + for i = 1:length(session.epochs) + + if session.epochs{i}.startTime ~= MergePoints.timestamps(i, 1) || ... + session.epochs{i}.stopTime ~= MergePoints.timestamps(i, 2) + allMatch = false; + break; + end + + end + + if allMatch + disp(' Epochs already match MergePoints, skipping update'); + shouldUpdate = false; + end + + end + + end + + if shouldUpdate + disp('Creating epochs from MergePoints...'); + % Clear old epochs + session.epochs = {}; + + for i = 1:size(MergePoints.foldernames, 2) + session.epochs{i}.name = MergePoints.foldernames{i}; + session.epochs{i}.startTime = MergePoints.timestamps(i, 1); + session.epochs{i}.stopTime = MergePoints.timestamps(i, 2); + + % Add behavioral context if available + session.epochs{i}.behavioralParadigm = 'Unknown'; + session.epochs{i}.environment = 'Unknown'; + session.epochs{i}.manipulation = 'None'; + + % Try to infer from folder name + folderName = lower(MergePoints.foldernames{i}); + + if contains(folderName, 'sleep') || contains(folderName, 'rest') + session.epochs{i}.behavioralParadigm = 'Sleep'; + session.epochs{i}.environment = 'Sleep box'; + elseif contains(folderName, 'maze') || contains(folderName, 'task') || contains(folderName, 'session') + session.epochs{i}.behavioralParadigm = 'Maze'; + session.epochs{i}.environment = 'Maze'; + elseif contains(folderName, 'open') || contains(folderName, 'field') + session.epochs{i}.behavioralParadigm = 'Open field'; + session.epochs{i}.environment = 'Open field'; + end + end + + disp([' Created ', num2str(length(session.epochs)), ' epochs from MergePoints']); + end + + return + + catch ME + warning(['Could not load MergePoints: ', ME.message]); + end + + end + + % If no MergePoints, try to detect subsession folders + subsessionFolders = dir(fullfile(basepath, 'session*')); + subsessionFolders = subsessionFolders([subsessionFolders.isdir]); + + if length(subsessionFolders) >= 2 + % Multiple subsessions found - create epochs from them + disp(['Found ', num2str(length(subsessionFolders)), ' subsession folders']); + disp('WARNING: MergePoints.events.mat not found. Please run concatenateDats first!'); + disp('Creating placeholder epochs from subsession folders...'); + + session.epochs = {}; + cumulativeTime = 0; + + for i = 1:length(subsessionFolders) + session.epochs{i}.name = subsessionFolders(i).name; + session.epochs{i}.startTime = cumulativeTime; + + % Try to get duration from digitalin.dat in subsession folder + subDatFile = fullfile(basepath, subsessionFolders(i).name, 'digitalin.dat'); + + if exist(subDatFile, 'file') + fileInfo = dir(subDatFile); + % digitalin.dat: uint16, 2 bytes per sample + nSamples = fileInfo.bytes / 2; + duration = nSamples / session.extracellular.sr; + session.epochs{i}.stopTime = cumulativeTime + duration; + cumulativeTime = cumulativeTime + duration; + else + session.epochs{i}.stopTime = inf; + end + + % Add behavioral context + session.epochs{i}.behavioralParadigm = 'Unknown'; + session.epochs{i}.environment = 'Unknown'; + session.epochs{i}.manipulation = 'None'; + + % Infer from folder name + folderName = lower(subsessionFolders(i).name); + + if contains(folderName, 'sleep') || contains(folderName, 'rest') + session.epochs{i}.behavioralParadigm = 'Sleep'; + session.epochs{i}.environment = 'Sleep box'; + elseif contains(folderName, 'maze') || contains(folderName, 'task') || contains(folderName, 'session') + session.epochs{i}.behavioralParadigm = 'Maze'; + session.epochs{i}.environment = 'Maze'; + elseif contains(folderName, 'open') || contains(folderName, 'field') + session.epochs{i}.behavioralParadigm = 'Open field'; + session.epochs{i}.environment = 'Open field'; + end + + end + + disp(' NOTE: Epoch times are approximate. Run concatenateDats for accurate timestamps.'); + return + end + + % No MergePoints and no subsessions - check if we need a default epoch + if ~isfield(session, 'epochs') || isempty(session.epochs) + session = createDefaultEpoch(session, basename); + else + disp(' Using existing epochs (no MergePoints or subsessions found)'); + end + +end + +function session = createDefaultEpoch(session, basename) + % createDefaultEpoch - Create single epoch spanning entire session + % + % Used when no MergePoints or subsession folders exist. Creates one + % epoch covering the full recording duration with placeholder metadata. + % + % INPUTS: + % session - Partially filled session struct + % basename - Base name for session (used as epoch name) + % + % OUTPUTS: + % session - Session struct with single epoch: + % - name: basename + % - startTime: 0 + % - stopTime: session.general.duration (or inf if unknown) + % - behavioralParadigm: 'Unknown' + % - environment: 'Unknown' + % - manipulation: 'None' + % + % NOTES: + % - Used as fallback when no epoch information available + % - stopTime set to inf if duration not yet calculated + + session.epochs{1}.name = basename; + session.epochs{1}.startTime = 0; + + if isfield(session.general, 'duration') && session.general.duration > 0 + session.epochs{1}.stopTime = session.general.duration; + else + session.epochs{1}.stopTime = inf; % Will be updated when duration is known + end + + session.epochs{1}.behavioralParadigm = 'Unknown'; + session.epochs{1}.environment = 'Unknown'; + session.epochs{1}.manipulation = 'None'; + + disp(' Created default epoch for entire session'); + +end + +function session = extractSessionDateTime(session, basepath, ~) + % extractSessionDateTime - Determine recording date and time + % + % Attempts to extract recording timestamp from folder naming conventions + % or file metadata. Supports common Intan and OpenEphys naming patterns. + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % (unused) - Placeholder for basename (not needed) + % + % OUTPUTS: + % session - Session struct with populated: + % - general.date (format: 'YYYY-MM-DD') + % - general.time (format: 'HH:MM:SS') + % + % SUPPORTED FORMATS: + % Intan: basename_YYMMDD_HHMMSS (e.g., session_210315_143022) + % OpenEphys: basename_YYYY-MM-DD_HH-MM-SS (e.g., session_2021-03-15_14-30-22) + % + % FALLBACK ORDER: + % 1. Folder name pattern matching + % 2. .rhd file modification timestamp + % 3. Current date (last resort) + % + % NOTES: + % - Uses regex to parse common date/time patterns + % - Assumes YY dates are 2000+ + % - Warns when falling back to current date + + % Try to extract from folder name first (common formats) + pathParts = strsplit(basepath, filesep); + sessionFolder = pathParts{end}; + + % Common date/time patterns + % Intan: basename_YYMMDD_HHMMSS + % OpenEphys: basename_YYYY-MM-DD_HH-MM-SS + % General: basename_date + + % Try Intan format: _YYMMDD_HHMMSS + intanPattern = '_(\d{6})_(\d{6})'; + tokens = regexp(sessionFolder, intanPattern, 'tokens'); + + if ~isempty(tokens) + dateStr = tokens{1}{1}; % YYMMDD + timeStr = tokens{1}{2}; % HHMMSS + + % Parse date + year = 2000 + str2double(dateStr(1:2)); + month = str2double(dateStr(3:4)); + day = str2double(dateStr(5:6)); + + % Parse time + hour = str2double(timeStr(1:2)); + minute = str2double(timeStr(3:4)); + second = str2double(timeStr(5:6)); + + session.general.date = sprintf('%04d-%02d-%02d', year, month, day); + session.general.time = sprintf('%02d:%02d:%02d', hour, minute, second); + + disp([' Extracted date: ', session.general.date, ' ', session.general.time]); + return + end + + % Try OpenEphys format: _YYYY-MM-DD_HH-MM-SS + oePattern = '_(\d{4})-(\d{2})-(\d{2})_(\d{2})-(\d{2})-(\d{2})'; + tokens = regexp(sessionFolder, oePattern, 'tokens'); + + if ~isempty(tokens) + year = str2double(tokens{1}{1}); + month = str2double(tokens{1}{2}); + day = str2double(tokens{1}{3}); + hour = str2double(tokens{1}{4}); + minute = str2double(tokens{1}{5}); + second = str2double(tokens{1}{6}); + + session.general.date = sprintf('%04d-%02d-%02d', year, month, day); + session.general.time = sprintf('%02d:%02d:%02d', hour, minute, second); + + disp([' Extracted date: ', session.general.date, ' ', session.general.time]); + return + end + + % Fall back to file modification time + rhdFiles = dir(fullfile(basepath, '**', '*.rhd')); + + if ~isempty(rhdFiles) + fileDate = datetime(rhdFiles(1).date, 'InputFormat', 'dd-MMM-yyyy HH:mm:ss'); + session.general.date = datestr(fileDate, 'yyyy-mm-dd'); + session.general.time = datestr(fileDate, 'HH:MM:SS'); + disp([' Using file timestamp: ', session.general.date, ' ', session.general.time]); + return + end + + % Last resort: use current date + session.general.date = datestr(now, 'yyyy-mm-dd'); + disp([' Could not extract date from folder or files. Using current date: ', session.general.date]); + +end + +function session = addBehavioralPaths(session, basepath, basename) + % addBehavioralPaths - Detect and catalog behavioral data files + % + % Searches for tracking, video, DeepLabCut, and Optitrack files and + % adds their paths to the session struct. Useful for automatic discovery + % of processed behavioral data. + % + % INPUTS: + % session - Partially filled session struct + % basepath - Path to session directory + % basename - Base name for session files + % + % OUTPUTS: + % session - Session struct with populated fields: + % - behavioralTracking: Cell array of tracking file names + % - videos: Cell array of video file paths (relative) + % - deepLabCut: Cell array of DLC CSV file paths + % - optitrack: Cell array of Optitrack CSV file paths + % + % DETECTED FILES: + % Tracking: *.position.behavior.mat, *.tracking.behavior.mat, *.animal.behavior.mat + % Videos: *.mp4, *.avi, *.mov, *.mkv (recursive search) + % DeepLabCut: *DLC*.csv (recursive search) + % Optitrack: *optitrack*.csv (recursive search) + % + % NOTES: + % - Stores relative paths (from basepath) for portability + % - Searches recursively in subdirectories + % - Only called if 'addBehaviorFiles' parameter is true + % - Displays summary of found files + + % Check for tracking files + trackingFiles = { + [basename, '.position.behavior.mat'] + [basename, '.tracking.behavior.mat'] + [basename, '.animal.behavior.mat'] + }; + + for i = 1:length(trackingFiles) + + if exist(fullfile(basepath, trackingFiles{i}), 'file') + + if ~isfield(session, 'behavioralTracking') + session.behavioralTracking = {}; + end + + [~, name, ~] = fileparts(trackingFiles{i}); + session.behavioralTracking{end + 1} = name; + disp([' Found tracking file: ', trackingFiles{i}]); + end + + end + + % Check for video files + videoExtensions = {'*.mp4', '*.avi', '*.mov', '*.mkv'}; + + for i = 1:length(videoExtensions) + videoFiles = dir(fullfile(basepath, '**', videoExtensions{i})); + + if ~isempty(videoFiles) + + if ~isfield(session, 'videos') + session.videos = {}; + end + + for j = 1:length(videoFiles) + videoPath = fullfile(videoFiles(j).folder, videoFiles(j).name); + % Store relative path + session.videos{end + 1} = strrep(videoPath, [basepath, filesep], ''); + end + + disp([' Found ', num2str(length(videoFiles)), ' video files (', videoExtensions{i}, ')']); + end + + end + + % Check for DeepLabCut output + dlcFiles = dir(fullfile(basepath, '**', '*DLC*.csv')); + + if ~isempty(dlcFiles) + + if ~isfield(session, 'deepLabCut') + session.deepLabCut = {}; + end + + for i = 1:length(dlcFiles) + dlcPath = fullfile(dlcFiles(i).folder, dlcFiles(i).name); + session.deepLabCut{end + 1} = strrep(dlcPath, [basepath, filesep], ''); + end + + disp([' Found ', num2str(length(dlcFiles)), ' DeepLabCut files']); + end + + % Check for Optitrack CSV + optitrackFiles = dir(fullfile(basepath, '**', '*optitrack*.csv')); + + if ~isempty(optitrackFiles) + + if ~isfield(session, 'optitrack') + session.optitrack = {}; + end + + for i = 1:length(optitrackFiles) + optiPath = fullfile(optitrackFiles(i).folder, optitrackFiles(i).name); + session.optitrack{end + 1} = strrep(optiPath, [basepath, filesep], ''); + end + + disp([' Found ', num2str(length(optitrackFiles)), ' Optitrack files']); + end + +end