sets2specGND

 sets2specGND() - Create a MATLABmk specGND struct variable from a set of
                  EEGLAB *.set files.  specGND variables contain power 
                  spectra from individual participants and the grand
                  average of those spectra (as well as supportive 
                  information like channel locations). Power spectra is
                  estimated using Slepian multitapers or a Hann taper and
                  the discrete Fourier transfrom (DFT).  A Hann taper is
                  used if the number of tapers requested is 1, otherwise
                  Slepian tapers are used.

 Usage:
  >> specGND=sets2specGND(gui_infiles_or_tmplt,varargin);

 Required Inputs:
   gui_infiles_or_tmplt - ['gui', a cell array of strings, or a string template]
                          If 'gui', a GUI is created that allows you to select
                          which set files to import (this is probably the
                          easiest way to import files).  If a cell array of
                          of strings, each element of the cell array should
                          contain the filename of an EEGLAB set file (e.g.,
                          {'visodbl01.set','visodbl02.set'}.  Otherwise, this
                          input should be a filename template (i.e., a string with
                          # where the subject number should be--'visodbl#.set').
                          If you provide a template, you must use the option
                          'sub_ids' (see below) to specify the ID numbers of
                          each subject. Include the files' path unless the files
                          are in the current working directory.


 Optional Inputs:
   sub_ids          - [integer vector] A set of integers specifying the
                      subject ID numbers to include in the grand average.
                      Only necessary if a filename template is given as
                      the input to gui_infiles_or_tmplt.

   use_bins         - [integer vector] A set of integers specifying which
                      bins to import into MATLAB.  If not specified, all
                      bins will be imported. Note, if you import only a
                      subset of bins, the bin numberings will start at 1
                      and run to the number you've imported (i.e., they may
                      differ from the bin numbers in the set files.
                      {default: import all bins}

   exclude_chans    - A cell array of channel labels to exclude from the
                      importation (e.g., {'A2','lle','rhe'}). You cannot
                      use both this option and 'include_chans' (below).{default:
                      not used, import all channels}

   include_chans    - A cell array specifying a subset of channel labels to import
                      (e.g., {'Fz','Cz','Pz'}).  All other channels will
                      be ignored. You cannot use both this option and
                      'exclude_chans' (above). {default: not used, import
                      all channels}

   exp_name         - [string] Name of the experiment. {default: 'An
                      Experiment'}

   out_fname        - [string] Filename to save specGND variable to.  If empty
                      (i.e., not specified), a GUI will be created to prompt
                      you for a filename. If the string 'no save', the specGND
                      variable will not be saved to disk. If no file
                      extension is given '.specGND' will be added to the
                      filename.  If not specified, a GUI will be created
                      to ask for a filename. {default: not specified}

   time_wind        - [start stop] The peri-event time window to perform
                      the DFT on in milliseconds. For example if 
                      time_wind=[-100 900], all time points from 100 ms 
                      before to 900 ms after selected events will be
                      used for DFT computation.  For epoched EEGLAB set 
                      files, the default is to use the full epoch length.
                      For continuous EEGLAB set files, the default is 
                      [500 1500].

   n_tapers         - [integer] The number of Slepian tapers to use when 
                      estimating the power spectra. Increasing the number 
                      of tapers increases the reliability of the estimate 
                      at the cost of reduced frequency resolution. The
                      frequency resolution half-bandwidth equals
                      (K+1)/(2*T), where K=# of tapers and T=the length of
                      the time window (see 'time_wind' argument) in seconds.
                      Minimum # of tapers is 1. Maximum # of tapers is the
                      biggest integer less than (T*sampling rate)-1.Note, 
                      'n_tapers' and 'half_bandwidth' optional inputs 
                      cannot BOTH be used. Moreover, if only 1 taper is 
                      requested, a Hann taper is used instead of a Slepian
                      taper. {default: 1}

   half_bandwidth   - [scalar] The frequency resolution "half bandwidth" 
                      (in Hz).  Multiply this by two to get the effective 
                      minimal frequency resolution of your estimate. Note, 
                      'n_tapers' and 'half_bandwidth' optional inputs cannot 
                      BOTH be used.{default: (K+1)/(2*T)<-see n_tapers 
                      optional input above}

   pad              - [integer] The amount of zero padding to add to each 
                      window of data before computing the Fourier
                      transform of the data. Zero padding effectively 
                      gives you spectral estimates at more frequencies by 
                      interpolating between frequencies. It can also speed 
                      up the Fourier transform by making the number of time
                      points a multiple of two. pad can take values
                      -1,0,1,2,etc... -1 corresponds to no padding, 0 
                      corresponds to padding to the next highest power of 
                      2. 1 corresponds to padding to the second next 
                      highest power of two etc... For example, for N = 500,
                      if pad=-1, we do not pad; if pad=0, we pad the DFT
                         to 512 points, if pad=1, we pad to 1024 points etc.
                         {default: 0}
   rm_mean          - ['yes' or 'no'] If 'yes', the mean of each epoch of
                      data is removed before computing power spectra.
                      {default: 'yes'}

   verblevel        - An integer specifiying the amount of information you want
                      this function to provide about what it is 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}
                        3 - stuff that might help you debug (show all
                            reports)

 Output:
   specGND  - Struct variable containing grand average power spectra, 
              individual subject power spectra, mass univariate t-test 
              results and more.


 Global Variables:
   VERBLEVEL - MATLABmk level of verbosity (i.e., tells 
               functions how much to report about what they're doing during
               runtime) set by the optional function argument 'verblevel'

 Notes:
 -The spectral estimation in this function is performed by the pmtm.m 
 function that is part of Mathworks's Signal Processing Toolbox

 -If you specify a desired number of tapers, n, this function chooses a
 frequency resolution bandwidth such that you can derive n+1 tapers and
 ignore the last taper (which is apparently not that useful).

 -The mean of each epoch of data is removed to dampen drift and other low
 frequency activity before performing the DFT.

 -This function expects the EEG variable to have information about bins
 stored in it.  More specifically, there needs to be (1) an EEG.bindesc
 field that contains a text description of each bin and (2) epochs of data
 labeled as bin types (e.g., 'bin1').  Bin types are stored in the field
 EEG.event(#).type for both continuous and epoched data. For epoched data, 
 bin types are also stored in the field EEG.epoch(#).eventtype. Use the 
 function bin_info2EEG.m to add bin information to an EEG variable.

 -ICs labeled as artifacts will be removed from the data before DFT is
 performed.  This is done by the function remove_artifact_ics.m.  IC labeling
 can be done via EEGLAB or MATLABmk functions/conventions.  If both
 conventions appear in the EEG variable, only MATLABmk IC labels will be used.

 -This function is not yet compatible with ERPLAB's conventions for
 storing bin information.

 -The specGND field indiv_bin_raw_ct and cals are specific to Kutaslab data.
 Other labs should be able to ignore them. 

 Author:
 David Groppe
 Kutaslab, 12/2010

0001 % sets2specGND() - Create a MATLABmk specGND struct variable from a set of
0002 %                  EEGLAB *.set files.  specGND variables contain power
0003 %                  spectra from individual participants and the grand
0004 %                  average of those spectra (as well as supportive
0005 %                  information like channel locations). Power spectra is
0006 %                  estimated using Slepian multitapers or a Hann taper and
0007 %                  the discrete Fourier transfrom (DFT).  A Hann taper is
0008 %                  used if the number of tapers requested is 1, otherwise
0009 %                  Slepian tapers are used.
0010 %
0011 % Usage:
0012 %  >> specGND=sets2specGND(gui_infiles_or_tmplt,varargin);
0013 %
0014 % Required Inputs:
0015 %   gui_infiles_or_tmplt - ['gui', a cell array of strings, or a string template]
0016 %                          If 'gui', a GUI is created that allows you to select
0017 %                          which set files to import (this is probably the
0018 %                          easiest way to import files).  If a cell array of
0019 %                          of strings, each element of the cell array should
0020 %                          contain the filename of an EEGLAB set file (e.g.,
0021 %                          {'visodbl01.set','visodbl02.set'}.  Otherwise, this
0022 %                          input should be a filename template (i.e., a string with
0023 %                          # where the subject number should be--'visodbl#.set').
0024 %                          If you provide a template, you must use the option
0025 %                          'sub_ids' (see below) to specify the ID numbers of
0026 %                          each subject. Include the files' path unless the files
0027 %                          are in the current working directory.
0028 %
0029 %
0030 % Optional Inputs:
0031 %   sub_ids          - [integer vector] A set of integers specifying the
0032 %                      subject ID numbers to include in the grand average.
0033 %                      Only necessary if a filename template is given as
0034 %                      the input to gui_infiles_or_tmplt.
0035 %
0036 %   use_bins         - [integer vector] A set of integers specifying which
0037 %                      bins to import into MATLAB.  If not specified, all
0038 %                      bins will be imported. Note, if you import only a
0039 %                      subset of bins, the bin numberings will start at 1
0040 %                      and run to the number you've imported (i.e., they may
0041 %                      differ from the bin numbers in the set files.
0042 %                      {default: import all bins}
0043 %
0044 %   exclude_chans    - A cell array of channel labels to exclude from the
0045 %                      importation (e.g., {'A2','lle','rhe'}). You cannot
0046 %                      use both this option and 'include_chans' (below).{default:
0047 %                      not used, import all channels}
0048 %
0049 %   include_chans    - A cell array specifying a subset of channel labels to import
0050 %                      (e.g., {'Fz','Cz','Pz'}).  All other channels will
0051 %                      be ignored. You cannot use both this option and
0052 %                      'exclude_chans' (above). {default: not used, import
0053 %                      all channels}
0054 %
0055 %   exp_name         - [string] Name of the experiment. {default: 'An
0056 %                      Experiment'}
0057 %
0058 %   out_fname        - [string] Filename to save specGND variable to.  If empty
0059 %                      (i.e., not specified), a GUI will be created to prompt
0060 %                      you for a filename. If the string 'no save', the specGND
0061 %                      variable will not be saved to disk. If no file
0062 %                      extension is given '.specGND' will be added to the
0063 %                      filename.  If not specified, a GUI will be created
0064 %                      to ask for a filename. {default: not specified}
0065 %
0066 %   time_wind        - [start stop] The peri-event time window to perform
0067 %                      the DFT on in milliseconds. For example if
0068 %                      time_wind=[-100 900], all time points from 100 ms
0069 %                      before to 900 ms after selected events will be
0070 %                      used for DFT computation.  For epoched EEGLAB set
0071 %                      files, the default is to use the full epoch length.
0072 %                      For continuous EEGLAB set files, the default is
0073 %                      [500 1500].
0074 %
0075 %   n_tapers         - [integer] The number of Slepian tapers to use when
0076 %                      estimating the power spectra. Increasing the number
0077 %                      of tapers increases the reliability of the estimate
0078 %                      at the cost of reduced frequency resolution. The
0079 %                      frequency resolution half-bandwidth equals
0080 %                      (K+1)/(2*T), where K=# of tapers and T=the length of
0081 %                      the time window (see 'time_wind' argument) in seconds.
0082 %                      Minimum # of tapers is 1. Maximum # of tapers is the
0083 %                      biggest integer less than (T*sampling rate)-1.Note,
0084 %                      'n_tapers' and 'half_bandwidth' optional inputs
0085 %                      cannot BOTH be used. Moreover, if only 1 taper is
0086 %                      requested, a Hann taper is used instead of a Slepian
0087 %                      taper. {default: 1}
0088 %
0089 %   half_bandwidth   - [scalar] The frequency resolution "half bandwidth"
0090 %                      (in Hz).  Multiply this by two to get the effective
0091 %                      minimal frequency resolution of your estimate. Note,
0092 %                      'n_tapers' and 'half_bandwidth' optional inputs cannot
0093 %                      BOTH be used.{default: (K+1)/(2*T)<-see n_tapers
0094 %                      optional input above}
0095 %
0096 %   pad              - [integer] The amount of zero padding to add to each
0097 %                      window of data before computing the Fourier
0098 %                      transform of the data. Zero padding effectively
0099 %                      gives you spectral estimates at more frequencies by
0100 %                      interpolating between frequencies. It can also speed
0101 %                      up the Fourier transform by making the number of time
0102 %                      points a multiple of two. pad can take values
0103 %                      -1,0,1,2,etc... -1 corresponds to no padding, 0
0104 %                      corresponds to padding to the next highest power of
0105 %                      2. 1 corresponds to padding to the second next
0106 %                      highest power of two etc... For example, for N = 500,
0107 %                      if pad=-1, we do not pad; if pad=0, we pad the DFT
0108 %                         to 512 points, if pad=1, we pad to 1024 points etc.
0109 %                         {default: 0}
0110 %   rm_mean          - ['yes' or 'no'] If 'yes', the mean of each epoch of
0111 %                      data is removed before computing power spectra.
0112 %                      {default: 'yes'}
0113 %
0114 %   verblevel        - An integer specifiying the amount of information you want
0115 %                      this function to provide about what it is doing during runtime.
0116 %                       Options are:
0117 %                        0 - quiet, only show errors, warnings, and EEGLAB reports
0118 %                        1 - stuff anyone should probably know
0119 %                        2 - stuff you should know the first time you start working
0120 %                            with a data set {default value}
0121 %                        3 - stuff that might help you debug (show all
0122 %                            reports)
0123 %
0124 % Output:
0125 %   specGND  - Struct variable containing grand average power spectra,
0126 %              individual subject power spectra, mass univariate t-test
0127 %              results and more.
0128 %
0129 %
0130 % Global Variables:
0131 %   VERBLEVEL - MATLABmk level of verbosity (i.e., tells
0132 %               functions how much to report about what they're doing during
0133 %               runtime) set by the optional function argument 'verblevel'
0134 %
0135 % Notes:
0136 % -The spectral estimation in this function is performed by the pmtm.m
0137 % function that is part of Mathworks's Signal Processing Toolbox
0138 %
0139 % -If you specify a desired number of tapers, n, this function chooses a
0140 % frequency resolution bandwidth such that you can derive n+1 tapers and
0141 % ignore the last taper (which is apparently not that useful).
0142 %
0143 % -The mean of each epoch of data is removed to dampen drift and other low
0144 % frequency activity before performing the DFT.
0145 %
0146 % -This function expects the EEG variable to have information about bins
0147 % stored in it.  More specifically, there needs to be (1) an EEG.bindesc
0148 % field that contains a text description of each bin and (2) epochs of data
0149 % labeled as bin types (e.g., 'bin1').  Bin types are stored in the field
0150 % EEG.event(#).type for both continuous and epoched data. For epoched data,
0151 % bin types are also stored in the field EEG.epoch(#).eventtype. Use the
0152 % function bin_info2EEG.m to add bin information to an EEG variable.
0153 %
0154 % -ICs labeled as artifacts will be removed from the data before DFT is
0155 % performed.  This is done by the function remove_artifact_ics.m.  IC labeling
0156 % can be done via EEGLAB or MATLABmk functions/conventions.  If both
0157 % conventions appear in the EEG variable, only MATLABmk IC labels will be used.
0158 %
0159 % -This function is not yet compatible with ERPLAB's conventions for
0160 % storing bin information.
0161 %
0162 % -The specGND field indiv_bin_raw_ct and cals are specific to Kutaslab data.
0163 % Other labs should be able to ignore them.
0164 %
0165 % Author:
0166 % David Groppe
0167 % Kutaslab, 12/2010
0168 
0169 %%%%%%%%%%%%%%%% Revision History %%%%%%%%%%%%%%%%%
0170 % 5/1/2011: Pre-allocated ize of cal pulse ERPs is taken from cal pulse
0171 % ERPs stored in EEG variable as opposed to size of EEG epochs (i.e.
0172 % size(EEG.data))
0173 
0174 %%%%%%%%%%%%%%%% Future Work %%%%%%%%%%%%%%%%%
0175 % -add flexible high pass filter option?
0176 
0177 
0178 function specGND=sets2specGND(gui_infiles_or_tmplt,varargin)
0179 %Note, condesc is ignored since we don't have to write to avg files
0180 
0181 p=inputParser;
0182 p.addRequired('gui_infiles_or_tmplt',@(x) ischar(x) || iscell(x));
0183 p.addParamValue('sub_ids',[],@isnumeric);
0184 p.addParamValue('exp_name','An Experiment',@ischar);
0185 p.addParamValue('use_bins',[],@isnumeric);
0186 p.addParamValue('exclude_chans',[],@(x) ischar(x) || iscell(x));
0187 p.addParamValue('include_chans',[],@(x) ischar(x) || iscell(x));
0188 p.addParamValue('out_fname',[],@ischar);
0189 p.addParamValue('rm_mean','yes',@ischar);
0190 p.addParamValue('n_tapers',[],@(x) isnumeric(x) && (length(x)==1));
0191 p.addParamValue('half_bandwidth',[],@(x) isnumeric(x) && (length(x)==1)); %in Hz
0192 p.addParamValue('pad',0,@(x) isnumeric(x) && (length(x)==1));
0193 p.addParamValue('time_wind',[],@(x) isnumeric(x) && (length(x)==2));
0194 p.addParamValue('verblevel',[],@(x) isnumeric(x) && (length(x)==1));
0195 
0196 p.parse(gui_infiles_or_tmplt,varargin{:});
0197 
0198 
0199 global EEG;
0200 global VERBLEVEL;
0201 
0202 if isempty(p.Results.verblevel),
0203     if isempty(VERBLEVEL),
0204         VERBLEVEL=2;
0205     end
0206 else
0207     VERBLEVEL=p.Results.verblevel;
0208 end
0209 
0210 
0211 %% Figure out which channels to ignore if any
0212 %But first make sure exclude & include options were not both used.
0213 if ~isempty(p.Results.include_chans) && ~isempty(p.Results.exclude_chans)
0214     error('You cannot use BOTH ''include_chans'' and ''exclude_chans'' options.');
0215 end
0216 if ischar(p.Results.exclude_chans),
0217     exclude_chans{1}=p.Results.exclude_chans;
0218 elseif isempty(p.Results.exclude_chans)
0219     exclude_chans=[];
0220 else
0221     exclude_chans=p.Results.exclude_chans;
0222 end
0223 if ischar(p.Results.include_chans),
0224     include_chans{1}=p.Results.include_chans;
0225 elseif isempty(p.Results.include_chans)
0226     include_chans=[];
0227 else
0228     include_chans=p.Results.include_chans;
0229 end
0230 
0231 
0232 %%  Select files for loading
0233 if strcmpi(gui_infiles_or_tmplt,'GUI'),
0234     loading=1;
0235     infiles=[];
0236     while loading,
0237         [neofname, inpath]=uigetfile({'*.set','*.set files'; ...
0238             '*.*','All Files (*.*)'},'EEGLAB set Files to Import','MultiSelect','on');
0239         if ischar(neofname),
0240             clear infname;
0241             infname{1}=neofname; %make it a cell array for consistent syntax below
0242         else
0243             infname=neofname;
0244         end
0245         if ~inpath,
0246             if isempty(infiles),
0247                 fprintf('File selection cancelled.  Aborting sets2specGND.\n');
0248                 specGND=[];
0249                 return;
0250             else
0251                 loading=0;
0252             end
0253         else
0254             if isempty(infiles),
0255                 infiles=cell(1,length(infname)); %preallocate mem
0256                 for a=1:length(infname),
0257                     infiles{a}=[inpath infname{a}];
0258                 end
0259             else
0260                 n_files=length(infiles);
0261                 ct=0;
0262                 for a=(n_files+1):(n_files+length(infname)),
0263                     ct=ct+1;
0264                     infiles{a}=[inpath infname{ct}];
0265                 end
0266             end
0267             % trigger GUI to see if user wants to load more
0268             resp=questdlg('Do you want to load more files?',...
0269                 'OUTPUT','Yes','No','Yes');
0270             if strcmpi(resp,'No'),
0271                 loading=0;
0272             end
0273         end
0274     end
0275 elseif iscell(p.Results.gui_infiles_or_tmplt)
0276     infiles=p.Results.gui_infiles_or_tmplt;
0277 else
0278     if isempty(p.Results.sub_ids)
0279         error('If gui_infiles_or_tmplt is a string filename template, you must specify the subject numbers with the argument ''sub_ids''.');
0280     elseif ~ismember('#',p.Results.gui_infiles_or_tmplt),
0281         error('The filename template %s needs to contain a # to indicate where the subject numbers go (e.g., odbl#.nrm).', ...
0282             p.Results.gui_infiles_or_tmplt);
0283     else
0284         lb_id=find(p.Results.gui_infiles_or_tmplt=='#');
0285         prefix=p.Results.gui_infiles_or_tmplt(1:(lb_id(1)-1)); %indexing lb_id by 1 in case there are multiple pound signs
0286         postfix=p.Results.gui_infiles_or_tmplt((lb_id(1)+1):end);
0287         n_infiles=length(p.Results.sub_ids);
0288         infiles=cell(1,n_infiles);
0289         for s=1:n_infiles,
0290             no_pad=[prefix num2str(p.Results.sub_ids(s)) postfix];
0291             if p.Results.sub_ids(s)<10,
0292                 padded=[prefix num2str(p.Results.sub_ids(s),'%.2d') postfix];
0293                 [sP wP]=unix(['ls ' padded]); %if sp==0, then the file exists, need wP argument so that it isn't automatically displayed in command window
0294                 [sNP wNP]=unix(['ls ' no_pad]); %Need wP argument so that it isn't automatically displayed in command window
0295                 if ~sP
0296                     if ~sNP,
0297                         error('You have a file named %s and one named %s.  I don''t know which one to load!\n', ...
0298                             padded,no_pad);
0299                     else
0300                         infiles{s}=padded;
0301                     end
0302                 else
0303                     infiles{s}=no_pad;
0304                 end
0305             else
0306                 infiles{s}=no_pad;
0307             end
0308         end
0309     end
0310 end
0311 
0312 %Get rid of any redundant input files
0313 infiles=unique(infiles);
0314 n_infiles=length(infiles);
0315 
0316 
0317 %% Load first data set and draft specGND fields
0318 EEG=pop_loadset(infiles{1});
0319 EEG_var_check(EEG,infiles{1});
0320 EEG=epoch_if_necessary(EEG,p);
0321 [n_EEG_chans, n_tpts, n_epochs]=size(EEG.data);
0322 
0323 %Get indices of usable channels
0324 if ~isempty(exclude_chans),
0325     use_chans=zeros(1,n_EEG_chans);
0326     ex_ct=0;
0327     for a=1:n_EEG_chans,
0328         if ~ismember_ci(EEG.chanlocs(a).labels,exclude_chans)
0329             use_chans(a)=1;
0330         else
0331             ex_ct=ex_ct+1;
0332             ex_labels{ex_ct}=EEG.chanlocs(a).labels;
0333         end
0334     end
0335     n_chans=sum(use_chans);
0336     use_chans=find(use_chans==1);
0337     
0338     if VERBLEVEL>1,
0339         missed=setdiff_ci(exclude_chans,ex_labels);
0340         n_missed=length(missed);
0341         if n_missed,
0342             if n_missed==1,
0343                 msg=sprintf('I attempted to exclude the following channel, but it was not found:');
0344             else
0345                 msg=sprintf('I attempted to exclude the following channels, but they were not found:');
0346             end
0347             for a=1:n_missed,
0348                 msg=[msg ' ' missed{a}];
0349             end
0350             watchit(msg);
0351         end
0352     end
0353 elseif ~isempty(include_chans),
0354     use_chans=zeros(1,n_EEG_chans);
0355     in_ct=0;
0356     for a=1:n_EEG_chans,
0357         if ismember_ci(EEG.chanlocs(a).labels,include_chans)
0358             use_chans(a)=1;
0359             in_ct=in_ct+1;
0360             in_labels{in_ct}=EEG.chanlocs(a).labels;
0361         end
0362     end
0363     n_chans=sum(use_chans);
0364     use_chans=find(use_chans==1);
0365     
0366     if VERBLEVEL>1,
0367         missed=setdiff_ci(include_chans,in_labels);
0368         n_missed=length(missed);
0369         if n_missed,
0370             if n_missed==1,
0371                 msg=sprintf('I attempted to include the following channel, but it was not found:');
0372             else
0373                 msg=sprintf('I attempted to include the following channels, but they were not found:');
0374             end
0375             for a=1:n_missed,
0376                 msg=[msg ' ' missed{a}];
0377             end
0378             watchit(msg);
0379         end
0380     end
0381 else
0382     n_chans=n_EEG_chans;
0383     use_chans=1:n_chans;
0384 end
0385 
0386 n_EEG_bins=length(EEG.bindesc);
0387 if isempty(p.Results.use_bins),
0388     use_bins=1:n_EEG_bins;
0389     n_bins=n_EEG_bins;
0390 else
0391     %modicum of error checking
0392     if min(p.Results.use_bins)<1,
0393         error('You cannot request to import bins with indices less than 1 (e.g., "Bin -1").');
0394     elseif max(p.Results.use_bins)>n_EEG_bins,
0395         error('File %s only has %d bins, but you''ve requested to import Bin %d.\n', ...
0396             infiles{1},n_EEG_bins,max(p.Results.use_bins));
0397     end
0398     use_bins=unique(p.Results.use_bins);
0399     n_bins=length(use_bins);
0400 end
0401 %initialize cell array for keeping track of used urevents
0402 used_urevents=cell(n_infiles,n_bins); %at most there will be n_infiles subjects
0403 
0404 specGND.exp_desc=p.Results.exp_name;
0405 specGND.filename=[];
0406 specGND.filepath=[];
0407 specGND.saved='no';
0408 if isempty(p.Results.time_wind),
0409     use_tpts=1:n_tpts;
0410 else
0411     tpt1=find_tpt(p.Results.time_wind(1),EEG.times);
0412     tpt2=find_tpt(p.Results.time_wind(2),EEG.times);
0413     use_tpts=tpt1:tpt2;
0414     if length(use_tpts)<2,
0415         error('Your time window is too small!  It only contains one time point.');
0416     end
0417 end
0418 n_use_tpts=length(use_tpts);
0419 if VERBLEVEL>1,
0420     fprintf('Using all time points between %d and %d ms (i.e., from time point %d to %d).\n', ...
0421         EEG.times(use_tpts(1)),EEG.times(use_tpts(end)),use_tpts(1),use_tpts(end));
0422 end
0423 
0424 %Set mtspectrum.m/pwr_spec.m parameters (the functions that perform the DFT)
0425 %params.Fs=EEG.srate; %Old Chronux code
0426 T=n_use_tpts/EEG.srate; %length of epoch in seconds
0427 if ~isempty(p.Results.n_tapers),
0428     K=p.Results.n_tapers;
0429     if K>=(T*EEG.srate)-1,
0430         error('Too many tapers.  The maximum number of possible tapers for this size time window and sampling rate is: %d', ...
0431             (T*EEG.srate)-2);
0432     end
0433     delta_f=(K+1)/(2*T);
0434     if ~isempty(p.Results.half_bandwidth) && VERBLEVEL,
0435         watchit('''bandwidth'' and ''n_tapers'' arguments to sets2specGND.m both used.  Only one is needed and I will use ''n_tapers''.');
0436     end
0437     if VERBLEVEL>1
0438         fprintf('Frequency resolution half-bandwidth is: %f Hz\n',delta_f);
0439     end
0440 elseif ~isempty(p.Results.half_bandwidth),
0441     delta_f=p.Results.half_bandwidth;
0442     K=2*delta_f*T-1;
0443     K=round(delta_f*2*T-1);
0444     if K<1,
0445         error('Your frequency resolution half-bandwidth is too small.  The minimum for this sized time window is: %g.', ...
0446             1/T);
0447     else
0448         delta_f=(K+1)/(2*T);
0449         fprintf('Desired half bandwidth (%f Hz) will be approximated by %f Hz.\n', ...
0450             p.Results.half_bandwidth,delta_f);
0451     end
0452     if K>=(T*EEG.srate)-1,
0453         error('Too many tapers.  The maximum number of possible tapers for this size time window and sampling rate is: %d.\nUse a smaller half bandwidth.', ...
0454             (T*EEG.srate)-2);
0455     end
0456 else
0457     K=1; 
0458     delta_f=(K+1)/(2*T); %delta_f=1/T (i.e., the Raleigh frequency)
0459     if VERBLEVEL>1
0460         fprintf('Frequency resolution half-bandwidth is: %f Hz\n',delta_f);
0461     end
0462 end
0463 TW=T*delta_f; %time-bandwidth product
0464 if K<=2,
0465     DropLastTaper=false; %this is required by pmtm.m, the function that computes the DFT
0466 else
0467     DropLastTaper=true;
0468 end
0469 
0470 %params.tapers=[T*delta_f K];
0471 %params.pad=p.Results.pad;
0472 %params.trialave=0; %because we're going to keep a running sum and then average ourselves (just in case a subject's data is split among multiple set files)
0473 if VERBLEVEL>1,
0474     fprintf('Number of tapers=%d\n',K);
0475 end
0476 
0477 
0478 %figure out how many time points are going to be fed into DFT
0479 n_dft_tpts=max(2^(nextpow2(n_use_tpts)+p.Results.pad),n_use_tpts);
0480 dft_T=n_dft_tpts/EEG.srate;
0481 freq_intrvl=1/dft_T;
0482 specGND.freqs=0:freq_intrvl:(EEG.srate/2);
0483 n_freqs=length(specGND.freqs);
0484 specGND.grands_pow_dB=zeros(n_chans,n_freqs,n_bins)*NaN;
0485 specGND.grands_pow_dB_stder=specGND.grands_pow_dB;
0486 specGND.grands_pow_dB_t=specGND.grands_pow_dB;
0487 specGND.dif=zeros(1,n_bins); %1 if a bin represents a difference between two conditions, 0 otherwise
0488 specGND.sub_ct=zeros(1,n_bins);
0489 specGND.chanlocs=EEG.chanlocs(use_chans);
0490 for b=1:n_bins,
0491     %Note, we're ignoring Kutaslab condition codes
0492     specGND.bindesc{b}=EEG.bindesc{use_bins(b)};
0493 end
0494 specGND.srate=EEG.srate;
0495 specGND.time_pts=EEG.times(use_tpts);
0496 orig_file_times=EEG.times; %for comparison with future files
0497 specGND.n_tapers=K;
0498 specGND.half_bandwidth=delta_f;
0499 specGND.indiv_fnames=infiles;
0500 if isempty(EEG.subject),
0501     if VERBLEVEL,
0502         fprintf('Set file, %s, does not have a subject name. I will use the filename as the subject''s name.\n', ...
0503             infiles{1});
0504     end
0505     specGND.indiv_subnames{1}=infiles{1};
0506 else
0507     specGND.indiv_subnames{1}=EEG.subject;
0508 end
0509 specGND.indiv_traits=[];
0510 specGND.indiv_bin_ct=zeros(n_infiles,n_bins);
0511 if isfield(EEG,'rawtrials_per_bin'),
0512     specGND.indiv_bin_raw_ct=zeros(n_infiles,n_bins); %This field might not be that useful, but just in case
0513 else
0514     specGND.indiv_bin_raw_ct=zeros(n_infiles,n_bins)*NaN; %This field is only relevant to Kutaslab data.  It keeps track of how many trials there were before artifact rejection.
0515 end
0516 specGND.indiv_pow_dB=zeros(n_chans,n_freqs,n_bins,n_infiles);
0517 specGND.indiv_art_ics=[];
0518 if isfield(EEG,'cal_info'),
0519     %This field is only relevant to Kutaslab data.
0520     n_cal_tpts=size(EEG.cal_info.erps,2);
0521     specGND.cals.indiv_cals=zeros(n_chans,n_cal_tpts);
0522     specGND.cals.indiv_cal_ct=zeros(1,n_infiles);
0523     specGND.cals.grand_cals=zeros(n_chans,n_cal_tpts);
0524     cal_info_present=1;
0525 else
0526     cal_info_present=0;
0527     specGND.cals=[];
0528 end
0529 specGND.history=[];
0530 specGND.t_tests=[];
0531 
0532 %% Loop through set files
0533 sub_ct=1;
0534 for filenum=1:n_infiles,
0535     new_sub=1; %Assume data from this subject has not been already loaded, until we learn otherwise (This is to deal with the fact that data from the same participant may be distributed among multiple set files)
0536     if filenum>1,
0537         global EEG; %this needs to be done to be compatible with EEGLAB (which stores EEG globally)
0538         if VERBLEVEL,
0539             fprintf('\n\n');
0540         end
0541         EEG=pop_loadset(infiles{filenum});
0542         EEG_var_check(EEG,infiles{filenum});
0543         EEG=epoch_if_necessary(EEG,p);
0544         
0545         %collect sub names
0546         if isempty(EEG.subject),
0547             if VERBLEVEL,
0548                 fprintf('Set file, %s, does not have a subject name. I will use the filename as the subject''s name.\n', ...
0549                     infiles{filenum});
0550             end
0551             sub_name=infiles{filenum};
0552         else
0553             sub_name=EEG.subject;
0554         end
0555         ur_sub_id=0;
0556         for old_sub=1:sub_ct,
0557             if strcmpi(specGND.indiv_subnames{old_sub},sub_name),
0558                 ur_sub_id=old_sub;
0559                 if VERBLEVEL
0560                     fprintf('Set file, %s, will be appended to data already loaded from Subject %s.\n', ...
0561                         infiles{filenum},specGND.indiv_subnames{old_sub});
0562                 end
0563                 new_sub=0;
0564                 break;
0565             end
0566         end
0567         if ~ur_sub_id,
0568             sub_ct=sub_ct+1;
0569             sub=sub_ct;
0570             specGND.indiv_subnames{sub}=sub_name;
0571         else
0572             sub=ur_sub_id;
0573         end
0574         
0575         EEG_var_check(EEG,infiles{filenum});
0576         
0577         [n_EEG_chans2, n_tpts2, n_epochs]=size(EEG.data);
0578         n_EEG_bins2=length(EEG.bindesc);
0579         
0580         %Figure out which channels to use
0581         if ~isempty(exclude_chans),
0582             use_chans=zeros(1,n_EEG_chans2);
0583             ex_ct=0;
0584             ex_labels=[];
0585             for a=1:n_EEG_chans2,
0586                 if ~ismember_ci(EEG.chanlocs(a).labels,exclude_chans)
0587                     use_chans(a)=1;
0588                 else
0589                     ex_ct=ex_ct+1;
0590                     ex_labels{ex_ct}=EEG.chanlocs(a).labels;
0591                 end
0592             end
0593             n_chans2=sum(use_chans);
0594             use_chans=find(use_chans==1);
0595             
0596             if VERBLEVEL>1,
0597                 missed=setdiff_ci(exclude_chans,ex_labels);
0598                 n_missed=length(missed);
0599                 if n_missed,
0600                     if n_missed==1,
0601                         msg=sprintf('I attempted to exclude the following channel, but it was not found:');
0602                     else
0603                         msg=sprintf('I attempted to exclude the following channels, but they were not found:');
0604                     end
0605                     for a=1:n_missed,
0606                         msg=[msg ' ' missed{a}];
0607                     end
0608                     watchit(msg);
0609                 end
0610             end
0611         elseif ~isempty(include_chans),
0612             use_chans=zeros(1,n_EEG_chans2);
0613             in_ct=0;
0614             in_labels=[];
0615             for a=1:n_EEG_chans2,
0616                 if ismember_ci(EEG.chanlocs(a).labels,include_chans)
0617                     use_chans(a)=1;
0618                     in_ct=in_ct+1;
0619                     in_labels{in_ct}=EEG.chanlocs(a).labels;
0620                 end
0621             end
0622             n_chans2=sum(use_chans);
0623             use_chans=find(use_chans==1);
0624             
0625             if VERBLEVEL>1,
0626                 missed=setdiff_ci(include_chans,in_labels);
0627                 n_missed=length(missed);
0628                 if n_missed,
0629                     if n_missed==1,
0630                         msg=sprintf('I attempted to include the following channel, but it was not found:');
0631                     else
0632                         msg=sprintf('I attempted to include the following channels, but they were not found:');
0633                     end
0634                     for a=1:n_missed,
0635                         msg=[msg ' ' missed{a}];
0636                     end
0637                     watchit(msg);
0638                 end
0639             end
0640         else
0641             n_chans2=n_EEG_chans2;
0642             use_chans=1:n_chans2;
0643         end
0644         
0645         %Check for consistency with specGND
0646         if EEG.srate~=specGND.srate,
0647             error('File %s has a different sampling rate than file %s.\n', ...
0648                 infiles{1},infiles{filenum});
0649         end
0650         if n_chans2~=n_chans,
0651             error('File %s has a different number of channels to import than file %s.\n', ...
0652                 infiles{1},infiles{filenum});
0653         end
0654         if n_tpts2~=n_tpts,
0655             error('File %s has a different number of time points than file %s.\n', ...
0656                 infiles{1},infiles{filenum});
0657         end
0658         if ~min(EEG.times==orig_file_times),
0659             error('The epochs in file %s begin and end at different times than file %s.\n', ...
0660                 infiles{1},infiles{filenum});
0661         end
0662         if isempty(p.Results.use_bins),
0663             if n_EEG_bins2~=n_bins,
0664                 error('File %s has a different number of bins than file %s.\n', ...
0665                     infiles{1},infiles{filenum});
0666             end
0667         else
0668             if max(p.Results.use_bins)>n_EEG_bins2,
0669                 error('File %s only has %d bins, but you''ve requested to import Bin %d.\n', ...
0670                     infiles{filenum},n_EEG_bins2,max(p.Results.use_bins));
0671             end
0672         end
0673         
0674         bin_ct=0;
0675         for b=use_bins,
0676             bin_ct=bin_ct+1;
0677             if ~strcmpi(specGND.bindesc{bin_ct},EEG.bindesc{b}),
0678                 error('The #%d imported bin in file %s is different than that of file %s.\n', ...
0679                     b,infiles{1},infiles{filenum});
0680             end
0681         end
0682         try
0683             [fs1, fs2, er]=comp_struct_quiet(specGND.chanlocs,EEG.chanlocs(use_chans));
0684         catch
0685             error('File %s''s imported channel location information differs from that of file %s.\n', ...
0686                 infiles{1},infiles{filenum});
0687         end
0688         if ~isempty(er),
0689             error('File %s''s imported channel location information differs from that of file %s.\n', ...
0690                 infiles{1},infiles{filenum});
0691         end
0692     else
0693         sub=sub_ct; %for first set file, sub=1 (since sub_ct=1)
0694     end
0695     
0696     %Remove ICs labeled as artifacts
0697     fldnames=fieldnames(EEG);
0698     art_ics=[];
0699     ics_removed=0;
0700     if sum(EEG.reject.gcompreject),
0701         ics_removed=1;
0702     else
0703         for fn=1:length(fldnames),
0704             if strcmpi(fldnames{fn},'iclabels'),
0705                 ics_removed=1;
0706                 break;
0707             end
0708         end
0709     end
0710     if ics_removed,
0711         if strcmpi(p.Results.rm_mean,'yes'),
0712             art_ics=remove_artifact_ics([EEG.times(use_tpts(1)) EEG.times(use_tpts(end))]);
0713         else
0714             art_ics=remove_artifact_ics([]);
0715         end
0716     end
0717     specGND.indiv_art_ics{filenum}=art_ics;
0718     
0719     %remove epoch means to dampen drift
0720     if strcmpi(p.Results.rm_mean,'yes'),
0721         if ~ics_removed,
0722             if VERBLEVEL>=2,
0723                 fprintf('Baselining data by removing mean of selected time window (%d to %d).\n', ...
0724                     round(EEG.times(use_tpts(1))),round(EEG.times(use_tpts(end))));
0725             end
0726             %Note, you could probably speed this up by just baselining epochs of
0727             %interest ??
0728             EEG.data=rmbase(EEG.data,n_tpts,use_tpts);
0729             EEG.data=reshape(EEG.data,n_EEG_chans,n_tpts,n_epochs);
0730         end
0731     end
0732      
0733     epoch_in_bin=zeros(n_bins,n_epochs); %could replace n_bins with n_use_bins to save memory
0734     %Collect epoch ids for all trials that fall in a particular bin
0735     for a=1:n_epochs,
0736         for b=1:length(EEG.epoch(a).eventtype),
0737             if (length(EEG.epoch(a).eventtype{b})>2) && strcmpi(EEG.epoch(a).eventtype{b}(1:3),'bin')
0738                 bin=str2num(EEG.epoch(a).eventtype{b}(4:end));
0739                 if ismember(bin,use_bins) && ~index_cellarray_or_vector(EEG.epoch(a).eventlatency,b,VERBLEVEL),
0740                     %only consider event types with latency 0
0741                  
0742                     binind=find(use_bins==bin); %new bin index (in case only a subset of bins selected
0743                     
0744                     new_event=1;
0745                     urevent=index_cellarray_or_vector(EEG.epoch(a).eventurevent,b,VERBLEVEL);
0746                     if ~new_sub,
0747                         %check to see if event has already been used for
0748                         %this bin
0749                         if ismember(urevent,used_urevents{sub,binind}), % TPU
0750                             if VERBLEVEL>1,
0751                                 fprintf('Epoch %d has already been added to Bin %d for Subject %s.\n', ...
0752                                     a,bin,specGND.indiv_subnames{sub});
0753                             end
0754                             new_event=0;
0755                         end
0756                     end
0757                     if new_event,
0758                         epoch_in_bin(binind,a)=1;
0759                         used_urevents{sub,binind}=[used_urevents{sub,binind} urevent];
0760                     end
0761                 end
0762             end
0763         end
0764     end
0765     
0766     for b=use_bins,
0767         binind=find(use_bins==b);  % index to the vector of used bins
0768         bin_ids=find(epoch_in_bin(binind,:));
0769         specGND.indiv_bin_ct(sub,binind)=specGND.indiv_bin_ct(sub,binind)+length(bin_ids);
0770         chan_ct=0;
0771         for c=use_chans,
0772             chan_ct=chan_ct+1;
0773             %perform DFT on epochs that fall in that bin and add them to a
0774             %running sum
0775             %if you ever want to use another taper, uncomment the code
0776             %below
0777             %[S,f]=pwr_spec(squeeze(EEG.data(c,use_tpts,bin_ids)),params.Fs,'hann',params.pad,1);
0778             %specGND.indiv_pow_dB(c,:,binind,sub)=specGND.indiv_pow_dB(c,:,binind,sub)+sum(10*log10(S),2)';
0779             
0780             ep_ct=0;
0781             S=zeros(n_freqs,length(bin_ids));
0782             for ep=bin_ids,
0783                 ep_ct=ep_ct+1;
0784                 [S(:,ep_ct), f]=pmtm(squeeze(EEG.data(c,use_tpts,bin_ids(ep_ct))),TW,n_dft_tpts,EEG.srate,'DropLastTaper',DropLastTaper,'unity');        
0785             end
0786             %Chronux based code (note, Chronux power values are half of the
0787             %pmtm values save for the DC component.  This is because pmtm
0788             %doubles the power of frequencies that are redundant (due the
0789             %symmetric nature of the DFT of a real valued time series)
0790             %[S2,f2]=mtspectrumc(squeeze(EEG.data(c,use_tpts,bin_ids)),params); %time x trials
0791             specGND.indiv_pow_dB(chan_ct,:,binind,sub)=specGND.indiv_pow_dB(chan_ct,:,binind,sub)+sum(10*log10(S),2)';
0792         end
0793         if (binind==1) && (chan_ct==1),
0794             specGND.freqs=f;
0795         end
0796         % Code for error checking (make sure we're getting the right trials)
0797         % dat=rmbase(EEG.data(c,use_tpts,bin_ids),256,1:25);
0798         % indiv_erps(c,:,binind,sub)=mean(dat,3);
0799     end
0800 
0801     if isfield(EEG,'rawtrials_per_bin'),
0802         %Again, this field is only relevant to Kutaslab data.  It keeps
0803         %track of how many trials there were before artifact rejection.
0804         specGND.indiv_bin_raw_ct(sub,:)=specGND.indiv_bin_raw_ct(sub,:)+EEG.rawtrials_per_bin(use_bins);
0805     end
0806     
0807     if isfield(EEG,'cal_info')
0808         %Again, this field is only relevant to Kutaslab data.
0809         %Note, I assume that even if a single participant's data are
0810         %distributed among multiple set files, EEG.cal_info is the same in
0811         %all the set files.
0812         specGND.cals.indiv_cals(:,:,sub)=EEG.cal_info.erps(use_chans,:);
0813         specGND.cals.indiv_cal_ct(sub)=EEG.cal_info.npulse_used;
0814         specGND.cals.caldesc='cal pulses';
0815         cal_info_present=1;
0816     end
0817     
0818     clear global EEG; %make sure it's not saved since ICs might have been removed
0819 end %filenum loop
0820 
0821 
0822 %remove unused elements of specGND variable (if any), because there are fewer
0823 %subjects than set files
0824 specGND.indiv_bin_ct=specGND.indiv_bin_ct(1:sub_ct,:);
0825 specGND.indiv_bin_raw_ct=specGND.indiv_bin_raw_ct(1:sub_ct,:);
0826 specGND.indiv_pow_dB=specGND.indiv_pow_dB(:,:,:,1:sub_ct);
0827 if cal_info_present,
0828     specGND.cals.indiv_cals(:,:,1:sub_ct)=specGND.cals.indiv_cals(:,:,1:sub_ct);
0829     specGND.cals.indiv_cal_ct=specGND.cals.indiv_cal_ct(1:sub_ct);
0830 end
0831 
0832 %turn running sums into means
0833 for sub=1:sub_ct,
0834     if VERBLEVEL>1,
0835         fprintf('\nTrials per bin for Subject %s:\n',specGND.indiv_subnames{sub});
0836     end
0837     for b=1:n_bins,
0838         bin_ct=specGND.indiv_bin_ct(sub,b);
0839         if bin_ct,
0840             if VERBLEVEL>1,
0841                 fprintf('Bin %d (%s): %d trials\n',b,specGND.bindesc{b},bin_ct);
0842             end
0843             specGND.indiv_pow_dB(:,:,b,sub)=specGND.indiv_pow_dB(:,:,b,sub)/bin_ct;
0844             specGND.sub_ct(b)=specGND.sub_ct(b)+1;
0845         else
0846             watchit(sprintf('Subject %s has no epochs that fall into bin %d.',specGND.indiv_subnames{sub},b));
0847             specGND.indiv_pow_dB(:,:,b,sub)=specGND.indiv_pow_dB(:,:,b,sub)*NaN;
0848         end
0849     end
0850 end
0851 
0852 %% Compute grands
0853 if VERBLEVEL>1,
0854     fprintf('\n\n'); %add a couple lines between between counts and baselining info
0855 end
0856 
0857 %Compute grand average power spectra
0858 for b=1:n_bins,
0859     bin_subs=find(specGND.indiv_bin_ct(:,b));
0860     if specGND.sub_ct(b),        
0861         specGND.grands_pow_dB(:,:,b)=mean(specGND.indiv_pow_dB(:,:,b,bin_subs),4);
0862         specGND.grands_pow_dB_stder(:,:,b)=std(specGND.indiv_pow_dB(:,:,b,bin_subs),0,4)/sqrt(specGND.sub_ct(b));
0863         specGND.grands_pow_dB_t(:,:,b)=specGND.grands_pow_dB(:,:,b)./specGND.grands_pow_dB_stder(:,:,b);
0864     else
0865         watchit(sprintf('No average files contribute to bin %d.',b));
0866     end
0867 end
0868 
0869 %Compute grand average cal pulses (if individual subject cal pulses present)
0870 if cal_info_present,
0871     specGND.cals.grand_cals=mean(specGND.cals.indiv_cals,3);
0872 end
0873 
0874 %% Save specGND variable
0875 if isempty(p.Results.out_fname),
0876     %Create GUI
0877     [jname, jpath]=uiputfile({'*.specGND','*.specGND files'; ...
0878         '*','All files'},'Save specGND variable as:','untitled.specGND');
0879     if ~jpath,
0880         fprintf('Output filename selection cancelled.  specGND variable NOT saved to disk.\n');
0881     else
0882         %test to make sure file can be created
0883         isW=isWriteable([jpath jname]);
0884         if isW,
0885             specGND=save_matmk(specGND,jname,jpath,1); % 1 means that user won't be asked again about saving file
0886         else
0887             fprintf('specGND file could not be saved, but should still exist in MATLAB workspace.');
0888         end
0889     end
0890 elseif ~strcmpi(p.Results.out_fname,'no save'),
0891     [jpath, jname]=pathNname(p.Results.out_fname);
0892     %Add .specGND extension if no extension given
0893     if ~ismember('.',jname),
0894         jname=[jname '.specGND'];
0895     end
0896     specGND=save_matmk(specGND,jname,jpath);
0897 end
0898 
0899 %%%%%%%% END OF MAIN FUNCTION %%%%%%%
0900 
0901 
0902 function EEG_var_check(EEG,fname)
0903 
0904 if ~isfield(EEG,'bindesc'),
0905     error('The variable stored in file %s does not contain a bin descriptor. You need to apply bin_info2EEG.m to it.',fname);
0906 end
0907 
0908 function value=index_cellarray_or_vector(cellarray_or_vector,index,VERBLEVEL)
0909 
0910 if index>length(cellarray_or_vector),
0911     %This if statement is necessary for compatibility with Kutaslab .set
0912     %files, which can have multiple types per epoch but typically only a
0913     %single latency value
0914     index=length(cellarray_or_vector);
0915     if VERBLEVEL>2,
0916         fprintf('EEG.epoch(#).eventtype index exceeds EEG.epoch(#).latency index.  I will simply use the last element of EEG.epoch(#).latency.\n');
0917     end
0918 end
0919 
0920 if iscell(cellarray_or_vector),
0921     value=cellarray_or_vector{index};
0922 else
0923     value=cellarray_or_vector(index);
0924 end
0925 
0926 
0927 function EEG=epoch_if_necessary(EEG,p)
0928 %Note, when epoching continuous data using EEGLAB's pop_epoch.m function
0929 %you need to include at least one data point before the event of interest
0930 %to ensure that information about the event of interest is included in
0931 %the EEG.epoch field.  When you compute the DFT you can extract only the
0932 %post-event window of interest.
0933 
0934 %Are data epoched or continuous?
0935 if length(size(EEG.data))==2
0936     %epoch the data
0937     time_wind=p.Results.time_wind/1000; %p.Results.time wind should be in milliseconds
0938     if isempty(time_wind),
0939        time_wind=[-.1 1.5]; %default
0940     elseif time_wind(1)>=0,
0941        time_wind(1)=-2/EEG.srate; %start at two samples before 0 to ensure that events of interest are included in the epoch
0942     end
0943     use_bins=p.Results.use_bins;
0944     if isempty(use_bins),
0945         %go through EEG.event and collect all the types of bins
0946         n_ev=length(EEG.event);
0947         for a=1:n_ev,
0948             if length(EEG.event(a).type)>3,
0949                 if strcmpi(EEG.event(a).type(1:3),'bin'),
0950                     ev_bin=str2double(EEG.event(a).type(4:end));
0951                     if ~ismember(ev_bin,use_bins)
0952                         use_bins=[use_bins ev_bin];
0953                     end
0954                 end
0955             end
0956         end
0957     end
0958     n_use_bins=length(use_bins);
0959     if ~n_use_bins,
0960         error('This set file does not contain any events of type ''bin#''.  You need to run bin_info2EEG.m on it.');
0961     end
0962     use_types=cell(1,n_use_bins);
0963     for a=1:n_use_bins,
0964         use_types{a}=['bin' int2str(use_bins(a))];
0965     end
0966     EEG=pop_epoch(EEG,use_types,time_wind, 'newname',[EEG.setname '_epoched'], 'epochinfo', 'yes');
0967 end
0968 
0969 
0970 function isW = isWriteable(inFile)
0971 % Function checks to make sure "inFile" can be written to disk. "inFile" is
0972 % not modified by the function.
0973 %
0974 % Author: Tom Urbach
0975 %
0976 
0977 isW = 0;
0978 [fid, message] = fopen(inFile,'a+');
0979 if (fid ~= -1)
0980     isW = 1;
0981     fclose(fid);
0982 else
0983     fprintf('Error opening %s: %s\n', inFile, message);
0984 end
0985 
0986 
0987 
0988 function yesno=ismember_ci(str,str_array)
0989 % function yesno=ismember_ci(str,str_array)
0990 % A case insensitive version of ismember.m but just for strings
0991 %
0992 % Inputs:
0993 %  str       - a string
0994 %  str_array - a cell array of strings
0995 %
0996 % Outputs:
0997 %  yesno - 1 if str is a member of str_array.  0 otherwise.  Comparison is
0998 %          not case sensitive.
0999 
1000 yesno=0;
1001 n_str=length(str_array);
1002 
1003 for a=1:n_str,
1004     if strcmpi(str,str_array{a})
1005         yesno=1;
1006         break;
1007     end
1008 end
1009 
1010 
1011 
1012 function dif_str=setdiff_ci(superset,subset)
1013 %function dif_str=setdiff_ci(superset,subset)
1014 %
1015 % Inputs:
1016 %   superset - an cell array of strings
1017 %   subset   - an cell array of strings
1018 %
1019 % Outputs:
1020 %   dif_str - a cell array of the strings that are in superset but NOT
1021 %             subset. Comparison is not case sensitive.
1022 %
1023 
1024 n_super=length(superset);
1025 n_sub=length(subset);
1026 dif_ct=0;
1027 dif_str=[];
1028 for a=1:n_super,
1029     found=0;
1030     for b=1:n_sub,
1031         if strcmpi(superset{a},subset{b}),
1032             found=1;
1033             break
1034         end
1035     end
1036     if ~found,
1037         dif_ct=dif_ct+1;
1038         dif_str{dif_ct}=superset{a};
1039     end
1040 end