Home > utils > physio > read_biopac.m

read_biopac

PURPOSE ^

Reads physiological data measured by the BioPac system, epochs data

SYNOPSIS ^

function ri = read_biopac(fname,params)

DESCRIPTION ^

 Reads physiological data measured by the BioPac system, epochs data
 
   ri = read_biopac(fname,params)
 
 The BioPac system has been used by the Janata Lab at the 3T
 scanner at the UC Davis Imaging Research Center, to collect physiological
 measures while participants were being scanned. This system, though,
 could conceivably be used to collect data outside of the scanner.
 Therefore, this script does not assume that physiological data it
 analyzes was collected in an NMR environment, and it does not assume that
 a TR is available for epoching purposes, but it does allow for a TR to be
 used for epoching, and it also allows for a more general presentation
 logfile-based epoching.
 
 REQUIRES
   - reading data -
   fname - path to single file to be processed
   params.channels(:) - struct array containing information about the
     channels to be extracted and returned
     .name - channel name, that will be converted to a fieldname in ri
       i.e. if a channel name is "cardiac", then that data will be found
       in ri.signal.cardiac in the output data
     .comment_text - name given to the data stream in AcqKnowledge/BioPac,
       as reflected in .hdr.per_chan_data.comment_text in the .acq file
     NOTE: special channel name 'trigger' treated as scanner pulses,
     assumed to be TTL pulses time-locked to volume acquisition onset
   params.signal_bounds - 2-length vector containing start and end times
       in seconds to initially bound the raw signal. this can be used in
       the case that one physio file contains data from multiple runs.
       Since data here are handled on a run-by-run basis, this can be used
       to limit the initial data that is extracted to the samples
       pertinent for the given run. a value of 0 in the first position
       will be replaced with a 1, while a value of 0 in the second
       position will be replaced with the final sample in the data set

   - epoching -
   params.bound_by_trigger - if this is true, and if there is a channel in
     params.channels(:) with name 'trigger', and if this data is not empty
     or a flat line, triggers will be used to define the epoch for this run
     ... the epoch will be first trigger through last trigger+TR*srate ...
     NOTE: if params.bound_by_trigger is set, and if a trigger channel
     exists, params.link2pres will be ignored
     NOTE: this requires params.TR ... if no params.TR is specified, it
     will be estimated as median(diff(trigger_onsets_in_ms))/1000
   params.baseline - the amount of time, in seconds, to allow at the
     beginning of the signal for a baseline calculation.
   params.TR - length (in seconds) between repetitions, for fMRI analyses.
     see also: params.bound_by_trigger
   params.link2pres - structure with settings to link up to a presentation
     response file ... MAKE SURE that you specify filt criteria that will
     return a set of items that then have a 1 to 1 correlation with event
     onsets in a channel within the acq file. If the return from your
     filtering criteria exactly matches the channel you've specified to
     link to, then all signals from the acq file will be aligned with the
     presentation data, and data falling before run_rel_time = 0 will be
     thrown out
     NOTE: if params.bound_by_trigger is set, and if a trigger channel
     exists, link2pres will be ignored
     NOTE: link2pres ASSUMES that the run begins AFTER physio data
     collection begins. Even if presentation is started before physio data
     collection, it will adjust all times in run_rel_time to what it
     thinks is the beginning of the run, and all signals processed by
     link2pres rely on this.
     .presfname - path and filename of the presentation file to link to
     .filt - filtering criteria applied to presfname data via
       ensemble_filter ... you may want to use this struct to filter by
       run, if you are calling read_biopac (or proc_physio_data) from
       within a loop over runs, for instance
     .channelname - ri.signal.(.channelname) will be compared to
       filtered results from presfname. If a one-to-one match is found
       between presfname filtered data and .channelname data, then epochs
       will be defined by the beginning of the presentation run_rel_time,
       and the last event returned by the presfname filtered data
     .diff_tol - a diff will be done for all onsets given for presentation
       data as well as link signal data, and then a diff will be done
       between those diffs. If the diffs don't line up in length, or if
       there is any diff between diffs greater than .diff_tol, link2pres
       will no be applied
     .bound_by_last_pulse - if set to 1, and if EVENT_TYPE={'Pulse'},
       EVENT_CODE = {'255'} filtering criteria returns any rows, all data
       after sample = ((lastPulseTime+TR*1000)/stime) will be dropped. If
       no params.TR is specified, TR=median(diff(pulse_onsets_in_ms))/1000
     .estimate_trigger - if a 'trigger' channel can not be found, but
       presentation data is provided and contains events in the form of
       EVENT_TYPE = {'Pulse'} && EVENT_CODE = 255, and if link2pres is
       specified, after temporal linking is achieved, the pulse events
       will be extracted and converted to triggers. Missing triggers will
       be interpolated, either using the specified TR, or using a TR
       estimated from the median pulse period.

 OUTPUT
   ri - 'runinfo' struct, in a form similar to that of read_keithley.m and
     read_mate.m, containing fields for data streams, events, and onsets,
     also containing a .meta field, with .meta.srate (sampling rate),
     .meta.stime (sampling time, or # ms per sample), and
     .meta.baseline_samps (# of baseline samples). This differs from that
     of read_keithley and read_mate in that the actual signal fields are
     found in ri.signal, rather than within ri
       NOTE: ri.trigger values are sample indices, NOT times in ms or sec

 FB 2008.10.12
 FB 2009.06.24 - added baseline calculation

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function ri = read_biopac(fname,params)
0002 
0003 % Reads physiological data measured by the BioPac system, epochs data
0004 %
0005 %   ri = read_biopac(fname,params)
0006 %
0007 % The BioPac system has been used by the Janata Lab at the 3T
0008 % scanner at the UC Davis Imaging Research Center, to collect physiological
0009 % measures while participants were being scanned. This system, though,
0010 % could conceivably be used to collect data outside of the scanner.
0011 % Therefore, this script does not assume that physiological data it
0012 % analyzes was collected in an NMR environment, and it does not assume that
0013 % a TR is available for epoching purposes, but it does allow for a TR to be
0014 % used for epoching, and it also allows for a more general presentation
0015 % logfile-based epoching.
0016 %
0017 % REQUIRES
0018 %   - reading data -
0019 %   fname - path to single file to be processed
0020 %   params.channels(:) - struct array containing information about the
0021 %     channels to be extracted and returned
0022 %     .name - channel name, that will be converted to a fieldname in ri
0023 %       i.e. if a channel name is "cardiac", then that data will be found
0024 %       in ri.signal.cardiac in the output data
0025 %     .comment_text - name given to the data stream in AcqKnowledge/BioPac,
0026 %       as reflected in .hdr.per_chan_data.comment_text in the .acq file
0027 %     NOTE: special channel name 'trigger' treated as scanner pulses,
0028 %     assumed to be TTL pulses time-locked to volume acquisition onset
0029 %   params.signal_bounds - 2-length vector containing start and end times
0030 %       in seconds to initially bound the raw signal. this can be used in
0031 %       the case that one physio file contains data from multiple runs.
0032 %       Since data here are handled on a run-by-run basis, this can be used
0033 %       to limit the initial data that is extracted to the samples
0034 %       pertinent for the given run. a value of 0 in the first position
0035 %       will be replaced with a 1, while a value of 0 in the second
0036 %       position will be replaced with the final sample in the data set
0037 %
0038 %   - epoching -
0039 %   params.bound_by_trigger - if this is true, and if there is a channel in
0040 %     params.channels(:) with name 'trigger', and if this data is not empty
0041 %     or a flat line, triggers will be used to define the epoch for this run
0042 %     ... the epoch will be first trigger through last trigger+TR*srate ...
0043 %     NOTE: if params.bound_by_trigger is set, and if a trigger channel
0044 %     exists, params.link2pres will be ignored
0045 %     NOTE: this requires params.TR ... if no params.TR is specified, it
0046 %     will be estimated as median(diff(trigger_onsets_in_ms))/1000
0047 %   params.baseline - the amount of time, in seconds, to allow at the
0048 %     beginning of the signal for a baseline calculation.
0049 %   params.TR - length (in seconds) between repetitions, for fMRI analyses.
0050 %     see also: params.bound_by_trigger
0051 %   params.link2pres - structure with settings to link up to a presentation
0052 %     response file ... MAKE SURE that you specify filt criteria that will
0053 %     return a set of items that then have a 1 to 1 correlation with event
0054 %     onsets in a channel within the acq file. If the return from your
0055 %     filtering criteria exactly matches the channel you've specified to
0056 %     link to, then all signals from the acq file will be aligned with the
0057 %     presentation data, and data falling before run_rel_time = 0 will be
0058 %     thrown out
0059 %     NOTE: if params.bound_by_trigger is set, and if a trigger channel
0060 %     exists, link2pres will be ignored
0061 %     NOTE: link2pres ASSUMES that the run begins AFTER physio data
0062 %     collection begins. Even if presentation is started before physio data
0063 %     collection, it will adjust all times in run_rel_time to what it
0064 %     thinks is the beginning of the run, and all signals processed by
0065 %     link2pres rely on this.
0066 %     .presfname - path and filename of the presentation file to link to
0067 %     .filt - filtering criteria applied to presfname data via
0068 %       ensemble_filter ... you may want to use this struct to filter by
0069 %       run, if you are calling read_biopac (or proc_physio_data) from
0070 %       within a loop over runs, for instance
0071 %     .channelname - ri.signal.(.channelname) will be compared to
0072 %       filtered results from presfname. If a one-to-one match is found
0073 %       between presfname filtered data and .channelname data, then epochs
0074 %       will be defined by the beginning of the presentation run_rel_time,
0075 %       and the last event returned by the presfname filtered data
0076 %     .diff_tol - a diff will be done for all onsets given for presentation
0077 %       data as well as link signal data, and then a diff will be done
0078 %       between those diffs. If the diffs don't line up in length, or if
0079 %       there is any diff between diffs greater than .diff_tol, link2pres
0080 %       will no be applied
0081 %     .bound_by_last_pulse - if set to 1, and if EVENT_TYPE={'Pulse'},
0082 %       EVENT_CODE = {'255'} filtering criteria returns any rows, all data
0083 %       after sample = ((lastPulseTime+TR*1000)/stime) will be dropped. If
0084 %       no params.TR is specified, TR=median(diff(pulse_onsets_in_ms))/1000
0085 %     .estimate_trigger - if a 'trigger' channel can not be found, but
0086 %       presentation data is provided and contains events in the form of
0087 %       EVENT_TYPE = {'Pulse'} && EVENT_CODE = 255, and if link2pres is
0088 %       specified, after temporal linking is achieved, the pulse events
0089 %       will be extracted and converted to triggers. Missing triggers will
0090 %       be interpolated, either using the specified TR, or using a TR
0091 %       estimated from the median pulse period.
0092 %
0093 % OUTPUT
0094 %   ri - 'runinfo' struct, in a form similar to that of read_keithley.m and
0095 %     read_mate.m, containing fields for data streams, events, and onsets,
0096 %     also containing a .meta field, with .meta.srate (sampling rate),
0097 %     .meta.stime (sampling time, or # ms per sample), and
0098 %     .meta.baseline_samps (# of baseline samples). This differs from that
0099 %     of read_keithley and read_mate in that the actual signal fields are
0100 %     found in ri.signal, rather than within ri
0101 %       NOTE: ri.trigger values are sample indices, NOT times in ms or sec
0102 %
0103 % FB 2008.10.12
0104 % FB 2009.06.24 - added baseline calculation
0105 
0106 ri = struct('signal',struct(),'meta',struct());
0107 global r
0108 
0109 r = init_results_struct;
0110 r.type = 'physio_biopac';  % Identify the type of this reporting instance
0111 r.report_on_fly = 1;
0112 
0113 if ~isempty(fname) && exist(fname,'file')
0114   % read data
0115   pdata = load_acq(fname);
0116   chans = {pdata.hdr.per_chan_data(:).comment_text};
0117 else
0118   % data can't be found, skip this file
0119   msg = sprintf('couldn''t find biopac file %s, SKIPPING\n',fname);
0120   r = update_report(r,msg);
0121   return
0122 end
0123 
0124 % get sampling rate
0125 stime = pdata.hdr.graph.sample_time; % msec per sample
0126 srate = 1000/stime; % 1000 msec / stime = sampling rate
0127 ri.meta.srate = srate;
0128 ri.meta.stime = stime;
0129 
0130 % get TR
0131 if isfield(params,'TR')
0132   TR = params.TR;
0133 else
0134   TR = 0;
0135   warning('no TR found in the params struct\n');
0136 end
0137 
0138 % signal bounds?
0139 if isfield(params,'signal_bounds') && length(params.signal_bounds) == 2
0140   % convert from params.signal_bounds in secs to sig_start/end in samples
0141   sig_start = params.signal_bounds(1)*srate;
0142   sig_end = params.signal_bounds(2)*srate;
0143   if ~isnumeric(sig_start) || sig_start < 1
0144     sig_start = 1;
0145   end
0146   if ~isnumeric(sig_end) || sig_end < 1 || sig_end > size(pdata.data,1)
0147     sig_end = size(pdata.data,1);
0148   end
0149 else
0150   sig_start = 1;
0151   sig_end = size(pdata.data,1);
0152 end
0153 
0154 % find channels, get data
0155 ctxts = {pdata.hdr.per_chan_data.comment_text};
0156 if isfield(params,'channels')
0157   nc = length(params.channels);
0158   for ic=1:nc
0159     ctxt = params.channels(ic).comment_text;
0160     name = params.channels(ic).name;
0161     cidx = strmatch(ctxt,ctxts);
0162     if ~isempty(cidx)
0163       ri.signal.(name) = pdata.data(sig_start:sig_end,cidx);
0164     end
0165   end
0166 end
0167 
0168 % find trigger onsets, if they exist
0169 if isfield(ri.signal,'trigger')
0170   % onsets and offsets are in samples, NOT time in msec or sec
0171   [ri.trigger.onsets,ri.trigger.offsets] = ...
0172       find_thresh_cross(ri.signal.trigger,params);
0173 end
0174 
0175 % extract run-length epochs
0176 % epoch_start and epoch_end are expressed in samples, NOT in time
0177 if isfield(ri,'trigger') && isfield(ri.trigger,'onsets') ...
0178       && ~isempty(ri.trigger.onsets) && isfield(params,'bound_by_trigger') ...
0179       && params.bound_by_trigger
0180   % check for TR
0181   if ~exist('TR','var') || ~TR
0182     warning('estimating TR from trigger onsets');
0183     TR = median(diff(ri.trigger.onsets*stime))/1000;
0184     warning('TR estimated to be %ds',TR);
0185   end
0186 
0187   % bound by triggers
0188   epoch_start = ri.trigger.onsets(1);
0189   epoch_end = ri.trigger.onsets(end)+TR*srate; % add one TR at the end
0190 
0191 elseif isfield(params,'link2pres') ...
0192         && isfield(params.link2pres,'channelname')...
0193         && ischar(params.link2pres.channelname) ...
0194         && isfield(ri.signal,params.link2pres.channelname)
0195     
0196   % extract signal
0197   linksig = ri.signal.(params.link2pres.channelname);
0198   [lsons,lsoffs] = find_thresh_cross(linksig,params);
0199   lsdiff_ms = diff(lsons)*stime;
0200   nlsdiff = length(lsdiff_ms);
0201   
0202   % link to presentation file, bound by presentation data
0203   pfname = params.link2pres.presfname;
0204   presdata = load(pfname);
0205   pcol = set_var_col_const(presdata.vars);
0206   
0207   if isfield(params.link2pres,'filt')
0208     pfilt = params.link2pres.filt;
0209     linkdata = ensemble_filter(presdata,pfilt);
0210   else
0211     linkdata = presdata;
0212   end
0213   pdiff_ms = diff(linkdata.data{pcol.RUN_REL_TIME});
0214   npdiff = length(pdiff_ms);
0215 
0216   if isfield(params.link2pres,'diff_tol')
0217     diff_tol = params.link2pres.diff_tol;
0218   else
0219     diff_tol = 100; % 100ms tolerance
0220   end
0221   nlsd = length(lsdiff_ms);
0222   npd  = length(pdiff_ms);
0223 
0224   % the first two conditionals within this if/else train are a hack to
0225   % account for an additional leading or trailing signal in either linksig
0226   % or presdata, that doesn't appear in the other ... essentially, if
0227   % linksig has items 1 2 3 and pdata has 2 3 or 1 2 (i.e. pdata is missing
0228   % a leading or trailing element, but otherwise is intact), or vice versa
0229   % (linksig is missing an element), we can still link the two measurements
0230   % together
0231   
0232 %   FIXME: right now, d3 = d1>d2, and the sum of # of indices where one is
0233 %   greater than the other is used to figure out where 'less difference'
0234 %   exists, but this is not really a metric of 'less difference', rather
0235 %   'indices with less difference' ... might another formulation work
0236 %   better? how about d4 = d1-d2, then compare abs(sum(d4(d3))) to
0237 %   abs(sum(d4(~de))) ... ? if problems arise, I will try this strategy
0238   if (nlsd - npd) == 1
0239     % one more event in the link signal than in the presentation file
0240     
0241     % assume the first pres signal is missing
0242     d1 = abs(diff([lsdiff_ms'; [lsdiff_ms(1); pdiff_ms]']));
0243     
0244     % assume the last pres signal is missing
0245     d2 = abs(diff([lsdiff_ms'; [pdiff_ms; lsdiff_ms(end)]']));
0246     
0247     % greater diff for first signal compared to last signal assumption?
0248     d3 = d1 > d2;
0249     
0250     if (sum(d3) < nlsd/2)
0251       % smaller difference between link sig and pres times when assuming
0252       % the first pres time is missing
0253       % in other words, linksig(2) = prestime(1)
0254       if ~any(d1 > diff_tol)
0255         % no differences greater than diff_tol, safe to assume the first
0256         % pres time is missing
0257         pbegin = linkdata.data{pcol.RUN_REL_TIME}(1); % 1st pres time
0258         % check to see that the first epoch time is not before the
0259         % beginning of the run as presentation has calculated it
0260         if lsons(2)*stime - lsons(1)*stime > pbegin
0261           msg = sprintf(['possible error calculating the beginning of '...
0262               'the epoch. There are more link signal events %d than '...
0263               'presentation linking events %d, and it appears that the '...
0264               'first presentation event is missing, but if this is the '...
0265               'case, the first link signal event is before the start '...
0266               'of the run, when comparing event times to presentation '...
0267               'run-relative times'],nlsd,npd);
0268           r = update_report(r,msg);
0269           ri = struct();
0270           return
0271         else
0272           % start of the epoch is the onset of the second link signal
0273           % event, minus the time of the first (but really, the second)
0274           % presentation linking event, converted to link signal samples
0275           epoch_start = round((lsons(2)*stime - pbegin)/stime);
0276         end
0277       else
0278         % differences greater than diff_tol, likely a middle time that is
0279         % missing
0280         msg = sprintf('missing presentation link, can not compute');
0281         r = update_report(r,msg);
0282         ri = struct();
0283         return
0284       end
0285     elseif (sum(d3) > nlsd/2)
0286       % smaller difference between link sig and pres times when assuming
0287       % the last pres time is missing
0288       % in other words, linksig(1) = prestime(1)
0289       if ~any(d2 > diff_tol)
0290         % no differences greater than diff_tol, safe to assume the last
0291         % pres time is missing
0292         pbegin = linkdata.data{pcol.RUN_REL_TIME}(1); % 1st pres time
0293         % start of the epoch is the onset of the first link signal event,
0294         % minus the time of the first presentation linking event, converted
0295         % to link signal samples
0296         epoch_start = round((lsons(1)*stime - pbegin)/stime);
0297       else
0298         % differences greater than diff_tol, likely a middle time that is
0299         % missing
0300         msg = sprintf('missing presentation link, can not compute');
0301         r = update_report(r,msg);
0302         ri = struct();
0303         return
0304       end
0305     else
0306       msg = sprintf('missing presentation link, can not compute');
0307       r = update_report(r,msg);
0308       ri = struct();
0309       return
0310     end % if (sum(d3) < nlsd/2
0311   elseif (npd - nlsd) == 1
0312     % one more event in the presentation file than in the link signal
0313     
0314     % assume the first link signal is missing
0315     d1 = abs(diff([pdiff_ms'; [pdiff_ms(1); lsdiff_ms]']));
0316     
0317     % assume the last link signal is missing
0318     d2 = abs(diff([pdiff_ms'; [lsdiff_ms; pdiff_ms(end)]']));
0319     
0320     % greater diff for first signal compared to last signal assumption?
0321     d3 = d1 > d2;
0322     
0323     if (sum(d3) < npd/2)
0324       % smaller difference between link sig and pres times when assuming
0325       % the first link signal is missing
0326       % in other words, prestime(2) == linksig(1)
0327       if ~any(d1 > diff_tol)
0328         % no differences greater than diff_tol, safe to assume the first
0329         % link signal is missing
0330         pbegin = linkdata.data{pcol.RUN_REL_TIME}(2); % 2nd pres time
0331         % check to see that the beginning of physio data collection is
0332         % before the beginning of the run as presentation has logged it
0333         if lsons(1)*stime < pbegin
0334           msg = sprintf(['possible error calculating the beginning of '...
0335               'the epoch. There are more presentation linking events %d'...
0336               ' than link signal events %d, and it appears that the '...
0337               'first link signal event is missing, but if this is the '...
0338               'case, the first presentation linkin event is before the '...
0339               'start of physio data collection, when comparing physio '...
0340               'data to presentation run-relative times'],npd,nlsd);
0341           r = update_report(r,msg);
0342           ri = struct();
0343           return
0344         else
0345           % start of the epoch is the onset of the first link signal
0346           % event, minus the time of the second presentation linking event,
0347           % converted to link signal samples
0348           epoch_start = round((lsons(1)*stime - pbegin)/stime);
0349         end
0350       else
0351         % differences greater than diff_tol, likely a middle time that is
0352         % missing
0353         msg = sprintf('missing link signal, can not compute');
0354         r = update_report(r,msg);
0355         ri = struct();
0356         return
0357       end
0358     elseif (sum(d3) > npd/2)
0359       % smaller difference between link sig and pres times when assuming
0360       % the last link signal is missing
0361       % in other words, linksig(1) = prestime(1)
0362       if ~any(d2 > diff_tol)
0363         % no differences greater than diff_tol, safe to assume the last
0364         % link signal is missing
0365         pbegin = linkdata.data{pcol.RUN_REL_TIME}(1); % 1st pres time
0366         % start of the epoch is the onset of the first link signal event,
0367         % minus the time of the first presentation linking event, converted
0368         % to link signal samples
0369         epoch_start = round((lsons(1)*stime - pbegin)/stime);
0370       else
0371         % differences greater than diff_tol, likely a middle time that is
0372         % missing
0373         msg = sprintf('missing link signal, can not compute');
0374         r = update_report(r,msg);
0375         ri = struct();
0376         return
0377       end
0378     else
0379       msg = sprintf('missing link signal, can not compute');
0380       r = update_report(r,msg);
0381       ri = struct();
0382       return
0383     end
0384   elseif npd ~= nlsd
0385     msg = sprintf(['link signal size (%d) is not equal in length to '...
0386         'presentation data size (%d), therefore link2pres skipped\n'],...
0387         nlsdiff,npdiff);
0388     r = update_report(r,msg);
0389     ri = struct();
0390     return
0391   elseif any(diff([lsdiff_ms'; pdiff_ms']) > diff_tol)
0392     msg = sprintf(['more than %d ms diff btwn link signal onset diffs '...
0393         'and presentation data onset diffs, therefore link2pres skipped\n'],...
0394         diff_tol);
0395     r = update_report(r,msg);
0396     ri = struct();
0397     return
0398   else
0399     pbegin = linkdata.data{pcol.RUN_REL_TIME}(1);
0400     epoch_start = round((lsons(1)*stime - pbegin)/stime);
0401   end
0402 
0403   epoch_end = length(linksig);
0404   
0405   if isfield(params.link2pres,'bound_by_last_pulse') || ...
0406           isfield(params.link2pres,'estimate_trigger')
0407     lpfilt.include.all.EVENT_TYPE={'Pulse'};
0408     lpfilt.include.all.EVENT_CODE={'255'};
0409     lpfilt.include.all.RUN=params.run;
0410     lpdata = ensemble_filter(presdata,lpfilt);
0411     rrtimes = lpdata.data{pcol.RUN_REL_TIME};
0412     
0413     if ~isempty(lpdata.data{1})
0414       if ~exist('TR','var') || ~TR
0415         warning('estimating TR from pulse data');
0416         TR = median(diff(rrtimes))/1000;
0417         warning('TR estimated to be %ds',TR);
0418       end
0419       
0420       if isfield(params.link2pres,'estimate_trigger')
0421         cpp = check_pulse_periods(rrtimes);
0422         pptl = find(cpp.data.pulse_period_too_long_mask);
0423         for ipp=1:length(pptl)
0424           intpp = rrtimes(pptl(ipp)):(TR*1000):rrtimes(pptl(ipp)+1);
0425           rrtimes = union(rrtimes,intpp);
0426         end
0427         samples = rrtimes/stime;
0428         ri.trigger.onsets = samples;
0429         ri.trigger.offsets = samples+1;
0430       end
0431       
0432       if isfield(params.link2pres,'bound_by_last_pulse')
0433         % the end of the epoch = run_rel_time/stime + epoch_start
0434         pend = ceil((lpdata.data{pcol.RUN_REL_TIME}(end)+TR*1000)/stime);
0435         epoch_end = pend + epoch_start;
0436       end
0437     else
0438       warning('no pulses found, no trigger estimation or bounding');
0439     end
0440   end
0441 
0442 else
0443   warning('could not generate run-length epochs');
0444   return
0445 end
0446 
0447 % add baseline period?
0448 if isfield(params,'baseline')
0449     blstart = epoch_start - params.baseline*srate;
0450     if blstart < 1, blstart = 1; end
0451     ri.meta.baseline_samps = epoch_start - blstart;
0452     epoch_start = blstart;
0453 else
0454     ri.meta.baseline_samps = 0;
0455 end
0456 
0457 % apply bounds
0458 flds = fieldnames(ri.signal);
0459 nflds = length(flds);
0460 for ifld = 1:nflds
0461     fld = flds{ifld};
0462     nsamp = length(ri.signal.(fld));
0463     if epoch_start > nsamp
0464       continue
0465     end
0466     if epoch_end > nsamp
0467       epoch_end = nsamp;
0468     end
0469     ri.signal.(fld) = ri.signal.(fld)(epoch_start:epoch_end);
0470 end
0471 
0472 % adjust triggers?
0473 if isfield(ri,'trigger') && isfield(ri.trigger,'onsets') ...
0474       && ~isempty(ri.trigger.onsets) && isfield(params,'bound_by_trigger') ...
0475       && params.bound_by_trigger
0476   ri.trigger.onsets = ri.trigger.onsets + ri.meta.baseline_samps;
0477   ri.trigger.offsets = ri.trigger.offsets + ri.meta.baseline_samps;
0478 end

Generated on Wed 20-Sep-2023 04:00:50 by m2html © 2003