Add EMG-BIDS import/export support with multiple coordinate systems

.gitignore

@@ -2,3 +2,8 @@ EEG-BIDS_testcases
 *.asv
 *~
 *.asv
+
+CLAUDE.md
+.DS_Store
+.context/
+.rules/

bids_check_regular_sampling.m

@@ -0,0 +1,67 @@
+% BIDS_CHECK_REGULAR_SAMPLING - Check if EEG data has regular sampling
+%
+% Usage:
+%   [isRegular, avgFreq] = bids_check_regular_sampling(EEG)
+%   [isRegular, avgFreq] = bids_check_regular_sampling(EEG, tolerance)
+%
+% Inputs:
+%   EEG       - [struct] EEGLAB dataset structure
+%   tolerance - [float] acceptable deviation from regular sampling (default: 0.0001 = 0.01%)
+%
+% Outputs:
+%   isRegular - [boolean] true if sampling is regular within tolerance
+%   avgFreq   - [float] average sampling frequency in Hz
+%
+% Note:
+%   EDF and BDF formats require perfectly regular sampling. This function
+%   checks if data has irregular timestamps and calculates the average
+%   frequency for potential resampling.
+%
+% Authors: Seyed Yahya Shirazi, 2025
+
+function [isRegular, avgFreq] = bids_check_regular_sampling(EEG, tolerance)
+
+if nargin < 1
+    help bids_check_regular_sampling;
+    return;
+end
+
+if nargin < 2
+    tolerance = 0.01;  % 1% tolerance for regular sampling detection
+end
+
+if isempty(EEG.data)
+    error('EEG.data is empty');
+end
+
+if EEG.trials > 1
+    isRegular = true;
+    avgFreq = EEG.srate;
+    return;
+end
+
+if isfield(EEG, 'times') && length(EEG.times) > 1
+    intervals = diff(EEG.times);
+
+    if length(unique(intervals)) == 1
+        isRegular = true;
+        avgFreq = EEG.srate;
+        return;
+    end
+
+    avgInterval = mean(intervals);
+    maxDeviation = max(abs(intervals - avgInterval)) / avgInterval;
+
+    isRegular = maxDeviation < tolerance;
+    % Check if times are in seconds (continuous) or milliseconds (epoched)
+    if EEG.trials == 1 && avgInterval < 1
+        % Continuous data: times in seconds
+        avgFreq = 1 / avgInterval;
+    else
+        % Epoched data: times in milliseconds
+        avgFreq = 1000 / avgInterval;
+    end
+else
+    isRegular = true;
+    avgFreq = EEG.srate;
+end

bids_export.m

@@ -347,7 +347,7 @@ function bids_export(files, varargin)
     'rmtempfiles'  'string'  {'on' 'off'}    'on';
     'exportformat' 'string'  {'same' 'eeglab' 'edf' 'bdf'}    'eeglab';
     'individualEventsJson' 'string'  {'on' 'off'}    'off';
-    'modality'     'string'  {'ieeg' 'meg' 'eeg' 'auto'}       'auto';
+    'modality'     'string'  {'ieeg' 'meg' 'eeg' 'emg' 'auto'}       'auto';
     'README'       'string'  {}    '';
     'CHANGES'      'string'  {}    '' ;
     'copydata'     'integer' {}    [0 1]; % legacy, does nothing now
@@ -371,6 +371,13 @@ function bids_export(files, varargin)
     opt.SourceDatasets = opt.sourceDatasets;
 end
 
+if strcmpi(opt.modality, 'emg')
+    if strcmpi(opt.exportformat, 'eeglab')
+        opt.exportformat = 'bdf';
+        fprintf('EMG data detected: changing export format to bdf\n');
+    end
+end
+
 % deleting folder
 % ---------------
 fprintf('Exporting data to %s...\n', opt.targetdir);
@@ -872,6 +879,13 @@ function copy_data_bids_eeg(sIn, subjectStr, sess, fileStr, opt)
 [pathIn,fileInNoExt,ext] = fileparts(fileIn);
 fprintf('Processing file %s\n', fileIn);
 [EEG,opt.modality] = eeg_import(fileIn, 'modality', opt.modality, 'noevents', opt.noevents, 'importfunc', opt.importfunc);
