Authored by Anthony Boyer

Re-work of spikedetection so it can deal with inconsistent sampling frequencies.

@@ -2,8 +2,8 @@ function ImaGIN_SpikesDetection(S) @@ -2,8 +2,8 @@ function ImaGIN_SpikesDetection(S)
2 % This function extracts and concatenates the baselines from cropped, qualitatively controlled, stimulation files. 2 % This function extracts and concatenates the baselines from cropped, qualitatively controlled, stimulation files.
3 % It runs Delphos' spike detection on extracted baselines after bipolar montage computation. 3 % It runs Delphos' spike detection on extracted baselines after bipolar montage computation.
4 % It then computes corrected spike rates for each bipolar recording channel. 4 % It then computes corrected spike rates for each bipolar recording channel.
5 -% /!\ Cropped files CAN have a variable number of channels and thus the monopolar <-> bipolar mapping has to be established for each cropped file  
6 -% and taken into account during baseline concatenation for consistency. /!\ 5 +% /!\ Cropped files CAN have a variable number of channels AND different sampling rates, the monopolar <-> bipolar mapping has to be established for each cropped file
  6 +% and taken into account during baseline concatenation for consistency AND cropped files with the same sampling rate have to be concatenated together. /!\
7 7
8 %% Vars parsing 8 %% Vars parsing
9 stims_files = S.stims_files; 9 stims_files = S.stims_files;
@@ -88,58 +88,134 @@ function ImaGIN_SpikesDetection(S) @@ -88,58 +88,134 @@ function ImaGIN_SpikesDetection(S)
88 end 88 end
89 end 89 end
90 90
91 - %% Create the new monopolar baseline object  
92 - for nn=1:numel(global_montage)  
93 - channel{nn} = global_montage{nn};  
94 - timeseries(nn,:) = [channels.(global_montage{nn}).timeserie{:}];  
95 - zeroed(nn,:) = [channels.(global_montage{nn}).zeroed{:}];  
96 - all_coordinates = cell2mat(channels.(global_montage{nn}).coordinates(:));  
97 - [lines cols] = find(~isnan(all_coordinates));  
98 - coordinates_with_values = all_coordinates(unique(lines),:);  
99 - if isempty(coordinates_with_values)  
100 - coordinates(nn,:) = NaN(1,3);  
101 - elseif sum(all(coordinates_with_values == coordinates_with_values(1,:))) == 3  
102 - coordinates(nn,:) = coordinates_with_values(1,:); 91 + %% Create monopolar baseline objects
  92 + unique_srates = unique(sampling);
  93 +
  94 + for ss=1:numel(unique_srates)
  95 + clear srate labels timeseries zeroed coordinates;
  96 + srate = unique_srates(ss);
  97 + for nn=1:numel(global_montage)
  98 + labels{nn} = global_montage{nn};
  99 + timeseries(nn,:) = [channels.(global_montage{nn}).timeserie{find(sampling==srate)}];
  100 + zeroed(nn,:) = [channels.(global_montage{nn}).zeroed{find(sampling==srate)}];
  101 + all_coordinates = cell2mat(channels.(global_montage{nn}).coordinates(:));
  102 + [lines cols] = find(~isnan(all_coordinates));
  103 + coordinates_with_values = all_coordinates(unique(lines),:);
  104 + if isempty(coordinates_with_values)
  105 + coordinates(nn,:) = NaN(1,3);
  106 + elseif sum(all(coordinates_with_values == coordinates_with_values(1,:))) == 3
  107 + coordinates(nn,:) = coordinates_with_values(1,:);
  108 + else
  109 + error('Inconsistent electrodes coordinates between cropped files.');
  110 + end
  111 + end
  112 + rates_tables{ss} = spikes_from_monopolar(srate,labels,timeseries,zeroed,coordinates,path_out);
  113 + end
  114 +
  115 + rates_labels = {};
  116 + rates_scores = [];
  117 + for rr=1:numel(rates_tables)
  118 + rates_labels = [rates_labels rates_tables{rr}{:,1}];
  119 + current_rates = rates_tables{rr}{:,2};
  120 + current_rates(isinf(current_rates)) = NaN;
  121 + rates_scores = [rates_scores current_rates];
  122 + end
  123 +
  124 + for ll=1:size(rates_labels,1)
  125 + if all(strcmp(rates_labels(ll,:),rates_labels{ll,1})==1)
  126 + bipolar_labels{ll} = rates_labels{ll,1};
