Authored by Anthony Boyer

New implementation of SpikeDetection (RAM usage) with new parameters.

... ... @@ -697,7 +697,7 @@ N = size(tf_stat,1);
if N > 16000;
decimate = floor(linspace(Fs,N-Fs,15000));
else
decimate = Fs:N-Fs;
decimate = floor(Fs:N-Fs);
end
tf_stat = tf_stat(decimate,:);
K = 1.5;
... ...
... ... @@ -38,7 +38,7 @@ elec= sensors(D,'eeg'); % add channels without positions into bad channels
pos = elec.elecpos;
Sens = elec.label;
chanLbs = D.chanlabels;
idxNaN = [] % find(isnan(pos(:,1))); We now assume channels without coordinates are not bad (CRFs older than 2009 do not have T1pre coordinates)
idxNaN = []; % find(isnan(pos(:,1))); We now assume channels without coordinates are not bad (CRFs older than 2009 do not have T1pre coordinates)
if ~isempty(idxNaN)
bIdx = [bIdx(:);idxNaN(:)];
bIdx = sort(unique(bIdx));
... ...
... ... @@ -18,8 +18,6 @@ function ImaGIN_SpikesDetection(S)
D = spm_eeg_load(stims_files{ff});
% Get the samping frequency
sampling(ff) = fsample(D);
% Get the coordinates
sensors_info = sensors(D,'EEG');
% Get the list of bad monopolar channels
all_channels_idxes = indchantype(D,'ALL');
good_eeg_channels_idxes = indchantype(D,'EEG','GOOD');
... ... @@ -48,6 +46,8 @@ function ImaGIN_SpikesDetection(S)
file_zeroed(all_bad_channels_idxes,:,1) = zeros(numel(all_bad_channels_idxes),numel(negative_time_idxes));
% Update the channels struct
eeg_chans_idxes = indchantype(D,'EEG');
% Sensor object to get coordinates
sensors_info = sensors(D,'EEG');
for cc=1:numel(eeg_chans_idxes)
channel_name = chanlabels(D,eeg_chans_idxes(cc));
if isvarname(channel_name{:})
... ... @@ -92,10 +92,10 @@ function ImaGIN_SpikesDetection(S)
%% Create monopolar baseline objects
unique_srates = unique(sampling);
for ss=1:numel(unique_srates)
clear srate labels timeseries zeroed coordinates;
srate = unique_srates(ss);
chunk_table = {};
for nn=1:numel(global_montage)
labels{nn} = global_montage{nn};
timeseries(nn,:) = [channels.(global_montage{nn}).timeserie{find(sampling==srate)}];
... ... @@ -111,14 +111,45 @@ function ImaGIN_SpikesDetection(S)
error('Inconsistent electrodes coordinates between cropped files.');
end
end
rates_tables{ss} = spikes_from_monopolar(srate,labels,timeseries,zeroed,coordinates,path_out);
% Split time series into chunks to solve performance (RAM) issues.
chunk_length = 300*srate; % 5 min chunks.
chunk_nb = fix(size(timeseries,2)/chunk_length);
if chunk_nb>0
for kk=1:chunk_nb
idx_start = (kk-1)*chunk_length;
idx_end = (kk)*chunk_length;
chunk_timeseries = timeseries(:,idx_start:idx_end);
chunk_zeroed = zeroed(:,idx_start:idx_end);
chunk_table{end+1} = spikes_from_monopolar(srate,labels,chunk_timeseries,chunk_zeroed,coordinates,path_out);
end
if idx_end<size(timeseries,2) && (size(timeseries,2)-idx_end)>srate
chunk_timeseries = timeseries(:,idx_end:end);
chunk_zeroed = zeroed(:,idx_end:end);
chunk_table{end+1} = spikes_from_monopolar(srate,labels,chunk_timeseries,chunk_zeroed,coordinates,path_out);
end
else
chunk_table{end+1} = spikes_from_monopolar(srate,labels,timeseries,zeroed,coordinates,path_out);
end
% Add the spikes from all chunks
tot_spikes = zeros(numel(chunk_table{1}{:,1}),1);
tot_time = zeros(numel(chunk_table{1}{:,1}),1);
for rr=1:numel(chunk_table)
current_spike_count = chunk_table{rr}{:,2};
current_recording_time = chunk_table{rr}{:,3};
current_spike_count(isinf(current_spike_count)) = 0;
current_spike_count(isnan(current_spike_count)) = 0;
tot_spikes = tot_spikes+current_spike_count;
tot_time = tot_time+current_recording_time;
end
rate_table=table(chunk_table{1}{:,1},round(tot_spikes./tot_time,1));
rates_tables{ss} = rate_table;
end
rates_labels = {};
rates_scores = [];
for rr=1:numel(rates_tables)
rates_labels = [rates_labels rates_tables{rr}{:,1}];
current_rates = rates_tables{rr}{:,2};
for ss=1:numel(rates_tables)
rates_labels = [rates_labels rates_tables{ss}{:,1}];
current_rates = rates_tables{ss}{:,2};
current_rates(isinf(current_rates)) = NaN;
rates_scores = [rates_scores current_rates];
end
... ... @@ -133,7 +164,7 @@ function ImaGIN_SpikesDetection(S)
mean_rates = mean(rates_scores,2,'omitnan');
% Store rate table in an individual files for easy extraction
rate_table_file = 'SPK_bBaseline_rates.mat';
rate_table_file = 'SPK_bBaseline_rates_updated_parameters.mat';
spike_rates = [bipolar_labels',num2cell(mean_rates)];
save(fullfile(path_out,rate_table_file), 'spike_rates');
... ... @@ -143,41 +174,29 @@ function ImaGIN_SpikesDetection(S)
end
function channels_map = monopolar_to_bipolar_mapping(monopolar_file,path_out)
function stim_mapping = monopolar_to_bipolar_mapping(monopolar_file,path_out)
mono_stimulation_obj = spm_eeg_load(monopolar_file);
monopolar_sensors = sensors(mono_stimulation_obj,'EEG');
[mono_path,mono_name,mono_ext] = fileparts(monopolar_file);
[~,mono_name,~] = fileparts(monopolar_file);
Sbp.Fname = monopolar_file;
Sbp.FileOut = fullfile(path_out,['b' mono_name]);
bp_stimulation_obj = ImaGIN_BipolarMontage(Sbp);
for cc=1:nchannels(bp_stimulation_obj)
bipolar_chan = chanlabels(bp_stimulation_obj,cc);
mono_chans = regexp(bipolar_chan,'-','split');
channels_map(cc).bp_chan = bipolar_chan;
channels_map(cc).bp_idx = cc;
channels_map(cc).mono_chan1 = mono_chans{1}{1};
channels_map(cc).mono_chan2 = mono_chans{1}{2};
stim_mapping(cc).bp_chan = bipolar_chan;
stim_mapping(cc).bp_idx = cc;
stim_mapping(cc).mono_chan1 = mono_chans{1}{1};
stim_mapping(cc).mono_chan2 = mono_chans{1}{2};
% Get the monopolar channels idxes.
chan1_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{1}));
chan2_idx = find(strcmp(monopolar_sensors.label,mono_chans{1}{2}));
% Test and correct for duplicated labels in the monopolar object
if numel(chan1_idx) == 0
error('Could not find the monopolar channel from the bipolar montage.');
elseif numel(chan1_idx) == 1 % If 1 match, as expected, do nothing
%pass
else
error('2 or more channels with the same label')
% Test if exactly 1 match in the monopolar object
if numel(chan1_idx) ~=1 || numel(chan2_idx) ~= 1
error('Could notidentifyh monopolar channeslused forr bipolar montage computation.');
end
% Same logic as for chan1_idx
if numel(chan2_idx) == 0
error('Could not find the monopolar channel from the bipolar montage.');
elseif numel(chan2_idx) == 1
%pass
else
error('2 or more channels with the same label')
end
channels_map(cc).mono_idx1 = chan1_idx;
channels_map(cc).mono_idx2 = chan2_idx;
stim_mapping(cc).mono_idx1 = chan1_idx;
stim_mapping(cc).mono_idx2 = chan2_idx;
end
delete([Sbp.FileOut '.mat']);delete([Sbp.FileOut '.dat']);delete([Sbp.FileOut '_log.txt']);
end
... ... @@ -187,23 +206,26 @@ function rates_table = spikes_from_monopolar(srate,labels,timeseries,zeroed,coor
% Create the new sensor
new_sensors_info.balance.current = 'none';
new_sensors_info.chanpos = coordinates;
new_sensors_info.elecpos = coordinates;
new_sensors_info.chantype = cell(1,numel(labels)); new_sensors_info.chantype(:) = {'eeg'};
new_sensors_info.chanunit = cell(1,numel(labels)); new_sensors_info.chanunit(:) = {'V'};
new_sensors_info.elecpos = coordinates;
new_sensors_info.label = labels;
new_sensors_info.type = 'ctf';
new_sensors_info.unit = 'mm';
% Create the new files
file_header = fullfile(path_out,['monopolar_baseline_' num2str(round(srate)) 'Hz.mat']);
file_data = fullfile(path_out,['monopolar_baseline_' num2str(round(srate)) 'Hz.dat']);
file_data = fullfile(path_out,['monopolar_baseline_' num2str(round(srate)) 'Hz.dat']);
file_log = fullfile(path_out,['monopolar_baseline_' num2str(round(srate)) 'Hz_log.txt']);
if exist(file_header, 'file') == 2
delete(file_header);
end
if exist(file_data, 'file') == 2
delete(file_data);
end
end
if exist(file_log, 'file') == 2
delete(file_data);
end
% Create the new object
D = meeg(size(timeseries,1),size(timeseries,2),1);
D = fname(D,['monopolar_baseline_' num2str(round(srate)) 'Hz']);
... ... @@ -230,7 +252,10 @@ function rates_table = spikes_from_monopolar(srate,labels,timeseries,zeroed,coor
% Run Delphos' spike detection on baselines
fs = bp_baseline_obj.fsample;
results = Delphos_detector(bp_baseline_obj,chanlabels(bp_baseline_obj),'SEEG',fs,{'Spk'},[],[],40,[]);
%old_param_thr = [1.4,10,1.3,11];
%results = Delphos_detector(bp_baseline_obj,chanlabels(bp_baseline_obj),'SEEG',fs,{'Spk'},[],[],40,[]); % Old parameters
new_param_thr = [1.4,10,1.5,9];
results = Delphos_detector(bp_baseline_obj,chanlabels(bp_baseline_obj),'SEEG',fs,{'Spk'},[],[],200,new_param_thr);
% Calculate and store corrected spike rates in a table
total_samples = nsamples(bp_baseline_obj);
... ... @@ -240,23 +265,10 @@ function rates_table = spikes_from_monopolar(srate,labels,timeseries,zeroed,coor
total_times(cc) = (total_samples - zeroed_samples) /fs/60; % Total recording time in minutes
end
chans = [results.labels]';
rates = round(results.n_Spk./total_times',1); % Round to nearest decimal
rates_table = table(chans,rates);
% Store spike features in a table
spks_positions = [results.markers.position]';
spks_type = {results.markers.label}';
spks_channels = {results.markers.channels}';
spks_values = {results.markers.value}';
spks_table = table(spks_channels,spks_type,spks_positions,spks_values);
% Store all tables in a .mat
[~,bname,~] = fileparts(bp_baseline_obj.fname);
spkFileName = ['Results_' bname '.mat'];
spkFile.spikeRate = rates_table;
spkFile.markers = spks_table;
spkFile.baseline = ['Total duration of ' num2str((nsamples(bp_baseline_obj)/fsample(bp_baseline_obj))/60) ' minutes'];
spkFile.comment = 'Spike rate is the number of spikes per minute for each bipolar channel';
results_file = fullfile(path_out,spkFileName);
save(results_file, 'spkFile');
spike_count = results.n_Spk;
recording_time = total_times';
rates_table = table(chans,spike_count,recording_time);
% Clean intermediate files
delete([Sbp.FileOut '.mat']);delete([Sbp.FileOut '.dat']);delete([Sbp.FileOut '_log.txt']);
delete(file_header);delete(file_data);delete(file_log);
end
... ...