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_summ0001 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