+
+% Set default export format to EDF for EMG data (after modality is determined)
+if strcmpi(opt.modality, 'emg') && strcmpi(opt.exportformat, 'eeglab')
+    opt.exportformat = 'edf';
+    fprintf('Note: EMG data will be exported to EDF format (default for EMG)\n');
+end
+
 if contains(fileIn, '_bids_tmp_') && strcmpi(opt.rmtempfiles, 'on')
     fprintf('Deleting temporary file %s\n', fileIn);
     delete(fileIn);
@@ -892,15 +906,37 @@ function copy_data_bids_eeg(sIn, subjectStr, sess, fileStr, opt)
 fileOut = [fileBase '_' opt.modality ext];
 
 % select data subset
-EEG = eeg_selectsegment(EEG, 'eventtype', eventtype, 'eventindex', eventindex, 'timeoffset', timeoffset );        
+EEG = eeg_selectsegment(EEG, 'eventtype', eventtype, 'eventindex', eventindex, 'timeoffset', timeoffset );
 
 if ~isequal(opt.exportformat, 'same') || isequal(ext, '.set')
     % export data if necessary
     [filePathTmp,fileOutNoExt,~] = fileparts(fileOut);
     if isequal(opt.exportformat, 'eeglab')
         pop_saveset(EEG, 'filename', [ fileOutNoExt '.set' ], 'filepath', filePathTmp);
     else
-        pop_writeeeg(EEG, fullfile(filePathTmp, [ fileOutNoExt '.' opt.exportformat]), 'TYPE',upper(opt.exportformat));
+        % Check for regular sampling when exporting to EDF/BDF
+        if strcmpi(opt.exportformat, 'edf') || strcmpi(opt.exportformat, 'bdf')
+            [isRegular, avgFreq] = bids_check_regular_sampling(EEG);
+            if ~isRegular
+                error(['EDF/BDF export requires regular sampling. Your data has irregular sampling (avg %.2f Hz).\n' ...
+                       'Please resample your data before exporting:\n' ...
+                       '  EEG = pop_resample(EEG, %.0f);\n' ...
+                       'Then re-run bids_export.'], avgFreq, round(avgFreq));
+            end
+        end
+
+        % For EMG (and other modalities), save without events in EDF
+        % Events are saved separately in events.tsv
+        if strcmpi(opt.modality, 'emg') && (strcmpi(opt.exportformat, 'edf') || strcmpi(opt.exportformat, 'bdf'))
+            % Temporarily remove events before writing EDF
+            tmpEvents = EEG.event;
+            EEG.event = [];
+            pop_writeeeg(EEG, fullfile(filePathTmp, [ fileOutNoExt '.' opt.exportformat]), 'TYPE',upper(opt.exportformat));
+            % Restore events for events.tsv export
+            EEG.event = tmpEvents;
+        else
+            pop_writeeeg(EEG, fullfile(filePathTmp, [ fileOutNoExt '.' opt.exportformat]), 'TYPE',upper(opt.exportformat));
+        end
     end
 else
     % copy the file
@@ -966,18 +1002,30 @@ function copy_data_bids_eeg(sIn, subjectStr, sess, fileStr, opt)
     EEG=pop_chanedit(EEG, 'lookup', opt.chanlookup);
 end
 
+% Set datatype before writing channel/electrode files
+if strcmpi(opt.modality, 'eeg')
+    EEG.etc.datatype = 'eeg';
+elseif strcmpi(opt.modality, 'ieeg')
+    EEG.etc.datatype = 'ieeg';
+elseif strcmpi(opt.modality, 'meg')
+    EEG.etc.datatype = 'meg';
+elseif strcmpi(opt.modality, 'emg')
+    EEG.etc.datatype = 'emg';
+end
+
 % Write electrode file information (electrodes.tsv and coordsystem.json)
 bids_writechanfile(EEG, fileOutRed);
-bids_writeelectrodefile(EEG, fileOutRed, 'export', opt.elecexport);
+bids_writeelectrodefile(EEG, fileOutRed, 'export', opt.elecexport, 'rootdir', opt.targetdir);
+
+% Write modality-specific info files
 if strcmpi(opt.modality, 'eeg')
     bids_writetinfofile(EEG, tInfo, notes, fileOutRed);
-    EEG.etc.datatype = 'eeg';
 elseif strcmpi(opt.modality, 'ieeg')
     bids_writeieegtinfofile(EEG, tInfo, notes, fileOutRed);
