Home > matlabmk > ic_distortion.m

ic_distortion

PURPOSE ^

ic_distortion() - Quantifies how much ICA distorts artifact free ERPs,

SYNOPSIS ^

function [dstrt_ratio_per_chan, mn_dstrt_ratio, cln_erp_ct, drty_erp_ct]=ic_distortion(EEG,varargin)

DESCRIPTION ^

 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/2009

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 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

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005