diff --git a/dependencies/csc_events_to_hypnogram.m b/dependencies/csc_events_to_hypnogram.m index 32f60ee088bc2fffaa4e8e9e4818131812b7f39f..3622bf877cfac64caf8aeb7685e6fa500e885ffd 100644 --- a/dependencies/csc_events_to_hypnogram.m +++ b/dependencies/csc_events_to_hypnogram.m @@ -1,32 +1,84 @@ -function [EEG, table_data, handles] = csc_events_to_hypnogram(EEG, flag_plot) +function [EEG, table_data, handles] = csc_events_to_hypnogram(EEG, flag_plot, flag_mode) % turns event data from the csc_eeg_plotter into sleep stages and plots data -handles = []; - % Notes % ^^^^^ -% ignores event 4 labels +% flag_mode denotes whether its classical (0) or transitional scoring (1) +% the essential distinction is how to treat "event 4": in classical scoring this +% is used to mark arousals (one at the start and end of the arousal), in +% transitional scoring arousals are marked simply as brief wake periods and +% "event 4" is used to denote signal artefacts (and ignored for scoring) -% pre-allocate to wake -stages = single(zeros(1, EEG.pnts)); +% custom options +if nargin < 2 + flag_plot = false; + flag_mode = 0; % default to classical sleep scoring +elseif nargin < 3 + flag_mode = 0; +end -% convert event timing from seconds to samples -event_timing = floor([EEG.csc_event_data{:, 2}] * EEG.srate); -% loop over each event -for n = 1 : length(EEG.csc_event_data) +% remove all event 4 and save in separate arousal/artefact table +% '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +% find event 4s +event4_logical = [EEG.csc_event_data{:, 3}] == 4; +event4_ind = find(event4_logical); + +% create new table of only event 4s +event4_table = [EEG.csc_event_data(event4_logical, :)]; + +% pre-allocate arousal logical +EEG.swa_scoring.arousals = false(EEG.pnts, 1); + +% loop over each event 4 +if flag_mode == 0 + % there should be an even number of event 4s since each marks start and end + % of an arousal + if mod(sum(event4_logical), 2) + error('There is an odd number of event 4s, check your events') + end - % ignore "event 4" (artefacts) for basic sleep analysis - if EEG.csc_event_data{n, 3} == 4 - continue + for n = 1 : 2 : sum(event4_logical) + % find samples + start_sample = floor(event4_table{n, 2} * EEG.srate); + end_sample = floor(event4_table{n + 1, 2} * EEG.srate); + % mark the interval samples as true + EEG.swa_scoring.arousals(start_sample : end_sample) = true; end - if n ~= length(EEG.csc_event_data) +elseif flag_mode == 1 + % each artefact's end is simply marked by the next stage + for n = 1 : sum(event4_logical) + % find samples + start_sample = floor(event4_table{n, 2} * EEG.srate); + end_sample = floor(EEG.csc_event_data{event4_ind(n) + 1, 2} * EEG.srate); + % mark the interval samples as true + EEG.swa_scoring.arousals(start_sample : end_sample) = true; + end +end + + +% produce the sleep scoring stages % +% '''''''''''''''''''''''''''''''' % +% create a temporary table without event 4 +tmp_events = EEG.csc_event_data(~event4_logical, :); + +% convert event timing from seconds to samples +event_timing = floor([tmp_events{:, 2}] * EEG.srate); + +% pre-allocate to wake +stages = int8(ones(1, EEG.pnts) * -1); + +% loop over each scoring event +for n = 1 : length(tmp_events) + + % check for last event + if n ~= length(tmp_events) stages(event_timing(n) : event_timing(n + 1)) = ... - EEG.csc_event_data{n, 3}; + tmp_events{n, 3}; else stages(event_timing(n) : end) = ... - EEG.csc_event_data{n, 3}; + tmp_events{n, 3}; end end @@ -38,8 +90,7 @@ stages(stages == 6) = 0; EEG.swa_scoring.stages = stages; % get the sleep stats -table_data = swa_sleep_statistics(EEG, 0, 'deutsch'); - +table_data = swa_sleep_statistics(EEG, 0, 'deutsch', flag_mode); if flag_plot @@ -51,11 +102,14 @@ if flag_plot % convert time to hours for plot time_range = EEG.times / 1000 / 60 / 60; - REM_range = stages; REM_range(REM_range ~= 1) = NaN; + + % create a new range to colour REM in red + REM_range = double(stages); REM_range(REM_range ~= 1) = NaN; % plot hypnogram handles.fig = figure('color', 'w'); handles.ax = axes('nextplot', 'add', ... + 'xlim', [0, time_range(end)], ... 'ytick', 0:4, ... 'yticklabel', {'WACH', 'REM', 'N1', 'N2', 'N3'}, ... 'yDir', 'reverse', ... @@ -65,14 +119,20 @@ if flag_plot handles.REM_plot = plot(time_range, REM_range, ... 'lineWidth', 4, ... 'color', 'r'); + + % plot arousals/artefacts + arousal_range = double(EEG.swa_scoring.arousals); + arousal_range(arousal_range == 0) = NaN; + handles.arousals = plot(time_range, arousal_range - 1.2, ... + 'lineWidth', 3, ... + 'color', 'k'); + % label axes xlabel('time (hours)') ylabel('sleep stage') % make a pie chart pie_data = cell2mat(table_data(3:7, 4)); -% pie_labels = {'wake', 'N1', 'N2', 'N3', 'REM'}; figure('color', 'w'); -% handles.pie = pie(pie_data, pie_labels); handles.pie = pie(pie_data, {'', '', '', '', ''}); end \ No newline at end of file diff --git a/dependencies/csc_subplot.m b/dependencies/csc_subplot.m new file mode 100644 index 0000000000000000000000000000000000000000..27ca2a4436a8c2be6af0302b020c1c08f39b8367 --- /dev/null +++ b/dependencies/csc_subplot.m @@ -0,0 +1,29 @@ +function handles = csc_subplot(n_rows, n_columns, axes_space) + +if nargin < 3 + axes_space = 0.05; +end + +% open a figure +handles.fig = figure('color', 'w', 'units', 'norm'); + +% determine axes spacing (normalised units) +axes_space = 0.05; +axes_width = 1 / n_columns - axes_space * 1.5; +axes_height = 1 / n_rows - axes_space * 1.5; + +% determine axes positions +axes_x_s = linspace(axes_width + axes_space, 1 - axes_space, n_columns) - axes_width; +axes_y_s = linspace(axes_height + axes_space, 1 - axes_space, n_rows) - axes_height; + +% create all axes +for h = 1 : n_columns + for v = 1 : n_rows + + handles.axes(v, h) = axes(... + 'position', [... + axes_x_s(h), axes_y_s(n_rows - v + 1), ... + axes_width, axes_height], ... + 'nextplot', 'add'); + end +end \ No newline at end of file diff --git a/io/mff_scripts_internal/mff_convert_to_EEGLAB.m b/io/mff_scripts_internal/mff_convert_to_EEGLAB.m index 924d593fc75c5d1e4cbd85c7adb31201e8735970..f2514ed7018443fc082e7e3df3e6321bc410f70a 100644 --- a/io/mff_scripts_internal/mff_convert_to_EEGLAB.m +++ b/io/mff_scripts_internal/mff_convert_to_EEGLAB.m @@ -3,6 +3,7 @@ function EEG = mff_convert_to_EEGLAB(fileName, save_name, recording_system) FLAG_TO_STRUCT = true; +FLAG_DOWNSAMPLE = false; % if the filename is not specified open a user dialog if nargin < 1 @@ -53,22 +54,57 @@ EEG.chanlocs(257:end)=[]; % either write directly to file or to struct if FLAG_TO_STRUCT - % calculate total size and pre-allocate - EEG.data = zeros(EEG.nbchan, mffData.signal_binaries(recording_system).channels.num_samples(1), 'single'); - % open a progress bar waitHandle = waitbar(0,'Please wait...', 'Name', 'Importing Channels'); - % fill the EEG.data - for current_channel = 1 : EEG.nbchan + if ~FLAG_DOWNSAMPLE + % calculate total size and pre-allocate + EEG.data = zeros(EEG.nbchan, mffData.signal_binaries(recording_system).channels.num_samples(1), 'single'); - % update the waitbar - waitbar(current_channel/EEG.nbchan, waitHandle, sprintf('Channel %d of %d', current_channel, EEG.nbchan)) + % fill the EEG.data + for current_channel = 1 : EEG.nbchan + + % update the waitbar + waitbar(current_channel/EEG.nbchan, waitHandle, sprintf('Channel %d of %d', current_channel, EEG.nbchan)) + + % get all the data + temp_data = mff_import_signal_binary(mffData.signal_binaries(recording_system), current_channel, 'all'); + EEG.data(current_channel, :) = temp_data.samples; + + end + + else - % get all the data - temp_data = mff_import_signal_binary(mffData.signal_binaries(recording_system), current_channel, 'all'); - EEG.data(current_channel, :) = temp_data.samples; + % get the first channel + temp_data = mff_import_signal_binary(mffData.signal_binaries(recording_system), 1, 'all'); + % NOTE: decimate function takes care of filtering +% % design filter (cheby 1 _ 10th order _ nyquist) +% filter_design = designfilt('lowpassiir', 'FilterOrder', 10, ... +% 'PassbandFrequency', temp_data.sampling_rate / 4, 'PassbandRipple', 1, 'SampleRate', temp_data.sampling_rate); +% +% % filter the channel +% temp_data.samples = filtfilt(filter_design, double(temp_data.samples)); + + % downsample the channel + temp_data.samples = single(decimate(double(temp_data.samples), 2)); + + % pre-allocate EEG to remaining size + EEG.data = zeros(EEG.nbchan, length(temp_data.samples), 'single'); + EEG.data(1, :) = temp_data.samples; + + % loop for rest of channels + for current_channel = 2 : EEG.nbchan + + % update the waitbar + waitbar(current_channel/EEG.nbchan, waitHandle, sprintf('Channel %d of %d', current_channel, EEG.nbchan)) + + % get all the data + temp_data = mff_import_signal_binary(mffData.signal_binaries(recording_system), current_channel, 'all'); + temp_data.samples = single(decimate(double(temp_data.samples), 2)); + EEG.data(current_channel, :) = temp_data.samples; + + end end % delete the progress bar diff --git a/preprocessing/artifact_rejection/csc_artifact_detection_fft.m b/preprocessing/artifact_rejection/csc_artifact_detection_fft.m index bb42f66c55ed0fc1ef23d21911ea33310e850819..99ae4b44db412c1c53d92ff255cf94b0abff632e 100644 --- a/preprocessing/artifact_rejection/csc_artifact_detection_fft.m +++ b/preprocessing/artifact_rejection/csc_artifact_detection_fft.m @@ -1,22 +1,23 @@ -function bad_epochs = csc_artifact_detection_fft(fft_bands, bands_of_interest, method) +function bad_epochs = csc_artifact_detection_fft(fft_all, freq_range, options) % plot the fft power of each channel for each epoch and find samples above % the manually set threshold -% specify parameters - %TODO: make these parameters options -default_percentile = 99; - -% only select the bands of interest from the fft -fft_bands = fft_bands(:,:, bands_of_interest); - -% calculate the default percentiles over epochs (dim=2) of each band -channel_thresholds = squeeze(prctile(fft_bands, default_percentile, 2)); - % plot the channels and bands of interest -if strcmp(method, 'semi_automatic') +if strcmp(options.method, 'semi_automatic') + + [channel_thresholds, fft_bands] = csc_artifact_rejection_fft_gui(fft_all, freq_range, options); + +else + + % concatenate the ffts + fft_bands = csc_calculate_freq_bands(fft_all, freq_range, options); + + % only select the bands of interest from the fft + fft_bands = fft_bands(:,:, options.bands_of_interest); + + % calculate the default percentiles over epochs (dim=2) of each band + channel_thresholds = squeeze(prctile(fft_bands, options.default_percentile, 2)); - channel_thresholds = csc_artifact_rejection_fft_gui(fft_bands, channel_thresholds); - end % create threshold matrix @@ -24,9 +25,9 @@ threshold_matrix = repmat(channel_thresholds, [size(fft_bands, 2), 1]); threshold_matrix = reshape(threshold_matrix, size(fft_bands)); % find the above threshold samples -bad_channel_epoch_bands = fft_bands > threshold_matrix; +bad_epochs = fft_bands > threshold_matrix; -% find samples with minimum number above threshold epochs - %TODO: provide minimum number of channels -bad_epochs = any(sum(bad_channel_epoch_bands, 3), 1); +% % find samples with minimum number above threshold epochs +% % TODO: provide minimum number of channels +% bad_epochs = any(sum(bad_epochs, 3), 1); \ No newline at end of file diff --git a/preprocessing/artifact_rejection/csc_artifact_rejection.m b/preprocessing/artifact_rejection/csc_artifact_rejection.m index a5fcc8bc4511f386d259aad443e549e3c443778a..bd02c70cc167913c21aedf1c9ceb9329829f396a 100644 --- a/preprocessing/artifact_rejection/csc_artifact_rejection.m +++ b/preprocessing/artifact_rejection/csc_artifact_rejection.m @@ -25,22 +25,19 @@ switch method % calculate the fft % ~~~~~~~~~~~~~~~~~ [fft_all, freq_range] = csc_average_reference_and_FFT(EEG, options); - - % concatenate the ffts - %TODO: have the bands of interest as options - fft_bands = csc_calculate_freq_bands(fft_all, freq_range, options); - + % run the artifact detection % ~~~~~~~~~~~~~~~~~~~~~~~~~~ % calculate the bad epochs using thresholding of the band of % interest (2nd parameters) - bad_epochs = csc_artifact_detection_fft(fft_bands, [1, 5], 'semi_automatic'); - - % calculate the bad epochs in time - bad_epochs_time = (find(bad_epochs) * options.epoch_length) - options.epoch_length; + % TODO: bands of interest should be an option as well + EEG.bad_epochs = csc_artifact_detection_fft(fft_all, freq_range, options); - % rearrange the time stamps into regions - EEG.bad_regions = [bad_epochs_time', bad_epochs_time' + options.epoch_length]; +% % calculate the bad epochs in time +% bad_epochs_time = (find(bad_epochs) * options.epoch_length) - options.epoch_length; +% +% % rearrange the time stamps into regions +% EEG.bad_regions = [bad_epochs_time', bad_epochs_time' + options.epoch_length]; otherwise fprintf(1, 'Error: unrecognised option call: %s', method); @@ -58,12 +55,15 @@ saveName = [EEG.filename(1:end-4), '_fftANok.mat']; varargin = varargin{1}; options = struct(... - 'ave_ref', 0 ,... - 'bad_channels', [] ,... - 'save_file', 0 ,... - 'save_name', saveName,... - 'epoch_length', 6 ,... - 'freq_limit', 40 ); + 'ave_ref', 0, ... + 'bad_channels', [], ... + 'bands_of_interest', [1, 5], ... + 'default_percentile', 99, ... + 'method', 'semi_automatic', ... + 'save_file', 0, ... + 'save_name', saveName, ... + 'epoch_length', 6, ... + 'freq_limit', 40 ); % read the acceptable names optionNames = fieldnames(options); diff --git a/preprocessing/artifact_rejection/csc_artifact_rejection_fft_gui.m b/preprocessing/artifact_rejection/csc_artifact_rejection_fft_gui.m index 8dbfc4f908401c846946d2b95767b829e9b0bbf7..fa05f0add43d50bf484d374e85459fdf6d2d3d8e 100644 --- a/preprocessing/artifact_rejection/csc_artifact_rejection_fft_gui.m +++ b/preprocessing/artifact_rejection/csc_artifact_rejection_fft_gui.m @@ -1,7 +1,19 @@ -function channel_thresholds = csc_artifact_rejection_fft_gui(fft_bands, thresholds) +function [channel_thresholds, fft_bands] = csc_artifact_rejection_fft_gui(fft_all, freq_range, options) % set some default options -options.ylimitmax = 0; +options.ylimitmax = 1; + + +% calculate the intial bands +% concatenate the ffts +fft_bands = csc_calculate_freq_bands(fft_all, freq_range, options); + +% only select the bands of interest from the fft +fft_bands = fft_bands(:,:, options.bands_of_interest); + +% calculate the default percentiles over epochs (dim=2) of each band +thresholds = squeeze(prctile(fft_bands, options.default_percentile, 2)); + % create the figure handles.fig = figure(... @@ -14,6 +26,8 @@ handles.fig = figure(... set(handles.fig,... 'KeyPressFcn', {@cb_KeyPressed}); +set(handles.fig, ... + 'CloseRequestFcn', {@pb_accept}); % slider for channel navigation @@ -52,10 +66,10 @@ band_minimum = squeeze(min(min(fft_bands, [], 1), [], 2)); band_maximum = squeeze(max(max(fft_bands, [], 1), [], 2)); % draw axes -for a = 1:no_bands +for a = 1 : no_bands handles.axes(a) = axes(... 'parent', handles.fig ,... - 'position', [axes_pos_x(a) 0.1, axes_width, 0.75] ,... + 'position', [axes_pos_x(a) 0.35, axes_width, 0.5] ,... 'nextPlot', 'add' ,... 'color', [0.2, 0.2, 0.2] ,... 'xcolor', [0.9, 0.9, 0.9] ,... @@ -66,7 +80,7 @@ for a = 1:no_bands % set the y limits of the axes to be minimal and maximal possible if options.ylimitmax set(handles.axes(a),... - 'ylim', [band_minimum(a), band_maximum(a)]); + 'ylim', [0, band_maximum(a)]); end % set the callback @@ -76,13 +90,36 @@ end % set the xlimits of the axes set(handles.axes, 'xlim', [0, size(fft_bands, 2)]); -% push button to finish editing thresholds +% draw spectral axes +handles.spectral_ax = axes(... + 'parent', handles.fig ,... + 'position', [0.075, 0.05, 0.85, 0.2] ,... + 'nextPlot', 'add' ,... + 'color', [0.2, 0.2, 0.2] ,... + 'xcolor', [0.9, 0.9, 0.9] ,... + 'ycolor', [0.9, 0.9, 0.9] ,... + 'fontName', 'Century Gothic' ,... + 'fontSize', 8 ); + + +% push button to reset channel threshold handles.pb_accept = uicontrol(... 'Parent', handles.fig,... 'Style', 'pushbutton',... - 'String', 'accept',... + 'String', 'reset',... + 'Units', 'normalized',... + 'Position', [0.1 .875 0.1 0.05],... + 'FontName', 'Century Gothic',... + 'FontSize', 11); +set(handles.pb_accept, 'Callback', {@pb_accept}) + +% push button to finish editing thresholds +handles.pb_accept = uicontrol(...fft_bands + 'Parent', handles.fig,... + 'Style', 'pushbutton',... + 'String', 'done',... 'Units', 'normalized',... - 'Position', [0.8 .875 0.05 0.05],... + 'Position', [0.8 .875 0.1 0.05],... 'FontName', 'Century Gothic',... 'FontSize', 11); set(handles.pb_accept, 'Callback', {@pb_accept}) @@ -90,6 +127,8 @@ set(handles.pb_accept, 'Callback', {@pb_accept}) % set the data into the figure structure guidata(handles.fig, handles); +setappdata(handles.fig, 'fft_all', fft_all); +setappdata(handles.fig, 'freq_range', freq_range); setappdata(handles.fig, 'fft_bands', fft_bands); setappdata(handles.fig, 'thresholds', thresholds); @@ -101,6 +140,7 @@ uiwait(handles.fig); % once uiwait resumes get the most current thresholds channel_thresholds = getappdata(handles.fig, 'thresholds'); +fft_bands = getappdata(handles.fig, 'fft_bands'); % and delete the figure delete (handles.fig); @@ -115,6 +155,8 @@ handles = guidata(object); % get the data fft_bands = getappdata(handles.fig, 'fft_bands'); thresholds = getappdata(handles.fig, 'thresholds'); +fft_all = getappdata(handles.fig, 'fft_all'); +freq_range = getappdata(handles.fig, 'freq_range'); % always plot the first channel in the initial plot nCh = 1; @@ -122,7 +164,7 @@ nCh = 1; % loop plot the channel for each band on the appropriate axes for b = 1:size(fft_bands, 3); % plot the data - handles.plot(b) = plot(fft_bands(nCh,:,b),... + handles.plot(b) = plot(fft_bands(nCh,:,b),...channel_thresholds 'parent', handles.axes(b), ... 'color', [0.8, 0.8, 0.8]); @@ -147,6 +189,29 @@ set(handles.labels,... 'interpreter', 'latex'); +% calculate the original spectrum +original_spectra = sqrt(mean(fft_all(nCh, :, :) .^2 , 3)); + +% calculate the "thresholded" spectrum +currently_good_epochs = ... + fft_bands(nCh, :, 1) < thresholds(nCh, 1) | ... + fft_bands(nCh, :, 2) < thresholds(nCh, 2); +thresholded_spectra = sqrt(mean(fft_all(nCh, :, currently_good_epochs) .^2 , 3)); + +% plot the overall spectrum +flag_log = true; +if flag_log + handles.original_spectra = plot(freq_range, log(original_spectra), ... + 'parent', handles.spectral_ax, ... + 'color', [0.4, 0.4, 0.4], ... + 'lineWidth', 6); + + handles.thresholded_spectra = plot(freq_range, log(thresholded_spectra), ... + 'parent', handles.spectral_ax, ... + 'color', [0.8, 0.8, 0.8], ... + 'lineWidth', 1); +end + % set the handles structure guidata(handles.fig, handles); @@ -165,7 +230,7 @@ thresholds = getappdata(handles.fig, 'thresholds'); nCh = handles.jslider.getValue(); % loop for replotting -for b = 1:size(fft_bands, 3); +for b = 1:size(fft_bands, 3) set(handles.plot(b),... 'ydata', fft_bands(nCh,:,b)); @@ -177,6 +242,36 @@ end % update the channel indicator set(handles.channel_indicator, 'string', num2str(nCh)); +function fcn_update_spectra(object) +% fast update of the thresholded ydata in the spectrum + +% get the GUI figure handles +handles = guidata(object); + +% get the data +fft_bands = getappdata(handles.fig, 'fft_bands'); +thresholds = getappdata(handles.fig, 'thresholds'); +fft_all = getappdata(handles.fig, 'fft_all'); + +% get the current channel +nCh = handles.jslider.getValue(); + +% calculate the original spectrum +original_spectra = sqrt(mean(fft_all(nCh, :, :) .^2 , 3)); + +% calculate the "thresholded" spectrum +currently_good_epochs = not(... + fft_bands(nCh, :, 1) > thresholds(nCh, 1) | ... + fft_bands(nCh, :, 2) > thresholds(nCh, 2)); +thresholded_spectra = sqrt(mean(fft_all(nCh, :, currently_good_epochs) .^2 , 3)); + +% replace the spectrum +flag_log = true; +if flag_log + set(handles.original_spectra, 'ydata', log(original_spectra)); + set(handles.thresholded_spectra, 'ydata', log(thresholded_spectra)); +end + function cb_KeyPressed(object, eventdata) % callback function when the keyboard is pressed @@ -195,7 +290,9 @@ switch eventdata.Key handles.jslider.setValue(nCh-1); end + % update the plots fcn_update_plots(object); + fcn_update_spectra(object) case 'rightarrow' @@ -204,14 +301,16 @@ switch eventdata.Key handles.jslider.setValue(nCh+1); end + % update the plots fcn_update_plots(object); + fcn_update_spectra(object) case 'uparrow' current_limits = get(handles.axes, 'ylim'); % set each axes individually for n = 1:length(current_limits) - set(handles.axes(n), 'ylim', current_limits{n}/1.3); + set(handles.axes(n), 'ylim', [0, current_limits{n}(2) / 1.3]); end case 'downarrow' @@ -219,7 +318,7 @@ switch eventdata.Key % set each axes individually for n = 1:length(current_limits) - set(handles.axes(n), 'ylim', current_limits{n}*1.3); + set(handles.axes(n), 'ylim', [0, current_limits{n}(2) * 1.3]); end end @@ -251,6 +350,7 @@ setappdata(handles.fig, 'thresholds', thresholds); % update the axes fcn_update_plots(object); +fcn_update_spectra(object); function pb_accept(~, ~) diff --git a/preprocessing/csc_quick_spectral_analysis.m b/preprocessing/csc_quick_spectral_analysis.m index 630d117afd067f013b05e3cb45fde2a252a2d266..d7715309af7bf6443ae279ff4b8e4bdc15d9b8c9 100644 --- a/preprocessing/csc_quick_spectral_analysis.m +++ b/preprocessing/csc_quick_spectral_analysis.m @@ -1,4 +1,10 @@ -function [spectral_data, spectral_range] = csc_quick_spectral_analysis(EEG) +function [spectral_data, spectral_range] = csc_quick_spectral_analysis(EEG, flag_plot) +% performs 1 second window spectral analysis using pwelch on EEG structure + +% check inputs +if nargin < 2 + flag_plot = true; +end % flag parameters % TODO: put as input arguments @@ -26,10 +32,17 @@ end % remove bad channels if flag_remove_channels if isfield (EEG, 'bad_channels') + % bad channels should be a cell array of lists fprintf(1, 'Removing bad channels...\n'); EEG.data(EEG.bad_channels{1}, :) = []; EEG.chanlocs(EEG.bad_channels{1}) = []; EEG.nbchan = size(EEG.data, 1); + elseif isfield (EEG, 'good_channels') + % good channels should be a logical vector + fprintf(1, 'Removing bad channels...\n'); + EEG.data(~EEG.good_channels, :) = []; + EEG.chanlocs(~EEG.good_channels) = []; + EEG.nbchan = size(EEG.data, 1); end end @@ -66,18 +79,20 @@ else end end -% define ranges of interest -delta_range = spectral_range >= 0.9 & spectral_range <= 4.1; -spindle_range = spectral_range >= 10.9 & spectral_range <= 16.1; - % topography % '''''''''' -% calculate average power -delta_power = double(sqrt(mean(spectral_data(delta_range, :) .^ 2))); -spindle_power = double(sqrt(mean(spectral_data(spindle_range, :) .^ 2))); -delta_median = median(delta_power); -spindle_median = median(spindle_power); - -% plot topographies -csc_Topoplot(log(delta_power), EEG.chanlocs); -csc_Topoplot(log(spindle_power), EEG.chanlocs); +if flag_plot + % define ranges of interest + delta_range = spectral_range >= 0.9 & spectral_range <= 4.1; + spindle_range = spectral_range >= 10.9 & spectral_range <= 16.1; + + % calculate average power + delta_power = double(sqrt(mean(spectral_data(delta_range, :) .^ 2))); + spindle_power = double(sqrt(mean(spectral_data(spindle_range, :) .^ 2))); + delta_median = median(delta_power); + spindle_median = median(spindle_power); + + % plot topographies + csc_Topoplot(log(delta_power), EEG.chanlocs); + csc_Topoplot(log(spindle_power), EEG.chanlocs); +end \ No newline at end of file diff --git a/tools/csc_phase_locking_factor.m b/tools/csc_phase_locking_factor.m new file mode 100644 index 0000000000000000000000000000000000000000..2b46bd8d8a644ae2ad9c47424773452cbcf3db2e --- /dev/null +++ b/tools/csc_phase_locking_factor.m @@ -0,0 +1,41 @@ +function [PLF, threshold] = csc_phase_locking_factor(time_series, srate, flag_filter) + +% filter data if necessary +if flag_filter + fprintf(1, 'Filtering data at > 8Hz... \n'); + fc = 8; % cut off frequency + fn = floor(srate / 2); % nyquivst frequency = sample frequency/2; + order = 3; % 3rd order filter, high pass + [b a] = butter(order, (fc/fn), 'high'); + + % reshape time series to long + num_trials = size(time_series, 2); + time_series = time_series(:); + time_series = filtfilt(b, a, double(time_series)); + + % reshape back to trials + time_series = reshape(time_series, [], num_trials); +end + +for k= 1 : size(time_series,2) + hilberT(:,k) = hilbert((time_series(:,k))); %for every trial calculate hilbert transform + hilberT(:,k) = hilberT(:,k)./abs(hilberT(:,k)); %divide for the absolute value to obtain the instantaneous phase +end +PLF=(abs(mean(hilberT,2))); + +% statistical significance +tendbaseline=726; % end of baseline, in sample +alpha=0.01; %alpha value +xaxe=0:0.001:1; +meanbase=mean(PLF(1:tendbaseline)); +xrayle=raylpdf(xaxe,meanbase./sqrt(pi/2)); %build the rayle dirtribution +acch=0; +is=1; +while acch==0 + if sum(xrayle(1:is))/sum(xrayle)>(1-alpha) + acch=1; + else + is=is+1; + end +end +threshold = xaxe(1, is); \ No newline at end of file diff --git a/visualisation/csc_component_plot.m b/visualisation/csc_component_plot.m index b9e19f051a04ada94b514aecebdd567e15779760..769379c4cf5e9bac2eb13e6a7bb496e40eaf7c2a 100644 --- a/visualisation/csc_component_plot.m +++ b/visualisation/csc_component_plot.m @@ -4,6 +4,16 @@ function component_list = csc_component_plot(EEG) % make the figure handles = define_interface(); +% check for good channels field in EEG +if ~isfield(EEG, 'good_channels') + EEG.good_channels = true(EEG.nbchan, 1); +end + +% check for marked bad_data +if ~isfield(EEG, 'bad_data') + EEG.bad_data = false(1, EEG.pnts); +end + % check for EEG.icaact if isempty(EEG.icaact) fprintf(1, 'recalculating EEG.icaact...\n') @@ -15,6 +25,9 @@ if isempty(EEG.icaact) end end +% save the EEG to the figure +setappdata(handles.fig, 'EEG', EEG); + % allocate the component list from scratch % TODO: look for previously run good_components and put on component list @@ -25,21 +38,23 @@ else handles.component_list = true(number_components, 1); end -% check for marked bad_data -if ~isfield(EEG, 'bad_data') - EEG.bad_data = false(1, EEG.pnts); -end - % set some specific properties from the data set([handles.ax_erp_time, handles.ax_erp_image],... 'xlim', [EEG.times(1), EEG.times(end)] / 1000); +% check if trialed data +if EEG.trials > 1 + ica_data = reshape(EEG.icaact, size(EEG.icaact, 1), []); +else + ica_data = EEG.icaact; +end + % run the frequency analysis on all components at once % use pwelch fprintf(1, 'Computing frequency spectrum for all components\n'); window_length = 5 * EEG.srate; [spectral_data, spectral_range] = pwelch(... - EEG.icaact(:, ~EEG.bad_data)' ,... % data (transposed to channels are columns) + ica_data' ,... % data (transposed to channels are columns) hanning(window_length) ,... % window length with hanning windowing floor(window_length / 2) , ... % overlap window_length ,... % points in calculation (window length) @@ -51,7 +66,6 @@ handles.spectral_data = spectral_data(spectral_range < 45, :, :)'; % update the figure handles guidata(handles.fig, handles) -setappdata(handles.fig, 'EEG', EEG); % initial plot initial_plots(handles.fig); @@ -264,7 +278,12 @@ handles.plots.image = ... % plot the evoked potential % % ------------------------- % data_to_plot = mean(EEG.icaact(no_comp, : , :), 3)'; -data_to_plot(EEG.bad_data) = nan; +if isfield(EEG, 'bad_data') + data_to_plot(EEG.bad_data) = nan; +elseif isfield(EEG, 'good_data') + data_to_plot(~EEG.good_data) = nan; +end + handles.plots.erp_time = ... plot(handles.ax_erp_time,... EEG.times / 1000, data_to_plot,... @@ -290,7 +309,7 @@ handles.plots.spectra = ... % plot the topography % % ------------------- % handles.plots.topo = ... - csc_Topoplot(EEG.icawinv(:, no_comp), EEG.chanlocs ,... + csc_Topoplot(EEG.icawinv(EEG.good_channels, no_comp), EEG.chanlocs(EEG.good_channels) ,... 'axes', handles.ax_topoplot ,... 'plotChannels', false); @@ -341,8 +360,11 @@ else end % re-set the image of all trials % +trial_images = squeeze(EEG.icaact(current_component, :, :))'; set(handles.plots.image, ... - 'cData', squeeze(EEG.icaact(current_component, :, :))'); + 'cData', trial_images); +% re-adjust the axes limits to match image percentiles +set(handles.ax_erp_image, 'CLim', [prctile(trial_images(:), 2), prctile(trial_images(:), 98)]); % re-set the evoked potential % data_to_plot = mean(EEG.icaact(current_component, : , :), 3)'; @@ -361,11 +383,15 @@ set(handles.plots.spectra, ... % ------------------- % % plot the topography % % ------------------- % +current_topo_data = EEG.icawinv(EEG.good_channels, current_component); handles.plots.topo = ... - csc_Topoplot(EEG.icawinv(:, current_component), EEG.chanlocs ,... + csc_Topoplot(current_topo_data, EEG.chanlocs(EEG.good_channels),... 'axes', handles.ax_topoplot ,... 'plotChannels', false); +% adjust the color limits +set(handles.plots.topo.CurrentAxes, 'CLim', [min(current_topo_data), max(current_topo_data)]); + guidata(handles.fig, handles) diff --git a/visualisation/csc_eeg_plotter.m b/visualisation/csc_eeg_plotter.m index 7a77ed410533680bf3d250f1db7e65fc9fbc1301..8ee6aa1599978e92b0772a2599ed229ed2eb09b2 100644 --- a/visualisation/csc_eeg_plotter.m +++ b/visualisation/csc_eeg_plotter.m @@ -12,7 +12,7 @@ function EEG = csc_eeg_plotter(varargin) % TODO: Fix this ugly default setting style (e.g. handles.options...) % declare defaults EPOCH_LENGTH = 30; -FILTER_OPTIONS = [0.3 40]; % default filter bandpass +FILTER_OPTIONS = [0.7 40; 10 40; 0.3 10; 0.1 10]; % default filter bandpass handles.n_disp_chans = 12; % number of initial channels to display handles.v_grid_spacing = 1; % vertical grid default spacing (s) handles.h_grid_spacing = 75; % horizontal grid default spacing (uV) @@ -127,7 +127,7 @@ handles.menu.filter_toggle = uimenu(handles.menu.options,... handles.menu.filter_settings = uimenu(handles.menu.options,... 'label', 'filter settings',... 'accelerator', 'f' ,... - 'callback', {@fcn_options, 'filter_settings'}); + 'callback', {@fcn_filter_settings}); handles.menu.icatoggle = uimenu(handles.menu.options,... 'label', 'toggle channels/components',... 'checked', 'off' ,... @@ -399,7 +399,7 @@ else % normal plotting of activity % recalculate data based on projections projection_data = EEG.icawinv(:, EEG.good_components) ... - * EEG.icaact(EEG.good_components, range); + * icaData(EEG.good_components, range); % subselect displayed channels data_to_plot = projection_data(... @@ -436,25 +436,38 @@ end if strcmp(get(handles.menu.filter_toggle, 'checked'), 'on') ... && handles.plotICA == 0 - % determine filtering parameters - % check for empty boxes for one-sided filters - if isnan(handles.filter_options(2)) + % calculate and filter separately for each channel type + for n = 1 : 4 - [filt_param_b, filt_param_a] = butter(2, ... - handles.filter_options(1)/(EEG.srate/2), 'high'); + % which channels to apply this to + channel_ind = handles.channel_types(handles.disp_chans) == n; - elseif isnan(handles.filter_options(1)) + % check if relevant + if sum(channel_ind) == 0 + continue + end - [filt_param_b, filt_param_a] = butter(2, ... - handles.filter_options(2)/(EEG.srate/2), 'low'); - else - [filt_param_b, filt_param_a] = ... - butter(2,[handles.filter_options(1)/(EEG.srate/2),... - handles.filter_options(2)/(EEG.srate/2)]); + % determine filtering parameters + % check for empty boxes for one-sided filters + if isnan(handles.filter_options(n, 2)) + + [filt_param_b, filt_param_a] = butter(2, ... + handles.filter_options(n, 1) / (EEG.srate / 2), 'high'); + + elseif isnan(handles.filter_options(1)) + + [filt_param_b, filt_param_a] = butter(2, ... + handles.filter_options(n, 2) / (EEG.srate / 2), 'low'); + else + [filt_param_b, filt_param_a] = ... + butter(2,[handles.filter_options(n, 1)/(EEG.srate/2),... + handles.filter_options(n, 2) / (EEG.srate / 2)]); + end + + % apply the filter to the data window + data_to_plot(channel_ind, :) = single(filtfilt(filt_param_b, filt_param_a, ... + double(data_to_plot(channel_ind, :)'))'); end - - % apply the filter to the data window - data_to_plot = single(filtfilt(filt_param_b, filt_param_a, double(data_to_plot'))'); end @@ -532,8 +545,12 @@ end % plot the channel data % ^^^^^^^^^^^^^^^^^^^^^ if flag_replot - % delete existing handles - if isfield(handles, 'plot_eeg'); delete(handles.plot_eeg); end + + % delete existing handles and lines + if isfield(handles, 'plot_eeg') + delete(handles.plot_eeg); + handles = rmfield(handles, 'plot_eeg'); + end % plot lines handles.plot_eeg = line(time, data_to_plot,... @@ -555,7 +572,10 @@ end % plot the labels in their own boxes % ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ % delete existing handles -if isfield(handles, 'labels'); delete(handles.labels); end +if isfield(handles, 'labels') + delete(handles.labels); + handles = rmfield(handles, 'labels'); +end handles.labels = zeros(handles.n_disp_chans, 1); % loop for each label for n = 1 : handles.n_disp_chans @@ -732,6 +752,7 @@ if ~isfield(EEG, 'csc_montage') EEG.csc_montage.label_channels = cell(EEG.nbchan, 1); for n = 1 : EEG.nbchan EEG.csc_montage.label_channels(n) = {num2str(n)}; + EEG.csc_montage.channel_type(n) = {'EEG'}; end EEG.csc_montage.channels(:, 1) = 1:EEG.nbchan; EEG.csc_montage.channels(:, 2) = EEG.nbchan; @@ -750,6 +771,21 @@ if ~isfield(EEG.csc_montage, 'scaling') EEG.csc_montage.scaling(:, 1) = ones(EEG.nbchan, 1); end +% check for channel type for backwards compatibility with old montages +if ~isfield(EEG.csc_montage, 'channel_type') + for n = 1 : length(EEG.csc_montage.label_channels) + EEG.csc_montage.channel_type{n, 1} = 'EEG'; + end +end + +% recalculate channel types +handles.channel_types = ones(length(EEG.csc_montage.channel_type), 1); +channel_types = {'EEG', 'EMG', 'EOG', 'Other'}; +for n = 1 : 4 + type_ind = cellfun(@(x) strcmp(x, channel_types{n}), EEG.csc_montage.channel_type); + handles.channel_types(type_ind) = n; +end + % check that the montage has enough channels to display if length(EEG.csc_montage.label_channels) < handles.n_disp_chans handles.n_disp_chans = length(EEG.csc_montage.label_channels); @@ -760,39 +796,29 @@ end % load ICA time courses if the information need to construct them is available. if isfield(EEG, 'icaweights') && isfield(EEG, 'icasphere') if ~isempty(EEG.icaweights) && ~isempty(EEG.icasphere) - % If we have the same number of components as channels... - if size(EEG.icaweights, 1) == size(eegData, 1) - % check for epoched data - ica_data = (EEG.icaweights * EEG.icasphere) ... - * EEG.data(EEG.icachansind, :); + % re-calculate the ICA activations + ica_data = (EEG.icaweights * EEG.icasphere) ... + * EEG.data(EEG.icachansind, :); + + % pad the remaining data in case not equal size + if ~[size(EEG.icaweights, 1) == size(eegData, 1)] + dimdiff = size(eegData, 1) - size(EEG.icaweights, 1); + pad = zeros(dimdiff, size(eegData, 2)); - setappdata(handles.fig, 'icaData', ica_data); - % If we have fewer components than channels (maybe you've already removed - % some of them), then pad the ICA weights with zeros and produce component - % activations as if you had the same number of components as channels. + % add to ica activations + ica_data = [ica_data; pad]; - elseif size(EEG.icaweights, 1) < size(eegData,1) - dimdiff = size(eegData, 1) - size(EEG.icaweights, 1); - pad = zeros(dimdiff, size(EEG.icaweights, 2)); - paddedweights = [EEG.icaweights ; pad]; - try - ica_data = paddedweights * EEG.icasphere * eegData; - catch - fprintf('%s. %s.\n',... - 'Data were reinterpolated after IC removal',... - 'Can no longer display IC activations'); - ica_data = zeros(size(eegData)); - end - setappdata(handles.fig, 'icaData', ica_data); - else - error('ICA unmixing matrix is too large for data'); end + + % save the ica data to the handles structure + setappdata(handles.fig, 'icaData', ica_data); + + % check for a "good_components" field + if ~isfield(EEG, 'good_components') + EEG.good_components = true(size(EEG.icaweights, 1), 1); + end + end - % check for a "good_components" field - if ~isfield(EEG, 'good_components') && ~isempty(EEG.icaweights) - EEG.good_components = true(size(ica_data, 1), 1); - end - end % adjust initially scaling to match the data @@ -809,20 +835,20 @@ if isfield(EEG, 'hidden_channels') handles.hidden_channels = EEG.hidden_channels; end +% look for bad trials +if isfield(EEG, 'marked_trials') + handles.trials = EEG.marked_trials; +else + % allocate marked trials + handles.trials = false(EEG.trials, 1); +end + % turn on the montage option set(handles.menu.montage, 'enable', 'on'); -% allocate marked trials -handles.trials = false(EEG.trials, 1); -guidata(handles.fig, handles) - % update the handles guidata(object, handles); -% reset the scrollbar values -handles.vertical_scroll.Max = -1; -handles.vertical_scroll.Min = -(EEG.nbchan - length(handles.disp_chans) + 1); - function fcn_close_window(object, ~) % just resume the ui if the figure is closed handles = guidata(object); @@ -869,7 +895,7 @@ handles.fig = figure(... handles.table = uitable(... 'parent', handles.fig ,... 'units', 'normalized' ,... - 'position', [0.05, 0.25, 0.9, 0.7] ,... + 'position', [0.05, 0.30, 0.9, 0.65] ,... 'backgroundcolor', handles.csc_plotter.colorscheme.bg_col_2 ,... 'foregroundcolor', handles.csc_plotter.colorscheme.fg_col_1 ,... 'columnName', {'label','time', 'type'}); @@ -880,7 +906,7 @@ handles.clear_button = uicontrol(... 'Style', 'pushbutton' ,... 'String', 'clear event(s)' ,... 'Units', 'normalized' ,... - 'Position', [0.05 0.15 0.9 0.04],... + 'Position', [0.05 0.20 0.9 0.04],... 'FontName', 'Century Gothic' ,... 'FontSize', 11,... 'tooltipString', 'delete above selected events'); @@ -892,11 +918,11 @@ handles.import_button = uicontrol(... 'Style', 'pushbutton' ,... 'String', 'import events' ,... 'Units', 'normalized' ,... - 'Position', [0.05 0.10 0.9 0.04] ,... + 'Position', [0.05 0.15 0.9 0.04] ,... 'FontName', 'Century Gothic',... 'FontSize', 11 ,... 'tooltipString', 'import events' ,... - 'enable', 'off'); + 'enable', 'on'); set(handles.import_button, 'callback', {@pb_event_option, 'import'}); % export to workspace @@ -905,12 +931,24 @@ handles.export_button = uicontrol(... 'Style', 'pushbutton',... 'String', 'export to workspace',... 'Units', 'normalized',... - 'Position', [0.05 0.05 0.9 0.04],... + 'Position', [0.05 0.10 0.9 0.04],... 'FontName', 'Century Gothic',... 'FontSize', 11,... 'tooltipString', 'export to workspace'); set(handles.export_button, 'callback', {@pb_event_option, 'export'}); +% export to file +handles.save_button = uicontrol(... + 'Parent', handles.fig,... + 'Style', 'pushbutton',... + 'String', 'save to file',... + 'Units', 'normalized',... + 'Position', [0.05 0.05 0.9 0.04],... + 'FontName', 'Century Gothic',... + 'FontSize', 11,... + 'tooltipString', 'export to file'); +set(handles.save_button, 'callback', {@pb_event_option, 'save'}); + % get the underlying java properties jscroll = findjobj(handles.table); jscroll.setVerticalScrollBarPolicy(jscroll.java.VERTICAL_SCROLLBAR_ALWAYS); @@ -1154,6 +1192,9 @@ handles.trial_borders = plot(x, y(1),... 'parent', handles.main_ax,... 'buttonDownFcn', {@bdf_mark_trial}); +% set previously marked trials color +set(handles.trial_borders(handles.trials), 'markerFaceColor', [0.9, 0.2, 0.2]); + % update the GUI handles guidata(handles.fig, handles) @@ -1200,12 +1241,39 @@ switch option fcn_redraw_events(handles.csc_plotter.fig, []); case 'import' - % TODO: import from workspace variable and EEG.event + % load an event file + [file_name, file_path] = uigetfile('*.mat'); + loaded_event_file = load(fullfile(file_path, file_name)); + + if isfield(loaded_event_file, 'event_data') + set(handles.table, 'data', loaded_event_file.event_data); + + % put the new table into the structure + % Get the EEG from the figure's appdata + EEG = getappdata(handles.csc_plotter.fig, 'EEG'); + EEG.csc_event_data = loaded_event_file.event_data; + setappdata(handles.csc_plotter.fig, 'EEG', EEG); + + % re-draw the event window in main + fcn_redraw_events(handles.csc_plotter.fig, []); + + else + % TODO: event error message + end case 'export' % assign events to base workspace event_data = fcn_compute_events(handles.csc_plotter); assignin('base', 'event_data', event_data); + + case 'save' + % save the events as a file + event_data = fcn_compute_events(handles.csc_plotter); + + % Ask where to put file... + [saveFile, savePath] = uiputfile('*.mat'); + save(fullfile(savePath, saveFile), 'event_data', '-mat'); + end @@ -1738,6 +1806,79 @@ elseif any(strcmp(event.Modifier, {'control', 'alt'})) end end +function fcn_filter_settings(object, event) +% get the original figure handles +handles.csc_plotter = guidata(object); +EEG = getappdata(handles.csc_plotter.fig, 'EEG'); + +% make a window +handles.fig = figure(... + 'name', 'csc filter settings',... + 'numberTitle', 'off',... + 'color', handles.csc_plotter.colorscheme.bg_col_1,... + 'menuBar', 'none',... + 'units', 'normalized',... + 'outerPosition', [0.4 0.4 0.2 0.2]); + +% montage table +handles.table = uitable(... + 'parent', handles.fig ,... + 'units', 'normalized' ,... + 'position', [0.05, 0.30, 0.9, 0.65] ,... + 'backgroundcolor', handles.csc_plotter.colorscheme.bg_col_2 ,... + 'foregroundcolor', handles.csc_plotter.colorscheme.fg_col_1 ,... + 'columnEditable', [false, true, true, false], ... + 'columnFormat', {'char', 'numeric', 'numeric', 'char'}, ... + 'columnName', {'type', 'high', 'low', 'color'}); + +% set the data +data = cell(4, 4); +data(:, 1) = {'EEG', 'EMG', 'EOG', 'Other'}; +data(:, [2, 3]) = num2cell(handles.csc_plotter.filter_options); +% TODO: make color options for individual channel types +data(:, 4) = {'d'}; +set(handles.table, 'Data', data); + +% automatically adjust the column width using java handle +jscroll = findjobj(handles.table); +jtable = jscroll.getViewport.getView; +jtable.setAutoResizeMode(jtable.AUTO_RESIZE_ALL_COLUMNS); + +% create the save button +handles.apply_filters = uicontrol(... + 'parent', handles.fig,... + 'style', 'push',... + 'string', 'apply',... + 'foregroundColor', 'k',... + 'units', 'normalized',... + 'position', [0.275 0.1 0.5 0.1],... + 'fontName', 'Century Gothic',... + 'fontWeight', 'bold',... + 'fontSize', 10); +set(handles.apply_filters, 'callback', {@fcn_apply_filters}); + +% save those handles +guidata(handles.fig, handles); + +function fcn_apply_filters(object, event) +% get filter figure handles +handles = guidata(object); +handles.csc_plotter = guidata(handles.csc_plotter.fig); + +% get relevant filter data +data = get(handles.table, 'Data'); +filter_data = cell2mat(data(:, [2, 3])); + +% put in plotter figure handles +handles.csc_plotter.filter_options = filter_data; + +% save those handles +guidata(handles.fig, handles); +guidata(handles.csc_plotter.fig, handles.csc_plotter); + +% update the plot +update_main_plot(handles.csc_plotter.fig, 1); + % Montage Functions % ^^^^^^^^^^^^^^^^^ @@ -1822,8 +1963,9 @@ handles.table = uitable(... 'position', [0.675, 0.1, 0.3, 0.8] ,... 'backgroundcolor', [0.1, 0.1, 0.1; 0.2, 0.2, 0.2],... 'foregroundcolor', [0.9, 0.9, 0.9] ,... - 'columnName', {'name','chan','ref', 'scale'},... - 'columnEditable', [true, true, true, true]); + 'columnName', {'name','chan','ref', 'scale', 'type'},... + 'columnFormat', {'char', 'numeric', 'numeric', 'numeric', {'EEG', 'EMG', 'EOG', 'Other'}}, ... + 'columnEditable', [true, true, true, true, true]); % automatically adjust the column width using java handle jscroll = findjobj(handles.table); @@ -1859,6 +2001,14 @@ handles.reference_list = uicontrol( ... 'selectionHighlight', 'off' ,... 'fontName', 'Century Gothic' ,... 'fontSize', 8); +% if montage exists set to reference already chosen +switch EEG.csc_montage.reference + case 'custom' + set(handles.reference_list, 'value', 2); + case 'average' + set(handles.reference_list, 'value', 3); +end +% set callback set(handles.reference_list, 'callback', {@fcn_montage_buttons, 'reference'}); % create the buttons @@ -1885,7 +2035,8 @@ handles.button_delete = uicontrol(... 'FontWeight', 'bold',... 'FontSize', 10); set(handles.button_delete, 'callback', {@fcn_button_delete}); - +guidata(handles.fig, handles); +guidata(handles.csc_plotter.fig, handles.csc_plotter); handles.button_reset = uicontrol(... 'Parent', handles.fig,... 'Style', 'push',... @@ -1916,6 +2067,7 @@ data = cell(length(EEG.csc_montage.label_channels), 3); data(:, 1) = deal(EEG.csc_montage.label_channels); data(:, [2,3]) = num2cell(EEG.csc_montage.channels); data(:, 4) = num2cell(EEG.csc_montage.scaling); +data(:, 5) = deal(EEG.csc_montage.channel_type); % put the data into the table set(handles.table, 'data', data); @@ -2082,6 +2234,7 @@ end EEG.csc_montage.label_channels = data(:,1); EEG.csc_montage.channels = cell2mat(data(:,[2,3])); EEG.csc_montage.scaling = cell2mat(data(:, 4)); +EEG.csc_montage.channel_type = data(:, 5); % compatibility with older matlab versions (handles dot notation). if verLessThan('matlab', '8.4') @@ -2102,10 +2255,6 @@ end % Reset hidden channels handles.csc_plotter.hidden_chans = []; -% update the slider to reflect new montage -handles.csc_plotter.vertical_scroll.Value = -1; -handles.csc_plotter.vertical_scroll.Min = -(EEG.nbchan - length(handles.csc_plotter.disp_chans)); - % change montage name % compatibility with older matlab versions (handles dot notation). tmp_list = get(handles.montage_list, 'string'); @@ -2115,6 +2264,14 @@ else EEG.csc_montage.name = tmp_list{get(handles.montage_list, 'Value')}; end +% recalculate channel type (faster plotting if calculated once here) +handles.csc_plotter.channel_types = ones(length(EEG.csc_montage.channel_type), 1); +channel_types = {'EEG', 'EMG', 'EOG', 'Other'}; +for n = 1 : 4 + type_ind = cellfun(@(x) strcmp(x, channel_types{n}), EEG.csc_montage.channel_type); + handles.csc_plotter.channel_types(type_ind) = n; +end + % update the handle structures guidata(handles.fig, handles); guidata(handles.csc_plotter.fig, handles.csc_plotter); @@ -2123,7 +2280,6 @@ setappdata(handles.csc_plotter.fig, 'EEG', EEG); % update the plot update_main_plot(handles.csc_plotter.fig); - function fcn_montage_buttons(object, ~, event_type) % get the montage handles handles = guidata(object); @@ -2132,7 +2288,7 @@ switch event_type case 'add' % add a row of ones to the table old_data = handles.table.Data; - new_row = {'', 1, 0, 1}; + new_row = {'', 1, 0, 1, 'EEG'}; set(handles.table, 'Data', [old_data; new_row]); case 'reset' @@ -2148,6 +2304,7 @@ switch event_type [montage(:, 2)] = deal(num2cell(1 : EEG.nbchan)); % actual number montage(:, 3) = {EEG.nbchan}; % reference montage(:, 4) = {1}; % scaling + [montage(:, 5)] = deal({'EEG'}); % assign to table set(handles.table, 'Data', montage); @@ -2176,7 +2333,15 @@ EEG.csc_montage.name = montage_name; if ~isempty(montage_name) && ~strcmp(montage_name, 'original') montage = load(fullfile(montage_dir, montage_name), '-mat'); if isfield(montage, 'data') + + % check for missing channel types + if size(montage.data, 2) < 5 + montage.data(:, 5) = deal({'EEG'}); + end + + % set into the table set(handles.table, 'data', montage.data); + else fprintf(1, 'Warning: could not find montage data in the file.\n'); end