diff --git a/code_review/ab/pop_subcomp.m b/code_review/ab/pop_subcomp.m new file mode 100755 index 0000000..822d482 --- /dev/null +++ b/code_review/ab/pop_subcomp.m @@ -0,0 +1,236 @@ +% POP_SUBCOMP - remove specified components from an EEG dataset. +% and subtract their activities from the data. Else, +% remove components already marked for rejection. When used +% with the options 'keepcomp', the function will retain [1] +% or reject[0] the components provided as input. +% Usage: +% >> OUTEEG = pop_subcomp( INEEG ); % pop-up window mode +% >> OUTEEG = pop_subcomp( INEEG, components, plotag); +% >> OUTEEG = pop_subcomp( INEEG, components, plotag, keepcomp); +% +% Pop-up window interface: +% "Component(s) to remove ..." - [edit box] Array of components to +% remove from the data. Sets the 'components' parameter +% in the command line call (see below). +% "Component(s) to retain ..." - [edit box] Array of components to +% to retain in the data. Sets the 'components' parameter in +% the command line call. Then, comp_to_remove = ... +% setdiff([1:size(EEG.icaweights,1)], comp_to_keep). See +% option 'keepcomp' for command line call. +% Overwrites "Component(s) to remove" (above). +% Command line inputs: +% INEEG - Input EEG dataset. +% components - Array of components to remove from the data. If empty, +% remove components previously marked for rejection (e.g., +% EEG.reject.gcompreject). +% plotag - [0|1] Display the difference between original and processed +% dataset. 1 = Ask for confirmation. 0 = Do not ask. {Default: 0} +% keepcomp - [0|1] If [1] will retain the components provided by the +% input variable 'components', [0] will reject them. Option +% intended to be used only with the components provided in +% 'components', and not with components marked for rejection +% {Default: 0} +% Outputs: +% OUTEEG - output dataset. +% +% Author: Arnaud Delorme, CNL / Salk Institute, 2001 +% +% See also: COMPVAR + +% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu +% +% This file is part of EEGLAB, see http://www.eeglab.org +% for the documentation and details. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: +% +% 1. Redistributions of source code must retain the above copyright notice, +% this list of conditions and the following disclaimer. +% +% 2. Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +% THE POSSIBILITY OF SUCH DAMAGE. + +% 01-25-02 reformated help & license -ad +% 02-15-02 propagate ica weight matrix -ad sm jorn + +function [EEG, com] = pop_subcomp( EEG, components, plotag, keepflag) + +com=''; +if nargin < 1 + help pop_subcomp; + return; +end +if nargin < 3 + plotag = 0; +end +if nargin < 4 + keepflag = 0; +end + +if nargin < 2 + components = []; + % popup window parameters + % ----------------------- + res = 'Manual rej.'; + if any(cellfun(@(x)any(x.gcompreject), { EEG.reject })) + if length(EEG) == 1 + compStr = sprintf('%d,', find(EEG.reject.gcompreject == 1)); + msg = sprintf('Components [%s] flagged for rejection.', compStr(1:end-1)); + res = questdlg2(strvcat(msg, 'Do you want to remove these components?', ' ', 'Note: we recommend removing components in STUDY instead', '(in the STUDY channel pre-computing interface, select', '"Remove ICA artifactual components pre-tagged in each dataset")'), 'Remove components from data', 'Cancel', 'Manual rej.', 'Yes', 'Cancel'); + else + msg = 'Components flagged for rejection detected in some datasets.'; + res = questdlg2(strvcat(msg, 'Do you want to remove these components?', ' ', 'Note: we recommend removing components in STUDY instead', '(in the STUDY channel pre-computing interface, select', '"Remove ICA artifactual components pre-tagged in each dataset")'), 'Remove components from data', 'Cancel', 'Yes', 'Cancel'); + end + if strcmpi(res, 'Cancel') + return; + end + if strcmpi(res, 'Yes') + components = ''; + end + else + if length(EEG) > 1 + warndlg2(strvcat('You have multiple datasets selected and no components', ... + 'flagged for rejection. Flag components first.')); + return; + end + end + + if strcmpi(res, 'Manual rej.') + if ~isempty(EEG.reject.gcompreject) + components = find(EEG.reject.gcompreject == 1); + components = components(:)'; + %promptstr = { ['Components to subtract from data' 10 '(default: pre-labeled components to reject):'] }; + else + components = []; + end + uilist = { { 'style' 'text' 'string' 'Note: for group level analysis, remove components in STUDY' } ... + { 'style' 'text' 'string' 'List of component(s) to remove from data' } ... + { 'style' 'edit' 'string' int2str(components) } ... + { 'style' 'text' 'string' 'Or list of component(s) to retain' } ... + { 'style' 'edit' 'string' '' } ... + }; + geom = { 1 [2 0.7] [2 0.7] }; + result = inputgui( 'uilist', uilist, 'geometry', geom, 'helpcom', 'pophelp(''pop_subcomp'')', ... + 'title', 'Remove components from data -- pop_subcomp()'); + if isempty(result) + return; + end + components = eval( [ '[' result{1} ']' ] ); + if ~isempty(result{2}) + componentsOld = components; + components = eval( [ '[' result{2} ']' ] ); + if isequal(components, componentsOld) + components = []; + end + keepflag = 1; %components = setdiff_bc([1:size(EEG.icaweights,1)], components); + end + end +end + +% process multiple datasets +% ------------------------- +if length(EEG) > 1 + if nargin < 2 + [ EEG, com ] = eeg_eval( 'pop_subcomp', EEG, 'params', { components, plotag, keepflag }, 'warning', 'on' ); + else + [ EEG, com ] = eeg_eval( 'pop_subcomp', EEG, 'params', { components, plotag, keepflag } ); + end + if isempty( components ) + com = [ com ' % [] or '' means removing components flagged for rejection' ]; + end + return; +end + +componentsOri = components; +if isempty(components) + if ~isempty(EEG.reject.gcompreject) + components = find(EEG.reject.gcompreject == 1); + else + fprintf('Warning: no components specified, no rejection performed\n'); + return; + end +else + if keepflag == 1; components = setdiff_bc([1:size(EEG.icaweights,1)], components); end + if (max(components) > size(EEG.icaweights,1)) || min(components) < 1 + error('Component index out of range'); + end +end + +fprintf('Computing projection and removing %d components ....\n', length(components)); +component_keep = setdiff_bc(1:size(EEG.icaweights,1), components); +compproj = EEG.icawinv(:, component_keep)*eeg_getdatact(EEG, 'component', component_keep, 'reshape', '2d'); +compproj = reshape(compproj, size(compproj,1), EEG.pnts, EEG.trials); + +% fprintf( 'The ICA projection accounts for %2.2f percent of the data\n', 100*varegg); + +if nargin < 2 || plotag ~= 0 + + ButtonName = 'continue'; + while ~strcmpi(ButtonName, 'Cancel') && ~strcmpi(ButtonName, 'Accept') + ButtonName=questdlg2( [ 'Please confirm your choice. Are you sure you want to remove the selected components from the data?' ], ... + 'Confirmation', 'Cancel', 'Plot ERPs', 'Plot single trials', 'Accept', 'Accept'); + if strcmpi(ButtonName, 'Plot ERPs') + if EEG.trials > 1 + tracing = [ squeeze(mean(compproj,3)) squeeze(mean(EEG.data(EEG.icachansind,:,:),3))]; + figure; + plotdata(tracing, EEG.pnts, [EEG.xmin*1000 EEG.xmax*1000 0 0], ... + 'Trial ERPs (blue) with and (red) without these components'); + else + warndlg2('Cannot plot ERPs for continuous data'); + end + elseif strcmpi(ButtonName, 'Plot single trials') + eegplot( EEG.data(EEG.icachansind,:,:), 'srate', EEG.srate, 'title', 'Black = channel before rejection; red = after rejection -- eegplot()', ... + 'limits', [EEG.xmin EEG.xmax]*1000, 'data2', compproj); + end + end + switch ButtonName + case 'Cancel' + disp('Operation cancelled'); + return; + case 'Accept' + disp('Components removed'); + end % switch +end +EEG.data(EEG.icachansind,:,:) = compproj; +EEG.setname = [ EEG.setname ' pruned with ICA']; +EEG.icaact = []; +goodinds = setdiff_bc(1:size(EEG.icaweights,1), components); +EEG.icawinv = EEG.icawinv(:,goodinds); +EEG.icaweights = EEG.icaweights(goodinds,:); +EEG.specicaact = []; +EEG.specdata = []; +EEG.reject = []; + +if isfield(EEG.etc, 'ic_classification') + if isfield(EEG.etc.ic_classification, 'ICLabel') + if isfield(EEG.etc.ic_classification.ICLabel, 'classifications') + if ~isempty(EEG.etc.ic_classification.ICLabel.classifications) + EEG.etc.ic_classification.ICLabel.classifications = EEG.etc.ic_classification.ICLabel.classifications(goodinds,:); + end + end + end +end + +try + EEG.dipfit.model = EEG.dipfit.model(goodinds); +catch, end + +com = sprintf('EEG = pop_subcomp( EEG, [%s], %d);', int2str(componentsOri(:)'), plotag); +if isempty( components ) + com = [ com ' % [] means removing components flagged for rejection' ]; +end +return; diff --git a/code_review/ab/preprocessing_pipeline.m b/code_review/ab/preprocessing_pipeline.m new file mode 100755 index 0000000..0b69937 --- /dev/null +++ b/code_review/ab/preprocessing_pipeline.m @@ -0,0 +1,320 @@ +function [] = preprocessing_pipeline(inputfile, outpath, lowBP, topBP, FLAG, condition, task, varargin) + +% what file are we using +if ~exist(inputfile,'file'), error('inputfile "%s" does not exist!', inputfile), end +[d, currentName, ext ] = fileparts(inputfile); + +%% cap locations +eeglabpath = fileparts(which('eeglab')); +cap_location = fullfile(eeglabpath,'/plugins/dipfit/standard_BESA/standard-10-5-cap385.elp'); +if ~exist(cap_location, 'file'), error('cannot find file for 128 channel cap: %s', cap_location), end +correction_cap_location = hera('Projects/7TBrainMech/scripts/eeg/Shane/resources/ChanLocMod128to64.ced'); +if ~exist(correction_cap_location, 'file'), error('cannot find file for correction 128 channel cap: %s', correction_cap_location), end + +%% Files +subj_files = file_locs(inputfile, outpath, task); + +% to know how far your script is with running +fprintf('==========\n%s:\n\t Initial Preprocessing(%s,%f,%f,%s)\n',... + currentName, inputfile, lowBP, topBP, outpath) + + +% where to save things +filter_folder = 'filtered'; +chanrj_folder = 'channels_rejected'; +epoch_folder = 'epoched'; +icawholein_folder ='rerefwhole'; +epoch_rj_marked_folder = 'marked_epochs'; +epochrj_folder = 'rejected_epochs'; +icaout = fullfile(outpath, 'ICA'); +icawholeout = fullfile(outpath, 'ICAwhole'); +icawholeouthomog = fullfile(outpath, 'AfterWhole/ICAwholeClean_homogenize'); + +% and what files will we create +rerefwhole_name = [currentName '_rerefwhole']; +chrm_name = [currentName '_badchannelrj']; +epochrj_name = [currentName '_epochs_rj']; +% FIXME: these is not actually used!? but is recoreded in data_removed WF20190911 +datarm_name = [currentName '_baddatarj']; +epochrm_name = [currentName '_badepochrj']; +icawholeout_name = [currentName '_rerefwhole_ICA']; +% TODO: collect other pop_saveset filenames here + +commonPlus = {'AFz','C1','C2','C3','C4','C5','C6','CP1','CP2','CP3','CP4',... + 'CP5','CP6','CPz','Cz','F1','F2','F3','F4','F5','F6','F7','F8','FC1',... + 'FC2','FC3','FC4','FC5','FC6','FCz','Fp1','Fp2','FT10','FT9','Fz','I1',... + 'I2','O1','O2','Oz','P1','P10','P2','P3','P4','P5','P6','P7','P8','P9',... + 'PO10','PO3','PO4','PO7','PO8','PO9','POz','Pz','T7','T8',... + 'AF8','AF7','AF4','AF3'}; + +epochrj = fullfile(outpath, epochrj_folder, [epochrj_name '.set']); + +if condition == 1 + icawholeoutFile = fullfile(icawholeout, [icawholeout_name '.set']); + +else + icawholeoutFile = 'no'; +end + +if exist(icawholeoutFile, 'file') + warning('%s already complete (have "%s")! todo load from file', currentName, icawholeout_name) + return +end + +if condition == 1 + xEEG = load_if_exists(subj_files.filter); +end + +if isstruct(xEEG) + [ALLEEG EEG CURRENTSET] = pop_newset([], xEEG, 0,... + 'setname',currentName,... + 'gui','off'); +else + + %% load EEG set and re- referance + EEG = pop_loadset(inputfile); + + if size(EEG.data,1) < 100 + Flag128 = 0; + EEG = pop_reref(EEG, [65 66]); %does this "restore" the 40dB of "lost SNR" ? -was it actually lost? ...this is potentially undone by PREP + EEG = eeg_checkset(EEG); + % find the cap to use + else + %[129 130] are the mastoid externals for the 128 electrode + EEG = pop_reref(EEG, [129 130]); %does this "restore" the 40dB of "lost SNR" ? -was it actually lost? ...this is potentially undone by PREP + EEG = eeg_checkset(EEG); + end + + %stores EEG set in ALLEEG, give setname + ALLEEG = []; + [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 0,... + 'setname',currentName,... + 'gui','on'); + [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG); + + + EEG.subject = currentName(1:findstr(currentName,'mgs')-2); + EEG.condition = currentName(findstr(currentName,'mgs'):end); + + %% Filtering + + %band-pass filter between low and max boundle (or 1 and 90 Hz) + % TODO: params are + % EEG, locutoff, hicutoff, filtorder, revfilt, usefft, plotfreqz, minphase); + % why 3380, 0, [], 0 + % > Warning: Transition band is wider than maximum stop-band width. + % For better results a minimum filter order of 6760 is recommended. + % Reported might deviate from effective -6dB cutoff frequency. + % [FLAGexist] = checkdone(fullfile(outpath,filter_folder,[currentName '_filtered'])); + + EEG = pop_eegfiltnew(EEG, lowBP, topBP, 3380, 0, [], 0); + % filtorder = 3380 - filter order (filter length - 1). Mandatory even. performing 3381 point bandpass filtering. + % pop_eegfiltnew() - transition band width: 0.5 Hz + % pop_eegfiltnew() - passband edge(s): [0.5 90] Hz + % pop_eegfiltnew() - cutoff frequency(ies) (-6 dB): [0.25 90.25] Hz + % pop_eegfiltnew() - filtering the data (zero-phase) + + %give a new setname and overwrite unfiltered data + EEG = pop_editset(EEG,'setname',[currentName '_bandpass_filtered']); + [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG); + + + + + %% Resample data + + % Downsample the data to 150Hz using antialiasing filter + if task == "MGS" || task == "Resting_State" || task == "anti" + EEG = pop_resample(EEG, 150, 0.8, 0.4); %0.8 is fc and 0.4 is dc. Default is .9 and .2. We dont know why Alethia changed them + elseif task == "SNR" + EEG = pop_resample(EEG, 512, 0.8, 0.4); %0.8 is fc and 0.4 is dc. Default is .9 and .2. We dont know why Alethia changed them + end + + + EEG = eeg_checkset(EEG); + + %change setname + EEG = pop_editset(EEG,'setname',[currentName '_filtered']); + [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); + %save filtered data + EEG = pop_saveset( EEG, 'filename',[currentName '_filtered'], ... + 'filepath',fullfile(outpath,filter_folder)); + [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG); +end + +%% CHANNELS +% remove external channels +EEG = pop_select( EEG,'nochannel',{'EX3' 'EX4' 'EX5' 'EX6' 'EX7' 'EX8' 'EXG1' 'EXG2' 'EXG3' 'EXG4' 'EXG5' 'EXG6' 'EXG7' 'EXG8' 'GSR1' 'GSR2' 'Erg1' 'Erg2' 'Resp' 'Plet' 'Temp' 'FT7' 'FT8' 'TP7' 'TP8' 'TP9' 'TP10'}); + +[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); +%import channel locations + +%change this maybe +% eeglab should have already been added with addpath + +EEG=pop_chanedit(EEG, 'lookup', cap_location); + +if size(EEG.data,1) > 100 + EEG = pop_select( EEG,'channel',commonPlus); + EEG = pop_chanedit(EEG, 'load', {correction_cap_location 'filetype' 'autodetect'}); + % 128 'AF8' --> 64 'AF6' + % 128 'AF7' --> 64 'AF5' + % 128 'AF4' --> 64 'AF2' + % 128 'AF3' --> 64 'AF1' + Flag128 = 1; +end +[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); + +%% bad channel rejection + +%bad channels are in general the channels that you can't save even if you +%reject 5-10% of your datapoints + +%different options for channel rejection are displayed. Option 3 is ued. + +% %1.look at standard deviation in bar plots and remove the channels +% with big std's --> manual +% stdData = std(EEG.data,0,2); +% figure(idx); bar(stdData) + +% %2.kurtosis +% EEG = pop_rejchan(EEG, 'elec',[1:8] ,'threshold',5,'norm','on','measure','kurt'); + +%3.clean_rawdata +if condition == 1 + + xEEG = load_if_exists(subj_files.chanrj); +end + +originalEEG = EEG; +if isstruct(xEEG) + [ALLEEG EEG] = eeg_store(ALLEEG, xEEG, CURRENTSET); +else + EEG = clean_rawdata(EEG, 8, [0.25 0.75], 0.7, 5, 15, 0.3); % we can discuss that + % vis_artifacts(EEG,originalEEG,'NewColor','red','OldColor','black'); + %change setname + EEG = pop_editset(EEG,'setname', chrm_name); + [ALLEEG EEG] = eeg_store(ALLEEG, EEG); + + EEG = pop_saveset( EEG,'filename', chrm_name, ... + 'filepath',sprintf('%s/%s/',outpath,chanrj_folder)); +end + +if condition == 1 + xEEG = load_if_exists(subj_files.rerefwhole_name); +end + +if isstruct(xEEG) + [ALLEEG EEG] = eeg_store(ALLEEG, xEEG, CURRENTSET); +else + [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); + + if ~any(find(cellfun (@any,regexpi (fieldnames(EEG.etc), 'clean_channel_mask')))); + EEG.etc.clean_channel_mask=42; + else + end + + %save the channels that were rejected in a variable + channels_removed{1} = chrm_name; %setname + channels_removed{2} = setdiff({originalEEG.chanlocs.labels},{EEG.chanlocs.labels}, 'stable'); + channels_removed{3} = find(EEG.etc.clean_channel_mask==0); + %also save the channels that were rejected in the EEG struct + EEG.channels_rj = channels_removed{2}; + EEG.channels_rj_nr = length(EEG.channels_rj); + + %save the proportion of the dataset that were rejected in a variable + data_removed{1} = datarm_name; %setname + data_removed{2} = length(find(EEG.etc.clean_sample_mask==0))/EEG.srate;% + data_removed{3} = length(find(EEG.etc.clean_sample_mask==0))/length(EEG.etc.clean_sample_mask);% + %also save the data that were rejected in the EEG struct + EEG.data_rj = data_removed{2}; + EEG.data_rj_nr = data_removed{3}; + + %% interpolate channels + % POSSIBLE PROBLEMS + % - injecting extra channels ontop of expected 64 (n>64) + % - 128 missing expected labels, adding too few back (n<64) + if Flag128 == 1 + nchan = 64; + ngood = length(EEG.chanlocs); + % 128 cap doesn't have exactly the same postions as 64 + % remove 4 that are in the wrong place and reinterpret + % AND interp any bad channels + % do this by removing the 4 128weirdos + % from the already trimmed (no bad channels) in EEG.chanlocs + + need_128interp = [2 3 35 36 ]; + % get the names of those to remove + n128name = {originalEEG.chanlocs(need_128interp).labels}; + % should always be {'AF5','AF1','AF2','AF6'} ?? + + % find where they are in current EEG files (if they haven't already been removed) + n128here_idx = find(ismember({EEG.chanlocs.labels},n128name)); + % keep those that aren't the ones we matched + % remove from chanlocs, data and update nbcan + % WARNING -- who knows what else we should have changed to update the set info! + keep_idx = setdiff(1:ngood, n128here_idx); + EEG.chanlocs = EEG.chanlocs(keep_idx); + EEG.data = EEG.data(keep_idx,:); + EEG.nbchan = length(keep_idx); + + %EEG_i = pop_interp(EEG, interp_ch, 'spherical'); + fprintf('%d channels in orig; want to interpolate %d bad and move %d\n',... + originalEEG.nbchan, nchan - ngood, length(need_128interp)) + EEG_i = pop_interp(EEG, originalEEG.chanlocs, 'spherical'); + + % could swap these channels (they're close, but not the same) + % BUT WE DONT + % 128 'AF7' --> 64 'AF5' In this point channel 2 + % 128 'AF3' --> 64 'AF1' In this point channel 3 + % 128 'AF4' --> 64 'AF2' In this point channel 35 + % 128 'AF8' --> 64 'AF6' In this point channel 36 + % lines above modify channel information and pocition in data to make + % it the same for 64 and 128 cap + + % need to do destructive swapping. need a copy + EEG = EEG_i; + EEG.chanlocs(2) = EEG_i.chanlocs(3);%EEG.chanlocs(2) must by 'AF1' in 64 cap + EEG.chanlocs(3) = EEG_i.chanlocs(2);%EEG.chanlocs(3) must by 'AF5' in 64 cap + % EEG.chaninfo.filecontent(4,:) = EEG_i.chaninfo.filecontent(3,:); This + % is not necessary i think, but just in case... + % EEG.chaninfo.filecontent(4,1) = '3'; + % EEG.chaninfo.filecontent(3,:) = EEG_i.chaninfo.filecontent(4,:); + % EEG.chaninfo.filecontent(3,1) = '2'; + EEG.data(2,:) = EEG_i.data(3,:);% ALERT ALERT Lines latelly added + EEG.data(3,:) = EEG_i.data(2,:);% ATERT ALERT Lines latelly added + else + EEG = pop_interp(EEG, originalEEG.chanlocs, 'spherical'); + end + % eeglab redraw + %% re-reference: average reference + if FLAG + EEG = pop_reref( EEG, []); + [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 1,... + 'setname',[currentName '_avref'],... + 'gui','off'); + + else + end + + %save whole rereferenced data for ICA whole + %save epochs rejected EEG data + EEG = pop_editset(EEG, 'setname', rerefwhole_name); + [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); + + EEG = pop_saveset(EEG, 'filename', rerefwhole_name, 'filepath', fullfile(outpath,icawholein_folder)); + [ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); +end + +%% Whole data ICA run + +icawholein = fullfile(outpath, icawholein_folder, [rerefwhole_name '.set']); + +if ~exist(subj_files.icawhole, 'file') + runICAs(icawholein,icawholeout,task) +else + fprintf('have %s, not rerunning\n', subj_files.icawholeout) +end + + + diff --git a/code_review/ab/runICAs.m b/code_review/ab/runICAs.m new file mode 100755 index 0000000..ccad9bb --- /dev/null +++ b/code_review/ab/runICAs.m @@ -0,0 +1,114 @@ +function [] = runICAs(inpath, outpath, task, varargin) +%% this script runs ICA on all filtered/cleaned/epoched EEG dmt data. +% PCA is used to decrease the number of components +% inpath like ....rerefwhole/*_rerefwhole.set +% ouptath like .../ICAwhole/ +% varargin can contain 'redo' to ignore file already done + + +% could use files from file_locs instead of defining here +files = file_locs(inpath, outpath, task); + +% -- icawhole defined elswhere like: +% icawholeout = fullfile(outpath, 'ICAwhole'); + + +%% find file +% if no file, try searching with *.set +if exist(inpath,'file') + EEGfileNames = dir(inpath); +else + EEGfileNames = dir([inpath, '*.set']); +end +currentEEG = EEGfileNames(1).name; +[~, name, ext] = fileparts(currentEEG); + +% did we already run? +% 20191220 - never creates, ICA_SAS code commented out -- why?! +% - instead we'll check for the actual final file _ICA.set +finalout = files.icawhole; +if exist(finalout, 'file') && ~isempty(strncmp('redo', varargin, 4)) + warning('already created %s, not running. add "redo" to ICA call to redo', finalout) + return +end + +%% run ICA +eeglab + +% load data to get bad channels and check channel # +EEG = pop_loadset('filename', currentEEG, 'filepath', fileparts(inpath)); +% TODO: just use clean_channel_mask +schans = {name, EEG.channels_rj, find(EEG.etc.clean_channel_mask==0)}'; +badchans = schans{3,1}; + +% if size(EEG.data,1) > 100 +% error('not 64 channel, not sure positions, skipping') +% end + +disp(currentEEG); + +%control for rank defficiency: if you want to analyse components as +%(and not only use the components to remove noise) you should keep the +%nr of components contstant across datasets. But there should not be +%more components than channels. (So check how many channels are removed). + +% PCAnr = 15; %for example 10 components because in 1 dataset 14 channels were removed +% %PCAnr = EEG.nbchan - EEG.channels_rj_nr; %or if you are not interested +% %in component analysis but want to use ICA only for artifact removal +% %change the PCAnr for each dataset individually + +ALLEEG = []; +[ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, 0, 'setname', name, 'gui','off'); + +% %% SASICA Command +% pop_runica +% %Then the user would pull up ICA'd results, and run SAS: +% EEG = eeg_SASICA(EEG, ... +% 'MARA_enable',0, ... +% 'FASTER_enable',0,... +% 'FASTER_blinkchanname','EX5', ... +% 'ADJUST_enable',0,... +% 'chancorr_enable',1,... +% 'chancorr_channames','No channel',... +% 'chancorr_corthresh','auto 4',... +% 'EOGcorr_enable',0,... +% 'EOGcorr_Heogchannames','No channel', ... +% 'EOGcorr_corthreshH','auto 4',... +% 'EOGcorr_Veogchannames','No channel',... +% 'EOGcorr_corthreshV','auto 4',... +% 'resvar_enable',0,... +% 'resvar_thresh',15,... +% 'SNR_enable',1,... +% 'SNR_snrcut',1,...q +% 'SNR_snrBL',[-Inf 0] ,... +% 'SNR_snrPOI',[0 Inf] ,... +% 'trialfoc_enable',0,... +% 'trialfoc_focaltrialout','auto',... +% 'focalcomp_enable',1,... +% 'focalcomp_focalICAout','auto',... +% 'autocorr_enable',1,... +% 'autocorr_autocorrint',20,... +% 'autocorr_dropautocorr','auto',... +% 'opts_nocompute',0,... +% 'opts_FontSize',14, ... +% 'opts_noplot', 1); +% +% fprintf('saving %s\n', fullfile(outpath, [name '_SAS.set'])) +% EEG = pop_saveset( EEG, 'filename',[name,'_SAS.set'], 'filepath', outpath); %save final preprocessed output + + +EEG = pop_runica(EEG, 'extended',1,'interupt','on','PCA',30); + +%create file name +name = EEG.setname; +%name = EEG.filename(1:end-4); + +name = [name, '_ICA']; + +%change set name +fprintf('saving %s\n', fullfile(outpath, [name '.set'])) +EEG = pop_editset(EEG, 'setname', name); +[ALLEEG EEG] = eeg_store(ALLEEG, EEG, CURRENTSET); +EEG = pop_saveset(EEG, 'filename', [name '.set'], 'filepath', outpath); + +end diff --git a/code_review/ab/run_preprocessing_pipeline.m b/code_review/ab/run_preprocessing_pipeline.m new file mode 100755 index 0000000..cbd30e3 --- /dev/null +++ b/code_review/ab/run_preprocessing_pipeline.m @@ -0,0 +1,122 @@ +function [] = run_preprocessing_pipeline(task) + +% load in paths to feildtrip and eeglab +addpath(genpath('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/Preprocessing_Functions/')); +addpath(genpath('/resources/Euge/')) +% run('/Volumes/Hera/Projects/7TBrainMech/scripts/eeg/Shane/resources/eeglab2022.1/eeglab.m'); +run('/Volumes/Hera/Abby/Resources/eeglab_current/eeglab2024.2/eeglab.m') + +% Outpath +% maindir = hera('Projects/7TBrainMech/scripts/eeg/Shane/preprocessed_data'); +maindir = hera('Abby/preprocessed_data'); + +disp(["i am running" task]) + +taskdirectory = [maindir '/' task]; + +% initial values +lowBP = 0.5; +topBP = 70; +FLAG = 1; + +%% settings +only128 = 0; % 0==do all, 1==only 128 channel subjects +condition = 1; %0 - if you want to overwrite an already existing file; 1- if you want it to skip subjects who have already been run through singlesubject +dryrun = 0; + +remark(task, taskdirectory, dryrun); % change the trigger values to be single digit + + +% gather all file paths for all subjects remarked data +setfilesDir = [taskdirectory, '/remarked/1*_20*.set']; +setfiles = all_remarked_set(setfilesDir); + +n = size(setfiles,1); %number of EEG sets to preprocess + +% loop through every subject to preprocess +for i = 1:n + inputfile = setfiles{i}; + try + preprocessing_pipeline(inputfile, taskdirectory, lowBP, topBP, FLAG, condition, task) + catch e + fprintf('Error processing "%s": %s\n',inputfile, e.message) + for s=e.stack + disp(s) + end + end +end + +%% select ICA values to reject + +ICA_Path = [taskdirectory '/ICAwhole']; +CleanICApath = [taskdirectory '/AfterWhole/ICAwholeClean/']; + + +EEGfileNames = dir([ICA_Path '/*.set']); + +for fidx = 1:length(EEGfileNames) + filename = EEGfileNames(fidx).name; + locs = file_locs(fullfile(ICA_Path,filename), taskdirectory, task); + if exist(locs.ICAwholeClean, 'file') + fprintf('skipping; already created %s\n', locs.ICAwholeClean); + continue + end + + selectcompICA +end + + +%% Homogenize Chanloc +datapath = [taskdirectory '/AfterWhole/ICAwholeClean']; +savepath = [taskdirectory '/AfterWhole/ICAwholeClean_homogenize']; + +setfiles0 = dir([datapath,'/*icapru.set']); +setfiles = {}; +for epo = 1:length(setfiles0) + setfiles{epo,1} = fullfile(datapath, setfiles0(epo).name); % cell array with EEG file names +end + +correction_cap_location = hera('Projects/7TBrainMech/scripts/eeg/Shane/resources/ELchanLoc.ced'); +for i = 1:length(setfiles) + homogenizeChanLoc(setfiles{i},correction_cap_location,savepath, taskdirectory, task) +end + + +%% filter out the 60hz artifact from electronics +datapath = [taskdirectory '/AfterWhole/ICAwholeClean_homogenize']; + +setfiles0 = dir([datapath,'/*icapru.set']); +setfiles = {}; +for epo = 1:length(setfiles0) + setfiles{epo,1} = fullfile(datapath, setfiles0(epo).name); % cell array with EEG file names +end + +for i = 1:length(setfiles) + inputfile = setfiles{i}; + [filepath,filename ,ext] = fileparts((setfiles{i})); + + EEG = pop_loadset(inputfile); + EEG = pop_eegfiltnew(EEG, 59, 61, [], 1, [], 0); + EEG = pop_saveset(EEG, 'filename', filename, 'filepath', datapath); +end + + +if task == "MGS" + %% Clean epochs to remove + + epoch_path = [outpath '/AfterWhole/ICAwholeClean/']; + epoch_folder = [outpath '/AfterWhole/epoch/']; + epoch_rj_marked_folder = [outpath '/AfterWhole/epochclean/']; + + EEGfileNames = dir([path_data, '/*_icapru.set']); + + revisar = {}; + for currentEEG = 1:size(EEGfileNames,1) + filename = [EEGfileNames(currentEEG).name]; + inputfile = [epoch_path,filename]; + revisar{currentEEG} = epochlean(inputfile,epoch_folder,epoch_rj_marked_folder); + end +end + + + diff --git a/code_review/ab/selectcompICA.m b/code_review/ab/selectcompICA.m new file mode 100755 index 0000000..01d47d2 --- /dev/null +++ b/code_review/ab/selectcompICA.m @@ -0,0 +1,26 @@ + +% quick check that ICA_Path is valid +if ~ exist(ICA_Path,'dir'); error('no ICA_Path directory "%s"!', ICA_Path); end + +% prefix ='_prep_allepoch_ica.set' +EEG = pop_loadset('filename',filename,'filepath',ICA_Path); +% EEG.setname = ['UG_' suj_n, '_',prefix(1:end-4)]; +EEG = eeg_checkset( EEG ); + +%% DB point on 12 +[EEG, com] = pop_selectcomps(EEG, [1:30] );pop_eegplot(EEG, 0, 1, 1); +pop_saveset( EEG, 'filename',[filename(1:end-4),'_pesos.set'],'filepath',CleanICApath); +%pop_eegplot( EEG, 0, 1, 1); +eeglab redraw +lowBP +% TEXTO = 'Lista la sel de comp? (If not put [0])'; +% interp_user_def = input(TEXTO); +% +% if interp_user_def~= 0 +cmps = find(EEG.reject.gcompreject); +eeglab redraw +OUTEEG = pop_subcomp( EEG, [cmps], 1); +eeglab redraw +pop_saveset( OUTEEG, 'filename',[filename(1:end-4),'_icapru.set'],'filepath',CleanICApath); +% end +close all diff --git a/code_review/ab/spectopo.m b/code_review/ab/spectopo.m new file mode 100755 index 0000000..c8307c9 --- /dev/null +++ b/code_review/ab/spectopo.m @@ -0,0 +1,993 @@ +% SPECTOPO - Plot the power spectral density (PSD) of winsize length segments of data +% epochs at all channels as a bundle of traces. At specified frequencies, +% plot the relative topographic distribution of PSD. If available, uses +% PWELCH from the Matlab signal processing toolbox, else the EEGLAB SPEC +% function. Plots the mean spectrum for all of the supplied data, not just +% the pre-stimulus baseline. +% Usage: +% >> spectopo(data, frames, srate); +% >> [spectra,freqs,speccomp,contrib,specstd] = ... +% spectopo(data, frames, srate, 'key1','val1', 'key2','val2' ...); +% Inputs: +% data = If 2-D (nchans,time_points); % may be a continuous single epoch, +% else a set of concatenated data epochs, else a 3-D set of data +% epochs (nchans,frames,epochs) +% frames = frames per epoch {default|0 -> data length} +% srate = sampling rate per channel (Hz) +% +% Optional 'keyword',[argument] input pairs: +% 'freq' = [float vector (Hz)] vector of frequencies at which to plot power +% scalp maps, or else a single frequency at which to plot component +% contributions at a single channel (see also 'plotchan') +% 'chanlocs' = [electrode locations filename or EEG.chanlocs structure] +% For format, see >> topoplot example +% 'limits' = [xmin xmax ymin ymax cmin cmax] axis limits. Sets x, y, and color +% axis limits. May omit final values or use NaNs. +% Ex: [0 60 NaN NaN -10 10], [0 60], ... +% Default color limits are symmetric around 0 and are different +% for each scalp map {default|all NaN's: from the data limits} +% 'title' = [quoted string] plot title {default: none} +% 'freqfac' = [integer] ntimes to oversample (to adjust frequency resolution) {default: 1} +% 'nfft' = [integer] Data points to zero-pad data windows to (overwrites 'freqfac') +% 'winsize' = [integer] window size in data points {default: Sampling Rate} +% 'overlap' = [integer] window overlap in data points {default: 0} +% 'percent' = [float 0 to 100] percent of the data to sample for computing the +% spectra. Values < 100 speed up the computation. {default: 100} +% 'freqrange' = [min max] frequency range to plot. Changes x-axis limits {default: +% 1 Hz for the min and Nyquist (srate/2) for the max. If specified +% power distribution maps are plotted, the highest mapped frequency +% determines the max freq}. +% 'wintype' = ['hamming','blackmanharris'] Window type used on the power spectral +% density estimation. The Blackman-Harris windows offers better attenuation +% than Hamming windows, but lower spectral resolution. {default: 'hamming'} +% 'blckhn' = [integer] Parameter to scale the windows length when Blackman-Harris +% windows is used in 'wintype' {default: 2} +% 'reref' = ['averef'|'off'] convert data to average reference {default: 'off'} +% 'mapnorm' = [float vector] If 'data' contain the activity of an independent +% component, this parameter should contain its scalp map. In this case +% the spectrum amplitude will be scaled to component RMS scalp power. +% Useful for comparing component strengths {default: none} +% 'boundaries' = data point indices of discontinuities in the signal {default: none} +% 'plot' = ['on'|'off'] 'off' -> disable plotting {default: 'on'} +% 'rmdc' = ['on'|'off'] 'on' -> remove DC {default: 'off'} +% 'plotmean' = ['on'|'off'] 'on' -> plot the mean channel spectrum {default: 'off'} +% 'plotchans' = [integer array] plot only specific channels {default: all} +% 'verbose' = ['on'|'off'] 'on' shows information on command line {default: 'on'} +% +% Optionally plot component contributions: +% 'weights' = ICA unmixing matrix. Here, 'freq' (above) must be a single frequency. +% ICA maps of the N ('nicamaps') components that account for the most +% power at the selected frequency ('freq') are plotted along with +% the spectra of the selected channel ('plotchan') and components +% ('icacomps'). +% 'plotchan' = [integer] channel at which to compute independent conmponent +% contributions at the selected frequency ('freq'). If 0, plot RMS +% power at all channels. {defatul|[] -> channel with highest power +% at specified 'freq' (above)). Do not confuse with +% 'plotchans' which select channels for plotting. +% 'mapchans' = [int vector] channels to plot in topoplots {default: all} +% 'mapframes'= [int vector] frames to plot {default: all} +% 'nicamaps' = [integer] number of ICA component maps to plot {default: 4}. +% 'icacomps' = [integer array] indices of ICA component spectra to plot ([] -> all). +% 'icamode' = ['normal'|'sub'] in 'sub' mode, instead of computing the spectra of +% individual ICA components, the function computes the spectrum of +% the data minus their contributions {default: 'normal'} +% 'icamaps' = [integer array] force plotting of selected ICA component maps +% {default: [] = the 'nicamaps' largest contributing components}. +% 'icawinv' = [float array] inverse component weight or mixing matrix. Normally, +% this is computed by inverting the ICA unmixing matrix 'weights' (above). +% However, if any components were removed from the supplied 'weights'mapchans +% then the component maps will not be correctly drawn and the 'icawinv' +% matrix should be supplied here {default: from component 'weights'} +% 'memory' = ['low'|'high'] a 'low' setting will use less memory for computing +% component activities, will take longer {default: 'high'} +% +% Replotting options: +% 'specdata' = [freq x chan array ] spectral data +% 'freqdata' = [freq] array of frequencies +% +% Topoplot options: +% other 'key','val' options are propagated to TOPOPLOT for map display +% (See >> help topoplot) +% +% Outputs: +% spectra = (nchans,nfreqs) power spectra (mean power over epochs), in dB +% freqs = frequencies of spectra (Hz) +% speccomp = component spectra (ncomps,nfreqs). Warning, only the function +% only computes the spectrum of the components given as input using +% the 'icacomps' parameter. Other component spectrum are filled +% with 0. +% contrib = contribution of each component (if 'icamode' is 'normal', relative +% variance, if 'icamode' is 'sub', percent variance accounted for) +% specstd = spectrum standard deviation in dB +% +% Notes: The original input format is still functional for backward compatibility. +% PSD has been replaced by PWELCH (see Matlab note 24750 on their web site) +% +% Non-backward compatible change (Nov 15 2015): +% Default winsize was set to the sampling rate (giving a default window +% length of 1 sec). Also, the y-axis label in the plot was corrected +% to read, "Log Power Spectral Density 10*log_{10}(\muV^{2}/Hz)' +% Finally, when winsize is not a power of 2, it is no longer promoted to +% the next higher power of 2. Thanks to Andreas Widmann for his comments. +% +% +% Authors: Scott Makeig, Arnaud Delorme & Marissa Westerfield, +% SCCN/INC/UCSD, La Jolla, 3/01 +% +% See also: TIMTOPO, ENVTOPO, TFTOPO, TOPOPLOT + +% Copyright (C) 3/01 Scott Makeig & Arnaud Delorme & Marissa Westerfield, SCCN/INC/UCSD, +% scott@sccn.ucsd.edu +% +% This file is part of EEGLAB, see http://www.eeglab.org +% for the documentation and details. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: +% +% 1. Redistributions of source code must retain the above copyright notice, +% this list of conditions and the following disclaimer. +% +% 2. Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF +% THE POSSIBILITY OF SUCH DAMAGE. + +% 3-20-01 added limits arg -sm +% 01-25-02 reformated help & license -ad +% 02-15-02 scaling by epoch number line 108 - ad, sm & lf +% 03-15-02 add all topoplot options -ad +% 03-18-02 downsampling factor to speed up computation -ad +% 03-27-02 downsampling factor exact calculation -ad +% 04-03-02 added axcopy -sm + +% Uses: MATLAB PWELCH, CHANGEUNITS, TOPOPLOT, TEXTSC + +function [eegspecdB,freqs,compeegspecdB,resvar,specstd] = spectopo(data,frames,srate,varargin) + +% formerly: ... headfreqs,chanlocs,limits,titl,freqfac, percent, varargin) + +icadefs; +LOPLOTHZ = 1; % low Hz to plot +FREQFAC = 1; % approximate frequencies/Hz (default) +allcolors = { [0 0.7500 0.7500] + [1 0 0] + [0 0.5000 0] + [0 0 1] + [0.2500 0.2500 0.2500] + [0.7500 0.7500 0] + [0.7500 0 0.7500] }; % colors from real plots }; + +if nargin<3 + help spectopo + return +end +if nargin <= 3 || ischar(varargin{1}) + % 'key' 'val' sequence + fieldlist = { 'freq' 'real' [] [] ; + 'specdata' 'real' [] [] ; + 'freqdata' 'real' [] [] ; + 'chanlocs' '' [] [] ; + 'freqrange' 'real' [0 srate/2] [] ; + 'memory' 'string' {'low','high'} 'high' ; + 'plot' 'string' {'on','off'} 'on' ; + 'plotmean' 'string' {'on','off'} 'off' ; + 'title' 'string' [] ''; + 'limits' 'real' [] [nan nan nan nan nan nan]; + 'freqfac' 'integer' [] FREQFAC; + 'percent' 'real' [0 100] 100 ; + 'reref' 'string' { 'averef','off','no' } 'off' ; + 'boundaries' 'integer' [] [] ; + 'nfft' 'integer' [1 Inf] [] ; + 'winsize' 'integer' [1 Inf] [] ; + 'overlap' 'integer' [1 Inf] 0 ; + 'icamode' 'string' { 'normal','sub' } 'normal' ; + 'weights' 'real' [] [] ; + 'mapnorm' 'real' [] [] ; + 'plotchan' 'integer' [1:size(data,1)] [] ; + 'plotchans' 'integer' [1:size(data,1)] [] ; + 'nicamaps' 'integer' [] 4 ; + 'blckhn' 'integer' [] 2 ; + 'icawinv' 'real' [] [] ; + 'icacomps' 'integer' [] [] ; + 'icachansind' 'integer' [] [1:size(data,1)] ; + 'icamaps' 'integer' [] [] ; + 'rmdc' 'string' {'on','off'} 'off'; + 'verbose' 'string' {'on','off'} 'on'; + 'wintype' 'string' {} 'hamming'; + 'mapchans' 'integer' [1:size(data,1)] [] + 'mapframes' 'integer' [1:size(data,2)] []}; + + [g, varargin] = finputcheck( varargin, fieldlist, 'spectopo', 'ignore'); + if ischar(g), error(g); end + if ~isempty(g.freqrange), g.limits(1:2) = g.freqrange; end + if ~isempty(g.weights) + if isempty(g.freq) || length(g.freq) > 2 + if ~isempty(get(0,'currentfigure')) && strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end + error('spectopo(): for computing component contribution, one must specify a (single) frequency'); + end + end +else + if ~isnumeric(data) + error('spectopo(): Incorrect call format (see >> help spectopo).') + end + if ~isnumeric(frames) || round(frames) ~= frames + error('spectopo(): Incorrect call format (see >> help spectopo).') + end + if ~isnumeric(srate) % 3rd arg must be the sampling rate in Hz + error('spectopo(): Incorrect call format (see >> help spectopo).') + end + + if nargin > 3, g.freq = varargin{1}; + else g.freq = []; + end + if nargin > 4, g.chanlocs = varargin{2}; + else g.chanlocs = []; + end + if nargin > 5, g.limits = varargin{3}; + else g.limits = [nan nan nan nan nan nan]; + end + if nargin > 6, g.title = varargin{4}; + else g.title = ''; + end + if nargin > 7, g.freqfac = varargin{5}; + else g.freqfac = FREQFAC; + end + if nargin > 8, g.percent = varargin{6}; + else g.percent = 100; + end + if nargin > 10, g.reref = 'averef'; + else g.reref = 'off'; + end + g.weights = []; + g.icamaps = []; +end +if g.percent > 1 + g.percent = g.percent/100; % make it from 0 to 1 +end +if ~isempty(g.freq) && isempty(g.chanlocs) + error('spectopo(): needs channel location information'); +end +if isempty(g.weights) && ~isempty(g.plotchans) + data = data(g.plotchans,:); + if ~isempty(g.chanlocs) + g.chanlocs = g.chanlocs(g.plotchans); + end +end + +if strcmpi(g.rmdc, 'on') + data = data - repmat(mean(data,2), [ 1 size(data,2) 1]); +end +data = reshape(data, size(data,1), size(data,2)*size(data,3)); + +if frames == 0 + frames = size(data,2); % assume one epoch +end + +%if ~isempty(g.plotchan) & g.plotchan == 0 & strcmpi(g.icamode, 'sub') +% if ~isempty(get(0,'currentfigure')) & strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end +% error('Cannot plot data component at all channels (option not implemented)'); +%end + +if ~isempty(g.freq) && min(g.freq)<0 + if ~isempty(get(0,'currentfigure')) && strcmp(get(gcf, 'tag'), 'spectopo'), close(gcf); end + fprintf('spectopo(): freqs must be >=0 Hz\n'); + return +end + +g.chanlocs2 = g.chanlocs; +if ~isempty(g.specdata) + eegspecdB = g.specdata; + freqs = g.freqdata; +else + epochs = round(size(data,2)/frames); + if frames*epochs ~= size(data,2) + error('Spectopo: non-integer number of epochs'); + end + if ~isempty(g.weights) + ncompsori = size(g.weights,1); + if isempty(g.icawinv) + g.icawinv = pinv(g.weights); % maps + end + if ~isempty(g.icacomps) + g.weights = g.weights(g.icacomps, :); + g.icawinv = g.icawinv(:,g.icacomps); + else + g.icacomps = [1:size(g.weights,1)]; + end + end + compeegspecdB = []; + resvar = NaN; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % compute channel spectra using pwelch() + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + epoch_subset = ones(1,epochs); + if g.percent ~= 1 && epochs == 1 + fprintf('Selecting the first %2.1f%% of data for analysis...\n', g.percent*100); + frames = round(size(data,2)*g.percent); + data = data(:, 1:frames); + g.boundaries(find(g.boundaries > frames)) = []; + if ~isempty(g.boundaries) + g.boundaries(end+1) = frames; + end; + end + if g.percent ~= 1 && epochs > 1 + epoch_subset = zeros(1,epochs); + nb = ceil( g.percent*epochs); + while nb>0 + index = ceil(rand*epochs); + if ~epoch_subset(index) + epoch_subset(index) = 1; + nb = nb-1; + end + end; + epoch_subset = find(epoch_subset == 1); + fprintf('Randomly selecting %d of %d data epochs for analysis...\n', length(epoch_subset),epochs); + else + epoch_subset = find(epoch_subset == 1); + end + if isempty(g.weights) + %%%%%%%%%%%%%%%%%%%%%%%%%%% + % compute data spectra + %%%%%%%%%%%%%%%%%%%%%%%%%%% + myfprintf(g.verbose, 'Computing spectra') + [eegspecdB freqs specstd] = spectcomp( data, frames, srate, epoch_subset, g); + if ~isempty(g.mapnorm) % normalize by component map RMS power (if data contain 1 component + myfprintf(g.verbose, 'Scaling spectrum by component RMS of scalp map power\n'); + eegspecdB = sqrt(mean(g.mapnorm.^4)) * eegspecdB; + % the idea is to take the RMS of the component activity (compact) projected at each channel + % spec = sqrt( power(g.mapnorm(1)*compact).^2 + power(g.mapnorm(2)*compact).^2 + ...) + % spec = sqrt( g.mapnorm(1)^4*power(compact).^2 + g.mapnorm(1)^4*power(compact).^2 + ...) + % spec = sqrt( g.mapnorm(1)^4 + g.mapnorm(1)^4 + ... )*power(compact) + end + + tmpc = find(eegspecdB(:,1)); % > 0 power chans + tmpindices = find(eegspecdB(:,1) == 0); + if ~isempty(tmpindices) + zchans = int2str(tmpindices); % 0-power chans + else zchans = []; + end + if length(tmpc) ~= size(eegspecdB,1) + myfprintf(g.verbose, '\nWarning: channels [%s] have 0 values, so will be omitted from the display', ... + zchans); + eegspecdB = eegspecdB(tmpc,:); + if ~isempty(specstd), specstd = specstd(tmpc,:); end + if ~isempty(g.chanlocs) + g.chanlocs2 = g.chanlocs(tmpc); + end + end + eegspecdB = 10*log10(eegspecdB); + specstd = 10*log10(specstd); + myfprintf(g.verbose, '\n'); + else + % compute data spectrum + if isempty(g.plotchan) || g.plotchan == 0 + myfprintf(g.verbose, 'Computing spectra') + [eegspecdB freqs specstd] = spectcomp( data, frames, srate, epoch_subset, g); + myfprintf(g.verbose, '\n'); % log below + else + myfprintf(g.verbose, 'Computing spectra at specified channel') + g.reref = 'off'; + [eegspecdB freqs specstd] = spectcomp( data(g.plotchan,:), frames, srate, epoch_subset, g); + myfprintf(g.verbose, '\n'); % log below + end + g.reref = 'off'; + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % select channels and spectra + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if isempty(g.plotchan) % find channel of minimum power + [tmp indexfreq] = min(abs(g.freq-freqs)); + [tmp g.plotchan] = min(eegspecdB(:,indexfreq)); + myfprintf(g.verbose, 'Channel %d has maximum power at %g\n', g.plotchan, g.freq); + eegspecdBtoplot = eegspecdB(g.plotchan, :); + elseif g.plotchan == 0 + myfprintf(g.verbose, 'Computing RMS power at all channels\n'); + eegspecdBtoplot = sqrt(mean(eegspecdB.^2,1)); % RMS before log as for components + else + eegspecdBtoplot = eegspecdB; + end + tmpc = find(eegspecdB(:,1)); % > 0 power chans + zchans = int2str(find(eegspecdB(:,1) == 0)); % 0-power chans + if length(tmpc) ~= size(eegspecdB,1) + myfprintf(g.verbose, '\nWarning: channels [%s] have 0 values, so will be omitted from the display', ... + zchans); + eegspecdB = eegspecdB(tmpc,:); + if ~isempty(specstd), specstd = specstd(tmpc,:); end + if ~isempty(g.chanlocs) + g.chanlocs2 = g.chanlocs(tmpc); + end + end + specstd = 10*log10(specstd); + eegspecdB = 10*log10(eegspecdB); + eegspecdBtoplot = 10*log10(eegspecdBtoplot); + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + % compute component spectra + %%%%%%%%%%%%%%%%%%%%%%%%%%% + newweights = g.weights; + if strcmp(g.memory, 'high') && strcmp(g.icamode, 'normal') + myfprintf(g.verbose, 'Computing component spectra: ') + [compeegspecdB freqs] = spectcomp( newweights*data(:,:), frames, srate, epoch_subset, g); + else % in case out of memory error, multiply conmponent sequentially + if strcmp(g.icamode, 'sub') && ~isempty(g.plotchan) && g.plotchan == 0 + % scan all electrodes + myfprintf(g.verbose, 'Computing component spectra at each channel: ') + for index = 1:size(data,1) + g.plotchan = index; + [compeegspecdB(:,:,index) freqs] = spectcomp( data, frames, srate, epoch_subset, g, newweights); + end + g.plotchan = 0; + else + myfprintf(g.verbose, 'Computing component spectra: ') + [compeegspecdB freqs] = spectcomp( data, frames, srate, epoch_subset, g, newweights); + end + end + myfprintf(g.verbose, '\n'); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % rescale spectra with respect to projection (rms = root mean square) + % all channels: component_i power = rms(inverseweigths(component_i)^2)*power(activation_component_i); + % one channel: component_i power = inverseweigths(channel_j,component_i)^2)*power(activation_component_i); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if strcmpi(g.icamode, 'normal') + for index = 1:size(compeegspecdB,1) + if g.plotchan == 0 % normalize by component scalp map power + compeegspecdB(index,:) = 10*log10( sqrt(mean(g.icawinv(:,index).^4)) * compeegspecdB(index,:) ); + else + compeegspecdB(index,:) = 10*log10( g.icawinv(g.plotchan,index)^2 * compeegspecdB(index,:) ); + end + end + else % already spectrum of data-components + compeegspecdB = 10*log10( compeegspecdB ); + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + % select components to plot + %%%%%%%%%%%%%%%%%%%%%%%%%%% + if isempty(g.icamaps) + [tmp indexfreq] = min(abs(g.freq-freqs)); + g.icafreqsval = compeegspecdB(:, indexfreq); + [g.icafreqsval g.icamaps] = sort(g.icafreqsval); + if strcmp(g.icamode, 'normal') + g.icamaps = g.icamaps(end:-1:1); + g.icafreqsval = g.icafreqsval(end:-1:1); + end + if g.nicamaps < length(g.icamaps), g.icamaps = g.icamaps(1:g.nicamaps); end + else + [tmp indexfreq] = min(abs(g.freq-freqs)); + g.icafreqsval = compeegspecdB(g.icamaps, indexfreq); + end + end +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% compute axis and caxis g.limits +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if length(g.limits)<1 || isnan(g.limits(1)) + g.limits(1) = LOPLOTHZ; +end + +if ~isempty(g.freq) + if length(g.limits)<2 || isnan(g.limits(2)) + maxheadfreq = max(g.freq); + if rem(maxheadfreq,5) ~= 0 + g.limits(2) = 5*ceil(maxheadfreq/5); + else + g.limits(2) = maxheadfreq*1.1; + end + end + + g.freq = sort(g.freq); % Determine topoplot frequencies + freqidx = zeros(1,length(g.freq)); % Do not interpolate between freqs + for f=1:length(g.freq) + [tmp fi] = min(abs(freqs-g.freq(f))); + freqidx(f)=fi; + end + +else % no freq specified + + if isnan(g.limits(2)) + g.limits(2) = srate/2; + end +end + +[tmp maxfreqidx] = min(abs(g.limits(2)-freqs)); % adjust max frequency +[tmp minfreqidx] = min(abs(g.limits(1)-freqs)); % adjust min frequency + +if length(g.limits)<3|isnan(g.limits(3)) + reallimits(1) = min(min(eegspecdB(:,minfreqidx:maxfreqidx))); +else + reallimits(1) = g.limits(3); +end +if length(g.limits)<4|isnan(g.limits(4)) + reallimits(2) = max(max(eegspecdB(:,minfreqidx:maxfreqidx))); + dBrange = reallimits(2)-reallimits(1); % expand range a bit beyond data g.limits + reallimits(1) = reallimits(1)-dBrange/7; + reallimits(2) = reallimits(2)+dBrange/7; +else + reallimits(2) = g.limits(4); +end + +if length(g.limits)<5 % default caxis plotting g.limits + g.limits(5) = nan; +end +if length(g.limits)<6 + g.limits(6) = nan; +end + +if isnan(g.limits(5))+isnan(g.limits(6)) == 1 + fprintf('spectopo(): limits 5 and 6 must either be given or nan\n'); + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% plot spectrum of each channel +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if strcmpi(g.plot, 'on') + mainfig = gca; axis off; + if ~isempty(g.freq) + specaxes = sbplot(3,4,[5 12], 'ax', mainfig); + end + + if isempty(g.weights) + %pl=plot(freqs(1:maxfreqidx),eegspecdB(:,1:maxfreqidx)'); % old command + if strcmpi(g.plotmean, 'on'), specdata = mean(eegspecdB,1); % average channels + else specdata = eegspecdB; + end + for index = 1:size(specdata,1) % scan channels + tmpcol = allcolors{mod(index, length(allcolors))+1}; + command = [ 'disp(''Channel ' int2str(index) ''')' ]; + pl(index)=plot(freqs(1:maxfreqidx),specdata(index,1:maxfreqidx)', ... + 'color', tmpcol, 'ButtonDownFcn', command); hold on; + end + else + for index = 1:size(eegspecdBtoplot,1) + tmpcol = allcolors{mod(index, length(allcolors))+1}; + command = [ 'disp(''Channel ' int2str(g.plotchan(index)) ''')' ]; + pl(index)=plot(freqs(1:maxfreqidx),eegspecdBtoplot(index,1:maxfreqidx)', ... + 'color', tmpcol, 'ButtonDownFcn', command); hold on; + end + end + set(pl,'LineWidth',2); + set(gca,'TickLength',[0.02 0.02]); + try, + axis([freqs(minfreqidx) freqs(maxfreqidx) reallimits(1) reallimits(2)]); + catch, disp('Could not adjust axis'); end + xl=xlabel('Frequency (Hz)'); + set(xl,'fontsize',AXES_FONTSIZE_L); + % yl=ylabel('Rel. Power (dB)'); + yl=ylabel('Log Power Spectral Density 10*log_{10}(\muV^{2}/Hz)');%yl=ylabel('Power 10*log_{10}(\muV^{2}/Hz)'); + set(yl,'fontsize',AXES_FONTSIZE_L); + set(gca,'fontsize',AXES_FONTSIZE_L) + box off; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% plot component contribution % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +colrs = {'r','b','g','m','c'}; % component spectra trace colors +if ~isempty(g.weights) + if strcmpi(g.plot, 'on') + if strcmpi(g.icamode, 'sub') + set(pl,'LineWidth',5, 'color', 'k'); + else + set(pl, 'linewidth', 2, 'color', 'k'); + end + + hold on; + for f=1:length(g.icamaps) + colr = colrs{mod((f-1),5)+1}; + command = [ 'disp(''Component ' int2str(g.icamaps(f)) ''')' ]; + pl2(index)=plot(freqs(1:maxfreqidx),compeegspecdB(g.icamaps(f),1:maxfreqidx)', ... + 'color', colr, 'ButtonDownFcn', command); hold on; + end + othercomps = setdiff_bc(1:size(compeegspecdB,1), g.icamaps); + if ~isempty(othercomps) + for index = 1:length(othercomps) + tmpcol = allcolors{mod(index, length(allcolors))+1}; + command = [ 'disp(''Component ' int2str(othercomps(index)) ''')' ]; + pl(index)=plot(freqs(1:maxfreqidx),compeegspecdB(othercomps(index),1:maxfreqidx)', ... + 'color', tmpcol, 'ButtonDownFcn', command); hold on; + end + end + if length(g.limits)<3|isnan(g.limits(3)) + newaxis = axis; + newaxis(3) = min(newaxis(3), min(min(compeegspecdB(:,1:maxfreqidx)))); + newaxis(4) = max(newaxis(4), max(max(compeegspecdB(:,1:maxfreqidx)))); + axis(newaxis); + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % indicate component contribution % + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + maxdatadb = max(eegspecdBtoplot(:,freqidx(1))); + [tmp indexfreq] = min(abs(g.freq-freqs)); + for index = 1:length(g.icacomps) + if strcmp(g.icamode, 'normal') + % note: maxdatadb = eegspecdBtoplot (RMS power of data) + resvar(index) = 100*exp(-(maxdatadb-compeegspecdB(index, indexfreq))/10*log(10)); + myfprintf(g.verbose, 'Component %d percent relative variance: %6.2f\n', g.icacomps(index), resvar(index)); + else + if g.plotchan == 0 + resvartmp = []; + for chan = 1:size(eegspecdB,1) % scan channels + resvartmp(chan) = 100 - 100*exp(-(eegspecdB(chan,freqidx(1))-compeegspecdB(index, indexfreq, chan))/10*log(10)); + end + resvar(index) = mean(resvartmp); % mean contribution for all channels + stdvar(index) = std(resvartmp); + myfprintf(g.verbose, 'Component %d percent variance accounted for: %6.2f ± %3.2f\n', ... + g.icacomps(index), resvar(index), stdvar(index)); + else + resvar(index) = 100 - 100*exp(-(maxdatadb-compeegspecdB(index, indexfreq))/10*log(10)); + myfprintf(g.verbose, 'Component %d percent variance accounted for: %6.2f\n', g.icacomps(index), resvar(index)); + end + end + end + + % for icamode=sub and plotchan == 0 -> take the RMS power of all channels + % ----------------------------------------------------------------------- + if ndims(compeegspecdB) == 3 + compeegspecdB = exp( compeegspecdB/10*log(10) ); + compeegspecdB = sqrt(mean(compeegspecdB.^2,3)); % RMS before log (dim1=comps, dim2=freqs, dim3=chans) + compeegspecdB = 10*log10( compeegspecdB ); + end + +end + +if ~isempty(g.freq) && strcmpi(g.plot, 'on') + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % plot vertical lines through channel trace bundle at each headfreq + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if isempty(g.weights) + for f=1:length(g.freq) + hold on; + plot([freqs(freqidx(f)) freqs(freqidx(f))], ... + [min(eegspecdB(:,freqidx(f))) max(eegspecdB(:,freqidx(f)))],... + 'k','LineWidth',2.5); + end + else + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % plot vertical line at comp analysis freq + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + mincompdB = min([min(eegspecdB(:,freqidx(1))) min(compeegspecdB(:,freqidx(1)))]); + maxcompdB = max([max(eegspecdB(:,freqidx(1))) max(compeegspecdB(:,freqidx(1)))]); + hold on; + plot([freqs(freqidx(1)) freqs(freqidx(1))],[mincompdB maxcompdB],'k', 'LineWidth',2.5); + end + + %%%%%%%%%%%%%%%%%%%%%%%%%% + % create axis for topoplot + %%%%%%%%%%%%%%%%%%%%%%%%%% + tmpmainpos = get(gca, 'position'); + headax = zeros(1,length(g.freq)); + for f=1:length(g.freq)+length(g.icamaps) + headax(f) = sbplot(3,length(g.freq)+length(g.icamaps),f, 'ax', mainfig); + axis([-1 1 -1 1]); + + %axis x coords and use + tmppos = get(headax(f), 'position'); + allaxcoords(f) = tmppos(1); + allaxuse(f) = 0; + end + large = sbplot(1,1,1, 'ax', mainfig); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % compute relative positions on plot + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if ~isempty(g.weights) + freqnormpos = tmpmainpos(1) + tmpmainpos(3)*(freqs(freqidx(1))-g.limits(1))/(g.limits(2)-g.limits(1)); + for index = 1:length(g.icamaps)+1 + [realpos(index) allaxuse] = closestplot( freqnormpos, allaxcoords, allaxuse ); + end + + % put the channel plot a little bit higher + tmppos = get(headax(realpos(1)), 'position'); + tmppos(2) = tmppos(2)+0.04; + set(headax(realpos(1)), 'position', tmppos); + else + realpos = 1:length(g.freq); % indices giving order of plotting positions + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % plot connecting lines using changeunits() + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + for f=1:length(g.freq)+length(g.icamaps) + if ~isempty(g.weights) + if f>length(g.freq) % special case of ICA components + from = changeunits([freqs(freqidx(1)),g.icafreqsval(f-1)],specaxes,large); + %g.icafreqsval contains the sorted frequency values at the specified frequency + else + from = changeunits([freqs(freqidx(f)),maxcompdB],specaxes,large); + end + else + from = changeunits([freqs(freqidx(f)),max(eegspecdB(:,freqidx(f)))],specaxes,large); + end + pos = get(headax(realpos(f)),'position'); + to = changeunits([0,0],headax(realpos(f)),large)+[0 -min(pos(3:4))/2.5]; + hold on; + if f<=length(g.freq) + colr = 'k'; + else + colr = colrs{mod((f-2),5)+1}; + end + li(realpos(f)) = plot([from(1) to(1)],[from(2) to(2)],colr,'LineWidth',PLOT_LINEWIDTH_S); + axis([0 1 0 1]); + axis off; + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % plot selected channel head using topoplot() + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + myfprintf(g.verbose, 'Plotting scalp distributions: ') + for f=1:length(g.freq) + axes(headax(realpos(f))); + + topodata = eegspecdB(:,freqidx(f))-nan_mean(eegspecdB(:,freqidx(f))); + + if isnan(g.limits(5)), + maplimits = 'absmax'; + else + maplimits = [g.limits(5) g.limits(6)]; + end + + % + % If 1 channel in g.plotchan + % + + if ~isempty(g.plotchan) && g.plotchan ~= 0 + % if ~isempty(varargin) % if there are extra topoplot() flags + % topoplot(g.plotchan,g.chanlocs,'electrodes','off', ... + % 'style', 'blank', 'emarkersize1chan', 10, varargin{:}); + % else + topoplot(g.plotchan,g.chanlocs,'electrodes','off', ... + 'style', 'blank', 'emarkersize1chan', 10); + % end + if isstruct(g.chanlocs) + tl=title(g.chanlocs(g.plotchan).labels); + else + tl=title([ 'c' int2str(g.plotchan)]); + end + + else % plot all channels in g.plotchans + + if isempty(g.mapframes) || g.mapframes(1) == 0 + g.mapframes = 1:size(eegspecdB,1); % default to plotting all chans + end + if ~isempty(varargin) + topoplot(topodata(g.mapframes),g.chanlocs2,'maplimits',maplimits, varargin{:}); + else + topoplot(topodata(g.mapframes),g.chanlocs2,'maplimits',maplimits); + end + if f= 3 + tmp = compeegspecdB; + compeegspecdB = zeros(ncompsori, size(tmp,2)); + compeegspecdB(g.icacomps,:) = tmp; +end + +%%%%%%%%%%%%%%%% +% Turn on axcopy (disabled to allow to click on curves) +%%%%%%%%%%%%%%%% +if strcmpi(g.plot, 'on') + disp('Click on each trace for channel/component index'); + axcopy(gcf, 'if ~isempty(get(gca, ''''userdata'''')), eval(get(gca, ''''userdata'''')); end;'); + % will not erase the commands for the curves +end + +%%%%%%%%%%%%%%%% +% Plot color bar +%%%%%%%%%%%%%%%% +function plotcolbar(g) + cb=cbar; + pos = get(cb,'position'); + set(cb,'position',[pos(1) pos(2) 0.03 pos(4)]); + set(cb,'fontsize',12); + try + if isnan(g.limits(5)) + ticks = get(cb,'ytick'); + [tmp zi] = find(ticks == 0); + ticks = [ticks(1) ticks(zi) ticks(end)]; + set(cb,'ytick',ticks); + set(cb,'yticklabel',strvcat('-',' ','+')); + end + catch, end; % in a single channel is plotted + +%%%%%%%%%%%%%%%%%%%%%%% +% function closest plot +%%%%%%%%%%%%%%%%%%%%%%% +function [index, usedplots] = closestplot(xpos, xcentercoords, usedplots); + notused = find(usedplots == 0); + xcentercoords = xcentercoords(notused); + [tmp index] = min(abs(xcentercoords-xpos)); + index = notused(index); + usedplots(index) = 1; + return; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% function computing spectrum +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function [eegspecdB, freqs, specstd] = spectcomp( data, frames, srate, epoch_subset, g, newweights); + usepwelch = license('checkout','Signal_Toolbox') && exist('pwelch'); + + if exist('newweights') == 1 + nchans = size(newweights,1); + else + nchans = size(data,1); + end + %fftlength = 2^round(log(srate)/log(2))*g.freqfac; + if isempty(g.winsize) +% winlength = max(pow2(nextpow2(frames)-3),4); %*2 since diveded by 2 later +% winlength = min(winlength, 512); +% winlength = max(winlength, 256); + winlength = min(round(srate), frames); + else + winlength = g.winsize; + end + + if strcmpi(g.wintype,'blackmanharris') + if usepwelch + winlength = blackmanharris(round(winlength/g.blckhn)); + else + g.wintype = 'hamming'; + fprintf('\nSignal processing toolbox (SPT) absent: unable to use Blackman-Harris window\n'); + fprintf(' using pwelch function from Octave\n'); + end + end + + if isempty(g.nfft) && strcmp(g.wintype,'hamming') + %fftlength = 2^(nextpow2(winlength))*g.freqfac; + fftlength = winlength*g.freqfac; + elseif ~isempty(g.nfft) && strcmp(g.wintype,'hamming') + fftlength = g.nfft; + elseif strcmp(g.wintype,'blackmanharris') + fftlength = 2^(nextpow2(length(winlength)))*g.freqfac; + end + + if ~usepwelch + myfprintf(g.verbose, '\nSignal processing toolbox (SPT) absent: spectrum computed using the pwelch()\n'); + myfprintf(g.verbose, 'function from Octave which is supposedly 100%% compatible with the Matlab pwelch function\n'); + end + myfprintf(g.verbose,' (window length %d; fft length: %d; overlap %d):\n', winlength, fftlength, g.overlap); + + for c=1:nchans % scan channels or components + if exist('newweights') == 1 + if strcmp(g.icamode, 'normal') + tmpdata = newweights(c,:)*data; % component activity + else % data - component contribution + tmpdata = data(g.plotchan,:) - (g.icawinv(g.plotchan,c)*newweights(c,:))*data; + end + else + tmpdata = data(c,:); % channel activity + end + if strcmp(g.reref, 'averef') + tmpdata = averef(tmpdata); + end + for e=epoch_subset + if isempty(g.boundaries) + if usepwelch + [tmpspec,freqs] = pwelch(matsel(tmpdata,frames,0,1,e),... + winlength,g.overlap,fftlength,srate); + else + [tmpspec,freqs] = spec(matsel(tmpdata,frames,0,1,e),fftlength,srate,... + winlength,g.overlap); + end + %[tmpspec,freqs] = psd(matsel(tmpdata,frames,0,1,e),fftlength,srate,... + % winlength,g.overlap); + if c==1 && e==epoch_subset(1) + eegspec = zeros(nchans,length(freqs)); + specstd = zeros(nchans,length(freqs)); + end + eegspec(c,:) = eegspec(c,:) + tmpspec'; + specstd(c,:) = specstd(c,:) + tmpspec'.^2; + else + g.boundaries = round(g.boundaries); + for n=1:length(g.boundaries)-1 + if g.boundaries(n+1) - g.boundaries(n) >= winlength % ignore segments of less than winlength + if usepwelch + [tmpspec,freqs] = pwelch(tmpdata(e,g.boundaries(n)+1:g.boundaries(n+1)),... + winlength,g.overlap,fftlength,srate); + else + [tmpspec,freqs] = spec(tmpdata(e,g.boundaries(n)+1:g.boundaries(n+1)),... + fftlength,srate,winlength,g.overlap); + end + if exist('eegspec') ~= 1 + eegspec = zeros(nchans,length(freqs)); + specstd = zeros(nchans,length(freqs)); + end + eegspec(c,:) = eegspec(c,:) + tmpspec'* ... + ((g.boundaries(n+1)-g.boundaries(n)+1)/g.boundaries(end)); + specstd(c,:) = eegspec(c,:) + tmpspec'.^2 * ... + ((g.boundaries(n+1)-g.boundaries(n)+1)/g.boundaries(end)); + end + end + end + end + myfprintf(g.verbose,'.'); + end + + n = length(epoch_subset); + eegspecdB = eegspec/n; % normalize by the number of sections + if n>1 % normalize standard deviation by the number of sections + specstd = sqrt( (specstd + eegspec.^2/n)/(n-1) ); + else specstd = []; + end + return; + + function myfprintf(verbose, varargin) + if strcmpi(verbose, 'on') + fprintf(varargin{:}); + end +