103 else 127 else
104 - error('Inconsistent electrodes coordinates between cropped files.');  
105 - end  
106 - end 128 + error('Inconsistent labels between bipolar baselines');
  129 + end
  130 + end
  131 + mean_rates = mean(rates_scores,2,'omitnan');
  132 +
  133 + % Store rate table in an individual files for easy extraction
  134 + rate_table_file = 'SPK_bBaseline_rates.mat';
  135 + spike_rates = [bipolar_labels',num2cell(mean_rates)];
  136 + save(fullfile(path_out,rate_table_file), 'spike_rates');
  137 +
  138 + % Ending step
  139 + set_final_status('OK');
  140 + disp('Done');
  141 +
  142 +end
  143 +
  144 +function channels_map = monopolar_to_bipolar_mapping(monopolar_file,path_out)
  145 + mono_stimulation_obj = spm_eeg_load(monopolar_file);
  146 + monopolar_sensors = sensors(mono_stimulation_obj,'EEG');
  147 + [mono_path,mono_name,mono_ext] = fileparts(monopolar_file);
  148 + Sbp.Fname = monopolar_file;
  149 + Sbp.FileOut = fullfile(path_out,['b' mono_name]);
  150 + bp_stimulation_obj = ImaGIN_BipolarMontage(Sbp);
  151 + for cc=1:nchannels(bp_stimulation_obj)
  152 + bipolar_chan = chanlabels(bp_stimulation_obj,cc);
  153 + mono_chans = regexp(bipolar_chan,'-','split');
  154 + channels_map(cc).bp_chan = bipolar_chan;
  155 + channels_map(cc).bp_idx = cc;
  156 + channels_map(cc).mono_chan1 = mono_chans{1}{1};
  157 + channels_map(cc).mono_chan2 = mono_chans{1}{2};
  158 + % Get the monopolar channels idxes.
  159 + chan1_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{1}));
  160 + chan2_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{2}));
  161 + % Test and correct for duplicated labels in the monopolar object
  162 + if numel(chan1_idx) == 0
  163 + error('Could not find the monopolar channel from the bipolar montage.');
  164 + elseif numel(chan1_idx) == 1 % If 1 match, as expected, do nothing
  165 + %pass
  166 + else
  167 + error('2 or more channels with the same label')
  168 + end
  169 + % Same logic as for chan1_idx
  170 + if numel(chan2_idx) == 0
  171 + error('Could not find the monopolar channel from the bipolar montage.');
  172 + elseif numel(chan2_idx) == 1
  173 + %pass
  174 + else
  175 + error('2 or more channels with the same label')
  176 + end
  177 + channels_map(cc).mono_idx1 = chan1_idx;
  178 + channels_map(cc).mono_idx2 = chan2_idx;
  179 + end
  180 + delete([Sbp.FileOut '.mat']);delete([Sbp.FileOut '.dat']);delete([Sbp.FileOut '_log.txt']);
  181 +end
  182 +
  183 +function rates_table = spikes_from_monopolar(srate,labels,timeseries,zeroed,coordinates,path_out)
  184 +
107 % Create the new sensor 185 % Create the new sensor
108 new_sensors_info.balance.current = 'none'; 186 new_sensors_info.balance.current = 'none';
109 new_sensors_info.chanpos = coordinates; 187 new_sensors_info.chanpos = coordinates;
110 - new_sensors_info.chantype = cell(1,numel(global_montage)); new_sensors_info.chantype(:) = {'eeg'};  
111 - new_sensors_info.chanunit = cell(1,numel(global_montage)); new_sensors_info.chanunit(:) = {'V'}; 188 + new_sensors_info.chantype = cell(1,numel(labels)); new_sensors_info.chantype(:) = {'eeg'};
  189 + new_sensors_info.chanunit = cell(1,numel(labels)); new_sensors_info.chanunit(:) = {'V'};
112 new_sensors_info.elecpos = coordinates; 190 new_sensors_info.elecpos = coordinates;
113 - new_sensors_info.label = channel; 191 + new_sensors_info.label = labels;
114 new_sensors_info.type = 'ctf'; 192 new_sensors_info.type = 'ctf';
115 new_sensors_info.unit = 'mm'; 193 new_sensors_info.unit = 'mm';
  194 +
116 % Create the new files 195 % Create the new files
117 - file_header = fullfile(path_out,'monopolar_baseline.mat');  
118 - file_data = fullfile(path_out,'monopolar_baseline.dat'); 196 + file_header = fullfile(path_out,['monopolar_baseline_' num2str(srate) 'Hz.mat']);
  197 + file_data = fullfile(path_out,['monopolar_baseline_' num2str(srate) 'Hz.dat']);
119 if exist(file_header, 'file') == 2 198 if exist(file_header, 'file') == 2
120 delete(file_header); 199 delete(file_header);
121 end 200 end
122 if exist(file_data, 'file') == 2 201 if exist(file_data, 'file') == 2
123 delete(file_data); 202 delete(file_data);
124 end 203 end
  204 +
