Home > database > physio > ensemble_physio_summ.m

ensemble_physio_summ

PURPOSE ^

calculates summary statistics for physio data

SYNOPSIS ^

function outdata = ensemble_physio_summ(indata,defs)

DESCRIPTION ^

 calculates summary statistics for physio data
 
 This function calculates a variety of summary statistics for epochs of
 physiological data signals, matches those statistics with responses from
 presentation logfiles (optional), saves these data to disk (signals and
 peaks) in a form amenable to import into SAS (optional), and also saves
 waveform graphics that may be used in assessing the validity of the
 summary statistics that have been calculated, as well as the processing
 steps taken during those calculations.
 
 This script has the option of either completing everything in one step,
 or using a two-step process whereby epochs and peaks are extracted in the
 first step, and then a second step takes the epoched and peaked data and
 calculates summary statistics from them. This allows manual intervention
 between steps to visually inspect and remove bad peaks or waveform
 artifacts from the data.
 
 NOTE: this script assumes that pre-processed physio data have been saved
 into eeglab .set files, and that presentation log files are available to
 guide signal epoching
 
 REQUIRES
 
   indata - cell array of ensemble data structs
       'sinfo' - data struct, containing sinfo data for all (but only)
           those subjects that you wish to analyze
       'paths' - ensemble data struct, output of ensemble_fmri_init,
           containing 'physio_outdir', which is used as the output
           location of all figures that are constructed
       'physio_paths' - output of ensemble_physio, ensemble data
           struct containing one row per subject/session/run, with paths
           to pre-processed physio data files for each run
       'pres_paths' - output of ensemble_parse_present, an ensemble
           data struct containing one row per subject/session/run, with
           paths to pre-processed presentation response and event data
 
 LOADING PHYSIO DATA FILES
   
   physio files are loaded on a per-subject, per-session, per-run basis,
   from locations identified in 'physio_paths', for those subjects that
   exist in 'sinfo'. By default, ensemble_physio (which generates
   physio_paths) generates a non-filtered raw data file for each
   subject/session/run, and saves these to physio_paths, and given
   filtering criteria, also saves a filtered data file to disk in the same
   location as the raw data file, and saves the path to the filtered data
   file into physio_paths, with '_filtered' appended to the end of the
   filename, before the extension. If the following parameter is set to 1,
   ensemble_physio_summ will search for and attempt to load the
   filtered data file as opposed to the raw data file. If set to 0, it
   will attempt to load the raw, unfiltered data file.
 
   defs.physio_summ.use_filtered (default 1)
 
 EPOCHING
 
   ensemble_physio_summ is intended to calculate summary statistics
   on epochs of physiological data signals. These epochs may, or may not,
   have responses or experimental classes associated with them. Epochs can
   currently be defined in terms of events from a presentation log file,
   or vectors of epoch onsets and offsets in seconds (calculated from the
   beginning of the provided preprocessed physiological signal)
 
   defs.physio_summ.epoch_start (default 0) - if a struct, with sub-field
       'filter', this is treated as a filter criteria for the included
       presentation data, to filter for events to be treated as epoch
       onsets. If a numeric vector, treated as a vector of epochs onsets
       in seconds. If zero, no epoching will be conducted.
   defs.physio_summ.epoch_end (default 0) - if single-length numeric, this
       is treated as the epoch length in seconds. If a numeric vector of
       the same length as epoch_start or the onset vector returned by
       epoch_start, this is treated as a vector of epoch offsets in
       seconds. If a struct, with sub-field 'filter', this is treated as
       filter criteria for the included presentation data, to filter for
       events to be treated as epoch offsets. By default, epoch_starts
       will be matched with epoch_ends, and epoching will only continue if
       an equal number of starts and ends is found, with all matching ends
       occuring after their supposed matched start. By default, the max
       epoch length will be used if a vector is found as the result of a
       filter. If a filter criterion is used and the field 'minimum_end'
       exists and is set to = 1, the minimum epoch time (offset-onset)
       will be used for all epochs. If this variable is set to 0 and
       defs.physio_summ.epoch_start is not set to 0, epoch lengths will be
       set to the maximum interval between epoch starts. NOTE: this will
       cause overlapping signal between some epochs where there is a
       variance in the inter-start interval. If this is a problem, you
       should set 'minimum_end' = 1, so that the epoch length gets set to
       the minimum inter-start interval.
   defs.physio_summ.baseline (default 2) - pre-epoch baseline to be added to
       epoch when calculating and extracting epochs. this value is assumed
       to be in 'seconds'.
 
       FIXME: In the case that there is not sufficient room
       for the given baseline period between the first epoch onset and the
       very beginning of the signal acquisition (for instance, first epoch
       onset is at 1s, but the baseline period is 2s), the first epoch
       will be not be returned by eeglab. This is because eeglab
       requires an equal number of datapoints for all epochs extracted
       from a given signal, and this case creates an out-of-bounds epoch
       definition that is dropped by eeglab.
 
       FIXME: pad the signal with 0s to provide enough data points to
       allow for a baseline for the first epoch, and then adjust timings
       of all other events? this would necessarily require special
       calculation of the baseline of the first epoch, limited to
       real-data only, when subtracting the baseline from the signal of
       interest, so as not to bias the baseline with artifically added 0s.
 
 SUMMARY DATA CALCULATION
 
   The current analyses calculate summary statistics for both pulse and
   skin conductance data. For pulse, mean heart rate within an epoch, as
   well as heart rate variance, is calculated. For scr, a number of
   statistics are calculated to investigate both individual transient skin
   conductance responses (SCRs), and change in overall skin conductance
   level (SCL).
 
   For SCRs, individual responses are defined as the highest peak within a
   range of +/- 'proxlim' seconds (default 0.6) within an epoch, whose
   height is greater than height_thresh (default 0.02). A standard minimum
   SCR height, as defined by Dawson et al 2000, is 0.02 microsiemens,
   therefore if the incoming signal is defined in microsiemens (as is
   biopac data when both high pass filters on the GSR100 are set to DC),
   the default height_thresh should be used. Total number of SCRs within
   an epoch (npeaks), mean peak height, peak height variance, first peak
   height, and first peak time are calculated. If preprocessed data are
   provided in microsiemens, then heights and variance will be returned in
   microsiemens.
 
   For SCLs, epoched skin conductance data are standardized
   within-subject, across epochs, to control for individual differences in
   SCL variability and range (Dawson et al 2000, p209; also Ben-Shakhar
   1985). Specific integral and mean SCL over the expected period of a
   specific SCR (spc_integral.being, spc_integral.end, default 2s to 6s)
   are calculated, as well as integral and mean SCL over the entire epoch
   (not including the baseline period).
 
   NOTE: if entering these data into mixed effects multi-level regression
   models, certain effects such as different rates of habituation between
   individuals will be somewhat controlled for by pooling variance within
   subjects, but if not, additional measures should be considered if you
   want to control for individual differences in (at least) rate of
   habituation when using repeated stimuli or chronic stimuli. Different
   rates of habituation may effect mean peak height and peak height
   variance, and may also have some effect on integral.
 
   defs.physio_summ.channels - channel names, as they appear in your
       eeglab chanlocs labels, that you wish to calculate summary
       information for.
   defs.physio_summ.scr - settings for SCR/GSR data summary
       .spc_integral.begin (default 2) - time, in seconds, to start
           calculation of the specific SCR integral
       .spc_integral.end (default 6) - time, in seconds, to end
           calculation of the specific SCR integral
       .find_peaks - the following parameters are sent to find_peaks:
           .thresh (default 0) - proportion of total signal
               range below which peaks will not be considered. for
               example, for .thresh = 0.5, peaks will only be returned if
               their height is more than half of the total signal range.
               This is passed on to find_peaks, but is always set to 0,
               since we want to be able to find peaks across the entire
               range of signal.
           .peak_heights.calcFrom (default 'left') - which side of a peak
               to find the trough used to calculate peak height
           .peak_heights.transformation (default 'none') - the
               transformation, if any, performed on peak height values
               before returning them.
           .peak_height_window (default 0.6) - sliding window in
               which only the highest peak is retained. for each peak
               found, if there are other peaks within +/- seconds of that
               peak, all but the highest peak will be discarded.
           .peak_height_thresh (default 0.02) - peak heights below
               this value will be rejected
       .keep_baseline (default: 1) - keep the baseline period in the
           final signal that is returned (1), or remove it (0)? either
           way, the baseline period (if greater than 0) mean is removed
           from the signal of interest for each epoch, and peaks during
           the baseline period are removed after peak finding
   defs.physio_summ.pulse - settings for pulse ox. data summary
       .find_peaks.thresh (default 0.5) - proportion of total signal
           range below which peaks will not be considered. for example,
           for .thresh = 0.5, peaks will only be returned if their height 
           is more than half of the total signal range
 
 RESPONSE MATCHING FOR EPOCHS
 
   if there are responses expected and present within a presentation
   logfile that can be matched in a 1 to 1 relationship with epochs that
   are extracted or defined, these responses can be included in any data
   files that are exported, matched to each epoch and the summary
   statistics calculated for each epoch. Also, overlays of all epochs
   corresponding to each unique response within a response set, and
   average waveform for each level of each response, can be calculated and
   output. presentation data must be included, in the form of a pres_paths
   data struct (created by presentation_parse) containing paths to
   presentation data that has been parsed and saved to disk for the given
   participant/session.
 
   defs.physio_summ.responses (default 0) - struct array, each struct
       identifying filtering criteria to be applied to the given
       presentation data, to extract events containing response data that
       should be mapped to each trial. This response data will be exported
       along with summary information calculated for each trial.
       .type - string, name of the given response. this string will be
           used as the variable name for the given response data in the
           resulting output data struct
       .event - struct, assumed to include either an .include or .exclude
           fieldname, and subsequent ensemble_filter criteria for the
           given response class
 
 DATA EXTRACTION VS SUMMARY STATISTIC CALCULATION
 
   the default is to extract epochs, extract peaks, and calculate summary
   statistics in one step. You can control the execution of each step with
   the following parameters:
 
   defs.physio_summ.EXTRACT_DATA (default: 1) - extract epochs and peaks
       from provided waveforms. These data are returned in 'gsr_epochs'
       and 'cardiac_epochs' datasets. REQUIRES: physio_paths and
       pres_paths input datasets
   defs.physio_summ.CALCULATE_STATS (default: 1) - calculate summary
       statistics on extracted epochs and peaks. REQUIRES: either
       defs.physio_summ.EXTRACT_DATA, physio_paths, and pres_paths, OR
       'gsr_epochs' and 'cardiac_epochs' input data structures
 
   to use the two step process, you can first execute this script with
   EXTRACT_DATA=1, CALCULATE_STATS=0, then import the results from that
   step into another script (possibly manual_adjust_signal_data.m) to
   manually edit the waveforms and peak information, and then return those
   edited datasets to this script with EXTRACT_DATA=0, CALCULATE_STATS=1
   to complete the summary statistic calculation.
 
 DATA OUTPUT
 
   defs.physio_summ.GRAPH_EPOCHS (default 0)
   defs.physio_summ.EPOCHS_BY_RESP (default 0)
   defs.physio_summ.SAVE_EPOCHS (default 1)
   defs.physio_summ.SAVE_DATA (default 0)
   defs.physio_summ.export (default 0)
   defs.physio_summ.sas (default 0)
   defs.physio_summ.export_stub (default '')
 
 FIXME - should be re-written, to include a mechanism by which channels to
 analyze are identified by EEG.chanlocs name as sub-fields of the params
 struct, and under that, sub-fields of each channel should be named
 exactly as the names of functions used to either filter or report the
 data. Then, sub-fields of those fields should be params that get sent to
 those functions as the second argument, the first argument being the
 data output of the previous function (in may or all cases, an eeglab
 struct? ... OR a cell array of structs containing a series of processing
 jobs, params for those jobs, function handles/references, etc ...?
 
 CITATIONS
   Dawson ME, Schell AM, & Filion DL (2000). "The Electrodermal System."
       In JT Cacioppo, LG Tassinary & GG Berntson (Eds.) "Handbook of
       Psychophysiology" (pp 201-223). Cambridge: Cambridge Univ. Press.
   Ben-Shakhar  (1985). "Standardization within individuals: Simple method
       to neutralize individual differences in skin conductance."
       Psychophysiology, 22:292-9.
 
 FB 2009.04.21
 FB 2009.05.04 - adding documentation, now zscores signal for SCL
 calculation but leaves signal for SCR calculation in raw form.
 FB 2009.06.18 - changed "params.lead" to "params.baseline", and checked
 the script for hard-coded reliance on anything fmri-related, with aims to
 generalize from ensemble_fmri_physio_summ to ensemble_physio_summ

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function outdata = ensemble_physio_summ(indata,defs)
0002 
0003 % calculates summary statistics for physio data
0004 %
0005 % This function calculates a variety of summary statistics for epochs of
0006 % physiological data signals, matches those statistics with responses from
0007 % presentation logfiles (optional), saves these data to disk (signals and
0008 % peaks) in a form amenable to import into SAS (optional), and also saves
0009 % waveform graphics that may be used in assessing the validity of the
0010 % summary statistics that have been calculated, as well as the processing
0011 % steps taken during those calculations.
0012 %
0013 % This script has the option of either completing everything in one step,
0014 % or using a two-step process whereby epochs and peaks are extracted in the
0015 % first step, and then a second step takes the epoched and peaked data and
0016 % calculates summary statistics from them. This allows manual intervention
0017 % between steps to visually inspect and remove bad peaks or waveform
0018 % artifacts from the data.
0019 %
0020 % NOTE: this script assumes that pre-processed physio data have been saved
0021 % into eeglab .set files, and that presentation log files are available to
0022 % guide signal epoching
0023 %
0024 % REQUIRES
0025 %
0026 %   indata - cell array of ensemble data structs
0027 %       'sinfo' - data struct, containing sinfo data for all (but only)
0028 %           those subjects that you wish to analyze
0029 %       'paths' - ensemble data struct, output of ensemble_fmri_init,
0030 %           containing 'physio_outdir', which is used as the output
0031 %           location of all figures that are constructed
0032 %       'physio_paths' - output of ensemble_physio, ensemble data
0033 %           struct containing one row per subject/session/run, with paths
0034 %           to pre-processed physio data files for each run
0035 %       'pres_paths' - output of ensemble_parse_present, an ensemble
0036 %           data struct containing one row per subject/session/run, with
0037 %           paths to pre-processed presentation response and event data
0038 %
0039 % LOADING PHYSIO DATA FILES
0040 %
0041 %   physio files are loaded on a per-subject, per-session, per-run basis,
0042 %   from locations identified in 'physio_paths', for those subjects that
0043 %   exist in 'sinfo'. By default, ensemble_physio (which generates
0044 %   physio_paths) generates a non-filtered raw data file for each
0045 %   subject/session/run, and saves these to physio_paths, and given
0046 %   filtering criteria, also saves a filtered data file to disk in the same
0047 %   location as the raw data file, and saves the path to the filtered data
0048 %   file into physio_paths, with '_filtered' appended to the end of the
0049 %   filename, before the extension. If the following parameter is set to 1,
0050 %   ensemble_physio_summ will search for and attempt to load the
0051 %   filtered data file as opposed to the raw data file. If set to 0, it
0052 %   will attempt to load the raw, unfiltered data file.
0053 %
0054 %   defs.physio_summ.use_filtered (default 1)
0055 %
0056 % EPOCHING
0057 %
0058 %   ensemble_physio_summ is intended to calculate summary statistics
0059 %   on epochs of physiological data signals. These epochs may, or may not,
0060 %   have responses or experimental classes associated with them. Epochs can
0061 %   currently be defined in terms of events from a presentation log file,
0062 %   or vectors of epoch onsets and offsets in seconds (calculated from the
0063 %   beginning of the provided preprocessed physiological signal)
0064 %
0065 %   defs.physio_summ.epoch_start (default 0) - if a struct, with sub-field
0066 %       'filter', this is treated as a filter criteria for the included
0067 %       presentation data, to filter for events to be treated as epoch
0068 %       onsets. If a numeric vector, treated as a vector of epochs onsets
0069 %       in seconds. If zero, no epoching will be conducted.
0070 %   defs.physio_summ.epoch_end (default 0) - if single-length numeric, this
0071 %       is treated as the epoch length in seconds. If a numeric vector of
0072 %       the same length as epoch_start or the onset vector returned by
0073 %       epoch_start, this is treated as a vector of epoch offsets in
0074 %       seconds. If a struct, with sub-field 'filter', this is treated as
0075 %       filter criteria for the included presentation data, to filter for
0076 %       events to be treated as epoch offsets. By default, epoch_starts
0077 %       will be matched with epoch_ends, and epoching will only continue if
0078 %       an equal number of starts and ends is found, with all matching ends
0079 %       occuring after their supposed matched start. By default, the max
0080 %       epoch length will be used if a vector is found as the result of a
0081 %       filter. If a filter criterion is used and the field 'minimum_end'
0082 %       exists and is set to = 1, the minimum epoch time (offset-onset)
0083 %       will be used for all epochs. If this variable is set to 0 and
0084 %       defs.physio_summ.epoch_start is not set to 0, epoch lengths will be
0085 %       set to the maximum interval between epoch starts. NOTE: this will
0086 %       cause overlapping signal between some epochs where there is a
0087 %       variance in the inter-start interval. If this is a problem, you
0088 %       should set 'minimum_end' = 1, so that the epoch length gets set to
0089 %       the minimum inter-start interval.
0090 %   defs.physio_summ.baseline (default 2) - pre-epoch baseline to be added to
0091 %       epoch when calculating and extracting epochs. this value is assumed
0092 %       to be in 'seconds'.
0093 %
0094 %       FIXME: In the case that there is not sufficient room
0095 %       for the given baseline period between the first epoch onset and the
0096 %       very beginning of the signal acquisition (for instance, first epoch
0097 %       onset is at 1s, but the baseline period is 2s), the first epoch
0098 %       will be not be returned by eeglab. This is because eeglab
0099 %       requires an equal number of datapoints for all epochs extracted
0100 %       from a given signal, and this case creates an out-of-bounds epoch
0101 %       definition that is dropped by eeglab.
0102 %
0103 %       FIXME: pad the signal with 0s to provide enough data points to
0104 %       allow for a baseline for the first epoch, and then adjust timings
0105 %       of all other events? this would necessarily require special
0106 %       calculation of the baseline of the first epoch, limited to
0107 %       real-data only, when subtracting the baseline from the signal of
0108 %       interest, so as not to bias the baseline with artifically added 0s.
0109 %
0110 % SUMMARY DATA CALCULATION
0111 %
0112 %   The current analyses calculate summary statistics for both pulse and
0113 %   skin conductance data. For pulse, mean heart rate within an epoch, as
0114 %   well as heart rate variance, is calculated. For scr, a number of
0115 %   statistics are calculated to investigate both individual transient skin
0116 %   conductance responses (SCRs), and change in overall skin conductance
0117 %   level (SCL).
0118 %
0119 %   For SCRs, individual responses are defined as the highest peak within a
0120 %   range of +/- 'proxlim' seconds (default 0.6) within an epoch, whose
0121 %   height is greater than height_thresh (default 0.02). A standard minimum
0122 %   SCR height, as defined by Dawson et al 2000, is 0.02 microsiemens,
0123 %   therefore if the incoming signal is defined in microsiemens (as is
0124 %   biopac data when both high pass filters on the GSR100 are set to DC),
0125 %   the default height_thresh should be used. Total number of SCRs within
0126 %   an epoch (npeaks), mean peak height, peak height variance, first peak
0127 %   height, and first peak time are calculated. If preprocessed data are
0128 %   provided in microsiemens, then heights and variance will be returned in
0129 %   microsiemens.
0130 %
0131 %   For SCLs, epoched skin conductance data are standardized
0132 %   within-subject, across epochs, to control for individual differences in
0133 %   SCL variability and range (Dawson et al 2000, p209; also Ben-Shakhar
0134 %   1985). Specific integral and mean SCL over the expected period of a
0135 %   specific SCR (spc_integral.being, spc_integral.end, default 2s to 6s)
0136 %   are calculated, as well as integral and mean SCL over the entire epoch
0137 %   (not including the baseline period).
0138 %
0139 %   NOTE: if entering these data into mixed effects multi-level regression
0140 %   models, certain effects such as different rates of habituation between
0141 %   individuals will be somewhat controlled for by pooling variance within
0142 %   subjects, but if not, additional measures should be considered if you
0143 %   want to control for individual differences in (at least) rate of
0144 %   habituation when using repeated stimuli or chronic stimuli. Different
0145 %   rates of habituation may effect mean peak height and peak height
0146 %   variance, and may also have some effect on integral.
0147 %
0148 %   defs.physio_summ.channels - channel names, as they appear in your
0149 %       eeglab chanlocs labels, that you wish to calculate summary
0150 %       information for.
0151 %   defs.physio_summ.scr - settings for SCR/GSR data summary
0152 %       .spc_integral.begin (default 2) - time, in seconds, to start
0153 %           calculation of the specific SCR integral
0154 %       .spc_integral.end (default 6) - time, in seconds, to end
0155 %           calculation of the specific SCR integral
0156 %       .find_peaks - the following parameters are sent to find_peaks:
0157 %           .thresh (default 0) - proportion of total signal
0158 %               range below which peaks will not be considered. for
0159 %               example, for .thresh = 0.5, peaks will only be returned if
0160 %               their height is more than half of the total signal range.
0161 %               This is passed on to find_peaks, but is always set to 0,
0162 %               since we want to be able to find peaks across the entire
0163 %               range of signal.
0164 %           .peak_heights.calcFrom (default 'left') - which side of a peak
0165 %               to find the trough used to calculate peak height
0166 %           .peak_heights.transformation (default 'none') - the
0167 %               transformation, if any, performed on peak height values
0168 %               before returning them.
0169 %           .peak_height_window (default 0.6) - sliding window in
0170 %               which only the highest peak is retained. for each peak
0171 %               found, if there are other peaks within +/- seconds of that
0172 %               peak, all but the highest peak will be discarded.
0173 %           .peak_height_thresh (default 0.02) - peak heights below
0174 %               this value will be rejected
0175 %       .keep_baseline (default: 1) - keep the baseline period in the
0176 %           final signal that is returned (1), or remove it (0)? either
0177 %           way, the baseline period (if greater than 0) mean is removed
0178 %           from the signal of interest for each epoch, and peaks during
0179 %           the baseline period are removed after peak finding
0180 %   defs.physio_summ.pulse - settings for pulse ox. data summary
0181 %       .find_peaks.thresh (default 0.5) - proportion of total signal
0182 %           range below which peaks will not be considered. for example,
0183 %           for .thresh = 0.5, peaks will only be returned if their height
0184 %           is more than half of the total signal range
0185 %
0186 % RESPONSE MATCHING FOR EPOCHS
0187 %
0188 %   if there are responses expected and present within a presentation
0189 %   logfile that can be matched in a 1 to 1 relationship with epochs that
0190 %   are extracted or defined, these responses can be included in any data
0191 %   files that are exported, matched to each epoch and the summary
0192 %   statistics calculated for each epoch. Also, overlays of all epochs
0193 %   corresponding to each unique response within a response set, and
0194 %   average waveform for each level of each response, can be calculated and
0195 %   output. presentation data must be included, in the form of a pres_paths
0196 %   data struct (created by presentation_parse) containing paths to
0197 %   presentation data that has been parsed and saved to disk for the given
0198 %   participant/session.
0199 %
0200 %   defs.physio_summ.responses (default 0) - struct array, each struct
0201 %       identifying filtering criteria to be applied to the given
0202 %       presentation data, to extract events containing response data that
0203 %       should be mapped to each trial. This response data will be exported
0204 %       along with summary information calculated for each trial.
0205 %       .type - string, name of the given response. this string will be
0206 %           used as the variable name for the given response data in the
0207 %           resulting output data struct
0208 %       .event - struct, assumed to include either an .include or .exclude
0209 %           fieldname, and subsequent ensemble_filter criteria for the
0210 %           given response class
0211 %
0212 % DATA EXTRACTION VS SUMMARY STATISTIC CALCULATION
0213 %
0214 %   the default is to extract epochs, extract peaks, and calculate summary
0215 %   statistics in one step. You can control the execution of each step with
0216 %   the following parameters:
0217 %
0218 %   defs.physio_summ.EXTRACT_DATA (default: 1) - extract epochs and peaks
0219 %       from provided waveforms. These data are returned in 'gsr_epochs'
0220 %       and 'cardiac_epochs' datasets. REQUIRES: physio_paths and
0221 %       pres_paths input datasets
0222 %   defs.physio_summ.CALCULATE_STATS (default: 1) - calculate summary
0223 %       statistics on extracted epochs and peaks. REQUIRES: either
0224 %       defs.physio_summ.EXTRACT_DATA, physio_paths, and pres_paths, OR
0225 %       'gsr_epochs' and 'cardiac_epochs' input data structures
0226 %
0227 %   to use the two step process, you can first execute this script with
0228 %   EXTRACT_DATA=1, CALCULATE_STATS=0, then import the results from that
0229 %   step into another script (possibly manual_adjust_signal_data.m) to
0230 %   manually edit the waveforms and peak information, and then return those
0231 %   edited datasets to this script with EXTRACT_DATA=0, CALCULATE_STATS=1
0232 %   to complete the summary statistic calculation.
0233 %
0234 % DATA OUTPUT
0235 %
0236 %   defs.physio_summ.GRAPH_EPOCHS (default 0)
0237 %   defs.physio_summ.EPOCHS_BY_RESP (default 0)
0238 %   defs.physio_summ.SAVE_EPOCHS (default 1)
0239 %   defs.physio_summ.SAVE_DATA (default 0)
0240 %   defs.physio_summ.export (default 0)
0241 %   defs.physio_summ.sas (default 0)
0242 %   defs.physio_summ.export_stub (default '')
0243 %
0244 % FIXME - should be re-written, to include a mechanism by which channels to
0245 % analyze are identified by EEG.chanlocs name as sub-fields of the params
0246 % struct, and under that, sub-fields of each channel should be named
0247 % exactly as the names of functions used to either filter or report the
0248 % data. Then, sub-fields of those fields should be params that get sent to
0249 % those functions as the second argument, the first argument being the
0250 % data output of the previous function (in may or all cases, an eeglab
0251 % struct? ... OR a cell array of structs containing a series of processing
0252 % jobs, params for those jobs, function handles/references, etc ...?
0253 %
0254 % CITATIONS
0255 %   Dawson ME, Schell AM, & Filion DL (2000). "The Electrodermal System."
0256 %       In JT Cacioppo, LG Tassinary & GG Berntson (Eds.) "Handbook of
0257 %       Psychophysiology" (pp 201-223). Cambridge: Cambridge Univ. Press.
0258 %   Ben-Shakhar  (1985). "Standardization within individuals: Simple method
0259 %       to neutralize individual differences in skin conductance."
0260 %       Psychophysiology, 22:292-9.
0261 %
0262 % FB 2009.04.21
0263 % FB 2009.05.04 - adding documentation, now zscores signal for SCL
0264 % calculation but leaves signal for SCR calculation in raw form.
0265 % FB 2009.06.18 - changed "params.lead" to "params.baseline", and checked
0266 % the script for hard-coded reliance on anything fmri-related, with aims to
0267 % generalize from ensemble_fmri_physio_summ to ensemble_physio_summ
0268 
0269 outdata = ensemble_init_data_struct();
0270 
0271 global r
0272 
0273 r = init_results_struct;
0274 
0275 r.type = 'physio_summ';  % Identify the type of this reporting instance
0276 r.report_on_fly = 1;
0277 
0278 % % load colormap
0279 % load new_seismic;
0280 % colormap(new_seismic);
0281 
0282 % Parse out the input data
0283 for idata = 1:length(indata)
0284   if isfield(indata{idata},'type')
0285     switch indata{idata}.type
0286       case {'paths'}
0287         pathdata = indata{idata};
0288         pcol = set_var_col_const(pathdata.vars);
0289       case {'physio_paths'}
0290         physio_paths = indata{idata};
0291         physcol = set_var_col_const(physio_paths.vars);
0292       case {'pres_paths'}
0293         pres_paths = indata{idata};
0294         prescol = set_var_col_const(pres_paths.vars);
0295       case {'sinfo'}
0296         sinfo = indata{idata};
0297         sinfo = sinfo.data;
0298         proc_subs = {sinfo(:).id};
0299         nsub_proc = length(proc_subs);
0300       case {'gsr_epochs'}
0301         gsr_epochs = indata{idata};
0302         gecol = set_var_col_const(gsr_epochs.vars);
0303       case {'cardiac_epochs'}
0304         cardiac_epochs = indata{idata};
0305         cecol = set_var_col_const(cardiac_epochs.vars);
0306     end
0307   end
0308 end
0309 
0310 % check for required vars
0311 check_vars = {'sinfo','pathdata'};
0312 check_required_vars;
0313 
0314 if (iscell(indata) && ~isempty(indata) && isfield(indata{1},'task') && ...
0315         ~isempty(strmatch('return_outdir',indata{1}.task))) || ...
0316         (isstruct(indata) && isfield(indata,'task') && ...
0317         ~isempty(strmatch('return_outdir',indata.task)))
0318   if exist('pathdata','var') && ~isempty(pathdata.data{1})
0319     if exist('nsub_proc','var') && length(nsub_proc) == 1 && ...
0320             exist('proc_subs','var')
0321       pfilt = struct();
0322       pfilt.include.all.subject_id = proc_subs;
0323       lpathdata = ensemble_filter(pathdata,pfilt);
0324       if ~isempty(lpathdata.data{1})
0325         sfilt = pfilt;
0326         sfilt.include.all.path_type = {'physio_outdir'};
0327         spathdata = ensemble_filter(lpathdata,sfilt);
0328         if length(spathdata.data{1}) == 1
0329           % one epi outdir, save outdata = epi_outdir
0330           outdata = spathdata.data{pcol.path}{1};
0331         else
0332           sfilt = pfilt;
0333           sfilt.include.all.path_type = {'sess_outdir'};
0334           spathdata = ensemble_filter(lpathdata,sfilt);
0335           if length(spathdata.data{1}) == 1;
0336             outdata = spathdata.data{pcol.path}{1};
0337           else
0338             sfilt = pfilt;
0339             sfilt.include.all.path_type = {'sub_outdir'};
0340             spathdata = ensemble_filter(lpathdata,sfilt);
0341             if length(spathdata.data{1}) == 1;
0342               outdata = spathdata.data{pcol.path}{1};            
0343             end
0344           end
0345         end
0346       end
0347     end
0348   end
0349   if ~exist('outdata','var') || ~exist(outdata,'dir'), outdata = ''; end
0350   return
0351 end
0352 
0353 % initialize outdata struct for sinfo
0354 outdata.vars = [outdata.vars 'sinfo'];
0355 sinfo_idx = length(outdata.vars);
0356 outdata.data{sinfo_idx} = ensemble_init_data_struct();
0357 outdata.data{sinfo_idx}.type = 'sinfo';
0358 outdata.data{sinfo_idx}.data = sinfo;
0359 
0360 % check for physio_summ parameters within defs
0361 if isfield(defs,'physio_summ')
0362   ps = defs.physio_summ;
0363 else
0364   error('please provide params through defs.physio_summ\n');
0365 end
0366 
0367 % get parameters
0368 try chans_req = ps.channels; catch error('must provide channel names\n'); end
0369 try GRAPH_EPOCHS = ps.GRAPH_EPOCHS; catch GRAPH_EPOCHS = 0; end
0370 try EPOCH_BY_RESP = ps.EPOCH_BY_RESP; catch EPOCH_BY_RESP = 0; end
0371 try SAVE_EPOCHS = ps.SAVE_EPOCHS; catch SAVE_EPOCHS = 0; end
0372 try SAVE_DATA = ps.SAVE_DATA; catch SAVE_DATA = 0; end
0373 try EXTRACT_DATA = ps.EXTRACT_DATA; catch EXTRACT_DATA = 1; end
0374 try CALCULATE_STATS = ps.CALCULATE_STATS; catch CALCULATE_STATS = 1; end
0375 try export = ps.export; catch export.write2file = 0; end
0376 try sas = ps.sas; catch sas.libsave = 0; end
0377 try epoch_start = ps.epoch_start; catch epoch_start = ''; end
0378 try epoch_end = ps.epoch_end; catch epoch_end = ''; end
0379 try responses = ps.responses; catch responses = 0; end
0380 try baseline = ps.baseline; catch baseline = 2; end % in seconds
0381 try use_filtered = ps.use_filtered; catch use_filtered = 0; end
0382 try ep_per_fig = ps.ep_per_fig; catch ep_per_fig = 3; end
0383 try pargs = ps.pargs; catch pargs = {'-dpsc','-append'}; end
0384 try export_stub = ps.export_stub; catch export_stub = ''; end
0385 
0386 % check required datasets, based upon steps taken
0387 if EXTRACT_DATA
0388   check_vars = {'physio_paths','pres_paths'};
0389   check_required_vars;
0390 end
0391 if CALCULATE_STATS && ~EXTRACT_DATA
0392   check_vars = {{'gsr_epochs','cardiac_epochs'}};
0393   check_required_vars;
0394 end
0395 
0396 % get scr analysis parameters
0397 try scr.zscore = ps.scr.zscore; catch scr.zscore = ''; end
0398 try scr.int_begin = ps.scr.int_begin; catch scr.int_begin = 2; end % in sec
0399 try scr.int_end = ps.scr.int_end; catch scr.int_end = 6; end % in seconds
0400 try scr.pk.thresh = ps.scr.find_peaks.thresh; catch scr.pk.thresh = 0; end
0401 try scr.pk.peak_heights = ps.find_peaks.peak_heights;
0402   catch scr.pk.peak_heights.calcFrom = 'left';
0403         scr.pk.peak_heights.transformation = 'none';
0404 end
0405 try scr.pk.peak_height_window = ps.scr.find_peaks.peak_height_window;
0406   catch scr.pk.peak_height_window = 0.6; end % in seconds
0407 try scr.pk.peak_height_thresh = ps.scr.find_peaks.peak_height_thresh;
0408   catch scr.pk.peak_height_thresh = 0.02; end
0409 try keep_baseline = ps.scr.keep_baseline; catch keep_baseline = 1; end
0410   
0411 % verify integral beginning and ending points
0412 if scr.int_begin > scr.int_end
0413   error(['must provide valid integral beginning and ending points, in '...
0414       'seconds\n']);
0415 end
0416 
0417 % get pulse analysis parameters
0418 try pulse.pk = ps.pulse.find_peaks;
0419 catch pulse.pk.thresh = 0.5; pulse.pk.peak_height_window = 0.3; end
0420 
0421 % get types of responses to include in output data
0422 if isstruct(responses) && (isnumeric(epoch_start) && ...
0423         (length(epoch_start) < 1 || epoch_start)) || ...
0424         (isstruct(epoch_start) && isfield(epoch_start,'filter'))
0425   resp_names = {responses.type};
0426   nresp = length(resp_names);
0427 else
0428   resp_names = {};
0429   nresp = 0;
0430 end
0431 
0432 % check channels against what is supported by this script
0433 supported = {'scr','gsr','cardiac','pulse'};
0434 nchanreq = length(chans_req);
0435 channels = {};
0436 
0437 % weed out unsupported channels
0438 for ic=1:nchanreq
0439   if ~isempty(strmatch(chans_req{ic},supported))
0440     channels = [channels chans_req{ic}];
0441   else
0442     warning('channel %s not supported, not being used\n',chans_req{ic});
0443   end
0444 end
0445 
0446 % get number of supported channels that have been requested
0447 nchan = length(channels);
0448 if ~nchan, error('no supported channels provided'), end
0449 
0450 %eeglab('initpaths');
0451 % Add only the eeglab paths we need
0452 eeglab_root = fileparts(which('eeglab'));
0453 addpath(fullfile(eeglab_root,'functions/popfunc'));
0454 
0455 % init default vars used in extracted datasets
0456 xvars = {'subject_id','session','ensemble_id','run','trial',...
0457     'signal','srate','peakidxs'};
0458 
0459 if EXTRACT_DATA
0460   lvars = {xvars{:},resp_names{:}};
0461   
0462   if (~isempty(strmatch('gsr',channels)) || ...
0463           ~isempty(strmatch('scr',channels)))
0464     % initialize a gsr summary data output structure
0465     outdata.vars = [outdata.vars 'gsr_epochs'];
0466     gsre_idx = length(outdata.vars);
0467     outdata.data{gsre_idx} = ensemble_init_data_struct();
0468     outdata.data{gsre_idx}.type = 'gsr_epochs';
0469     outdata.data{gsre_idx}.vars = lvars;
0470     gsre_cols = set_var_col_const(outdata.data{gsre_idx}.vars);
0471     outdata.data{gsre_idx}.data{gsre_cols.subject_id} = {};
0472     outdata.data{gsre_idx}.data{gsre_cols.session} = [];
0473     outdata.data{gsre_idx}.data{gsre_cols.ensemble_id} = [];
0474     outdata.data{gsre_idx}.data{gsre_cols.run} = [];
0475     outdata.data{gsre_idx}.data{gsre_cols.trial} = [];
0476     outdata.data{gsre_idx}.data{gsre_cols.signal} = {};
0477     outdata.data{gsre_idx}.data{gsre_cols.srate} = [];
0478     outdata.data{gsre_idx}.data{gsre_cols.peakidxs} = {};
0479     for in=1:length(resp_names)
0480       outdata.data{gsre_idx}.data{gsre_cols.(resp_names{in})} = {};
0481     end
0482   end
0483 
0484   if ~isempty(strmatch('cardiac',channels)) || ...
0485           ~isempty(strmatch('pulse',channels))
0486     % initialize a cardiac summary data output structure
0487     outdata.vars = [outdata.vars 'cardiac_epochs'];
0488     carde_idx = length(outdata.vars);
0489     outdata.data{carde_idx} = ensemble_init_data_struct();
0490     outdata.data{carde_idx}.type = 'cardiac_epochs';
0491     outdata.data{carde_idx}.vars = lvars;
0492     carde_cols = set_var_col_const(outdata.data{carde_idx}.vars);
0493     outdata.data{carde_idx}.data{carde_cols.subject_id} = {};
0494     outdata.data{carde_idx}.data{carde_cols.session} = [];
0495     outdata.data{carde_idx}.data{carde_cols.ensemble_id} = [];
0496     outdata.data{carde_idx}.data{carde_cols.run} = [];
0497     outdata.data{carde_idx}.data{carde_cols.trial} = [];
0498     outdata.data{carde_idx}.data{carde_cols.signal} = {};
0499     outdata.data{carde_idx}.data{carde_cols.srate} = [];
0500     outdata.data{carde_idx}.data{carde_cols.peakidxs} = {};
0501     for in=1:length(resp_names)
0502       outdata.data{carde_idx}.data{carde_cols.(resp_names{in})} = {};
0503     end
0504   end
0505 end
0506 
0507 % % % extract data?
0508 if EXTRACT_DATA
0509 
0510   %
0511   % START OF THE SUBJECT LOOP
0512   %
0513 
0514   for isub=1:nsub_proc
0515     subid = sinfo(isub).id;
0516     msg = sprintf('\t\tPROCESSING SUBJECT (%d/%d): %s\n', isub, nsub_proc,subid);
0517     r = update_report(r,msg);
0518 
0519     % Determine number of sessions for this subject
0520     nsess = length(sinfo(isub).sessinfo);
0521   
0522     %
0523     % START OF THE SESSION LOOP
0524     %
0525   
0526     for isess = 1:nsess
0527       sess = sinfo(isub).sessinfo(isess);
0528     
0529       if isfield(sess,'use_session') && ~sess.use_session
0530         msg = sprintf('\t\t\tSkipping session %s\n',sess.id);
0531         r = update_report(r,msg);
0532         continue
0533       elseif ~isfield(sess,'physio') || ~isstruct(sess.physio)
0534         msg = sprintf(['\t\t\tNo Physio sinfo for session %s, '...
0535             'SKIPPING!\n'],sess.id);
0536         r = update_report(r,msg);
0537         continue
0538       end
0539 
0540       pfilt.include.all.subject_id = {subid};
0541       pfilt.include.all.session = {sess.id};
0542       pdata = ensemble_filter(pres_paths,pfilt);
0543     
0544       if isempty(pdata.data{prescol.path}) || ...
0545               ~exist(pdata.data{prescol.path}{1},'file')
0546         error('can not find presentation data for subject %s, session %s',...
0547             subid,sess.id);
0548       else
0549         presdata = load(pdata.data{prescol.path}{1});
0550         presdcol = set_var_col_const(presdata.vars);
0551       end
0552     
0553       nruns = length(sess.use_epi_runs);
0554       for irun=1:nruns
0555       
0556         rnum = sess.use_epi_runs(irun);
0557         msg = sprintf('\t\tPROCESSING RUN %d (%d/%d)\n',rnum,irun,nruns);
0558         r = update_report(r,msg);
0559       
0560         lpresfilt.include.all.SUB_ID = {subid};
0561         lpresfilt.include.all.RUN = rnum;
0562       
0563         % get presentation data for this run
0564         lpres = ensemble_filter(presdata,lpresfilt);
0565       
0566         lresp = cell(nresp,1);
0567         for ir=1:nresp
0568           lrespfilt = responses(ir).event;
0569           lrespfilt.include.all.RUN = rnum;
0570           lrespdata = ensemble_filter(lpres,lrespfilt);
0571         
0572           nlresp = length(lrespdata.data{presdcol.RESP_CODE});
0573           lresp{ir} = cell(nlresp,1);
0574           for iresp=1:nlresp
0575             if ir > length(lresp) || iresp > length(lresp{ir}) || ...
0576                     iresp > length(lrespdata.data{presdcol.RESP_CODE}) || ...
0577                     isempty(lrespdata.data{presdcol.RESP_CODE}{iresp})
0578               fprintf(1,'error?\n');
0579             end % if ir > length(lresp) || iresp > lresp{ir
0580             lresp{ir}{iresp} = lrespdata.data{presdcol.RESP_CODE}{iresp}{1};
0581           end
0582         end
0583 
0584         lphysfilt = pfilt;
0585         lphysfilt.include.all.run = rnum;
0586         if use_filtered
0587           lphysfilt.include.all.path = {'.*filtered.*'};
0588         end
0589       
0590         lphys = ensemble_filter(physio_paths,lphysfilt);
0591         if isempty(lphys.data{physcol.path}) || ...
0592                 ~exist(lphys.data{physcol.path}{1},'file')
0593           warning('no physio file found for %s, %s, run %d\n',...
0594               subid,sess.id,rnum);
0595           continue
0596         end
0597       
0598         % load data
0599         [fp,fn,fx] = fileparts(lphys.data{physcol.path}{1});
0600         EEG = pop_loadset('filename',sprintf('%s%s',fn,fx),'filepath',fp);
0601       
0602         % set up epochs using presentation data
0603         if (isnumeric(epoch_start) && (length(epoch_start) > 1 || ...
0604                 epoch_start)) || ...
0605                 (isstruct(epoch_start) && isfield(epoch_start,'filter'))
0606 
0607           %%%% find beginnings of epochs
0608           % FIXME: add ability to handle multiple epoch definitions, for
0609           % instance a def for rest period and a def for task period
0610           if isnumeric(epoch_start) && length(epoch_start) > 1
0611             % epoch times provided in a vector, assumed in seconds
0612             estimes = epoch_start;
0613             ns = length(estimes);
0614           elseif isstruct(epoch_start) && isfield(epoch_start,'filter')
0615             % filter criteria for presentation data provided
0616             % events matching these criteria will be taken as epoch starts
0617             esf = epoch_start.filter;
0618             esdata = ensemble_filter(lpres,esf);
0619             if isempty(esdata.data{presdcol.EVENT_TYPE})
0620               warning(['epoch start data missing, subject %s, session '...
0621                   '%s, run %d\n'],subid,sess.id,rnum);
0622               continue
0623             end
0624             estimes = esdata.data{presdcol.RUN_REL_TIME}/1000;
0625             ns = length(estimes);
0626           else
0627             warning(['no epoch start information, skipping run %d, sub '...
0628                 '%s, session %s\n'],rnum,subid,sess.id);
0629             continue
0630           end
0631         
0632           %%%% find ends of epochs
0633           if isnumeric(epoch_end) && length(epoch_end) == 1
0634             % epoch time provided, assumed in seconds
0635             eetimes = estimes + epoch_end;
0636             ne = length(eetimes);
0637           elseif isnumeric(epoch_end) && length(epoch_end) > 1
0638             % epoch end times provided in a vector, assumed in seconds
0639             % will later be compared to epoch start times
0640             eetimes = epoch_end;
0641             ne = length(eetimes);
0642           elseif isstruct(epoch_end) && isfield(epoch_end,'filter')
0643             % filter criteria for presentation data provided
0644             % events matching these criteria will be taken as epoch ends
0645             eef = epoch_end.filter;
0646             eedata = ensemble_filter(lpres,eef);
0647             if isempty(eedata.data{presdcol.EVENT_TYPE})
0648               warning(['epoch end data missing, subject %s, session '...
0649                   '%s, run %d\n'],subid,sess.id,rnum);
0650               continue
0651             end
0652             eetimes = eedata.data{presdcol.RUN_REL_TIME}/1000;
0653             ne = length(eetimes);
0654           else
0655             warning(['no epoch end information provided for sub %s, sess '...
0656                 '%s, run %d, calculating as minimum distance between '...
0657                 'epoch onsets\n'],subid,sess.id,rnum);
0658             etime = min(diff(estimes));
0659           end
0660           
0661           % account for raw signal baseline period
0662           estimes = estimes - EEG.xmin;
0663           eetimes = eetimes - EEG.xmin;
0664 
0665           if ~exist('etime','var') && ne ~= ns || any (estimes > eetimes)
0666             % end time not specified,
0667             warning(['outlier starting or ending events, or events not '...
0668                 'aligned, sub %s, session %s, run %d\n'],subid,sess.id,rnum);
0669             continue
0670             % FIXME: should have some logic here to at least attempt to
0671             % reconcile estimes and eetimes ... rules such as "no end time
0672             % before the first start, no start time before the first end", as
0673             % well as "line up each start/end, 1 to 1, pick out outliers, and
0674             % reconcile pairs whose end is before it's given start
0675           end
0676 
0677           % find epoch length, if not explicitly provided
0678           if ~exist('etime','var') 
0679             % make sure to remove epochs that exceed the samples in EEG
0680             ee_end_idxs = eetimes > (EEG.pnts/EEG.srate);
0681             if any(ee_end_idxs)
0682               warning('%d epochs out of bounds, removing\n',length(ee_end_idxs));
0683             end
0684             eetimes(ee_end_idxs) = []; % remove exceeding epochs from eetimes
0685             estimes(ee_end_idxs) = []; % remove exceeding epochs from estimes
0686             dtimes = eetimes - estimes;
0687             if isstruct(epoch_end) && isfield(epoch_end,'minimum_end') ...
0688                     && epoch_end.minimum_end
0689               etime = min(dtimes);
0690             else
0691               etime = max(dtimes);
0692             end
0693           end
0694 
0695           % double check that all epochs end before the end of the data set
0696           badtimes = (estimes+etime) > (EEG.pnts/EEG.srate);
0697           if any(badtimes)
0698             ee_end_idxs = find(badtimes);
0699             warning('%d epochs out of bounds, removing\n',length(ee_end_idxs));
0700             eetimes(ee_end_idxs) = []; % remove exceeding epochs from eetimes
0701             estimes(ee_end_idxs) = []; % remove exceeding epochs from estimes
0702           end
0703           ns = length(estimes);
0704           ne = length(eetimes);
0705 
0706           % import event onsets
0707           EEG = pop_importevent(EEG,'append','no','event',[ones(ns,1) ...
0708               estimes],'fields',{'type','latency'},'timeunit',1);
0709 
0710           % epoch the data based upon event onsets and the largest event duration
0711           EEG = pop_epoch(EEG,[],[-baseline etime],'eventindices',1:ns);
0712 
0713           if EEG.trials ~= ns
0714             error('one or more epochs was lost, sub %s, %s, run %d\n',...
0715                 subid,sess.id,rnum);
0716           end
0717           
0718         else
0719           ns = 1;
0720         end % if (isnumeric(epoch_start
0721         
0722         % set sampling rate
0723         scr.pk.Fs = EEG.srate;
0724         pulse.pk.Fs = EEG.srate;
0725 
0726         % set up a figure file for tsub/sess/run
0727         if GRAPH_EPOCHS && SAVE_EPOCHS
0728           figfname = fullfile(fp,sprintf('%s_epochs.ps',fn));
0729         end
0730         
0731         for ic=1:nchan
0732           c = channels{ic};
0733 
0734           baseline_end = (baseline*EEG.srate)+1; % next sample after end of baseline
0735 
0736           % extract channel location
0737           cidx = strmatch(c,{EEG.chanlocs.labels});
0738           if isempty(cidx)
0739             warning(['no channel label found for %s, skipping for '...
0740                 'sub %s, session %s, run %d'],c,subid,sess.id,rnum);
0741             continue
0742           end
0743 
0744           switch c
0745             case {'gsr','scr'}
0746             
0747               % extract all epochs for this channel
0748               epochs = squeeze(EEG.data(cidx,:,:));
0749               
0750               sepoch = size(epochs);
0751               if (sepoch(1) == ns && sepoch(2) ~= ns)
0752                 epochs = epochs';
0753                 sepoch = size(epochs);
0754               end
0755 
0756               % subtract baseline for each epoch
0757               for iep=1:ns
0758                 if baseline
0759                   mbase = mean(epochs(1:(baseline*EEG.srate),iep));
0760                   pksig = epochs(:,iep) - mbase;
0761                   if ~keep_baseline
0762                     pksig(1:baseline*EEG.srate) = [];
0763                   end
0764                 else
0765                   pksig = epochs(:,iep);
0766                 end
0767                 
0768                 % get peaks
0769                 [pidxs,heights,bidxs] = find_peaks(pksig,scr.pk);
0770 
0771                 if baseline && keep_baseline
0772                   % remove peaks whose troughs occur before 1 sec after
0773                   % stimulus onset
0774                   rmbase = bidxs < (baseline+1)*EEG.srate;
0775                   pidxs(rmbase) = [];
0776                   heights(rmbase) = [];
0777                   bidxs(rmbase) = [];
0778                 end
0779                 
0780                 % graph??
0781                 if GRAPH_EPOCHS
0782                   m = mod(iep-1,ep_per_fig);
0783                   if m == 0, figure(); end
0784                   subplot(ep_per_fig,1,m+1);
0785                   plot(pksig);
0786                   title(sprintf(['%s, %s, run %d, signal %s, epoch'...
0787                       ' %d'],subid,sess.id,rnum,c,iep));
0788                   hold on;
0789                   for ipk = 1:length(pidxs)
0790                     plot(pidxs(ipk),pksig(pidxs(ipk)),'g*');
0791                   end
0792                   hold off;
0793                   
0794                   if SAVE_EPOCHS && (m+1) == ep_per_fig
0795                     % save graphed epochs to a file
0796                     print(pargs{:},figfname);
0797                   end
0798                 end
0799                 
0800                 % save signal out to gsr_epochs
0801                 outdata.data{gsre_idx} = ensemble_add_data_struct_row(...
0802                     outdata.data{gsre_idx},'subject_id',subid,'session',...
0803                     sess.id,'ensemble_id',sess.ensemble_id,'run',rnum,...
0804                     'trial',iep,'signal',{pksig},'srate',EEG.srate,...
0805                     'peakidxs',pidxs);
0806               end % for iep=1:ns
0807               
0808               % save final figure to file if not already saved
0809               if GRAPH_EPOCHS && SAVE_EPOCHS && (m+1) ~= ep_per_fig
0810                 print(pargs{:},figfname);
0811               end
0812 
0813             case {'cardiac','pulse'}
0814                 
0815               % iterate over epochs
0816               for iep=1:ns
0817                 % find peaks
0818                 pksig = EEG.data(cidx,baseline_end:end,iep);
0819                 pidxs = find_peaks(pksig,pulse.pk);
0820 
0821                 % save signal out to cardiac_epochs
0822                 outdata.data{carde_idx} = ensemble_add_data_struct_row(...
0823                     outdata.data{carde_idx},'subject_id',subid,'session',...
0824                     sess.id,'ensemble_id',sess.ensemble_id,'run',rnum,...
0825                     'trial',iep,'signal',{pksig'},'srate',EEG.srate,...
0826                     'peakidxs',pidxs);
0827 
0828                 % graph??
0829                 if GRAPH_EPOCHS
0830                   m = mod(iep-1,ep_per_fig);
0831                   if m == 0, figure(); end
0832                   subplot(ep_per_fig,1,m+1);
0833                   plot(pksig);
0834                   title(sprintf(['%s, %s, run %d, signal %s, epoch'...
0835                       ' %d'],subid,sess.id,rnum,c,iep));
0836                   hold on;
0837                   for ipk = 1:length(pidxs)
0838                     plot(pidxs(ipk),pksig(pidxs(ipk)),'g*');
0839                   end
0840                   hold off;
0841 
0842                   if SAVE_EPOCHS && (m+1) == ep_per_fig
0843                     % save graphed epochs to a file
0844                     print(pargs{:},figfname);
0845                   end
0846                 end
0847               end % for iep=1:ns
0848               
0849               % save final figure to file if not already saved
0850               if GRAPH_EPOCHS && SAVE_EPOCHS && (m+1) ~= ep_per_fig
0851                 print(pargs{:},figfname);
0852               end
0853 
0854             otherwise
0855               warning('unknown channel %s',c);
0856               
0857           end % switch c
0858 
0859           % add response information
0860           for iresp = 1:length(resp_names)
0861             % get all responses
0862             rname = resp_names{iresp};
0863             respdata = lresp{iresp};
0864             nrespd = length(respdata);
0865             % check against length of events (ns)
0866             if nrespd < ns
0867               warning(['data for response %s of different length (%d) '...
0868                   'than number of epochs (%d), these are not being '...
0869                   'added to the output\n'],rname,nrespd,ns);
0870               continue
0871             elseif ns < nrespd
0872               warning(['more responses (%d) than epochs (%d), assume '...
0873                   'that epochs were dropped from the end, and dropping '...
0874                   'extra responses'],nrespd,ns);
0875               respdata(ns+1:end) = [];
0876             end
0877 
0878             % get indices for this channel and this response
0879             switch c
0880               case {'gsr','scr'}
0881                 oidx = gsre_idx;
0882                 ridx = gsre_cols.(rname);
0883               case {'cardiac','pulse'}
0884                 oidx = carde_idx;
0885                 ridx = carde_cols.(rname);
0886               otherwise
0887                 warning('unknown channel type (%s)\n',c);
0888                 continue
0889             end
0890 
0891             % save into outdata
0892             outdata.data{oidx}.data{ridx}(end-ns+1:end) = respdata;
0893             
0894           end % for iresp=1:length(resp_names
0895         end % for ic=1:nchan
0896       end % for irun=1:
0897     end % for isess = 1:
0898 
0899     if (GRAPH_EPOCHS && EPOCH_BY_RESP)
0900       pfilt = struct();
0901       pfilt.include.all.subject_id = {subid};
0902       pfilt.include.all.path_type = {'physio_outdir'};
0903       lpathdata = ensemble_filter(pathdata,pfilt);
0904       if ~isempty(lpathdata.data{1})
0905         physdir = lpathdata.data{pcol.path}{1};
0906         check_dir(physdir);
0907       else
0908         warning(['could not find physio_outdir in pathdata, not graphing '...
0909             'epochs by response for subject %s'],subid);
0910         continue
0911       end
0912 
0913       if GRAPH_EPOCHS && EPOCH_BY_RESP
0914         figfname = fullfile(physdir,sprintf('%s_epochs_by_resp.ps',subid));
0915         % graph epochs by response type
0916         for iresp=1:length(resp_names)
0917           rname = resp_names{iresp};
0918           for ic=1:nchan
0919             c = channels{ic};
0920             % get channel and response indexes
0921             switch c
0922               case {'gsr','scr'}
0923                 eidx = gsre_idx;
0924                 ridx = gsre_cols.(rname);
0925                 sidx = gsre_cols.signal;
0926               case {'cardiac','pulse'}
0927                 eidx = carde_idx;
0928                 ridx = carde_cols.(rname);
0929                 sidx = carde_cols.signal;
0930               otherwise
0931                 warning('unknown channel type (%s)\n',c);
0932                 continue
0933             end % switch c
0934 
0935             % get mask matrix
0936             if iscell(outdata.data{eidx}.data{ridx}) && ...
0937                     all(cellfun(@isempty,outdata.data{eidx}.data{ridx}))
0938               warning('no responses for response_name %s\n',rname);
0939               continue
0940             end
0941 
0942             [rm,ur] = make_mask_mtx(outdata.data{eidx}.data{ridx});
0943         
0944             for ir=1:length(ur)
0945               figure();
0946               allsig = [];
0947               hold on;
0948               lidxs = find(rm(:,ir));
0949               for i=1:length(lidxs)
0950                 li=lidxs(i);
0951                 sig=outdata.data{eidx}.data{sidx}{li};              
0952                 if isempty(sig), continue, end
0953                 if ~isempty(strmatch(c,{'cardiac','pulse'}))
0954                   pidxs=find_peaks(sig,pulse.pk);
0955                   if isempty(pidxs), continue, end
0956                   sig=60./(diff(pidxs)/EEG.srate);
0957                 end
0958                 if size(sig,1) > size(sig,2), sig = sig'; end
0959                 plot(sig,'g');
0960                 ssdiff = size(allsig,2) - size(sig,2);
0961                 if ssdiff > 0
0962                   % allsig is longer, extend sig with nans at the end
0963                   sig = [sig nan(1,ssdiff)];
0964                 elseif ssdiff < 0
0965                   % sig is longer, extend allsig with nans at the end
0966                   allsig = [allsig nan(size(allsig,1),-ssdiff)];
0967                 end
0968                 allsig = [allsig; sig];
0969               end % for i=1:length(lidxs
0970               plot(nanmean(allsig),'k');
0971               hold off;
0972               lr = ur{ir};
0973               if isnumeric(lresp), lresp = num2str(lresp); end
0974               title(sprintf('signal: %s, response: %s, %s',c,rname,lr));
0975           
0976               if SAVE_EPOCHS, print(pargs{:},figfname); end
0977             end % for ir=1:length(ur
0978           end % for iresp=1:length(resp_names
0979         end % for ic=1:nchan
0980       end % if GRAPH_EPOCHS && EPOCH_BY_RESP
0981     end % if (GRAPH_EPOCHS && EPOCH_BY_RESP) ||
0982 
0983   end % for isub=1:
0984 
0985 end % if EXTRACT_DATA
0986 
0987 if CALCULATE_STATS
0988 
0989   % get gsr data?
0990   if exist('gsr_epochs','var')
0991     gdata = gsr_epochs;
0992   elseif exist('gsre_idx','var')
0993     gdata = outdata.data{gsre_idx};
0994   else
0995     warning('no gsr epoched data found');
0996   end
0997   
0998   % get cardiac data?
0999   if exist('cardiac_epochs','var')
1000     cdata = cardiac_epochs;
1001   elseif exist('carde_idx','var')
1002     cdata = outdata.data{carde_idx};
1003   else
1004     warning('no cardiac epoched data found');
1005   end
1006   
1007   if exist('gdata','var')
1008       
1009     % add response information
1010     ridxs = ~ismember(gdata.vars,xvars);
1011     rvars = {gdata.vars{ridxs}};
1012     if ~iscell(rvars), rvars = {rvars}; end
1013 
1014     % initialize a gsr summary data output structure
1015     outdata.vars = [outdata.vars 'gsr_summ'];
1016     gsr_idx = length(outdata.vars);
1017     outdata.data{gsr_idx} = ensemble_init_data_struct();
1018     outdata.data{gsr_idx}.type = 'gsr_summ';
1019     outdata.data{gsr_idx}.vars = {'subject_id','session','ensemble_id',...
1020         'run','trial','integral','scl','sclidx',...
1021         'ttl_integral','ttl_scl','ttl_sclidx',...
1022         'npeaks','mean_peak_height','peak_height_variance',...
1023         'first_peak_time','first_peak_height',rvars{:}};
1024     gsr_cols = set_var_col_const(outdata.data{gsr_idx}.vars);
1025     outdata.data{gsr_idx}.data{gsr_cols.subject_id} = {};
1026     outdata.data{gsr_idx}.data{gsr_cols.session} = [];
1027     outdata.data{gsr_idx}.data{gsr_cols.ensemble_id} = [];
1028     outdata.data{gsr_idx}.data{gsr_cols.run} = [];
1029     outdata.data{gsr_idx}.data{gsr_cols.trial} = [];
1030     outdata.data{gsr_idx}.data{gsr_cols.integral} = [];
1031     outdata.data{gsr_idx}.data{gsr_cols.scl} = [];
1032     outdata.data{gsr_idx}.data{gsr_cols.sclidx} = [];
1033     outdata.data{gsr_idx}.data{gsr_cols.ttl_integral} = [];
1034     outdata.data{gsr_idx}.data{gsr_cols.ttl_scl} = [];
1035     outdata.data{gsr_idx}.data{gsr_cols.ttl_sclidx} = [];
1036     outdata.data{gsr_idx}.data{gsr_cols.npeaks} = [];
1037     outdata.data{gsr_idx}.data{gsr_cols.mean_peak_height} = [];
1038     outdata.data{gsr_idx}.data{gsr_cols.peak_height_variance} = [];
1039     outdata.data{gsr_idx}.data{gsr_cols.first_peak_time} = [];
1040     outdata.data{gsr_idx}.data{gsr_cols.first_peak_height} = [];
1041     for in=1:length(rvars)
1042       outdata.data{gsr_idx}.data{gsr_cols.(rvars{in})} = {};
1043     end
1044 
1045     % initialize gsr vars
1046     gcol = set_var_col_const(gdata.vars);
1047     ng = length(gdata.data{1});
1048     
1049     % calculate stats for each epoch
1050     for ie=1:ng
1051       sig   = gdata.data{gcol.signal}{ie};
1052       srate = gdata.data{gcol.srate}(ie);
1053       pidxs = gdata.data{gcol.peakidxs}{ie};
1054       if iscell(pidxs), pidxs = pidxs{1}; end
1055       if ~isempty(epoch_end) && isnumeric(epoch_end)
1056         end_samp = epoch_end*srate;
1057         sig = sig(1:end_samp);
1058         pidxs(pidxs > end_samp) = [];
1059       end
1060       if isempty(pidxs)
1061         pidxs = [];
1062         hghts = [];
1063         mpkhght = NaN;
1064         pkhghtv = NaN;
1065         fpktime = NaN;
1066         fpkhght = NaN;
1067       else
1068         hghts = peak_heights('signal',sig,'params',scr.pk.peak_heights,...
1069             'peakIdxs',pidxs);
1070         mpkhght = mean(hghts);
1071         pkhghtv = var(hghts);
1072         fpktime = pidxs(1)/srate;
1073         fpkhght = hghts(1);
1074       end
1075       npeaks = length(pidxs);
1076       
1077       baseline_end = baseline*srate;
1078 
1079       % calculate integral within integral window defined by
1080       % int_begin and int_end
1081       startint = baseline_end + scr.int_begin*srate + 1;
1082       endint   = min([(baseline_end + scr.int_end*srate) length(sig)]);
1083 
1084       integral = trapz(sig(startint:endint));
1085       mscl = mean(sig(startint:endint));
1086       sclidx = trapz(abs(sig(startint:endint)))/(endint-startint);
1087       ttl_int = trapz(sig(baseline_end+1:end));
1088       ttlmscl = mean(sig(baseline_end+1:end));
1089       ttlsclidx = trapz(abs(sig(baseline_end+1:end)))/(length(sig)-baseline_end);
1090 
1091       % save to outdata
1092       outdata.data{gsr_idx} = ensemble_add_data_struct_row(...
1093           outdata.data{gsr_idx},...
1094           'subject_id',gdata.data{gcol.subject_id}{ie},...
1095           'session',gdata.data{gcol.session}(ie),...
1096           'ensemble_id',gdata.data{gcol.ensemble_id}(ie),...
1097           'run',gdata.data{gcol.run}(ie),...
1098           'trial',gdata.data{gcol.trial}(ie),...
1099           'integral',integral(:),'scl',mscl(:),'sclidx',sclidx(:),...
1100           'ttl_integral',ttl_int(:),'ttl_scl',ttlmscl(:),...
1101           'ttl_sclidx',ttlsclidx(:),...
1102           'npeaks',npeaks,'mean_peak_height',mpkhght,...
1103           'peak_height_variance',pkhghtv,...
1104           'first_peak_time',fpktime,'first_peak_height',fpkhght);
1105 
1106     end % for ie=1:ng
1107     
1108     % pass on response information
1109     for ir = 1:length(rvars)
1110       % get column indices
1111       ri = gcol.(rvars{ir});
1112       ro = gsr_cols.(rvars{ir});
1113 
1114       % save into outdata
1115       outdata.data{gsr_idx}.data{ro} = gdata.data{ri};
1116     end
1117 
1118   end %  if exist('gdata
1119   
1120   if exist('cdata','var')
1121       
1122     % add response information
1123     ridxs = ~ismember(cdata.vars,xvars);
1124     rvars = {cdata.vars{ridxs}};
1125     if ~iscell(rvars), rvars = {rvars}; end
1126 
1127     % initialize a cardiac summary data output structure
1128     outdata.vars = [outdata.vars 'cardiac_summ'];
1129     card_idx = length(outdata.vars);
1130     outdata.data{card_idx} = ensemble_init_data_struct();
1131     outdata.data{card_idx}.type = 'cardiac_summ';
1132     outdata.data{card_idx}.vars = {'subject_id','session','ensemble_id',...
1133         'run','trial','heart_rate','hr_std','hr_rmssd','hr_slope',rvars{:}};
1134     card_cols = set_var_col_const(outdata.data{card_idx}.vars);
1135     outdata.data{card_idx}.data{card_cols.subject_id} = {};
1136     outdata.data{card_idx}.data{card_cols.session} = [];
1137     outdata.data{card_idx}.data{card_cols.ensemble_id} = [];
1138     outdata.data{card_idx}.data{card_cols.run} = [];
1139     outdata.data{card_idx}.data{card_cols.trial} = [];
1140     outdata.data{card_idx}.data{card_cols.heart_rate} = [];
1141     outdata.data{card_idx}.data{card_cols.hr_rmssd} = [];
1142     outdata.data{card_idx}.data{card_cols.hr_std} = [];
1143     outdata.data{card_idx}.data{card_cols.hr_slope} = [];
1144     for in=1:length(rvars)
1145       outdata.data{card_idx}.data{card_cols.(rvars{in})} = {};
1146     end
1147 
1148     % initialize cardiac vars
1149     ccol = set_var_col_const(cdata.vars);
1150     nc = length(cdata.data{1});
1151     
1152     % iterate over epochs
1153     for iep=1:nc
1154       srate = cdata.data{ccol.srate}(iep);
1155       pidxs = cdata.data{ccol.peakidxs}{iep};
1156       if iscell(pidxs), pidxs = pidxs{1}; end
1157       if ~isempty(epoch_end) && isnumeric(epoch_end)
1158         end_samp = epoch_end*srate;
1159         sig = sig(1:end_samp);
1160         pidxs(pidxs > end_samp) = [];
1161       end
1162       if isempty(pidxs), pidxs = nan; end
1163 
1164       % get beat-to-beat intervals in ms
1165       dtimes = diff(pidxs)/srate*1000;
1166       dbak = dtimes;
1167       
1168       % remove out-of-bounds values, such that 300ms <= dtimes <= 2000ms
1169       % Timonen et al 2006, in Li et al 2009 IntJPsychophysiology
1170       dbak(dbak < 300) = [];
1171       dbak(dbak > 2000) = [];
1172       dtimes(dtimes > mean(dtimes)+std(dtimes)*3) = [];
1173       dtimes(dtimes < mean(dtimes)-std(dtimes)*3) = [];
1174       dtimes(dtimes < 300) = [];
1175       dtimes(dtimes > 2000) = [];
1176 
1177       % heart rate (in beats per minute) and heart rate variance
1178       hr_bpm = 1/(mean(dtimes)/1000)*60;
1179       hr_std = std(dtimes);
1180       hr_rmssd = sqrt(mean(diff(dtimes).^2));
1181 
1182       % slope of change in heart rate
1183       zlength = zscore(1:length(dtimes))';
1184       ztimes = zscore(dtimes);
1185       hr_slope = regress(ztimes,zlength);
1186 
1187       % save to outdata
1188       outdata.data{card_idx} = ensemble_add_data_struct_row(...
1189           outdata.data{card_idx},...
1190           'subject_id', cdata.data{ccol.subject_id}{iep},...
1191           'session',    cdata.data{ccol.session}(iep),...
1192           'ensemble_id',cdata.data{ccol.ensemble_id}(iep),...
1193           'run',        cdata.data{ccol.run}(iep),...
1194           'trial',      cdata.data{ccol.trial}(iep),...
1195           'heart_rate',hr_bpm,'hr_rmssd',hr_rmssd,...
1196           'hr_std',hr_std,'hr_slope',hr_slope);
1197 
1198     end % for iep=1:nc
1199 
1200     
1201     % pass on response information
1202     for ir = 1:length(rvars)
1203       % get column indices
1204       ri = ccol.(rvars{ir});
1205       ro = card_cols.(rvars{ir});
1206 
1207       % save into outdata
1208       outdata.data{card_idx}.data{ro} = cdata.data{ri};
1209     end
1210     
1211   end % if exist('cdata
1212 
1213   if SAVE_DATA
1214     for ic=1:nchan
1215       c = channels{ic};
1216       % get channel idxs
1217       switch c
1218         case {'gsr','scr'}
1219           oidx = gsr_idx;
1220           ocol = gsr_cols;
1221         case {'cardiac','pulse'}
1222           oidx = card_idx;
1223           ocol = card_cols;
1224         otherwise
1225           warning('unknown channel type (%s)\n',c);
1226       end % switch c
1227 
1228       % find unique subjects in calculated data
1229       us = unique(outdata.data{oidx}.data{ocol.subject_id});
1230       ns = length(us);
1231       
1232       % iterate over subjects
1233       for is=1:ns
1234         subid = us{is};
1235         
1236         % get output data for this subject
1237         pfilt = struct();
1238         pfilt.include.all.subject_id = {subid};
1239         odata = ensemble_filter(outdata.data{oidx},pfilt);
1240         
1241         % get path data for this subject
1242         pfilt.include.all.path_type = {'physio_outdir'};
1243         lpathdata = ensemble_filter(pathdata,pfilt);
1244         if isempty(lpathdata.data{1})
1245           warning(['could not find physio_outdir in pathdata, not '...
1246               'saving data for subject %s'],subid);
1247         elseif isempty(odata.data{1})
1248           warning('could not find output data for subject %s, SKIPPING',...
1249               subid);
1250         else
1251           physdir = lpathdata.data{pcol.path}{1};
1252           check_dir(physdir);
1253 
1254           xsp.export = export;
1255           xsp.export.fname = fullfile(physdir,...
1256               sprintf('%s_%s_%sexport.txt',subid,odata.type,export_stub));
1257           xsp.sas = sas;
1258           xsp.sas.fname = fullfile(physdir,sprintf('%s_%s_%sexport.sas',...
1259               subid,odata.type,export_stub));
1260           xsp.sas.libname=sprintf('s%s_%s%s',subid,odata.type,export_stub);
1261           ensemble_export_sastxt(odata,xsp);
1262         end % if ~isempty(lpathdata.data{1
1263       end % for is=1:ns
1264     end % for ic=1:nchan
1265   end % if SAVE_DATA
1266 
1267 end % if CALCULATE_STATS

Generated on Sun 19-May-2019 04:00:51 by m2html © 2003