-    EEG.etc.datatype = 'ieeg';
 elseif strcmpi(opt.modality, 'meg')
     bids_writemegtinfofile(EEG, tInfo, notes, fileOutRed);
-    EEG.etc.datatype = 'meg';
+elseif strcmpi(opt.modality, 'emg')
+    bids_writeemgtinfofile(EEG, tInfo, notes, fileOutRed);
 end
 
 % write channel information

bids_getinfofromfolder.m

@@ -27,11 +27,12 @@
 % along with this program; if not, to the Free Software
 % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
-function [tasklist,sessions,runs] = bids_getinfofromfolder(bidsFolder)
+function [tasklist,sessions,runs,recordings] = bids_getinfofromfolder(bidsFolder)
 
 tasklist = {};
 sessions = {};
 runs     = {};
+recordings = {};
 files = dir(bidsFolder);
 [files(:).folder] = deal(bidsFolder);
 %fprintf('Scanning %s\n', bidsFolder);
@@ -41,10 +42,11 @@
             sessions = union(sessions, { files(iFile).name });
         end
 
-        [tasklistTmp,sessionTmp,runsTmp] = bids_getinfofromfolder(fullfile(files(iFile).folder, fullfile(files(iFile).name)));
+        [tasklistTmp,sessionTmp,runsTmp,recordingsTmp] = bids_getinfofromfolder(fullfile(files(iFile).folder, fullfile(files(iFile).name)));
         tasklist = union(tasklist, tasklistTmp);
         sessions = union(sessions, sessionTmp);
         runs     = union(runs    , runsTmp);
+        recordings = union(recordings, recordingsTmp);
      else
         if (~isempty(strfind(files(iFile).name, 'eeg')) || ~isempty(strfind(files(iFile).name, 'meg'))) && ~isempty(strfind(files(iFile).name, '_task'))
             pos = strfind(files(iFile).name, '_task');
@@ -53,12 +55,29 @@
             newTask = tmpStr(1:underS(1)-1);
             tasklist = union( tasklist, { newTask });
         end
-        if (~isempty(strfind(files(iFile).name, 'eeg')) || ~isempty(strfind(files(iFile).name, 'meg')))  && ~isempty(strfind(files(iFile).name, '_run'))
+        if (~isempty(strfind(files(iFile).name, 'eeg')) || ~isempty(strfind(files(iFile).name, 'meg')) || ~isempty(strfind(files(iFile).name, 'emg')))  && ~isempty(strfind(files(iFile).name, '_run'))
             pos = strfind(files(iFile).name, '_run');
             tmpStr = files(iFile).name(pos+5:end);
             underS = find(tmpStr == '_');
             newRun = tmpStr(1:underS(1)-1);
             runs = union( runs, { newRun } );
         end
+        if (~isempty(strfind(files(iFile).name, 'eeg')) || ~isempty(strfind(files(iFile).name, 'meg')) || ~isempty(strfind(files(iFile).name, 'emg')) || ~isempty(strfind(files(iFile).name, 'ieeg'))) && ~isempty(strfind(files(iFile).name, '_recording-'))
+            pos = strfind(files(iFile).name, '_recording-');
+            tmpStr = files(iFile).name(pos+11:end);
+            underS = find(tmpStr == '_');
+            if ~isempty(underS)
+                newRecording = tmpStr(1:underS(1)-1);
+            else
+                % recording is last entity before extension
+                dotS = find(tmpStr == '.');
+                if ~isempty(dotS)
+                    newRecording = tmpStr(1:dotS(1)-1);
+                else
+                    newRecording = tmpStr;
+                end
+            end
+            recordings = union( recordings, { newRecording } );
+        end
     end
 end

bids_importchanlocs.m

@@ -55,6 +55,10 @@
             chanlocs(iChan-1).X = elecData{iChan,2};
             chanlocs(iChan-1).Y = elecData{iChan,3};
             chanlocs(iChan-1).Z = elecData{iChan,4};
+            % Import coordinate_system (5th column) if present (EMG)
+            if size(elecData,2) >= 5 && ~isempty(elecData{iChan,5}) && ~strcmpi(elecData{iChan,5}, 'n/a')
+                chanlocs(iChan-1).coordinate_system = elecData{iChan,5};
+            end
         end
     end