% pop_micro_selectdata() - selects and aggregates data for microstate segmentation.
%
% pop_micro_selectdata selects data for microstate segmentation, in
% severeral different ways:
% * For spontaneous data: selection of GFP peaks.
% * For event related data: ERPs as average over epochs. If aggregating
% multiple datasets they can either be concatenated or averaged into
% single ERP.
% It is possible to aggregate data from multiple datasets loaded in ALLEEG.
% When aggregating data a new dataset 'NewEEG' will be created and added to
% ALLEEG.
%
% Usage:
% >> EEG = pop_micro_selectdata ( EEG ); % pop up window
% >> [EEG, ALLEEG, CURRENTSET, NewEEG] = pop_micro_selectdata ( EEG, ...
% ALLEEG ); % pop up window
% >> [EEG, ALLEEG, CURRENTSET, NewEEG] = pop_micro_selectdata ( EEG, ...
% ALLEEG, 'key1', 'val1', 'key2', 'val2' ... )
%
% Please cite this toolbox as:
% Poulsen, A. T., Pedroni, A., Langer, N., & Hansen, L. K. (2018).
% Microstate EEGlab toolbox: An introductionary guide. bioRxiv.
%
% Inputs:
% EEG - EEGlab EEG structure.
%
% Optional inputs:
% ALLEEG - EEGlab structure containing all datasets read into
% EEGlab as EEG structures.
% 'datatype' - 'spontaneous': Finds GFP peaks. 'ERPavg': Averages
% over epochs and aggregated datasets. 'ERPconc':
% Averages over epochs and concatenates aggregated
% datasets. Default is 'ERPconc'.
% multiple datasets they can either be concatenated or averaged into
% single ERP. . ----- -**********
% 'dataset_idx' - Indices for datasets in ALLEEG to aggragate data
% from. Leave empty to use current dataset (default).
% 'avgref' - Calculate average reference. 1 - yes (default),
% 0 - no.
% 'normalise' - Normalise each dataset with average channel std.
% 1 - yes, 0 - no (default).
% 'Markers' - String that indicates the marker to which data is
% epoched. !Not implemented yet!
% 'MinPeakDist' - Minimum Distance between GFP peaks in ms
% (default = 10 ms).
% 'Npeaks' - Number of GFP peaks per subject that enter the
% segmentation. Note that the maximum number of peaks
% is restricted to the minimum number of GFP peaks
% across subjects.
% 'GFPthresh' - Reject peaks over threshold (multiples of std(GFP)).
% Set to avoid extreme GFP peaks that are likely
% caused by artifacts. Set to zero to turn off
% (default).
% Outputs:
% EEG - EEGlab EEG structure with added microstate field.
% ALLEEG - EEGlab structure containing all datasets read into
% EEGlab as EEG structures.
% CURRENTSET - Workspace variable index of the current dataset. Only
% relevant when creating a new dataset through aggragating
% data.
% NewEEG - New EEGlab EEG structure containing data selected for
% microstate analysis and info in added microstate field.
%
% Authors:
% Andreas Trier Poulsen, atpo@dtu.dk
% Technical University of Denmark, DTU Compute, Cognitive systems.
%
% Andreas Pedroni, andreas.pedroni@uzh.ch
% University of Zürich, Psychologisches Institut, Methoden der
% Plastizitätsforschung.
%
% August 2017.
%
% See also: eeglab
% Copyright (C) 2017 Andreas Trier Poulsen, atpo@dtu.dk
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
function [EEG, ALLEEG, com, NewEEG] = pop_micro_selectdata(EEG, ALLEEG, varargin)
%% Error check and initialisation
com = '';
if nargin < 1
help pop_micro_selectdata;
return;
elseif nargin == 1
ALLEEG = EEG; % in case ALLEEG was not given as input.
end
if nargin < 3
% pop-up window in case no further input is given
settings = input_popup(ALLEEG);
if strcmp(settings,'cancel')
return
end
else
settings = check_settings(varargin);
end
%% Readying settings
dataset_idx = settings.dataset_idx;
if isempty(dataset_idx)
aggregate_data = 0;
Ndatasets = 1;
NewEEG = [];
else
if length(dataset_idx) == 1
error('Please select more than one dataset, when aggregating data.')
end
aggregate_data = 1;
Ndatasets = length(dataset_idx);
% check consistency of selected datasets
NewEEG = check_datasets_consistency(ALLEEG(dataset_idx));
end
%% do the aggregation for spontaneous or erp data
if strcmp(settings.datatype,'spontaneous')
%% Readying settings for countinuous data
Npeaks = settings.Npeaks;
MinPeakDist = settings.MinPeakDist;
GFPthresh = settings.GFPthresh;
% loop through the subjects
for i=1:Ndatasets
%% Load data
if aggregate_data
X = ALLEEG(dataset_idx(i)).data;
fs = ALLEEG(dataset_idx(i)).srate;
else
X = EEG.data;
fs = EEG.srate;
end
X = reshape(X ,size(X,1),size(X,2)*size(X,3));
%% Average reference
if settings.avgref
X = bsxfun(@minus, X, mean(X));
end
%% Normalise by average channel std.
if settings.normalise
X = X ./ mean(std(X,0,2));
end
%% Parse GFP peaks
% do the segmentation on the GFP peaks, because they are high in
% signal-to-noise
% calculate the GFP
GFP = double(std(X,[],1));
% Minimum distance in ms:
MinTFdistance = MinPeakDist * fs/1000; % assuming ms
[~, peakidx{i,1}] = findpeaks(GFP,'MinPeakDistance',MinTFdistance);
% OPTION, which may be included: to avoid extreme GFP peaks that are likely caused by artifacts
if GFPthresh > 0
noisepeak_idx = GFP(peakidx{i}) > (mean(GFP) + GFPthresh*std(GFP));
peakidx{i} = peakidx{i}(~noisepeak_idx);
end
end
%% Select the GPF peaks
% select a number of random peaks (if there are many peaks or you are
% impatient)
% this is the maximum number that can be selected of each subject
maxNpeaks = min(cellfun('length',peakidx));
if Npeaks > maxNpeaks
Npeaks = maxNpeaks;
end
GFPdata = [];
for i = 1:Ndatasets
% load data
if aggregate_data
X = ALLEEG(dataset_idx(i)).data;
else
X = EEG.data;
end
X = reshape(X ,size(X,1),size(X,2)*size(X,3));
% Average reference
if settings.avgref
X = bsxfun(@minus, X, mean(X));
end
% Normalise by average channel std.
if settings.normalise
X = X ./ mean(std(X,0,2));
end
% find a number of random peaks
selection = randperm(length(peakidx{i}));
GFPpeakidx{i} = peakidx{i}(selection(1:Npeaks));
GFPdata = [GFPdata, X(:,GFPpeakidx{i})];
end
%% Save chosen data in the EEG struct
if aggregate_data
% Add information to NewEEG struct
NewEEG.microstate.data = 'set_data';
NewEEG.microstate.GFPpeakidx = GFPpeakidx;
NewEEG.microstate.data_origin.ALLEEG_idx = dataset_idx;
dataset_names = {ALLEEG.setname};
NewEEG.microstate.data_origin.dataset_names = dataset_names(dataset_idx);
NewEEG.data = GFPdata;
NewEEG.setname = 'MicroGFPpeakData';
NewEEG.pnts = size(NewEEG.data,2);
NewEEG.chanlocs = ALLEEG(dataset_idx(1)).chanlocs; % assuming datasets have the same chanlocs
NewEEG.times = 0:1000/NewEEG.srate:(length(NewEEG.data)*1000/NewEEG.srate)-1;
NewEEG.event = [];
NewEEG.urevent = [];
NewEEG.eventdescription = [];
ALLEEG = eeg_store(ALLEEG, NewEEG);
else
% Save data in EEG struct given as input.
EEG.microstate.data = GFPdata;
EEG.microstate.GFPpeakidx = GFPpeakidx{1};
end
elseif sum(strcmp(settings.datatype,{'ERPavg','ERPconc'}))
ERP = [];
for i = 1:Ndatasets
if aggregate_data
X = ALLEEG(dataset_idx(i)).data;
else
X = EEG.data;
end
%% Average reference
if settings.avgref
X = bsxfun(@minus, X, mean(X));
end
%% Normalise by average channel std.
if settings.normalise
X = X ./ mean(std(X(:,:),0,2));
end
set_avg = mean(X,3);
% preparing datasets in case of aggregation of multiple datasets
if strcmp(settings.datatype,'ERPavg')
% grand average across datasets
ERP = cat(3, ERP, set_avg);
else
% concatenating dataset ERPs
ERP = cat(2, ERP, set_avg);
end
end
if strcmp(settings.datatype,'ERPavg')
% averaging across datasets in case of aggregation
ERP = mean(ERP, 3);
end
%% Save chosen data in the EEG struct
if aggregate_data
% Add information to NewEEG struct
NewEEG.microstate.data = 'set_data';
NewEEG.microstate.data_origin.ALLEEG_idx = dataset_idx;
dataset_names = {ALLEEG.setname};
NewEEG.microstate.data_origin.dataset_names = dataset_names(dataset_idx);
NewEEG.data = ERP;
NewEEG.setname = 'MicroERPdata';
NewEEG.pnts = size(NewEEG.data,2);
NewEEG.chanlocs = ALLEEG(dataset_idx(1)).chanlocs; % assuming datasets have the same chanlocs
NewEEG.times = 0:1000/NewEEG.srate:(length(NewEEG.data)*1000/NewEEG.srate)-1;
NewEEG.event = [];
NewEEG.urevent = [];
NewEEG.eventdescription = [];
ALLEEG = eeg_store(ALLEEG, NewEEG);
else
% Save data in EEG struct given as input.
EEG.microstate.data = ERP;
end
end
%% Define command string
com = sprintf('[%s, %s] = pop_micro_selectdata( %s, %s', inputname(1), ...
inputname(2),inputname(1), inputname(2));
com = settings_to_string(com,settings);
com = [com ' );'];
end
% -------------- Pop-ups-------------- %
function settings = input_popup(ALLEEG)
% Function for creating popup window to input algorithm settings
%
%% Create Inputs for popup
% Title string
info_str1 = 'Please note that this is an early version of the plugin. Bug-reports and suggestions';
info_str2 = 'are welcome at atpo@dtu.dk.';
line.info = { {'Style' 'text' 'string' info_str1} ...
{'Style' 'text' 'string' info_str2} {} };
geo.info = {1 1 1};
% Datatype
style.datatype = 'popupmenu';
dropdown_data = {'ERP - Concatenate datasets' 'ERP - Average datasets' ... % For dropdown menu
'Spontaneous - GFP peaks'};
popmenu.datatype = {'ERPconc' 'ERPavg' 'spontaneous'}; % Corresponding calls for pop-function
data_str = dropdown_data{1}; %string for popupmenu
for data_idx = 2:length(dropdown_data); data_str = [data_str '|' dropdown_data{data_idx}]; end;
line.datatype = { {'Style' 'text' 'string' 'Data type:'}, ...
{'Style' style.datatype 'string' data_str 'tag' 'datatype' 'value' 1} };
geo.datatype = {[1 1]};
% Aggregate data?
style.aggregate_data = 'checkbox';
line.aggregate_data = { {'Style' style.aggregate_data 'value' 0 'string' ...
'Select and aggregate data from other datasets (will open new window).'...
'tag' 'aggregate_data'} };
geo.aggregate_data = {1};
% Average reference?
style.avgref = 'checkbox';
line.avgref = { {'Style' style.avgref 'value' 1 'string' 'Calculate average reference.' ...
'tag' 'avgref'} };
geo.avgref = {1};
% Normalise dataset(s)?
style.normalise = 'checkbox';
norm_tipstr = 'Normalise each dataset with average channel std. ';
line.normalise = { {'Style' style.normalise 'value' 0 'string' 'Normalise dataset(s).' ...
'tooltipstring' norm_tipstr 'tag' 'normalise'} };
geo.normalise = {1};
% GFP info string
gfp_info_str1 = 'Settings for the extraction of GFP peak maps.';
gfp_info_str2 = '(Settings will be ignored if ERP is chosen).';
line.gfp_info = { {} {'Style' 'text' 'string' gfp_info_str1 'fontweight' 'bold'}...
{'Style' 'text' 'string' gfp_info_str2}};
geo.gfp_info = {1 1 1};
% MinPeakDist
style.MinPeakDist = 'edit';
line.MinPeakDist = { {'Style' 'text' 'string' 'Minimum peak distance (ms):'}, ...
{'Style' style.MinPeakDist 'string' ' 10 ' 'tag' 'MinPeakDist'} };
geo.MinPeakDist = {[1 .2]};
% Npeaks
style.Npeaks = 'edit';
npeaks_tipstr = ['Note that the maximum number of peaks is restricted to '...
'the minimum number of GFP peaks across subjects.'];
line.Npeaks = { {'Style' 'text' 'string' ...
'No. of GFP peaks per subject that enter the segmentation.:' ...
'tooltipstring' npeaks_tipstr}, ...
{'Style' style.Npeaks 'string' ' 1000 ' 'tag' 'Npeaks'} };
geo.Npeaks = {[1 .2]};
% GFPthresh
style.GFPthresh = 'edit';
thresh_tipstr = ['Set to avoid extreme GFP peaks that are likely caused by '...
'artifacts. Set to zero to turn off.'];
line.GFPthresh = { {'Style' 'text' 'string' 'Reject peaks over threshold (multiples of std(GFP)):' ...
'tooltipstring' thresh_tipstr}, ...
{'Style' style.GFPthresh 'string' ' 0 ' 'tag' 'GFPthresh'} };
geo.GFPthresh = {[1 .2]};
% Order inputs for GUI
geometry = [geo.info geo.datatype geo.aggregate_data geo.avgref ...
geo.normalise geo.gfp_info geo.MinPeakDist geo.Npeaks geo.GFPthresh];
uilist = [line.info line.datatype line.aggregate_data line.avgref ...
line.normalise line.gfp_info line.MinPeakDist line.Npeaks line.GFPthresh];
%% Create Popup
[~,~,~,pop_out] = inputgui( geometry, uilist, ...
'pophelp(''pop_micro_selectdata'');', ...
'Select and aggregates data for microstate segmentation -- pop_micro_selectdata()');
%% Interpret output from popup
if isstruct(pop_out)
settings = struct;
settings = interpret_popup(pop_out, settings, style, popmenu);
% select the datasets, which should be included in the analysis
if settings.aggregate_data
datasetNames = {ALLEEG.setname};
dataset_idx = listdlg2('ListString',datasetNames,'PromptString',...
{'Select datasets to aggregate:'});
else
dataset_idx = [];
end
settings.dataset_idx = dataset_idx;
settings = rmfield(settings,'aggregate_data');
% remove settings related spontaneous data
if ~strcmp(settings.datatype,'spontaneous')
settings = rmfield(settings,'MinPeakDist');
settings = rmfield(settings,'Npeaks');
settings = rmfield(settings,'GFPthresh');
end
else
settings = 'cancel';
end
end
% -------------- helper functions -------------- %
function settings = interpret_popup(pop_out, settings, style, popmenu)
% Interpret output from pop_up window, "pop_out", and arrange it in
% "settings" struct. The fields in "style" should be the same as in "pop_out"
% (defined as tags in inputgui.m). The struct popmenu is optional and only
% needed if popmenus are used in the pop_up window.
names = fieldnames(style);
for i = 1:length(names)
switch style.(names{i})
case 'edit'
if isempty(pop_out.(names{i})) % empty?
settings.(names{i}) = [];
else
settings.(names{i}) = eval(pop_out.(names{i}));
end
case 'checkbox'
settings.(names{i}) = pop_out.(names{i});
case 'popupmenu'
settings.(names{i}) = popmenu.(names{i}){pop_out.(names{i})};
end
end
end
function settings = check_settings(vargs)
%% check settings
% Checks optional inputs for pop_micro_selectdata() and enters them in
% settings struct. Undefined inputs is set to default values.
varg_check = { 'datatype' 'string' [] 'ERPconc';
'avgref' 'integer' [] 1;
'normalise' 'integer' [] 0;
'dataset_idx' 'real' [] []};
% Checking if datatype is defined
dat_ind = find(strcmp(vargs, 'datatype'));
if ~isempty(dat_ind)
if strcmp(lower(vargs{dat_ind+1}), 'spontaneous')
% Adding check for spontaneous data related settings
varg_check = [varg_check;
{'MinPeakDist' 'real' [] 10;
'Npeaks' 'integer' [] 1000;
'GFPthresh' 'real' [] 0 } ];
end
end
settings = finputcheck( vargs, varg_check);
if ischar(settings), error(settings); end; % check for error
end
function NewEEG = check_datasets_consistency(ALLEEG)
%% checks if basic parameters of EEG data are consistent between datasets.
%% and creates new EEG structure that will later be filled depending on the type of data (spontaneous, ERP)
if isequal(ALLEEG(:).nbchan) && isequal(ALLEEG(:).srate); % number of channels and sampling rate
% Create a clean EEG structure that is used to store data.
NewEEG = eeg_emptyset();
NewEEG.srate = ALLEEG(1).srate;
NewEEG.nbchan = ALLEEG(1).nbchan;
else
error('The datasets differ in the number of channels and/or the sampling rate') ;
end
end
function com = settings_to_string(com,settings)
% Adds settings struct to existing com string in the form 'key1', 'val1',
% 'key2', 'val2' ... .
% Can handle structs, strings, vectors and scalars. I.e. not matrices.
names = fieldnames(settings);
for i = 1:length(names)
if isstruct(settings.(names{i})) % struct?
com = settings_to_string(com,settings.(names{i}));
elseif isempty(settings.(names{i})) % empty?
com = [ com sprintf(', ''%s'', []', names{i}) ];
elseif ischar(settings.(names{i})) % string?
com = [ com sprintf(', ''%s'', ''%s''', names{i}, settings.(names{i})) ];
elseif length(settings.(names{i})) > 1 % vector?
N_elements = length(settings.(names{i}));
range = max(settings.(names{i})) - min(settings.(names{i})) + 1;
if N_elements == range % write vetor as 'min_value:max_value'
com = [ com sprintf(', ''%s'', %g:%g', names{i}, ...
min(settings.(names{i})), max(settings.(names{i}))) ];
else % write vector with individual elements
com = [ com sprintf(', ''%s'', [%g', names{i}, ...
settings.(names{i})(1))];
for n = 2:N_elements
com = [ com sprintf(',%g', settings.(names{i})(n))];
end
com = [ com ']'];
end
else % scalar
com = [ com sprintf(', ''%s'', %g', names{i}, settings.(names{i})) ];
end
end
end