ic_distortion() - Quantifies how much ICA distorts artifact free ERPs,
when it is used to remove EEG artifacts (e.g., blinks,
muscle activity). Relies on the field EEG.iclabels to
identify ICs corresponding to EEG artifacts. See is_art.m
and ic_prop.m for more information about these labels.
ic_prop.m is an easy way to label ICs.
Usage:
>> [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, drty_erp_ct]=ic_distortion(EEG,varargin);
Required Input:
EEG - EEGLAB "EEG" struct variable (stored as .set files by EEGLAB)
Optional Inputs:
'bins' - [integer vector] Bins to extract trials from. [] means to
use all bins {default: use all bins}
'bsln' - [start_time stop_time] or 'allpre': The time window whose
mean will be removed from each epoch to baseline the data.
Times are in units of milliseconds (e.g., [-100 0]). If
'allpre' all time points up to and including time 0 are used.
{default: 'allpre'}
'time_wind' - [start_time stop_time] The time window in which ERP
distortion is measured. Times are in units of milliseconds
(e.g., [250 550]). {default: all time points}
'plot' - [1 | 0] If 1, the a figure illustrating the topography of
ERP distortion will be produced. Useful for seeing which
electrodes are most distorted. {default: 0}
'fig_id' - [integer] ID number of Matlab figure window in which to create
the distortion topography. If an ID is not specified, the
plot will be generated in a new figure window. {default:
new figure window}
'verblevel' - an integer specifiying the amount of information you want
functions to provide about what they are doing during runtime.
Options are:
0 - quiet, only show errors, warnings, and EEGLAB reports
1 - stuff anyone should probably know
2 - stuff you should know the first time you start working
with a data set {default value if not already globally
specified}
3 - stuff that might help you debug (show all reports)
Outputs:
dstrt_ratio_per_chan - A unitless measure of how much ICA-artifact
correction is distorting the clean ERPs at each
electrode. This derived by (1) computing the
ERPs (cln_erps) using just artifact-free epochs,
(2) filtering the clean ERPs with the ICA artifact
correction filter to produce cln_erpsICA, (3)
taking the root mean squared difference between
cln_erps and cln_erpsICA at each channel (rmsd),
(4) dividing rmsd by the root mean squared amplitude
of the clean erps, rmsa. A dstrt_ratio_per_chan
value of 0 means no distortion. A dstrt_ratio_per_chan
value of .5 means that the magnitude of the
distortion is about half the magnitude of the
ERP at that channel. A dstrt_ratio_per_chan
value of 1 means that the magnitude of
distortion is about equal to that of the ERP
at that channel. Conceptually, this is analagous
to the "rmsen" blinkcp diagnostic, but the two
values cannot be directly compared.
mn_dstrt_ratio - The mean of dstrt_ratio_per_chan (i.e., the mean
distortion across all electrodes).
cln_erp_ct - The number of artifact free trials/epochs
included in the selected bins.
drty_erp_ct - The number of polluted trials/epochs
included in the selected bins.
Global Variables:
VERBLEVEL = level of verbosity (i.e., tells functions how much
how much to report about what they're doing during
runtime) set by the optional function argument 'verblevel'
Example:
>> load('/homes/dgroppe/SANDBOX/VISODBL/yngvob05_cmprbcp.set','-MAT');
>> [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, ...
drty_erp_ct]=ic_distortion(EEG,'bins',[7 8],'plot',1,'bsln',[-100 0]);
Notes:
-The ith epoch is labelled as polluted or clean via the field
EEG.epoch(i).eventlogflag. A non-zero value means that the epoch is
polluted. These are automatically set from a Kutaslab arf file when the
data are imported into MATLAB via crw2set.m.
-The accuracy of this function's estimates of ICA distortion depends upon
the accuracy with which the epochs have been labelled as "clean" or
"polluted". If polluted epochs are mistakenly labelled as "clean"
ic_distortion.m will OVERestimate the amount of distortion. The function
ic_prop.m can be used to visualize single epochs, sorted by artifact
status (i.e., "clean" or "polluted"), and can alert you to inaccurately
labelled epohcs.
-This function cannot be used to estimate the amount distortion caused by
ICA artifact correction of EEG artifacts that appear in all or almost all
epochs (e.g., heartbeat), since you can't get a good estimate of clean
ERPs. If you correct for these artifacts with ICA, ic_distortion will
OVERestimate the amount of distortion.
Author:
David Groppe
Kutaslab, 12/20090001 function [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, drty_erp_ct]=ic_distortion(EEG,varargin) 0002 % ic_distortion() - Quantifies how much ICA distorts artifact free ERPs, 0003 % when it is used to remove EEG artifacts (e.g., blinks, 0004 % muscle activity). Relies on the field EEG.iclabels to 0005 % identify ICs corresponding to EEG artifacts. See is_art.m 0006 % and ic_prop.m for more information about these labels. 0007 % ic_prop.m is an easy way to label ICs. 0008 % 0009 % Usage: 0010 % >> [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, drty_erp_ct]=ic_distortion(EEG,varargin); 0011 % 0012 % 0013 % Required Input: 0014 % EEG - EEGLAB "EEG" struct variable (stored as .set files by EEGLAB) 0015 % 0016 % Optional Inputs: 0017 % 'bins' - [integer vector] Bins to extract trials from. [] means to 0018 % use all bins {default: use all bins} 0019 % 'bsln' - [start_time stop_time] or 'allpre': The time window whose 0020 % mean will be removed from each epoch to baseline the data. 0021 % Times are in units of milliseconds (e.g., [-100 0]). If 0022 % 'allpre' all time points up to and including time 0 are used. 0023 % {default: 'allpre'} 0024 % 'time_wind' - [start_time stop_time] The time window in which ERP 0025 % distortion is measured. Times are in units of milliseconds 0026 % (e.g., [250 550]). {default: all time points} 0027 % 'plot' - [1 | 0] If 1, the a figure illustrating the topography of 0028 % ERP distortion will be produced. Useful for seeing which 0029 % electrodes are most distorted. {default: 0} 0030 % 'fig_id' - [integer] ID number of Matlab figure window in which to create 0031 % the distortion topography. If an ID is not specified, the 0032 % plot will be generated in a new figure window. {default: 0033 % new figure window} 0034 % 'verblevel' - an integer specifiying the amount of information you want 0035 % functions to provide about what they are doing during runtime. 0036 % Options are: 0037 % 0 - quiet, only show errors, warnings, and EEGLAB reports 0038 % 1 - stuff anyone should probably know 0039 % 2 - stuff you should know the first time you start working 0040 % with a data set {default value if not already globally 0041 % specified} 0042 % 3 - stuff that might help you debug (show all reports) 0043 % 0044 % 0045 % Outputs: 0046 % dstrt_ratio_per_chan - A unitless measure of how much ICA-artifact 0047 % correction is distorting the clean ERPs at each 0048 % electrode. This derived by (1) computing the 0049 % ERPs (cln_erps) using just artifact-free epochs, 0050 % (2) filtering the clean ERPs with the ICA artifact 0051 % correction filter to produce cln_erpsICA, (3) 0052 % taking the root mean squared difference between 0053 % cln_erps and cln_erpsICA at each channel (rmsd), 0054 % (4) dividing rmsd by the root mean squared amplitude 0055 % of the clean erps, rmsa. A dstrt_ratio_per_chan 0056 % value of 0 means no distortion. A dstrt_ratio_per_chan 0057 % value of .5 means that the magnitude of the 0058 % distortion is about half the magnitude of the 0059 % ERP at that channel. A dstrt_ratio_per_chan 0060 % value of 1 means that the magnitude of 0061 % distortion is about equal to that of the ERP 0062 % at that channel. Conceptually, this is analagous 0063 % to the "rmsen" blinkcp diagnostic, but the two 0064 % values cannot be directly compared. 0065 % mn_dstrt_ratio - The mean of dstrt_ratio_per_chan (i.e., the mean 0066 % distortion across all electrodes). 0067 % cln_erp_ct - The number of artifact free trials/epochs 0068 % included in the selected bins. 0069 % drty_erp_ct - The number of polluted trials/epochs 0070 % included in the selected bins. 0071 % 0072 % Global Variables: 0073 % VERBLEVEL = level of verbosity (i.e., tells functions how much 0074 % how much to report about what they're doing during 0075 % runtime) set by the optional function argument 'verblevel' 0076 % 0077 % 0078 % Example: 0079 % >> load('/homes/dgroppe/SANDBOX/VISODBL/yngvob05_cmprbcp.set','-MAT'); 0080 % >> [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, ... 0081 % drty_erp_ct]=ic_distortion(EEG,'bins',[7 8],'plot',1,'bsln',[-100 0]); 0082 % 0083 % 0084 % Notes: 0085 % -The ith epoch is labelled as polluted or clean via the field 0086 % EEG.epoch(i).eventlogflag. A non-zero value means that the epoch is 0087 % polluted. These are automatically set from a Kutaslab arf file when the 0088 % data are imported into MATLAB via crw2set.m. 0089 % 0090 % -The accuracy of this function's estimates of ICA distortion depends upon 0091 % the accuracy with which the epochs have been labelled as "clean" or 0092 % "polluted". If polluted epochs are mistakenly labelled as "clean" 0093 % ic_distortion.m will OVERestimate the amount of distortion. The function 0094 % ic_prop.m can be used to visualize single epochs, sorted by artifact 0095 % status (i.e., "clean" or "polluted"), and can alert you to inaccurately 0096 % labelled epohcs. 0097 % 0098 % -This function cannot be used to estimate the amount distortion caused by 0099 % ICA artifact correction of EEG artifacts that appear in all or almost all 0100 % epochs (e.g., heartbeat), since you can't get a good estimate of clean 0101 % ERPs. If you correct for these artifacts with ICA, ic_distortion will 0102 % OVERestimate the amount of distortion. 0103 % 0104 % Author: 0105 % David Groppe 0106 % Kutaslab, 12/2009 0107 % 0108 0109 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0110 % 0111 0112 %%%%%%%%%%%%%%%% POSSIBLE FUTURE DEVELOPMENT %%%%%%%%%%%%%%%% 0113 % 0114 %-Add typical values of mn_dstrt_ratio 0115 % 0116 0117 0118 0119 %Input Parser 0120 p=inputParser; 0121 %Note: the order of the required arguments needs to match precisely their 0122 %order in the function definition (which is also the order used by p.parse 0123 %below) 0124 p.addRequired('EEG',@isstruct); 0125 p.addParamValue('bins',[],@isnumeric); %default is to use all bins 0126 p.addParamValue('bsln','allpre',@(x) isnumeric(x) || strcmpi(x,'allpre')); %default is to use all prestim time points 0127 p.addParamValue('time_wind',[],@isnumeric); %default is to use all time points 0128 p.addParamValue('plot',0,@isnumeric); %default is to use all time points 0129 p.addParamValue('verblevel',[],@(x) x>=0); 0130 p.addParamValue('fig_id',[],@(x) x>=0); 0131 0132 p.parse(EEG,varargin{:}); 0133 0134 0135 global VERBLEVEL 0136 0137 if isempty(VERBLEVEL) && isempty(p.Results.verblevel) 0138 VERBLEVEL=2; 0139 elseif ~isempty(p.Results.verblevel) 0140 VERBLEVEL=p.Results.verblevel; 0141 end 0142 0143 0144 if isempty(p.Results.bins), 0145 n_bin=length(EEG.bindesc); 0146 bins=1:n_bin; 0147 else 0148 bins=p.Results.bins; 0149 n_bin=length(bins); 0150 end 0151 0152 s=size(EEG.data); 0153 n_chan=s(1); 0154 n_tpt=s(2); 0155 n_epoch=s(3); 0156 0157 0158 %% Find time points corresponding to baseline time window 0159 if strcmpi(p.Results.bsln,'allpre'), 0160 %use all prestimulus time points 0161 bsln_tpt(1)=1; 0162 bsln_tpt(2)=find_tpt(0,EEG.times); 0163 else 0164 bsln_tpt(1)=find_tpt(p.Results.bsln(1),EEG.times); 0165 bsln_tpt(2)=find_tpt(p.Results.bsln(2),EEG.times); 0166 end 0167 bsln_win=bsln_tpt(1):bsln_tpt(2); 0168 if isempty(bsln_win), 0169 error('Baseline window from %.1f to %.1f ms includes no time points no time points in current data set. Please revise it.', ... 0170 bsln_tpt(1),bsln_tpt(2)); 0171 else 0172 VerbReport(sprintf('Removing mean EEG data from %.1f ms to %.1f ms for each epoch/trial.\n', ... 0173 EEG.times(bsln_tpt(1)),EEG.times(bsln_tpt(2))),2,VERBLEVEL); 0174 end 0175 0176 0177 %% Find time points corresponding to baseline time window 0178 if isempty(p.Results.time_wind), 0179 use_tpts=1:n_tpt; 0180 else 0181 strt_tpt=find_tpt(p.Results.time_wind(1),EEG.times); 0182 end_tpt=find_tpt(p.Results.time_wind(2),EEG.times); 0183 use_tpts=strt_tpt:end_tpt; 0184 if isempty(use_tpts), 0185 error('Analysis time window from %.1f to %.1f ms includes no time points in current data set. Please revise it.', ... 0186 p.Results.time_wind(1),p.Results.time_wind(2)); 0187 else 0188 VerbReport(sprintf('Using ERPs from %.1f ms to %.1f ms to quantify ICA distortion.\n', ... 0189 EEG.times(use_tpts(1)),EEG.times(use_tpts(end))),1,VERBLEVEL); 0190 end 0191 end 0192 0193 0194 %% Check to see if some independent components have been labeled as artifacts 0195 fldnames=fieldnames(EEG); 0196 ic_arts=0; 0197 for a=1:length(fldnames), 0198 if strcmpi(fldnames{a},'iclabels'), 0199 ic_arts=1; 0200 end 0201 end 0202 0203 if ic_arts, 0204 n=length(EEG.iclabels); 0205 art=[]; 0206 for dg=1:n, 0207 if is_art(EEG.iclabels{dg}), 0208 art=[art dg]; 0209 end 0210 end 0211 non_art=setdiff(1:n_chan,art); 0212 0213 unmix=EEG.icaweights*EEG.icasphere; 0214 mix=inv(unmix); %just in case unmixing matrix has been scaled differently than mixing matrix 0215 ic_filt=mix(:,non_art)*unmix(non_art,:); 0216 VerbReport(sprintf('Removing the contribution of %d ICs that have been labeled as EEG artifacts.',length(art)), ... 0217 1,VERBLEVEL); 0218 VerbReport('The ICs labeled by artifacts are:',1,VERBLEVEL); 0219 for dg=art, 0220 VerbReport(sprintf('IC %d: %s',dg,EEG.iclabels{dg}),1,VERBLEVEL); 0221 end 0222 eeglab_rej=find(EEG.reject.gcompreject==1); 0223 dif1=setdiff(eeglab_rej,art); 0224 dif2=setdiff(art,eeglab_rej); 0225 if ~isempty(dif1) || ~isempty(dif2) 0226 if isempty(eeglab_rej), 0227 rej_ic_str='None'; 0228 else 0229 rej_ic_str=int2str(eeglab_rej); 0230 end 0231 msg=['The ICs marked for rejection by EEGLAB in EEG.reject.gcompreject differ from those labeled as artifacts in EEG.iclabels.' ...); 0232 10 'Only the artifact labels in EEG.iclabels will be used.' 10 ... 0233 'The ICs labeled by artifacts by EEG.reject.gcompreject are: ' rej_ic_str]; 0234 watchit(msg); 0235 end 0236 else 0237 error('No ICs have been labeled as artifacts. Use ic_prop.m to label ICs.'); 0238 end 0239 0240 %% Make erps from clean VS dirty trials 0241 cln_erps=zeros(n_chan,n_tpt); 0242 cln_erpsICA=zeros(n_chan,n_tpt); 0243 cln_erp_ct=0; 0244 drty_erp_ct=0; 0245 VerbReport('Filtering data with ICA epoch by epoch.',1,VERBLEVEL); 0246 for trial=1:n_epoch, 0247 if ~rem(trial,100), 0248 VerbReport(sprintf('Now processing epoch: %d',trial),1,VERBLEVEL); 0249 % if VERBLEVEL>=1, 0250 % fprintf('Now processing epoch: %d\n',trial); %I don't use 0251 % VerbReport here to prevent extra line breaks 0252 % end 0253 end 0254 for types=1:length(EEG.epoch(trial).eventtype), 0255 if strcmpi(EEG.epoch(trial).eventtype{types}(1:3),'bin'), 0256 bin_num=str2double(EEG.epoch(trial).eventtype{types}(4:end)); 0257 bin_id=find(bin_num==bins); 0258 if ~isempty(bin_id), 0259 if ~EEG.epoch(trial).eventlogflag, 0260 trial_data=squeeze(EEG.data(:,:,trial)); %one trial of data, artifact ICs removed 0261 trial_data=rmbase(trial_data,n_tpt,bsln_tpt(1):bsln_tpt(2)); % re-baseline, just in case 0262 cln_erps=cln_erps+trial_data; 0263 cln_erp_ct=cln_erp_ct+1; 0264 0265 %filter clean trial erps with ICA 0266 trial_data=ic_filt*trial_data; %one trial of data, artifact ICs removed 0267 trial_data=rmbase(trial_data,n_tpt,bsln_tpt(1):bsln_tpt(2)); % re-baseline, just in case 0268 cln_erpsICA=cln_erpsICA+trial_data; 0269 else 0270 drty_erp_ct=drty_erp_ct+1; 0271 end 0272 break; 0273 end 0274 end 0275 end 0276 end 0277 % Error check for bins that catch zero trials 0278 if ~cln_erp_ct, 0279 error('No clean trials fall into bin(s) %s.',int2str(bins)); 0280 end 0281 0282 %Turn sums into means 0283 cln_erps=cln_erps/cln_erp_ct; 0284 cln_erpsICA=cln_erpsICA/cln_erp_ct; 0285 0286 0287 % Measure distortion of clean ERP with rms difference at each channel 0288 rms_dstrt=squeeze(sqrt(mean((cln_erps(:,use_tpts)-cln_erpsICA(:,use_tpts)).^2,2))); 0289 rms_cln=squeeze(sqrt(mean(cln_erps(:,use_tpts).^2,2))); 0290 0291 dstrt_ratio_per_chan=rms_dstrt./rms_cln; 0292 mn_dstrt_ratio=mean(dstrt_ratio_per_chan); 0293 0294 VerbReport(sprintf('Mean distortion ratio per electrode is %.3f.\n',mn_dstrt_ratio), ... 0295 1,VERBLEVEL); 0296 0297 if p.Results.plot 0298 if isempty(p.Results.fig_id) 0299 hf=figure; 0300 else 0301 hf=figure(p.Results.fig_id); clf; 0302 end 0303 topoplot(dstrt_ratio_per_chan,EEG.chanlocs,'maplimits',[0 max(dstrt_ratio_per_chan)]); 0304 if n_bin==1, 0305 ht=title(EEG.bindesc{bins}); 0306 elseif n_bin==length(EEG.bindesc) 0307 ht=title('Bins: All'); 0308 else 0309 ht=title(['Bins: ' int2str(bins)]); 0310 end 0311 set(ht,'fontsize',14,'fontweight','bold'); 0312 hc=colorbar; 0313 set(get(hc,'YLabel'),'String','Distortion Ratio','fontsize',14,'fontweight','bold'); 0314 set(hf,'name',[EEG.setname ', ICA distortion']); 0315 end 0316 0317 0318