125 % Create the new object 205 % Create the new object
126 D = meeg(size(timeseries,1),size(timeseries,2),1); 206 D = meeg(size(timeseries,1),size(timeseries,2),1);
127 - D = fname(D,'monopolar_baseline'); 207 + D = fname(D,['monopolar_baseline_' num2str(srate) 'Hz']);
128 D = path(D,path_out); 208 D = path(D,path_out);
129 - if all(sampling == sampling(1))  
130 - D = fsample(D,sampling(1));  
131 - else  
132 - error('Inconsistent sampling frequency across cropped files.');  
133 - end 209 + D = fsample(D,srate);
134 D = sensors(D,'EEG',new_sensors_info) 210 D = sensors(D,'EEG',new_sensors_info)
135 D = blank(D) ; 211 D = blank(D) ;
136 D = link(D,file_data); 212 D = link(D,file_data);
137 D(:,:,1) = timeseries; 213 D(:,:,1) = timeseries;
138 - save(D);  
139 - 214 + save(D);
  215 +
140 %% Compute bipolar montage on baseline object to get the bipolar baseline object 216 %% Compute bipolar montage on baseline object to get the bipolar baseline object
141 Sbp.Fname = file_header; 217 Sbp.Fname = file_header;
142 - Sbp.FileOut = fullfile(path_out,'bipolar_baseline'); 218 + Sbp.FileOut = fullfile(path_out,['bipolar_baseline_' num2str(srate) 'Hz']);
143 bp_baseline_obj = ImaGIN_BipolarMontage(Sbp); 219 bp_baseline_obj = ImaGIN_BipolarMontage(Sbp);
144 global_mapping = monopolar_to_bipolar_mapping(Sbp.Fname,path_out); 220 global_mapping = monopolar_to_bipolar_mapping(Sbp.Fname,path_out);
145 221
@@ -179,53 +255,6 @@ function ImaGIN_SpikesDetection(S) @@ -179,53 +255,6 @@ function ImaGIN_SpikesDetection(S)
179 spkFile.markers = spks_table; 255 spkFile.markers = spks_table;
180 spkFile.baseline = ['Total duration of ' num2str((nsamples(bp_baseline_obj)/fsample(bp_baseline_obj))/60) ' minutes']; 256 spkFile.baseline = ['Total duration of ' num2str((nsamples(bp_baseline_obj)/fsample(bp_baseline_obj))/60) ' minutes'];
181 spkFile.comment = 'Spike rate is number of spikes per minute whithin a bipolar channel'; 257 spkFile.comment = 'Spike rate is number of spikes per minute whithin a bipolar channel';
182 - save(fullfile(path_out,spkFileName), 'spkFile');  
183 -  
184 - % Store rate table in an individual files for easy extraction  
185 - rate_table_file = 'SPK_bBaseline_rates.mat';  
186 - spike_rates = [chans,num2cell(rates)];  
187 - save(fullfile(path_out,rate_table_file), 'spike_rates');  
188 -  
189 - % Ending step  
190 - set_final_status('OK');  
191 - disp('Done');  
192 -end  
193 -  
194 -function channels_map = monopolar_to_bipolar_mapping(monopolar_file,path_out)  
195 - mono_stimulation_obj = spm_eeg_load(monopolar_file);  
196 - monopolar_sensors = sensors(mono_stimulation_obj,'EEG');  
197 - [mono_path,mono_name,mono_ext] = fileparts(monopolar_file);  
198 - Sbp.Fname = monopolar_file;  
199 - Sbp.FileOut = fullfile(path_out,['b' mono_name]);  
200 - bp_stimulation_obj = ImaGIN_BipolarMontage(Sbp);  
201 - for cc=1:nchannels(bp_stimulation_obj)  
202 - bipolar_chan = chanlabels(bp_stimulation_obj,cc);  
203 - mono_chans = regexp(bipolar_chan,'-','split');  
204 - channels_map(cc).bp_chan = bipolar_chan;  
205 - channels_map(cc).bp_idx = cc;  
206 - channels_map(cc).mono_chan1 = mono_chans{1}{1};  
207 - channels_map(cc).mono_chan2 = mono_chans{1}{2};  
208 - % Get the monopolar channels idxes.  
209 - chan1_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{1}));  
210 - chan2_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{2}));  
211 - % Test and correct for duplicated labels in the monopolar object  
212 - if numel(chan1_idx) == 0  
213 - error('Could not find the monopolar channel from the bipolar montage.');  
214 - elseif numel(chan1_idx) == 1 % If 1 match, as expected, do nothing  
215 - %pass  
216 - else  
217 - error('2 or more channels with the same label')  
218 - end  
219 - % Same logic as for chan1_idx  
220 - if numel(chan2_idx) == 0  
221 - error('Could not find the monopolar channel from the bipolar montage.');  
222 - elseif numel(chan2_idx) == 1  
223 - %pass  
224 - else  
225 - error('2 or more channels with the same label')  
226 - end  
227 - channels_map(cc).mono_idx1 = chan1_idx;  
228 - channels_map(cc).mono_idx2 = chan2_idx;  
229 - end  
230 - delete([Sbp.FileOut '.mat']);delete([Sbp.FileOut '.dat']);delete([Sbp.FileOut '_log.txt']); 258 + results_file = fullfile(path_out,spkFileName);
  259 + save(results_file, 'spkFile');
231 end 260 end