Home > matlabmk > set2avg.m

set2avg

PURPOSE ^

set2avg() - Computes ERPs from an EEGLAB set file or EEG variable and writes

SYNOPSIS ^

function set2avg(EEG_or_setname,out_fname,hdr_fname,varargin)

DESCRIPTION ^

 set2avg() - Computes ERPs from an EEGLAB set file or EEG variable and writes
             them to a Kutaslab average file.  ERPs will be normalized
             via the cal pulse ERP and can be artifact corrected via ICA 
             (if artifact ICs have been labeled).  Thus the average file 
             produced is analogous to a Kutaslab ".nrm" file.
              
 Usage:
  >> set2avg(EEG_or_setname,out_fname,hdr_fname,varargin)

 Required Inputs:
  EEG_or_setname - [EEGLAB EEG variable | EEGLAB set filename] Identifies
                   the EEGLAB EEG variable in memory or the set file on disk
                   from which the ERPs will be generated.  If specifying a
                   filename, be sure to include the file's path unless the
                   file is in the current working directory.  IC's
                   corresponding to EEG artifacts should be labelled as
                   such in the field EEG.iclabels.  See the function
                   is_art.m for a list of IC labels recognized as artifacts. 
  out_fname      - [string] The root of the name of the file(s) that will be 
                   created (e.g., 'sbjct1').  The file extension ".nrm",
                   ".icnrm", ".iccav" will be automatically added. ".nrm"
                   is used for ERPs derived from only artifact free trials
                   (no ICA artifact correction).  ".icnrm" is the same as
                   ".nrm" but with ICA artifact correction (useful for
                   seeing how ICA might be distorting ERPs).  ".iccor" is
                   used for ERPs derived from all trials with ICA artifact
                   correction.
  hdr_fname      - [string] The name of an existing Kutaslab average file
                   (e.g., sbjct2.nrm) or a header file produced by the program 
                   mdh.  The header from this file will be used as a
                   template for the new file that is being created.  Any
                   average or header file should work (i.e., it does not
                   need to have the same channels or bins as the set
                   file).  The template file cannot have the same name as
                   one of the average files being produced by set2avg.

 Optional Inputs:
 'bsln'          - [start_time stop_time | 'allpre'] Time window (in milliseconds) 
                   that will be used to baseline each trial.  In other words, for 
                   a given trial and electrode the mean voltage in that time window 
                   will be removed from each time point in the trial. [] means no 
                   baselining (data will be averaged as it is in the set file after ICA
                   artifact correction).  'allpre' means that all pre-event time 
                   points (including Time 0) will be used for baselining. Note, if 
                   you perform ICA artifact correction, you'll want to baseline the
                   data here as any previous baselining might no longer hold. 
                   {Default: []} 
 'pp10uv'        - [scalar] Data resolution ("points per microvolt") at which the
                   averages will be produced. This parameter determines how the
                   data are scaled before writing to the average file.  Ideally,
                   you want to make the data as big as possible, while remaining
                   in the limits of numeric representation range.  Kutaslab default
                   for this is 1000. {Default: 1000}
 'ica_correct'   - [1 | 0] If 1, ICA will be used to correct for EEG artifacts
                   and the average file produced will have a ".iccor" 
                   extension. ICs are identified as artifacts via the field 
                   EEG.iclabels.  See is_art.m for a list of labels recognized 
                   as artifacts. If 0, ICA is not used to correct for artifacts
                   and the average file produced will have a ".nrm" extension.
                   {Default: 0}
 'just_cor_drty' - [1 | 0] If 1, ICA artifact correction will only be
                   applied to polluted trials and ERPs will be formed from both
                   artifact-free and corrected trials.  This argument has no 
                   effect if the argument 'ica_correct' is set to 0. {Default: 0}
 'cln_avgs'      - [1 | 0] If 1, two additional average files will be produced
                   with the same root filename as out_fname. The first of
                   these will have the extension ".nrm" and will be formed  
                   only from artifact free trials (artifact free according 
                   to the arf file used to import the data into EEGLAB). 
                   The second of these will have the extension ".icnrm" 
                   and will be composed of the exact same ERPs in the ".nrm" file 
                   after the ICA artifact correction filter has been applied.
                   Contrasting these two average files is useful for estimating 
                   how ICA artifact correction might be distorting the data.
                   {Default: 0}
 'overwrite'     - [1 | 0] If 1, existing average files that have the same 
                   name as the new files being produced will be automatically 
                   overwritten.  If 0, the user will be asked if s/he wants  
                   to overwrite them. {Default: 0}
 '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)
  query_ccode    - [1 | 0] If one, the user will be asked at the
                   command line for the condition code of bins with zero
                   trials (if any).  If zero, the function will simply
                   use the condition code for the previous bin. {default: 0}


 Outputs:
  An average file with a name specified by 'out_fname' is written to disk.
  Nothing is output to MATLAB.  If 'cln_avgs' option is used, two
  additional average files will be produced.


 Example: 
 >> set_fname='/homes/dgroppe/SANDBOX/PRSENT/prsent41.set';
 >> hdr_fname='/homes/dgroppe/SANDBOX/PRSENT/prsent.hdr';
 >> set2avg(set_fname,'prsent41',hdr_fname,'bsln',[-100 0],'cln_avgs',1,'ica_correct',1);

 Author: 
 David Groppe
 Kutaslab, 11/2009 

 Notes: 
 -This function relies on erpio to create the avg files.  Thus it only
 works on the Kutaslab network (unless you have a compiled copy of erpio
 on whatever machine you're using).

 -Code was written to minimize memory demands at the cost of being slower.
 If faster code is needed, the function could be revised.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function set2avg(EEG_or_setname,out_fname,hdr_fname,varargin)
0002 % set2avg() - Computes ERPs from an EEGLAB set file or EEG variable and writes
0003 %             them to a Kutaslab average file.  ERPs will be normalized
0004 %             via the cal pulse ERP and can be artifact corrected via ICA
0005 %             (if artifact ICs have been labeled).  Thus the average file
0006 %             produced is analogous to a Kutaslab ".nrm" file.
0007 %
0008 % Usage:
0009 %  >> set2avg(EEG_or_setname,out_fname,hdr_fname,varargin)
0010 %
0011 % Required Inputs:
0012 %  EEG_or_setname - [EEGLAB EEG variable | EEGLAB set filename] Identifies
0013 %                   the EEGLAB EEG variable in memory or the set file on disk
0014 %                   from which the ERPs will be generated.  If specifying a
0015 %                   filename, be sure to include the file's path unless the
0016 %                   file is in the current working directory.  IC's
0017 %                   corresponding to EEG artifacts should be labelled as
0018 %                   such in the field EEG.iclabels.  See the function
0019 %                   is_art.m for a list of IC labels recognized as artifacts.
0020 %  out_fname      - [string] The root of the name of the file(s) that will be
0021 %                   created (e.g., 'sbjct1').  The file extension ".nrm",
0022 %                   ".icnrm", ".iccav" will be automatically added. ".nrm"
0023 %                   is used for ERPs derived from only artifact free trials
0024 %                   (no ICA artifact correction).  ".icnrm" is the same as
0025 %                   ".nrm" but with ICA artifact correction (useful for
0026 %                   seeing how ICA might be distorting ERPs).  ".iccor" is
0027 %                   used for ERPs derived from all trials with ICA artifact
0028 %                   correction.
0029 %  hdr_fname      - [string] The name of an existing Kutaslab average file
0030 %                   (e.g., sbjct2.nrm) or a header file produced by the program
0031 %                   mdh.  The header from this file will be used as a
0032 %                   template for the new file that is being created.  Any
0033 %                   average or header file should work (i.e., it does not
0034 %                   need to have the same channels or bins as the set
0035 %                   file).  The template file cannot have the same name as
0036 %                   one of the average files being produced by set2avg.
0037 %
0038 % Optional Inputs:
0039 % 'bsln'          - [start_time stop_time | 'allpre'] Time window (in milliseconds)
0040 %                   that will be used to baseline each trial.  In other words, for
0041 %                   a given trial and electrode the mean voltage in that time window
0042 %                   will be removed from each time point in the trial. [] means no
0043 %                   baselining (data will be averaged as it is in the set file after ICA
0044 %                   artifact correction).  'allpre' means that all pre-event time
0045 %                   points (including Time 0) will be used for baselining. Note, if
0046 %                   you perform ICA artifact correction, you'll want to baseline the
0047 %                   data here as any previous baselining might no longer hold.
0048 %                   {Default: []}
0049 % 'pp10uv'        - [scalar] Data resolution ("points per microvolt") at which the
0050 %                   averages will be produced. This parameter determines how the
0051 %                   data are scaled before writing to the average file.  Ideally,
0052 %                   you want to make the data as big as possible, while remaining
0053 %                   in the limits of numeric representation range.  Kutaslab default
0054 %                   for this is 1000. {Default: 1000}
0055 % 'ica_correct'   - [1 | 0] If 1, ICA will be used to correct for EEG artifacts
0056 %                   and the average file produced will have a ".iccor"
0057 %                   extension. ICs are identified as artifacts via the field
0058 %                   EEG.iclabels.  See is_art.m for a list of labels recognized
0059 %                   as artifacts. If 0, ICA is not used to correct for artifacts
0060 %                   and the average file produced will have a ".nrm" extension.
0061 %                   {Default: 0}
0062 % 'just_cor_drty' - [1 | 0] If 1, ICA artifact correction will only be
0063 %                   applied to polluted trials and ERPs will be formed from both
0064 %                   artifact-free and corrected trials.  This argument has no
0065 %                   effect if the argument 'ica_correct' is set to 0. {Default: 0}
0066 % 'cln_avgs'      - [1 | 0] If 1, two additional average files will be produced
0067 %                   with the same root filename as out_fname. The first of
0068 %                   these will have the extension ".nrm" and will be formed
0069 %                   only from artifact free trials (artifact free according
0070 %                   to the arf file used to import the data into EEGLAB).
0071 %                   The second of these will have the extension ".icnrm"
0072 %                   and will be composed of the exact same ERPs in the ".nrm" file
0073 %                   after the ICA artifact correction filter has been applied.
0074 %                   Contrasting these two average files is useful for estimating
0075 %                   how ICA artifact correction might be distorting the data.
0076 %                   {Default: 0}
0077 % 'overwrite'     - [1 | 0] If 1, existing average files that have the same
0078 %                   name as the new files being produced will be automatically
0079 %                   overwritten.  If 0, the user will be asked if s/he wants
0080 %                   to overwrite them. {Default: 0}
0081 % 'verblevel'     - an integer specifiying the amount of information you want
0082 %                   functions to provide about what they are doing during runtime.
0083 %                    Options are:
0084 %                      0 - quiet, only show errors, warnings, and EEGLAB reports
0085 %                      1 - stuff anyone should probably know
0086 %                      2 - stuff you should know the first time you start working
0087 %                          with a data set {default value if not already globally
0088 %                          specified}
0089 %                      3 - stuff that might help you debug (show all reports)
0090 %  query_ccode    - [1 | 0] If one, the user will be asked at the
0091 %                   command line for the condition code of bins with zero
0092 %                   trials (if any).  If zero, the function will simply
0093 %                   use the condition code for the previous bin. {default: 0}
0094 %
0095 %
0096 % Outputs:
0097 %  An average file with a name specified by 'out_fname' is written to disk.
0098 %  Nothing is output to MATLAB.  If 'cln_avgs' option is used, two
0099 %  additional average files will be produced.
0100 %
0101 %
0102 % Example:
0103 % >> set_fname='/homes/dgroppe/SANDBOX/PRSENT/prsent41.set';
0104 % >> hdr_fname='/homes/dgroppe/SANDBOX/PRSENT/prsent.hdr';
0105 % >> set2avg(set_fname,'prsent41',hdr_fname,'bsln',[-100 0],'cln_avgs',1,'ica_correct',1);
0106 %
0107 % Author:
0108 % David Groppe
0109 % Kutaslab, 11/2009
0110 %
0111 % Notes:
0112 % -This function relies on erpio to create the avg files.  Thus it only
0113 % works on the Kutaslab network (unless you have a compiled copy of erpio
0114 % on whatever machine you're using).
0115 %
0116 % -Code was written to minimize memory demands at the cost of being slower.
0117 % If faster code is needed, the function could be revised.
0118 
0119 
0120 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%%
0121 %
0122 %12/11/09-Added ic_distortion option to measure degree of ERP distortion
0123 %caused by ICA artifact correction. Also, code now uses inputParser. -DG
0124 %
0125 %2/10/09-Made code able to deal with EEG.epoch(trial).eventccode being a cell
0126 %array (which it is if you run eeg_checkset.m on it). -DG
0127 %
0128 %7/29/10-Added query_ccode option and made it possible for set2avg to deal
0129 %with bins that have no trials. Also added a check to make sure that the
0130 %number of channels of cal pulse information matches the number of channels
0131 %in the data.-DG
0132 
0133 %Input Parser
0134 p=inputParser;
0135 %Note: the order of the required arguments needs to match precisely their
0136 %order in the function definition (which is also the order used by p.parse
0137 %below)
0138 p.addRequired('EEG_or_setname',@(x) isstruct(x) || ischar(x));
0139 p.addRequired('out_fname',@ischar);
0140 p.addRequired('hdr_fname',@ischar);
0141 p.addParamValue('est_ic_dstrt',0,@isnumeric);
0142 p.addParamValue('dstrt_bins',[],@isnumeric); %default is to use all bins
0143 p.addParamValue('bsln','allpre',@(x) isnumeric(x) || strcmpi(x,'allpre')); %default is to use all prestim time points
0144 p.addParamValue('pp10uv',1000,@isnumeric); %default is to use all time points
0145 p.addParamValue('verblevel',[],@(x) x>=0);
0146 p.addParamValue('cln_avgs',0,@isnumeric);
0147 p.addParamValue('just_cor_drty',0,@isnumeric);
0148 p.addParamValue('ica_correct',0,@isnumeric);
0149 p.addParamValue('overwrite',0,@isnumeric);
0150 p.addParamValue('query_ccode',0,@isnumeric);
0151 
0152 p.parse(EEG_or_setname,out_fname,hdr_fname,varargin{:});
0153 
0154 global VERBLEVEL
0155 
0156 if isempty(VERBLEVEL) && isempty(p.Results.verblevel)
0157     VERBLEVEL=2;
0158 elseif ~isempty(p.Results.verblevel)
0159     VERBLEVEL=p.Results.verblevel;
0160 end
0161 
0162 % Load/manage EEG struct
0163 if ischar(EEG_or_setname),
0164     %EEG_or_setname is a set filename
0165     EEG=pop_loadset(EEG_or_setname); %note this reports an informative error if the set file doesn't exist
0166     % It also reports on the command line what file it's attempting to load
0167 else
0168     %EEG_or_setname is an EEGLAB EEG struct
0169     EEG=EEG_or_setname; %note that this does not create a second copy of EEG
0170     %in memory since neither variable is yet modified.
0171     clear EEG_or_setname;
0172 end
0173 
0174 n_bin=length(EEG.bindesc);
0175 [n_chan, n_tpt, n_epoch]=size(EEG.data);
0176 
0177 %Check to make sure number of cal pulse channels equals the number of
0178 %EEG.data channels
0179 n_cal_pulse=size(EEG.cal_info.erps,1);
0180 if n_cal_pulse~=n_chan,
0181    error('You have %d channels of cal pulse averages (in EEG.cal_info.erps) but %d channels of EEG data (in EEG.data).\n', ...
0182        n_cal_pulse,n_chan);
0183 end
0184 n_cal_pulse=size(EEG.cal_info.usedcal_itemnum,1);
0185 if n_cal_pulse~=n_chan,
0186    error('You have %d channels of cal pulse item numbers (in EEG.cal_info.usedcal_itemnum) but %d channels of EEG data (in EEG.data).\n', ...
0187        n_cal_pulse,n_chan);
0188 end
0189 
0190 % Check to make sure header template exists
0191 [status, w]=unix(['ls ' hdr_fname]);
0192 if status,
0193     error('Header template file %s doesn''t exist.',hdr_fname);
0194 end
0195 
0196 % Check to make sure number of time points is a multiple of 256
0197 rmndr=rem(n_tpt,256);
0198 if rmndr,
0199     %fprintf('The number of time points per trial is not a multiple of 256.\n');
0200     n_pad=256-rmndr;
0201     padding=zeros(n_chan,n_pad);
0202     %fprintf('Adding %d time points of zeros to the end of each trial to correct this.\n',n_pad);
0203     VerbReport(sprintf(['The number of time points per trial is not a multiple of 256.\n' ...
0204         'Adding %d time points of zeros to the end of each trial to correct this.\n'],n_pad),1,VERBLEVEL);
0205 else
0206     n_pad=0;
0207     padding=[];
0208 end
0209 
0210 
0211 % Get baseline window
0212 if ~isempty(p.Results.bsln),
0213     if strcmpi(p.Results.bsln,'allpre'),
0214         %use all prestimulus time points
0215         bsln_tpt(1)=1;
0216         bsln_tpt(2)=find_tpt(0,EEG.times);
0217     else
0218         bsln_tpt(1)=find_tpt(p.Results.bsln(1),EEG.times);
0219         bsln_tpt(2)=find_tpt(p.Results.bsln(2),EEG.times);
0220     end
0221     bsln_win=bsln_tpt(1):bsln_tpt(2);
0222     if isempty(bsln_win),
0223         error('Baseline window from %.1f to %.1f ms includes no time points. Please revise it.', ...
0224             bsln_tpt(1),bsln_tpt(2));
0225     else
0226         VerbReport(sprintf('Removing mean EEG from %.1f ms to %.1f ms for each epoch/trial.\n', ...
0227             EEG.times(bsln_tpt(1)),EEG.times(bsln_tpt(2))),1,VERBLEVEL);
0228     end
0229 end
0230 
0231 
0232 
0233 
0234 %% Remove artifact ICs if they've been labeled
0235 
0236 % Check to see if some independent components have been labeled as artifacts
0237 fldnames=fieldnames(EEG);
0238 ic_arts=0;
0239 for a=1:length(fldnames),
0240     if strcmpi(fldnames{a},'iclabels'),
0241         ic_arts=1;
0242     end
0243 end
0244 
0245 if ic_arts,
0246     n=length(EEG.iclabels);
0247     art_ics=[];
0248     for dg=1:n,
0249         if is_art(EEG.iclabels{dg}),
0250             art_ics=[art_ics dg];
0251         end
0252     end
0253     
0254     unmix=EEG.icaweights*EEG.icasphere;
0255     mix=EEG.icawinv;
0256 
0257     if p.Results.est_ic_dstrt,
0258         crnt_level=VERBLEVEL;
0259         %Set VERBLEVEL 1 to suppress redudant report about baselining data
0260         ic_distortion(EEG,'bins',p.Results.dstrt_bins,'bsln',p.Results.bsln,'plot',1,'verblevel',1);
0261         VERBLEVEL=crnt_level; % Reset VERBLEVEL to what it was before
0262         drawnow;        
0263     else
0264         if VERBLEVEL>=1,
0265             fprintf('%d ICs have been labeled as EEG artifacts.\n',length(art_ics));
0266             if ~isempty(art_ics),
0267                 fprintf('The ICs labeled by artifacts are: \n');
0268                 for dg=art_ics,
0269                     fprintf('IC %d: %s\n',dg,EEG.iclabels{dg});
0270                 end
0271             end
0272         end
0273         eeglab_rej=find(EEG.reject.gcompreject==1);
0274         dif1=setdiff(eeglab_rej,art_ics);
0275         dif2=setdiff(art_ics,eeglab_rej);
0276         if ~isempty(dif1) || ~isempty(dif2)
0277             if isempty(eeglab_rej),
0278                 rej_ic_str='None';
0279             else
0280                 rej_ic_str=int2str(eeglab_rej);
0281             end
0282             msg=['The ICs marked for rejection by EEGLAB in EEG.reject.gcompreject differ from those labeled as artifacts in EEG.iclabels.' ...);
0283                 10 'Only the artifact labels in EEG.iclabels will be used.' 10 ...
0284                 'The ICs labeled by artifacts by EEG.reject.gcompreject are: ' rej_ic_str];
0285             watchit(msg);
0286         end
0287     end
0288     
0289     if p.Results.cln_avgs,
0290         % Create an avg file from just clean trials (i.e., clean according to arf
0291         % file)
0292         just_clean_trials=1;
0293         ica_correct=0;
0294         VerbReport('Making ERPs from just artifact free trials, without ICA artifact correction.',2,VERBLEVEL);
0295         write_avg(EEG,out_fname,hdr_fname,bsln_win,p.Results.pp10uv,just_clean_trials,ica_correct,p.Results.overwrite,n_pad,padding,art_ics,0,p.Results.query_ccode);
0296         
0297         % Filter the ERPs from just clean trials with ICA and write another average file.
0298         % Useful for seeing how much ICA might be distorting the data
0299         just_clean_trials=1;
0300         ica_correct=1;
0301         VerbReport('Making ERPs from just artifact free trials, with ICA artifact correction.',2,VERBLEVEL);
0302         write_avg(EEG,out_fname,hdr_fname,bsln_win,p.Results.pp10uv,just_clean_trials,ica_correct,p.Results.overwrite,n_pad,padding,art_ics,0,p.Results.query_ccode);
0303     end
0304     
0305     if p.Results.ica_correct,
0306         if p.Results.just_cor_drty
0307             VerbReport('Making ERPs from all trials, with ICA artifact correction applied only to polluted trials.',2,VERBLEVEL);
0308         else
0309             VerbReport('Making ERPs from all trials, with ICA artifact correction applied to all trials.',2,VERBLEVEL);
0310         end
0311         just_clean_trials=0;
0312         write_avg(EEG,out_fname,hdr_fname,bsln_win,p.Results.pp10uv,just_clean_trials,p.Results.ica_correct,p.Results.overwrite,n_pad,padding,art_ics,p.Results.just_cor_drty,p.Results.query_ccode);
0313     end
0314 else
0315     if p.Results.ica_correct,
0316         error(['This EEG variable does not have ICA components. You need to run ICA first before you can perform ICA ' ...
0317             'artifact correction.']);
0318     else
0319         just_clean_trials=1;
0320         VerbReport('Making ERPs from just clean trials, without ICA artifact correction.',2,VERBLEVEL);
0321         art_ics=[];
0322     end
0323     write_avg(EEG,out_fname,hdr_fname,bsln_win,p.Results.pp10uv,just_clean_trials,p.Results.ica_correct,p.Results.overwrite,n_pad,padding,art_ics,0,p.Results.query_ccode);
0324 end
0325 
0326 
0327 
0328 function write_avg(EEG,out_fname,hdr_fname,bsln_win,pp10uv,just_clean_trials,ica_correct,ovrwrite,n_pad,padding,art_ics,just_cor_drty,query_ccode)
0329 %function write_avg(EEG,out_fname,hdr_fname,bsln_win,pp10uv,just_clean_trials,ica_correct,ovrwrite,n_pad,padding,art_ics,just_cor_drty,query_ccode)
0330 %
0331 % Inputs:
0332 %  EEG                - EEGLAB EEG variable from which ERPs will be derived.
0333 %  out_fname          - Name of average file that will be produce.  The following
0334 %                       naming conventions are automatically tagged on:
0335 %                        *.drty=all trials, no ICA artifact correction
0336 %                        *.nrm=just clean trials (according to garv log flags), no
0337 %                              ICA artifact correction
0338 %                        *.iccor=all trials, yes ICA artifact correction
0339 %                        *.icnrm=just clean trials (according to garv log flags),
0340 %                              yes ICA artifact correction (useful for estimating
0341 %                      ERP distortion produced by ICA correction)
0342 %  hdr_fname          - Name of the header file or average file that will be used as
0343 %                       a template for the new average file being produced.  This file
0344 %                       will not be modified.
0345 %  bsln_win           - The time window (in milliseconds) that will be used to
0346 %                       baseline each trial
0347 %  pp10uv             - Points per 10 microvolts (determines scale of data in avg
0348 %                       file)
0349 %  just_clean_trials  - [1 | 0] If 1, only clean trials (according to garv
0350 %                       log flags) will be used to form ERPs.
0351 %  ica_correct        - [1 | 0] If 1, ICA will be used to remove EEG
0352 %                       artifacts from the EEG.
0353 %  ovrwrite           - [1 | 0] If 1, existing average files that have the
0354 %                       same name as the new files being produced will be
0355 %                       automatically overwritten.  If 0, user will be
0356 %                       asked if s/he wants to overwrite them.
0357 %  n_pad              - Number of 0 microvolt time points that will be
0358 %                       added to ERPs so that the number of time points are
0359 %                       a multiple of 256.
0360 %  art_ics            - [binary vector] If the ith element of art_ics is a
0361 %                       1, the ith IC corresponds to an EEG artifact.
0362 %  just_cor_drty      - [1 | 0] If 1, ICA artifact correction correction
0363 %                       is only applied to polluted trials (according to garv
0364 %                       log flags)
0365 %  query_ccode        - [scalar] If non-zero, the user will be asked at the
0366 %                       command line for the condition code of bins with zero
0367 %                       trials (if any).  If zero, the function will simply
0368 %                       use the condition code for the previous bin.
0369 
0370 global VERBLEVEL;
0371 
0372 % Data parameters
0373 s=size(EEG.data);
0374 n_tpt=s(2);
0375 n_chan=s(1);
0376 n_epoch=s(3);
0377 n_bin=length(EEG.bindesc);
0378 
0379 % Create avg file name
0380 if ~just_clean_trials && ica_correct,
0381     ext='.iccor';
0382 elseif just_clean_trials && ica_correct,
0383     ext='.icnrm';
0384 elseif ~just_clean_trials && ~ica_correct,
0385     ext='.drty'; % I don't know why would you ever want to make such an average,
0386     % but the option is included just in case.  It is not commented though in set2avg.m.
0387 else
0388     % just clean trials but no ica correction
0389     ext='.nrm';
0390 end
0391 
0392 if ~ismember('.',out_fname),
0393     %tack on auto-extension
0394     out_fname=[out_fname ext];
0395 else
0396     id=find(out_fname=='.');
0397     mx_id=max(id);
0398     if mx_id==1,
0399         error('Output filename %s cannot begin with a period.',out_fname);
0400     end
0401     %remove existing extension & tack on auto-extension
0402     out_fname=[out_fname(1:(mx_id-1)) ext];
0403 end
0404 
0405 
0406 % Check to make sure out_fname doesn't already exist
0407 [status, w]=unix(['ls ' out_fname]);
0408 if ~status,
0409     while ~ovrwrite
0410         rsp=input(sprintf('File %s already exists.  Do you want to overwrite it (Y or N)?',out_fname),'s');
0411         if strcmpi(rsp,'Y') || strcmpi(rsp,'yes'),
0412             ovrwrite=1;
0413         elseif strcmpi(rsp,'N') || strcmpi(rsp,'no'),
0414             ovrwrite=-1;
0415         end
0416     end
0417 
0418     if ovrwrite==1,
0419             % Remove current version of output file if it already exists
0420             VerbReport(sprintf('Removing old version of file %s.',out_fname),1,VERBLEVEL);
0421             [status, w]=unix(['rm -f ' out_fname]); % ?? add an error check?
0422             %Note, status and w provide no information about the success of rm -f
0423     elseif ovrwrite==-1,
0424         fprintf('\n'); % line break to make command line reports look better
0425         return;
0426     end
0427 end
0428 
0429 
0430 % Create avg file
0431 VerbReport(sprintf('Creating average file %s using header information from %s.', ...
0432     out_fname,hdr_fname),1,VERBLEVEL);
0433 fid=erpio('createavg',out_fname,hdr_fname);
0434 if (fid==-1)
0435     error('Could not create ERP file %s using erpio.',out_fname);
0436 end
0437 erpio('put_hdrvar',fid,'pp10uv',pp10uv);
0438 
0439 
0440 erpio('put_hdrvar',fid,'cprecis',(n_tpt+n_pad)/256);
0441 erpio('put_hdrvar',fid,'expdesc',EEG.comments);
0442 erpio('put_hdrvar',fid,'subdesc',EEG.subject);
0443 erpio('put_hdrvar',fid,'presampling',-EEG.times(1));
0444 %erpio('put_hdrvar',fid,'epoch',EEG.times(end)-EEG.times(1)+1000/EEG.srate);
0445 %--at least one of my set files had EEG.times in units of seconds (weird)
0446 erpio('put_hdrvar',fid,'epoch',1000*(n_pad+n_tpt)/EEG.srate);
0447 erpio('put_hdrvar',fid,'chans',n_chan);
0448 %erpio('put_hdrvar',fid,'prfdesc','average'); ?? makes matlab crash
0449 erpio('put_hdrvar',fid,'clktick',100000/EEG.srate);
0450 for c=1:n_chan,
0451     erpio('put_hdrvar',fid,'chndesc',c-1,EEG.chanlocs(c).labels); 
0452 end
0453 
0454 %pre-allocate memory
0455 erps=zeros(n_chan,n_tpt,n_bin);
0456 erp_ct=zeros(1,n_bin);
0457 ccode=-ones(1,n_bin);
0458 
0459 % create ICA filter (if requested)
0460 if ica_correct,
0461     non_art=setdiff(1:n_chan,art_ics);
0462     unmix=EEG.icaweights*EEG.icasphere;
0463     mix=inv(unmix); %just in case EEG.icawinv has been scaled differently than mixing matrix
0464     ic_filt=mix(:,non_art)*unmix(non_art,:);
0465 end
0466 
0467 %add up all trials
0468 for trial=1:n_epoch,
0469     if iscell(EEG.epoch(trial).eventlogflag)
0470         eventLogFlag=0;
0471         for lfLoop=1:length(EEG.epoch(trial).eventlogflag)
0472             eventLogFlag=eventLogFlag+EEG.epoch(trial).eventlogflag{lfLoop};
0473         end
0474     else
0475         eventLogFlag=EEG.epoch(trial).eventlogflag;
0476     end
0477     if ~just_clean_trials || ~eventLogFlag
0478         % clean trials have an eventlogflag of 0
0479         if ica_correct,
0480             if just_cor_drty && eventLogFlag,
0481                 sweep=ic_filt*squeeze(EEG.data(:,:,trial));
0482             elseif ~just_cor_drty,
0483                 sweep=ic_filt*squeeze(EEG.data(:,:,trial));
0484             else
0485                 sweep=squeeze(EEG.data(:,:,trial));
0486             end
0487         else
0488             sweep=squeeze(EEG.data(:,:,trial));
0489         end
0490         sweep=rmbase(sweep,n_tpt,bsln_win);
0491         for types=1:length(EEG.epoch(trial).eventtype),
0492             if strcmpi(EEG.epoch(trial).eventtype{types}(1:3),'bin'),
0493                 bin_num=str2num(EEG.epoch(trial).eventtype{types}(4:end));
0494                 erps(:,:,bin_num)=erps(:,:,bin_num)+sweep;
0495                 erp_ct(bin_num)=erp_ct(bin_num)+1;
0496                 if iscell(EEG.epoch(trial).eventccode), % if eeg_checkset has been run eventccode will be a cell array
0497                     if ccode(bin_num)==-1,
0498                         ccode(bin_num)=EEG.epoch(trial).eventccode{types};
0499                     elseif ccode(bin_num)~=EEG.epoch(trial).eventccode{types},
0500                         error('Trial %d has a Condition Code of %d and contributes to Bin %d, but at least one other trial in that bin has a Condition Code of %d.', ...
0501                             trial,EEG.epoch(trial).eventccode{types},bin_num,ccode(bin_num));
0502                     end
0503                 else
0504                     if ccode(bin_num)==-1,
0505                         ccode(bin_num)=EEG.epoch(trial).eventccode;
0506                     elseif ccode(bin_num)~=EEG.epoch(trial).eventccode,
0507                         error('Trial %d has a Condition Code of %d and contributes to Bin %d, but at least one other trial in that bin has a Condition Code of %d.', ...
0508                             trial,EEG.epoch(trial).eventccode,bin_num,ccode(bin_num));
0509                     end
0510                 end
0511             end
0512         end
0513     end
0514 end
0515 
0516 crnt_cond=0;
0517 for b=0:n_bin,
0518     if b==0,
0519         erpio('put_hdrvar',fid,'condcode',0);
0520         erpio('put_hdrvar',fid,'condesc','Amp Calibration');
0521         erpio('put_hdrvar',fid,'sums',EEG.cal_info.npulse_used);
0522         erpio('put_hdrvar',fid,'totrawrecs',EEG.cal_info.npulse_psbl);
0523         erpio('writebin',fid,[EEG.cal_info.erps padding]*pp10uv/10,'Calibration Pulses');
0524         %note erpio returns nothing for 'writebin' (i.e. no indication
0525         %of success)
0526     else
0527         if ccode(b)<0,
0528             if query_ccode,
0529                 fprintf('There are no trials for Bin %d thus I don''t know what its condition code is.\n',b);
0530                 ccode(b)=input(['Enter the condition code for Bin ' int2str(b) ': ']);
0531             else
0532                 if VERBLEVEL>=2,
0533                     fprintf('There are no trials for Bin %d. I will use the condition code for the previous bin (i.e., %d).\n',b,crnt_cond);
0534                 end
0535                 ccode(b)=crnt_cond;
0536             end
0537         end
0538         if crnt_cond~=ccode(b),
0539             crnt_cond=ccode(b);
0540             erpio('put_hdrvar',fid,'condcode',crnt_cond);
0541             erpio('put_hdrvar',fid,'condesc',EEG.condesc{crnt_cond});
0542         end
0543         if VERBLEVEL>=2,
0544             fprintf('Bin %d, CCode %d, Trials %d(%d): %s\n',b,ccode(b),erp_ct(b),EEG.rawtrials_per_bin(b),EEG.bindesc{b});
0545         end
0546         %write avg
0547         erpio('put_hdrvar',fid,'sums',erp_ct(b));
0548         erpio('put_hdrvar',fid,'totrawrecs',EEG.rawtrials_per_bin(b));
0549         if erp_ct(b),
0550             erpio('writebin',fid,[squeeze(erps(:,:,b)*pp10uv/(10*erp_ct(b))) padding],EEG.bindesc{b});
0551         else
0552             erpio('writebin',fid,[squeeze(erps(:,:,b)) padding],EEG.bindesc{b});
0553         end
0554     end
0555 end
0556 
0557 %close file
0558 erpio('close',fid);
0559 
0560 fprintf('\n'); % line break to make command line reports look better
0561 
0562

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