parse_event_wav

 [resp, events, slices] = parse_event_wav(fname);

 Reads in a .wav file that contains stimulus and magnet event information.
 The default arrangement assumes that left and right button presses are on the
 left channel, and left button signals are smaller in magnitude than right
 button signals.  Magnet trigger and event onset information, together with
 the receiver unblank (slice timing) signals are mixed on the right.  Trigger
 and onset information are larger in magnitude than the receiver unblank.

 Button press onsets are negative going at start

0001 function [resp, events, slices] = parse_event_wav(fname,p)
0002 %
0003 % [resp, events, slices] = parse_event_wav(fname);
0004 %
0005 % Reads in a .wav file that contains stimulus and magnet event information.
0006 % The default arrangement assumes that left and right button presses are on the
0007 % left channel, and left button signals are smaller in magnitude than right
0008 % button signals.  Magnet trigger and event onset information, together with
0009 % the receiver unblank (slice timing) signals are mixed on the right.  Trigger
0010 % and onset information are larger in magnitude than the receiver unblank.
0011 %
0012 % Button press onsets are negative going at start
0013 %
0014 
0015 % Originally written for prime_long
0016 
0017 % 02/28/01 PJ Started adapting to be useful for other experiments, e.g. schub3
0018 %             The function now only returns detected events on the two
0019 %             channels, without any further elimination of events that might
0020 %             occur on block boundaries.  That process is left to individual
0021 %             scripts.
0022 %
0023 %             Implemented correct handling of button and event polarity
0024 
0025 error(nargchk(1,2,nargin))
0026   
0027 LOAD_DATA = 1;
0028 
0029 if nargin < 2
0030   p = check_pinfo;                % check and fill pinfo
0031                                                 % structure
0032 else
0033   p = check_pinfo(p);
0034 end
0035 
0036 predicted_num_slices = p.slice_per_vol*(p.nvol+p.ndummy);
0037 
0038 % Types of information on the different channels
0039 info_type = {'button_press','stim_and_magnet'};
0040 
0041 violation = [];
0042 RESP = 1;
0043 EVENT = 2;
0044 
0045 LEFT = 1;
0046 RIGHT = 2;
0047 
0048 if LOAD_DATA
0049   disp(sprintf('\nLoading data from file: %s', fname))
0050   [data, Fs, nbits, opts] = wavread(fname);
0051 end
0052 
0053 nchan = size(data,2);
0054 
0055 min_bp_samps = p.min_bp_isi/1000*Fs;
0056 min_event_trig_samps = p.min_event_trig_isi/1000*Fs;
0057 
0058 interslice_samps = p.interslice_intrvl/1000*Fs;
0059 
0060 % Assume that the initial 1 sec contains no events, and use this to get
0061 % estimate of mean and std of signal
0062 
0063 basemean = mean(data(1:Fs,:));
0064 basestd = std(data(1:Fs,:));
0065 
0066 % Determine the maximum signal amplitude on each channel
0067 
0068 if p.verbosity > 1
0069   disp(sprintf('Getting data min and max'))
0070 end
0071 [chanmax, chanmax_samps] = max(data);
0072 [chanmin, chanmin_samps] = min(data);
0073 
0074 if p.verbosity > 1
0075   for ichan = 1:size(data,2)
0076     disp(sprintf('Chan(%d): Minimum: %2.2f; Maximum: %2.2f', ichan,chanmin(ichan), chanmax(ichan)))
0077   end
0078 end
0079 
0080 %
0081 %  HANDLE THE BUTTON PRESS CHANNEL
0082 %
0083 
0084 % Set thresholds for button presses
0085 if p.button_polarity == -1
0086   thresh(1).high = 0.66*chanmin(1);    % right
0087   thresh(1).low = 0.25*chanmin(1);    % left
0088 else
0089   if isfield(p,'button_threshhigh')
0090     thresh(1).high = p.button_threshhigh*chanmax(1);
0091   else
0092     thresh(1).high = 0.66*chanmax(1);    % right
0093   end
0094   if isfield(p,'button_threshlow')
0095     thresh(1).low = p.button_threshlow*chanmax(1);    % left
0096   else
0097     thresh(1).low = 0.25*chanmax(1);    % left
0098   end    
0099 end
0100 
0101 ichan = 1;
0102 
0103 if p.verbosity > 1
0104   disp(sprintf('Processing response data on channel %d', ichan))
0105   disp(sprintf('Finding samples that exceed high (%1.2f) and low (%1.2f) thresholds', thresh(ichan).high, thresh(ichan).low))
0106 end
0107 
0108 switch info_type{ichan}
0109   case 'button_press'
0110     min_isi_samps = min_bp_samps;
0111   otherwise
0112     disp(sprintf('Do not know how to handle information type: %s', info_type{ichan}))
0113 end
0114 
0115 % Find all samples that exceed right and left button press threshold
0116 if p.button_polarity == -1
0117   threshsamps_high = find(data(:,ichan) <= thresh(ichan).high);%right
0118   threshsamps_low = find(data(:,ichan) <= thresh(ichan).low); % left
0119 else
0120   threshsamps_high = find(data(:,ichan) >= thresh(ichan).high);
0121   threshsamps_low = find(data(:,ichan) >= thresh(ichan).low); % left
0122 end
0123 
0124 % Get the first sample of each threshold crossing
0125 cross_samps_high = threshsamps_high([1 (find(diff(threshsamps_high)>1)+1)']);
0126 
0127 % Find all samps that exceed left button press threshold
0128 cross_samps_low = threshsamps_low([1 (find(diff(threshsamps_low)>1)+1)']);
0129 
0130 % Find all samps in low threshold vector that correspond to events in high
0131 % threshold vector.  Allow several samples of slop in threshold determination
0132 
0133 search_range = repmat(cross_samps_high',length(-p.perisamps:1:p.perisamps),1) + ...
0134     repmat((-p.perisamps:1:p.perisamps)',1,length(cross_samps_high));
0135 
0136 search_range = search_range(:);
0137 resp.samps = cross_samps_low;
0138 
0139 resp.time = cross_samps_low/Fs;
0140 
0141 rt_idx = find(ismember(cross_samps_low,search_range));
0142 lft_idx = find(~ismember(cross_samps_low,search_range));
0143 resp.id(rt_idx) = RIGHT;
0144 resp.id(lft_idx) = LEFT;
0145 
0146 disp(sprintf('Found a total of %d threshold crossings: RIGHT: %d; LEFT: %d', length(resp.samps), length(rt_idx), length(lft_idx)))
0147 
0148 % Make sure that no putative events occur within the minimum button press ISI
0149 tmpdiff = diff(resp.samps);
0150 bad_idx = find(tmpdiff < min_bp_samps)+1;
0151 num_bad = length(bad_idx);  
0152 
0153 % Eliminate overlapping button press events
0154 if num_bad
0155   if p.verbosity
0156     disp(sprintf('Found %d overlapping button press events', num_bad))
0157     disp('Eliminating ...')
0158   end
0159   resp.samps(bad_idx) = [];
0160   resp.time(bad_idx) = [];
0161   resp.id(bad_idx) = [];
0162 end
0163 
0164 % Do some sanity checks
0165 if length(resp.samps) ~= p.expect_num_resp
0166   if p.verbosity
0167     disp(sprintf('     WARNING: Number of detected RESPONSES (%d) did not match expected number (%d)', length(resp.samps), p.expect_num_resp))  
0168   end
0169   violation(RESP) = 1;
0170 end
0171 
0172 if nargout == 1
0173   return
0174 end
0175 
0176 %
0177 %  HANDLE THE EVENT ONSET AND MAGNET TRIGGER CHANNEL
0178 %
0179 
0180 ichan = 2;
0181 
0182 if p.verbosity > 1
0183   disp(sprintf('\nProcessing stimulus event data on channel %d', ichan))
0184 end
0185 
0186 %
0187 %  Deal with the event signals first
0188 %
0189 
0190 min_isi_samps = min_event_trig_samps;
0191 
0192 % Set threshold for triggers
0193 
0194 if p.event_thresh & p.event_polarity
0195   thresh(ichan).high = p.event_thresh * chanmax(ichan) * p.event_polarity;
0196 else
0197   thresh(ichan).high = 0.5;
0198 end
0199 
0200 if p.verbosity > 1
0201   disp(sprintf('Finding samples that exceed high (%1.2f) threshold', thresh(ichan).high))
0202 end
0203 
0204 if p.event_polarity == -1
0205   threshsamps_high = find(data(:,ichan) <= thresh(ichan).high);
0206 else
0207   threshsamps_high = find(data(:,ichan) >= thresh(ichan).high);
0208 end
0209 
0210 cross_samps_high = threshsamps_high([1 (find(diff(threshsamps_high)>1)+1)']);
0211 bad_idx = find(diff(cross_samps_high(:)) < min_isi_samps);
0212 
0213 disp(sprintf('Found a total of %d threshold crossings', length(cross_samps_high)));
0214 
0215 num_bad = length(bad_idx);
0216 
0217 if num_bad
0218   disp(sprintf('Found %d bad event triggers.  Eliminating ...', num_bad))
0219   cross_samps_high(bad_idx) = [];
0220 end
0221 
0222 %
0223 % Realign everything to first event (assuming this is the magnet trigger
0224 %
0225 
0226 if ~isfield(p,'assume_first_is_magtrig')
0227   p.assume_first_is_magtrig = 1;
0228 end
0229 
0230 if p.assume_first_is_magtrig
0231   if p.verbosity > 1
0232     disp(['   Stripping first event (run trigger) and calculating event times' ...
0233       ' relative to it ...'])
0234   end
0235 
0236   run_startsamp = cross_samps_high(1);
0237   event_samps = cross_samps_high(2:end)-run_startsamp;
0238 
0239   % Remove sample offset due to dummy shots
0240   ndummy_samps = p.ndummy*2*Fs;
0241   event_samps = event_samps - ndummy_samps;
0242   events.samps = event_samps;
0243   events.time = events.samps/Fs;
0244 
0245   % Subtract start offset from reaction time data
0246   resp.samps = resp.samps-run_startsamp - ndummy_samps;
0247   resp.time = resp.samps/Fs;
0248 
0249 else
0250   events.samps = cross_samps_high;
0251   events.time = events.samps/Fs;
0252 end
0253 
0254 return
0255 
0256 %
0257 % DO NOT DO ANY FURTHER PROCESSING OF EVENT DATA.  NOW THAT EVENTS HAVE BEEN
0258 % IDENTIFIED, LET CALLING SCRIPT TO FURTHER PROCESSING.
0259 %
0260 
0261 
0262 % Eliminate the events corresponding to noise bursts and beginning and end of
0263 % block boundaries
0264 
0265 if p.verbosity > 1
0266   disp('   Eliminating noise burst markers ...')
0267 end
0268 tmp = find(diff(event_samps) > 20*Fs)+1;
0269 tmp = [tmp tmp-1];
0270 tmp = tmp(:);
0271 
0272 if max(tmp) ~= length(event_samps)
0273   event_samps([1 tmp' length(event_samps)]) = [];
0274 else
0275   event_samps([1 tmp']) = [];
0276 end
0277 
0278 events.samps = event_samps;
0279 events.time = events.samps/Fs;
0280 
0281 % Check to see if number of event samps matches the expected number
0282 if length(event_samps) ~= p.expect_num_events
0283   if p.verbosity
0284     disp(sprintf('     WARNING: Number of detected EVENTS (%d) did not match expected number (%d)', length(event_samps), p.expect_num_events))
0285   end
0286   violation(EVENT) = 1;
0287 end
0288 
0289 if nargout < 3
0290   return
0291 end
0292 
0293 %
0294 %  Deal with the slice timing info
0295 %
0296 
0297 thresh(ichan).low = 0.2 * chanmin(ichan); % threshold for slices
0298 
0299 threshsamps_low = find(data(:,ichan) <= thresh(ichan).low);
0300 cross_samps_low = threshsamps_low([1 (find(diff(threshsamps_low)>1)+1)']);
0301 
0302 % Get rid of first threshold crossing as this is the off edge of the start
0303 % trigger
0304 
0305 cross_samps_low(1) = [];
0306 
0307 % Toss out any negative crossings that occur earlier than known slice period
0308 
0309 tmp = [cross_samps_low cross_samps_low+interslice_samps-p.interslice_slop_samps]';
0310 tmpdiff = diff(tmp(:));
0311 
0312 badidx = find(tmpdiff(2:2:end) < 0)+1;
0313 longidx = find(tmpdiff(2:2:end) > (2*p.interslice_slop_samps)) ...
0314     + 1;
0315 
0316 % Some of the samples in badidx are actually good.  This can be determined by
0317 % the proximity of indices given in badidx.  If they are adjacent samples,
0318 % chances are very good that the first one represents the off-edge of an event
0319 % trigger and the following one represents the onset of the next slice.
0320 % Therefore, if two indices differ by only 1, then toss the first.
0321 % Sometimes the off-edge of an event trigger masks the onset of the next slice,
0322 % and the delay until the next slice onset is longer than the period between
0323 % slice onsets. In these cases, slice onsets have to be inferred and inserted
0324 % into the sequence.
0325 
0326 % Handle the case of missed slice onsets first.  The simplest fix is to replace
0327 % the sample numbers of the extra short intervals with corrected indices.
0328 
0329 if ~all(ismember(longidx-1,badidx))
0330   if p.verbosity
0331     disp('WARNING: Extra long duration not preceded by extra short duration')
0332   end
0333 end
0334 
0335 cross_samps_low(longidx-1) = mean(cross_samps_low([longidx-2 longidx]),2);
0336 badidx(find(ismember(badidx,longidx-1))) = [];
0337 
0338 badidx(find(diff(badidx)==1)+1) = [];
0339 
0340 cross_samps_low(badidx) = [];        % remove bad samples
0341 
0342 % Reiterate process to toss out any negative crossings that occur earlier than known slice period
0343 tmp = [cross_samps_low cross_samps_low+interslice_samps-p.interslice_slop_samps]';
0344 tmpdiff = diff(tmp(:));
0345 
0346 badidx = find(tmpdiff(2:2:end) < 0)+1;
0347 badidx(find(diff(badidx)==1)+1) = [];
0348 cross_samps_low(badidx) = [];        % remove bad samples
0349 
0350 if predicted_num_slices ~= length(cross_samps_low)
0351   error(sprintf('Predicted number of slices (%d) does not match number found (%d)', predicted_num_slices, length(cross_samps_low)))
0352 end
0353 
0354 function p = check_pinfo(p)
0355   if nargin < 1
0356     p = struct(...
0357     'verbosity', [], ...
0358     'min_bp_isi',[], ...
0359     'min_event_trig_isi',[], ...
0360     'perisamps', [], ...
0361     'interslice_slop_samps', [], ...
0362     'interslice_intrvl', [], ...
0363     'slice_per_vol', [], ...
0364     'nvol', [], ...
0365     'ndummy', [], ...
0366     'event_polarity', [], ...
0367     'event_thresh', [], ...
0368     'expect_num_events', [], ...
0369     'slice_polarity', [], ...
0370     'button_polarity' ...
0371     );
0372   end
0373   
0374   if isempty(p.verbosity)
0375     p.verbosity = 1;
0376   end
0377   
0378   if isempty(p.min_bp_isi)
0379     p.min_bp_isi = 2000;            % minimum duration between button presses (ms)
0380   end
0381   
0382   if isempty(p.min_event_trig_isi)
0383     p.min_event_trig_isi = 1000;
0384   end
0385   
0386   if isempty(p.perisamps)
0387     p.perisamps = 5;
0388   end
0389   
0390   if isempty(p.interslice_slop_samps)
0391     p.interslice_slop_samps = 10;
0392   end
0393 
0394   if isempty(p.interslice_intrvl)
0395     p.interslice_intrvl = 74;            % time between slice onsets (ms)
0396   end
0397                     
0398   if isempty(p.slice_per_vol)
0399     p.slice_per_vol = 27;
0400   end
0401   
0402   if isempty(p.nvol)
0403     p.nvol = 184;
0404   end
0405   
0406   if isempty(p.ndummy)
0407     p.ndummy = 4;
0408   end
0409   
0410   if isempty(p.event_polarity)
0411     p.event_polarity = 1;
0412   end
0413   
0414   if isempty(p.event_thresh)
0415     p.event_thresh = 0.4;
0416   end
0417   
0418   if isempty(p.expect_num_events)
0419     p.event_thresh = 0;
0420   end
0421   
0422   if isempty(p.slice_polarity)
0423     p.slice_polarity = -1;
0424   end
0425   
0426   if isempty(p.button_polarity)
0427     p.button_polarity = -1;
0428   end
0429   
0430  % end check_pinfo