check_pulse_periods

 [r] = check_pulse_periods(pulse_times,params)

 Given a vector of event onset times, e.g. MRI scanner pulses, that are
 expected to have a constant period, the script identifies aberrent pulse periods

 pulse_times - vector of pulse times (msec)

 r is the reporting structure which contains a data structure with the
 following fields:
 bad_pulse_period_mask - logical vector indicating which pulse periods are
                         aberrant.  The first element in this vector refers to
                         the period between the first and second elements in
                         the pulse_times vector.
 pulse_period_too_short_mask - those periods that were too short
 pulse_period_too_long_mask - those periods that were too long
 nmiss - estimate of number of missed pulses
 median_period - estimate of the median pulse period
 msglog - list of all messages that were generated

0001 function [r] = check_pulse_periods(pulse_times,params,r)
0002 % [r] = check_pulse_periods(pulse_times,params)
0003 %
0004 % Given a vector of event onset times, e.g. MRI scanner pulses, that are
0005 % expected to have a constant period, the script identifies aberrent pulse periods
0006 %
0007 % pulse_times - vector of pulse times (msec)
0008 %
0009 % r is the reporting structure which contains a data structure with the
0010 % following fields:
0011 % bad_pulse_period_mask - logical vector indicating which pulse periods are
0012 %                         aberrant.  The first element in this vector refers to
0013 %                         the period between the first and second elements in
0014 %                         the pulse_times vector.
0015 % pulse_period_too_short_mask - those periods that were too short
0016 % pulse_period_too_long_mask - those periods that were too long
0017 % nmiss - estimate of number of missed pulses
0018 % median_period - estimate of the median pulse period
0019 % msglog - list of all messages that were generated
0020 
0021 % 05/12/06 Petr Janata - updated code with additional error checking
0022 % 08/03/06 PJ - implemented within the emerging reporting framework
0023 % 12/11/07 PJ - fixed minor bug in which nmiss was not being created when
0024 %               no pulses were missed
0025 
0026 r.type = 'pulse_period_check';  % Identify the type of this reporting instance
0027 r.report_on_fly = 1;
0028 
0029 try max_pulse_dev = params.max_pulse_dev; catch
0030   max_pulse_dev = 20; % maximum deviation (in msec) of a pulse from the expected
0031                      % period.
0032 end
0033 
0034 try min_pulse_period = params.min_pulse_period; catch
0035   min_pulse_period = 50;  % won't be going faster than 20 slices/sec
0036 end
0037 
0038 % Look at variability of timing pulses, and check for skipped pulses
0039 pulse_periods = diff(pulse_times);  % determine all of the pulse periods
0040 median_period = median(pulse_periods);
0041 
0042 % Check to see if the median pulse period is too short
0043 pulse_period_too_short_mask = zeros(size(pulse_periods));
0044 if median_period < min_pulse_period
0045   msg = sprintf(['Median pulse period (%4.1f) is less than minimum allowed ' ...
0046     'period (%4.1f)'], median_period, min_pulse_period)
0047   r = update_report(r,msg);
0048  
0049   pulse_period_too_short_mask = pulse_periods < min_pulse_period;
0050   msg = sprintf('Detected %d pulse periods that were too short\n', sum(pulse_period_too_short_mask));
0051   r = update_report(r,msg);
0052   
0053   % Revise the median pulse period estimate
0054   median_period = median(pulse_periods(~pulse_period_too_short_mask));
0055 end
0056 
0057 % Check for pulse periods that are too long
0058 pulse_period_too_long_mask = pulse_periods > median_period+max_pulse_dev;
0059 bad_pulse_period_idxs = find(pulse_period_too_long_mask);
0060 
0061 nbad = length(bad_pulse_period_idxs);
0062 nmiss = 0;
0063 if nbad
0064   msg = sprintf('Detected %d pulse periods that were too long\n', nbad);
0065   r = update_report(r,msg);
0066 
0067   % Estimate the number of missed pulses
0068   for ibad = 1:nbad
0069     nmiss = nmiss + round(pulse_periods(bad_pulse_period_idxs(ibad))/median_period)-1;
0070   end % for ibad=
0071   msg = sprintf('Estimate of %d missed pulses\n', nmiss);
0072   r = update_report(r,msg);
0073 end
0074 
0075 % Pulse periods that are too long aren't technically bad in that there are only
0076 % pulses missing.  Short pulse periods are bad though, so add them to the bad_pulse_period_mask.
0077 r.data.bad_pulse_period_mask = pulse_period_too_short_mask; % | pulse_period_too_long_mask;
0078 r.data.pulse_period_too_short_mask = pulse_period_too_short_mask;
0079 r.data.pulse_period_too_long_mask = pulse_period_too_long_mask;
0080 r.data.nmiss = nmiss;
0081 r.data.median_period = median_period;