read_physmon

 [runinfo, params] = read_keithley(fname,params,data)

 Reads physiological data monitoring file.
 Note, these files do not contain information about the sampling rate, so this
 has to be specified by the user.
 The number of columns is equal to the number of channels, though the channel
 IDs are not specified.  

 Specific channel assignments are assumed:
    1:  Cardiac
    2:  Respiration
    3:  Slice timing (receiver unblank)
    4:  Button presses
    5:  Event markers

 Channel information should really be more dynamic and user-specifiable, but
 I'll deal with this when the need arises.

 If data has already been loaded once, it can be passed into the function
 without having to be loaded again.

 Note:  The data file can have multiple runs

0001 function [ri, p] = read_keithley(fname,p,data)
0002 % [runinfo, params] = read_keithley(fname,params,data)
0003 %
0004 % Reads physiological data monitoring file.
0005 % Note, these files do not contain information about the sampling rate, so this
0006 % has to be specified by the user.
0007 % The number of columns is equal to the number of channels, though the channel
0008 % IDs are not specified.
0009 %
0010 % Specific channel assignments are assumed:
0011 %    1:  Cardiac
0012 %    2:  Respiration
0013 %    3:  Slice timing (receiver unblank)
0014 %    4:  Button presses
0015 %    5:  Event markers
0016 %
0017 % Channel information should really be more dynamic and user-specifiable, but
0018 % I'll deal with this when the need arises.
0019 %
0020 % If data has already been loaded once, it can be passed into the function
0021 % without having to be loaded again.
0022 %
0023 % Note:  The data file can have multiple runs
0024 
0025 % 1/10/02 PJ
0026 % 11/10/02 PJ Accommodated files without behavioral or stim event data
0027 % 05/27/04 PJ Accommodating different polarity of event and button press
0028 %             channels
0029 % 08/31/04 PJ Added optional removal of DC offset from event channels. Added
0030 %             handling of negative key press events, including separate trigger
0031 %             threshold
0032 % 09/02/04 PJ Added option of removing events exceeding a threshold run duration
0033 
0034 ri = '';
0035 
0036 if nargin < 2
0037   error('read_keithley: Too few input arguments')
0038 end
0039 
0040 if p.bipolar
0041   error('Bipolar event marker option currently unsupported')
0042 end
0043 
0044 if p.loaddata | nargin < 3
0045   disp(sprintf('Loading data from file: %s', fname));
0046   data = load(fname);
0047   disp('Finished loading data')
0048 end
0049 
0050 % Get offset in file to start of first sequence
0051 
0052 if ~isempty(p.event_chan)
0053   if p.is_event_trig_pos
0054     start_samp = min(find(data(:,p.event_chan) > p.pos_trig_thresh));
0055   else
0056     start_samp = min(find(data(:,p.event_chan) < p.neg_trig_thresh));
0057   end
0058 else
0059   if p.is_event_trig_pos
0060     start_samp = min(find(data(:,p.slice_chan) > p.pos_trig_thresh));
0061   else
0062     start_samp = min(find(data(:,p.slice_chan) < p.neg_trig_thresh));
0063   end
0064 end
0065 
0066 % Eliminate earlier samps
0067 data(1:start_samp-2,:) = [];
0068 
0069 if p.remove_dc
0070   data(:,p.event_chan) = detrend(data(:,p.event_chan),0);
0071 end
0072 
0073 % Find key border onsets (pos/neg event sequences) and deviant marker offsets
0074 % (if bipolar)
0075 pos_idx = find(diff(data(:,p.event_chan) > p.pos_trig_thresh) == 1);
0076 
0077 % Find deviant marker onsets and key border offsets (if bipolar)
0078 neg_idx = find(diff(data(:,p.event_chan) < p.neg_trig_thresh) == 1);
0079 
0080 % Get all slice timing event onsets
0081 slice_idx = find(diff(data(:,p.slice_chan) > p.pos_trig_thresh) == 1);
0082 
0083 % Get all key press event onsets
0084 key_idx = [];
0085 if ~isempty(p.button_chan)
0086   if p.remove_dc
0087     data(:,p.button_chan) = detrend(data(:,p.button_chan),0);
0088   end
0089   if p.is_button_trig_pos
0090     key_idx = find(diff(data(:, p.button_chan) > p.pos_trig_thresh) == 1);
0091   else
0092     key_idx = find(diff(data(:,p.button_chan) < p.neg_trig_thresh_key) == 1);
0093   end
0094 end
0095 
0096 % Get all cardiac event onsets
0097 cardiac_idx = [];
0098 if ~isempty(p.cardiac_chan)
0099   cardiac_idx = find(diff(data(:,p.cardiac_chan) < p.neg_trig_thresh) == 1);
0100 end
0101 
0102 % Determine the number of runs and their starting locations in the data matrix
0103 if ~isempty(p.event_chan)
0104   start_idx = p.epi_start_to_first_stim - p.event_slop;
0105   stop_idx = p.epi_start_to_first_stim + p.event_slop;
0106   if p.is_event_trig_pos
0107     run_offset_idxs = ...
0108     pos_idx(find(ismember(diff(pos_idx),start_idx:stop_idx)));
0109   else
0110     run_offset_idxs = ...
0111     neg_idx(find(ismember(diff(neg_idx),start_idx:stop_idx)));
0112   end
0113 else
0114   run_offset_idxs = [1 slice_idx(find(diff(slice_idx) > ...
0115     p.slice_diff_thresh)+1)'];
0116 end
0117 
0118 nrun = length(run_offset_idxs);
0119 run_offset_idxs(end+1) = inf;
0120 for irun = 1:nrun
0121   % Find all slice onsets in this run
0122   curr_slice_idx = find((slice_idx >= run_offset_idxs(irun)) & (slice_idx < ...
0123       run_offset_idxs(irun+1)));
0124   
0125   % Eliminate the dummy shots
0126   num_dummy_slices = p.nslice_per_vol * p.ndummy_vol;
0127   dummy_s = p.ndummy_vol * p.TR;
0128   curr_slice_idx(1:num_dummy_slices) = [];
0129   
0130   ri(irun).slice_onsets = (slice_idx(curr_slice_idx)-run_offset_idxs(irun))/p.Fs - dummy_s;
0131   
0132   % Get indices of all positive and negative event markers for this run
0133   curr_pos_idx = find((pos_idx >= run_offset_idxs(irun)) & (pos_idx < run_offset_idxs(irun+1)));
0134   curr_neg_idx = find((neg_idx >= run_offset_idxs(irun)) & (neg_idx < run_offset_idxs(irun+1)));
0135   
0136   if p.bipolar
0137     %
0138     % NOTE: The bipolar part has not been debugged
0139     %
0140     
0141     % For each positive and negative marker pair (single event), determine which
0142     % came first
0143   
0144     % if value of idx_diff is negative, positive pulse came first
0145     idx_diff = pos_idx(curr_pos_idx) - neg_idx(curr_neg_idx);        
0146     ri(irun).pos_events = pos_idx(curr_pos_idx(find(sign(idx_diff) == -1)));
0147     ri(irun).neg_events = neg_idx(curr_neg_idx(find(sign(idx_diff) == 1)));
0148   else
0149     ri(irun).pos_events = (pos_idx(curr_pos_idx)-run_offset_idxs(irun))/p.Fs - dummy_s;
0150     ri(irun).neg_events = (neg_idx(curr_neg_idx)-run_offset_idxs(irun))/p.Fs - dummy_s;
0151   end
0152 
0153   %
0154   % Extract the key press times
0155   %
0156   curr_key_idx = find((key_idx >= run_offset_idxs(irun)) & (key_idx < ...
0157       run_offset_idxs(irun+1)));
0158   ri(irun).key_events = (key_idx(curr_key_idx)-run_offset_idxs(irun))/p.Fs - dummy_s;
0159   
0160   % Check to see if we should discard events that occur after some time from
0161   % run onsets
0162   if isfield(p,'max_run_dur') && ~isempty(p.max_run_dur)
0163     ri(irun).pos_events(ri(irun).pos_events > p.max_run_dur) = [];
0164     ri(irun).neg_events(ri(irun).neg_events > p.max_run_dur) = [];
0165     ri(irun).key_events(ri(irun).key_events > p.max_run_dur) = [];
0166   end
0167   
0168   %
0169   % Extract the onsets of heart beat markers. Negative going pulses
0170   %
0171   curr_cardiac_idx = find((cardiac_idx >= run_offset_idxs(irun)) & ...
0172       (cardiac_idx < run_offset_idxs(irun+1)));
0173   ri(irun).cardiac = (cardiac_idx(curr_cardiac_idx)-run_offset_idxs(irun))/p.Fs - dummy_s;
0174   
0175   % Remove cardiac markers occuring prior to start of actual image acquisition
0176   ri(irun).cardiac(find(ri(irun).cardiac < 0)) = [];
0177   
0178   run_start_samp = run_offset_idxs(irun) + p.Fs*dummy_s;
0179   run_stop_samp = slice_idx(curr_slice_idx(end));
0180 
0181   % Extract the respiratory data
0182   ri(irun).respir = data(run_start_samp:run_stop_samp,p.respir_chan);
0183   
0184   % Filter to get rid of higher frequency cardiac info
0185   if p.low_pass_cutoff
0186     [b,a] = butter(5,p.low_pass_cutoff/p.Fs);
0187     ri(irun).respir = filtfilt(b,a,ri(irun).respir);
0188   end
0189   
0190 end % for irun=
0191 
0192 if p.savedata
0193   [outpath, name] = fileparts(fname);
0194   outfname = fullfile(outpath, sprintf('%s%s.mat',name,p.outsuffix));
0195   fprintf('Saving data to %s ...\n', outfname);
0196   save(outfname,'ri')
0197 end