erpleeg() - continuous EEG analysi utility class: loads entire .crw, scales to microvolts, loads .log events
Usage:
>> e=erpleeg('ac08.crw');
These three steps provide calibrated continuous EEG and all the
event info available in the .crw and .log files.
>> e=erpleeg('ac08.crw');
>> e.calibrate('ac08.crw', 200, -16, 20, 10, 1, 2);
>> e.loadlog('ac08.log');
Use cases:
1. Slurp up .crw file into continous eeg data array (channels x sample),
along with a digheader structure and and timestamped event info (times,
event codes).
>> e = erpleeg('ac08.crw');
e =
erpleeg handle
Properties:
crwf: 'ac08.crw'
ceeg: [32x1009151 single]
digheader: [1x1 struct]
cevents: [5242x2 double]
cals: []
logf: ''
log: []
2. calibrate kutaslab EEG data using event-triggered calpulses. Like normerp but
for continous EEG. Calpulses are pattern-matched in the recordings
based on their wave shape so calibration doesn't need a log file or
condition code 0 to run.
USAGE: There are seven parameters, all obligatory, order fixed, ** NOT CHECKED **
.crw file w/ cals, pulse duration in ms, lo_cursor, high_cursor, cal
amplitude, polarity, number of points to average.
>> e.calibrate('ac08.crw', 200, -16, 20, 10, 1, 2)
e =
erpleeg handle
Properties:
crwf: 'ac08.crw'
ceeg: [32x1009151 single]
digheader: [1x1 struct]
cevents: [5242x2 double]
cals: [1x1 struct]
logf: ''
log: []
NOTE 1: some useful info about the calibration processing is tucked into the
erpleeg.cals structure.
>> e.cals
ans =
params: [1x1 struct]
stepfunc: [-5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 5 5 5 5 5 5 5 5 5 5 5 5]
pulsemodels: [156x32x3 double]
AD: [117x32x24 double]
ADmn: [32x24 double]
uVmn: [32x24 double]
3. Always check calibration parameters visually ...
>> e.plotcals()
4. Load log file info including artifact flags as dumped by logcat2
>> e.loadlog('ac08.log');
e =
Properties:
crwf: 'ac08.crw'
ceeg: [32x1009151 single]
digheader: [1x1 struct]
cevents: [5242x2 double]
cals: [1x1 struct]
logf: 'ac08.log'
log: [1x1 struct]
NOTE: the log structure includes column names and (strictly numeric)
log event information
e.log
ans =
names: {1x7 cell}
events: [5243x7 double]
logerrs: [0x9 double]
>> e.log.names
ans =
Columns 1 through 6
'evtno' 'evtcode' 'ccode' 'flags' 'clock_ticks' 'total_time'
Column 7
'delta_t'
Author: T. Urbach 04/20160001 % erpleeg() - continuous EEG analysi utility class: loads entire .crw, scales to microvolts, loads .log events 0002 % 0003 % Usage: 0004 % >> e=erpleeg('ac08.crw'); 0005 % 0006 % These three steps provide calibrated continuous EEG and all the 0007 % event info available in the .crw and .log files. 0008 % 0009 % >> e=erpleeg('ac08.crw'); 0010 % >> e.calibrate('ac08.crw', 200, -16, 20, 10, 1, 2); 0011 % >> e.loadlog('ac08.log'); 0012 % 0013 % Use cases: 0014 % 0015 % 1. Slurp up .crw file into continous eeg data array (channels x sample), 0016 % along with a digheader structure and and timestamped event info (times, 0017 % event codes). 0018 % 0019 % >> e = erpleeg('ac08.crw'); 0020 % 0021 % e = 0022 % 0023 % erpleeg handle 0024 % 0025 % Properties: 0026 % crwf: 'ac08.crw' 0027 % ceeg: [32x1009151 single] 0028 % digheader: [1x1 struct] 0029 % cevents: [5242x2 double] 0030 % cals: [] 0031 % logf: '' 0032 % log: [] 0033 % 0034 % 2. calibrate kutaslab EEG data using event-triggered calpulses. Like normerp but 0035 % for continous EEG. Calpulses are pattern-matched in the recordings 0036 % based on their wave shape so calibration doesn't need a log file or 0037 % condition code 0 to run. 0038 % 0039 % USAGE: There are seven parameters, all obligatory, order fixed, ** NOT CHECKED ** 0040 % 0041 % .crw file w/ cals, pulse duration in ms, lo_cursor, high_cursor, cal 0042 % amplitude, polarity, number of points to average. 0043 % 0044 % >> e.calibrate('ac08.crw', 200, -16, 20, 10, 1, 2) 0045 % e = 0046 % 0047 % erpleeg handle 0048 % 0049 % Properties: 0050 % crwf: 'ac08.crw' 0051 % ceeg: [32x1009151 single] 0052 % digheader: [1x1 struct] 0053 % cevents: [5242x2 double] 0054 % cals: [1x1 struct] 0055 % logf: '' 0056 % log: [] 0057 % 0058 % NOTE 1: some useful info about the calibration processing is tucked into the 0059 % erpleeg.cals structure. 0060 % 0061 % >> e.cals 0062 % 0063 % ans = 0064 % 0065 % params: [1x1 struct] 0066 % stepfunc: [-5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 5 5 5 5 5 5 5 5 5 5 5 5] 0067 % pulsemodels: [156x32x3 double] 0068 % AD: [117x32x24 double] 0069 % ADmn: [32x24 double] 0070 % uVmn: [32x24 double] 0071 % 0072 % 0073 % 3. Always check calibration parameters visually ... 0074 % 0075 % >> e.plotcals() 0076 % 0077 % 0078 % 4. Load log file info including artifact flags as dumped by logcat2 0079 % 0080 % >> e.loadlog('ac08.log'); 0081 % 0082 % e = 0083 % Properties: 0084 % crwf: 'ac08.crw' 0085 % ceeg: [32x1009151 single] 0086 % digheader: [1x1 struct] 0087 % cevents: [5242x2 double] 0088 % cals: [1x1 struct] 0089 % logf: 'ac08.log' 0090 % log: [1x1 struct] 0091 % 0092 % NOTE: the log structure includes column names and (strictly numeric) 0093 % log event information 0094 % 0095 % e.log 0096 % 0097 % ans = 0098 % 0099 % names: {1x7 cell} 0100 % events: [5243x7 double] 0101 % logerrs: [0x9 double] 0102 % 0103 % >> e.log.names 0104 % 0105 % ans = 0106 % 0107 % Columns 1 through 6 0108 % 0109 % 'evtno' 'evtcode' 'ccode' 'flags' 'clock_ticks' 'total_time' 0110 % 0111 % Column 7 0112 % 0113 % 'delta_t' 0114 % 0115 % 0116 % Author: T. Urbach 04/2016 0117 % 0118 0119 0120 % ------------------------------------------------------------ 0121 % BEGIN SOURCE 0122 % 0123 % Note: class inherits from handle so methods can change the property 0124 % objects w/out explict copy() for object property assignmnts. 0125 % ------------------------------------------------------------ 0126 0127 classdef erpleeg < handle 0128 properties (Access='public') 0129 % these are always assigned during construction 0130 crwf; % crw file 0131 ceeg; % continuous data matrix, chans x samples 0132 eventtrack; % events from crw, 0-filled at sampleswhen no events occured 0133 digheader; % headinfo: struct build from raw2asci -header output 0134 cevents; % row x column data table where columns 0135 % columns = timestamps, event code, cond code 0136 0137 % optional preprocessing 0138 cals = []; % structure w/ 2 sets of calibration snippets: before, after 0139 logf = ''; 0140 log = []; % rows x columns info from logcat2 output, checked against cevents 0141 0142 end 0143 0144 % public methods 0145 % - erpleeg: constructor, always loads crw/raw and events from the raw 0146 % - calibration: optionally calibrates 0147 % - loadlog: optionally load log file events incl. cond codes 0148 0149 methods (Static, Access='public') 0150 % these are auto-called by MATLAB save and load 0151 % stubs here in case default save/load needs to overridden 0152 function s = saveobj(obj) 0153 s = obj; 0154 end 0155 0156 function obj = loadobj(s) 0157 obj = s; 0158 end 0159 end 0160 0161 methods (Access='public') 0162 0163 function obj = erpleeg(crwf) 0164 % constructor sets the file name, loads the headinfo, the data, calibrates it 0165 % and computes some convenience variable 0166 0167 obj.crwf = crwf; 0168 obj.digheader = loaddigheader(obj); 0169 [obj.ceeg, obj.eventtrack, obj.cevents] = obj.loadcrw(obj.crwf); % this calls erpio2 to slurp data 0170 0171 % for convenience 0172 obj.digheader.srate = 100000 / obj.digheader.ctickt; % dig hardcodes 100K clocktick base 0173 0174 end 0175 0176 function loadlog(obj,logf) 0177 % optional function reads in logfile info via logcat2 0178 % no logcat2 on the cluster, using log2asci w/ a temp asci 0179 % log in /tmp 0180 0181 [status, uuid] = system('uuidgen'); 0182 if ( status ~= 0 ) 0183 error('system(uuidgen) failed ... not my problem') 0184 end 0185 0186 % nab event info w/ logcat2 0187 syscmd = sprintf('logcat2 %s -srate %d', ... 0188 logf, obj.digheader.srate ); 0189 0190 % rval has the labeled columns and data 0191 [status, rval] = system(syscmd); 0192 if (status ~= 0) 0193 fprintf('mklinux error\n %s', rval); 0194 return 0195 end 0196 0197 % initialize log events structure: names, events 0198 logevents = struct(); 0199 0200 % parse column digheader ... gotta be a better way 0201 for (c = [1:7]) 0202 [head, rval] = strtok(rval); 0203 logevents.names{c} = head; 0204 end 0205 0206 % scan and reshape the event data ... 7 cols w/ -srate option 0207 events = sscanf(rval, '%f'); 0208 events = reshape(events, 7, length(events)/7)'; 0209 logevents.events = events; 0210 0211 % sanity check against the events loaded from the crw/raw 0212 % atrocious hack: check one less log event b.c. 0213 % erpio2('get_cdata') drops last pause mark 0214 rawlogmismatch = find(events(1:end-1,2) ~= obj.cevents(:,2)); 0215 logevents.logerrs = [ 0216 obj.cevents(rawlogmismatch, :), ... 0217 events(rawlogmismatch,:) ]; 0218 0219 % if we make it here, set the file name and the log event structure 0220 obj.logf = logf; 0221 obj.log = logevents; 0222 0223 end 0224 0225 function calibrate(obj, cal_file, cal_len, ... 0226 lo_curs, hi_curs, cal_uV, cal_pol, ... 0227 avg_npts) 0228 % calibrate AD units to microvolts using event triggered calpuluses. 0229 % 0230 % 1. find calibration pulse signals in the EEG by linear model pattern match 0231 % to step function at (any) event codes. Looking around event codes 0232 % is more efficient but this generalizes to sweeping the entire 0233 % data file for running (non-triggered) calibration pulses as 0234 % well. 0235 % 0236 % 2. calculate calibration w/ normerp-like arguments except ... 0237 % "resolution" is not needed since cals are MATLAB floats 0238 % "cal_len" is added to give the cal pulse duty cycle. 0239 % 0240 % usage: arguments mandatory every time, no defaults 0241 % 0242 % cal_file: filename of .crw containing the pulses to use 0243 % cal_len: pulse length (=duty cycle) in millisecond, e.g., 200 0244 % lo_curs: time in ms to locate the bottom of the step relativ 0245 % hi_curs: time in ms to locate the top 0246 % cal_uV: amplitude of calibration pulse in microvolts, e.g., 10 0247 % cal_pol: pulse polarity, e.g. 1 or -1 0248 % avg_npts: average npts either side of cursor (total = 2*avgnpts+1) 0249 % 0250 0251 usagestr = ... 0252 ['USAGE: calibrate(cal_file, cal_len, lo_curs, hi_curs, cal_uV, ' ... 0253 'cal_pol, avg_npts) \n']; 0254 usagestr = [usagestr ... 0255 'cal_file: filename of .crw containing the pulses\n' ... 0256 'cal_len: pulse length (=duty cycle) in millisecond, e.g., 200\n' ... 0257 'lo_curs: time in ms before pulse to center measurement, e.g., -16\n' ... 0258 'hi_curs: time in ms after pulse to center measurement, e.g., 20\n' ... 0259 'cal_uV: amplitude of calibration pulse in microvolts, e.g., 10\n' ... 0260 'cal_pol: pulse polarity, e.g. 1 or -1\n' ... 0261 'avg_npts: average npts either side of cursor (total = 2*avgnpts+' ... 0262 '1)\n\n']; 0263 usagestr = [usagestr ... 0264 'EXAMPLE: e.calibrate(''ac08.crw'', 200, -16, 20, 10, 1, 2);']; 0265 0266 if ( ~exist('cal_file') | ~exist('cal_len') | ~exist('lo_curs') | ... 0267 ~exist('hi_curs') | ~exist('cal_uV') | ~exist('cal_pol') | ... 0268 ~exist('avg_npts') ) 0269 disp(sprintf(usagestr)) 0270 return 0271 end 0272 0273 % ------------------------------------------------------------ 0274 % Constants 0275 % ------------------------------------------------------------ 0276 0277 % hard coded pretty arbitrary R^2 > threshold for identifying a 0278 % "good enough" pulse on any channel. 0279 calR2_threshold = 0.90; 0280 0281 % make an ideal ca lpulse waveform 50% of the duty cycle 0282 % e.g. -50 to 50 ms step for a 200ms pulse 0283 % steplen samples, forced to even number 0284 steplen = round( (cal_len/(2.0 * 1000) * obj.digheader.srate )); 0285 steplen = steplen - mod(steplen,2); 0286 0287 % center step function 0288 stepfunc = [(-cal_uV/2.0) * ones(1,steplen/2), ... 0289 (cal_uV/2.0) * ones(1,steplen/2)]; 0290 0291 % replicate the step for each channel 0292 stepfuncs = repmat(stepfunc, obj.digheader.nchans,1); 0293 0294 % this will hold the pulse model fits 0295 pulsemodels = []; % number of cals x nchans x 3: b_0, b_1, R^2 0296 modelones = ones(length(stepfunc),1); % for regression intercept 0297 0298 % ------------------------------------------------------------ 0299 % Processing 0300 % ------------------------------------------------------------ 0301 0302 % 1. load up the cal pulses 0303 0304 % use original crw for cals, unless a different cal_file is given 0305 if ( strcmp(obj.crwf, cal_file) == 1 ) 0306 cal = obj; 0307 else 0308 % load up the different cal_file 0309 cal = erpleeg(cal_file); 0310 0311 % modicum of error checking 0312 if ( ~isequal(obj.digheader.chndes, cal.digheader.chndes )) 0313 error(sprintf('cal_file %s and EEG file %s have different channels', ... 0314 cal_file, obj.crwf)) 0315 end 0316 end 0317 0318 % loop through the events checking for cal pulses via good 0319 % linear fit to a square wave. Capture the linear model beta's 0320 % and R^2s for subsequent pruning. 0321 0322 % initialize an empty cal array: 0 cals x nchans x length stepfunc 0323 calcount = 0; 0324 calpulses = zeros(calcount, obj.digheader.nchans, length(stepfunc)); 0325 0326 for e = 1:length(cal.cevents(:,1)) 0327 0328 evtime = cal.cevents(e,1) + cal.digheader.odelay; % timestamp in ms .. 0329 evsamp = floor(evtime/(1000.0/cal.digheader.srate)); % = this sample 0330 0331 % nab all channels of snippet of data on either side of the event 0332 datasnippet = cal.ceeg(:, [evsamp-steplen/2: evsamp+(steplen/2)-1]); 0333 0334 % chans x 3: b_0, b_1, R^2 0335 pulsemodel = zeros(cal.digheader.nchans, 3); 0336 0337 % loop through channels calculating R^2 at each channel 0338 % if the fit is crummy, break out and try the next event 0339 % attrociously procedural but saves time since most events are not 0340 % cals ... first bad channel spoils the whole event and 0341 % the process moves on 0342 0343 for c = 1:cal.digheader.nchans 0344 0345 % nab this channel's pulse for processing 0346 thispulse = []; 0347 thispulse = datasnippet(c,:)'; 0348 0349 % drive a line through it, i.e., 0350 % this_pulse = b_0 + b_1 .* stepfunc' 0351 lmcal = []; 0352 lmcal = [modelones stepfunc'] \ thispulse; 0353 0354 % straight line estimate 0355 snippet_hat = [modelones stepfunc'] * lmcal; 0356 0357 % calculate R^2 for this channel, this snippet as 0358 % R2 = 1- sum square residuals over variance 0359 R2 = 1 - sum( (thispulse - snippet_hat).^2 ) / ... 0360 sum( (thispulse - mean(thispulse)).^2 ); 0361 0362 % no need to continue with if this channel if 0363 % regression slope is negative or fit is poor 0364 if (lmcal(2) < 0 | R2 < calR2_threshold) 0365 break 0366 end 0367 0368 % else nab the model and fit info 0369 pulsemodel(c, :) = [lmcal', R2]; 0370 0371 % AND. if we make it to the last channel without breaking out 0372 % of the loop ... 0373 % this is a reasonabl cal ... demean and capture it along w/ 0374 % the model info 0375 if ( c == cal.digheader.nchans ) 0376 calcount = calcount + 1; 0377 snipmn = repmat(mean(datasnippet,2),1,steplen); 0378 calpulses(calcount, :,:) = datasnippet-snipmn; 0379 pulsemodels(calcount, :, :) = pulsemodel; 0380 0381 % clf; plot(datasnippet'); % debugging 0382 % EEG.cevents(e,:) 0383 % dummy = input('Enter to continue'); 0384 end 0385 end 0386 end 0387 0388 % 2. Second pass, drop the bottom quartile of poorest fits 0389 0390 % ------------------------------------------------------------ 0391 % R2 based trimming ... keep top 3 quartiles 0392 % ------------------------------------------------------------ 0393 0394 % nab the R2s ... 3rd column in pulsemodels 0395 R2s = squeeze(pulsemodels(:,:,3)); % size(R2s) = ntrials x nchans 0396 0397 % for each trial, nab the smallest R2 = worst fit, at any channel 0398 % trials x 1 0399 minR2s = min(R2s, [],2); 0400 0401 % find the bottom quartile of worst fitting R^2s 0402 quartR2 = prctile(minR2s,25); 0403 0404 % trials x 1 column of 0s, 1, 0405 % droppedcaltrials = minR2s <= quartR2; 0406 goodR2trials = quartR2 < minR2s; 0407 0408 % ------------------------------------------------------------ 0409 % Linear slope based trimming ... trim slopes per Tukey 0410 % one bad channel and the whole trial is out 0411 % ------------------------------------------------------------ 0412 0413 % nab the linear model slopes ... 2nd column in pulsemodels: 0414 lmslopes = squeeze(pulsemodels(:,:,2)); % size(lmslopes) = ntrials x nchans 0415 0416 % nab quartiles across trials: 2 x chan 0417 slopepctiles = prctile(lmslopes, [25, 50, 75]); 0418 0419 % trim at median +/- 2.5 * interquartile range per Tukey 0420 trimlev = 2.5; 0421 slopetrimlo = slopepctiles(2,:) - (trimlev * (slopepctiles(3,:)-slopepctiles(1,:))); 0422 slopetrimhi = slopepctiles(2,:) + (trimlev * (slopepctiles(3,:)-slopepctiles(1,:))); 0423 0424 % matrix of cut off values ... 2 rows (hi, lo) x chan columns 0425 slopetrimvals = [slopetrimlo; slopetrimhi]; 0426 0427 % find the trials with *every channel* slope in the middle 2 quartiles 0428 % probably a way to vectorize 0429 goodslopetrials = []; 0430 for (i = [1: length(lmslopes(:,1))]) 0431 % check that slope at *all* channels is between the high and low trim vals 0432 goodslopetrials(i,1) = all(slopetrimvals(1,:) < lmslopes(i,:) & ... 0433 lmslopes(i,:) < slopetrimvals(2,:)); 0434 0435 end; 0436 0437 % update the keepers ... top 75% of the best model fits 0438 % calpulses = calpulses(find(minR2s > quartR2),:,:); 0439 % calpulses = calpulses(find(~droppedcaltrials),:,:); 0440 0441 % To be a good cal it must have a good slope with a good fit 0442 goodcaltrials = goodR2trials & goodslopetrials; 0443 calpulses = calpulses(find(goodcaltrials),:,:); 0444 0445 % ??? snapshot the bottom quartile for checking later 0446 % badcalpulses = calpulses(find(~goodcaltrials),:,:); 0447 0448 % fail if < 64 calibration pulses ... arbitrary choice 0449 caldims = size(calpulses); 0450 ncalmin = 60; 0451 if (caldims(1) < ncalmin) 0452 error(sprintf('Found %5d good cals after trimming, need at least %d', ... 0453 caldims(1), ncalmin)) 0454 end 0455 0456 % 3. calculate the time-domain average of the remaining good 0457 % calibration pulses. 0458 calpulses_mn = squeeze(mean(calpulses,1)); 0459 0460 % nab the mean cal points around specified cursors 0461 % .. normerp-like 0462 lcurspt = round(caldims(3)/2.0 + cal.ms2samp(lo_curs)); 0463 lcurspts = [lcurspt - avg_npts] + [1: (2*avg_npts) + 1]; 0464 lo_mn = mean(calpulses_mn(:,lcurspts),2); 0465 0466 hcurspt = round(caldims(3)/2.0 + cal.ms2samp(hi_curs)); 0467 hcurspts = [hcurspt - avg_npts] + [1: (2*avg_npts) + 1]; 0468 hi_mn = mean(calpulses_mn(:,hcurspts),2); 0469 0470 % calculate the cal step size at each channel: 32 x 1 0471 calscalars = (hi_mn - lo_mn)./cal_uV; 0472 0473 % ** FINALLY ** rescale single trial EEG to microvolts 0474 for c = [1:length(calscalars)] 0475 0476 % convert object eeg to uV 0477 obj.ceeg(c,:) = obj.ceeg(c,:)./calscalars(c); 0478 0479 % compute the mean uV cals and baseline on the low cursor samples 0480 % for visual check 0481 calsuVmn(c,:) = calpulses_mn(c,:)./calscalars(c); 0482 calsuVmn(c,:) = calsuVmn(c,:) - ... 0483 ( mean(calsuVmn(c,lcurspts)) - mean(stepfunc(lcurspts)) ); 0484 0485 0486 end 0487 0488 % 4. if we get here, update obj info on the way out ... 0489 calparams = struct(); 0490 calparams.cal_file = cal_file; 0491 calparams.cal_len = cal_len; 0492 calparams.lo_curs = lo_curs; 0493 calparams.hi_curs = hi_curs; 0494 calparams.cal_uV = cal_uV; 0495 calparams.cal_pol = cal_pol; 0496 calparams.avg_npts = avg_npts; 0497 calparams.cal_R2_threshold = calR2_threshold; 0498 calparams.lo_curs_pts = lcurspts; % points sampled before step 0499 calparams.hi_curs_pts = hcurspts; % point sampled after step 0500 0501 obj.cals.params = calparams; 0502 obj.cals.stepfunc = stepfunc; % 2-D array: chan x samp 0503 obj.cals.pulsemodels = pulsemodels; % 2-D array: b_0, b_1, R^2 0504 obj.cals.AD = calpulses; % 3-D array: trial x chan x samp 0505 obj.cals.ADmn = calpulses_mn; % 2-D array: chan x samp 0506 obj.cals.uVmn = calsuVmn; % 2-D array: chan x samp 0507 obj.cals.goodcaltrials = goodcaltrials; % 1-D array: trials x, 1==used for calibration 0508 0509 end % loadcals 0510 0511 0512 function plotcals(obj) 0513 % plots A/D before and microvolts after calibration 0514 % plots single trial A/D calpulses actually used in calibration by channel 0515 if ( ~isfield(obj.cals, 'uVmn') ) 0516 error('Run calibration function before plotcals') 0517 end 0518 0519 % ------------------------------------------------------------ 0520 % cals are identified in the continuous EEG by fitting against 0521 % a step function. This plots the linear model fits for all 0522 % the data snippets and the subset that remain after trimming. 0523 % ------------------------------------------------------------ 0524 figure 0525 subplot(2,2,1) 0526 hold on 0527 title('All Linear squarewave fit: Intercept'); 0528 boxplot(obj.cals.pulsemodels(:,:,1)) 0529 plot( [1:obj.digheader.nchans], obj.cals.pulsemodels(:,:,1),'b.') 0530 xlabel('Intercept') 0531 hold off 0532 0533 subplot(2,2,2) 0534 hold on 0535 title('All Linear squarewave fit: Slope'); 0536 boxplot(obj.cals.pulsemodels(:,:,2)) 0537 plot([1:obj.digheader.nchans], obj.cals.pulsemodels(:,:,2), 'r.') 0538 xlabel('Slope') 0539 hold off 0540 0541 subplot(2,2,3) 0542 hold on 0543 title('Used These Cals Linear squarewave fit: Intercept'); 0544 boxplot(obj.cals.pulsemodels(find(obj.cals.goodcaltrials),:,1)) 0545 plot( [1:obj.digheader.nchans], obj.cals.pulsemodels(find(obj.cals.goodcaltrials),:,1),'b.') 0546 xlabel('Intercept') 0547 hold off 0548 0549 subplot(2,2,4) 0550 hold on 0551 title('Used These Cals Linear squarewave fit: Slope'); 0552 boxplot(obj.cals.pulsemodels(find(obj.cals.goodcaltrials),:,2)) 0553 plot([1:obj.digheader.nchans], obj.cals.pulsemodels(find(obj.cals.goodcaltrials),:,2), 'r.') 0554 xlabel('Slope') 0555 hold off 0556 0557 % ------------------------------------------------------------ 0558 % over plot the single trial calibration trials actually used 0559 % after the trimming procedures. If wonky cals slip through, 0560 % they will show up here. 0561 % ------------------------------------------------------------ 0562 figure 0563 % four columns, however many rows it takes ... 0564 ncol = 4; 0565 nrow = ceil(obj.digheader.nchans/ncol); 0566 caldims = size(obj.cals.AD); 0567 pulsedims = size(obj.cals.pulsemodels); 0568 for (c = [1:obj.digheader.nchans] ) 0569 subplot(nrow, ncol, c); 0570 plot(squeeze(obj.cals.AD(:,c,:))'); 0571 xlabel(obj.digheader.chndes(c)); 0572 if (c==1) 0573 ylabel('Unscaled A/D units'); 0574 title(sprintf('Cal file %s: these %d used out of %d found', ... 0575 obj.cals.params.cal_file, ... 0576 caldims(1), pulsedims(1))); 0577 end 0578 end 0579 0580 % ------------------------------------------------------------ 0581 % Punchline: plot mean cal pulse for the good cals, before and after 0582 % converting to microvolts 0583 % ------------------------------------------------------------ 0584 figure; 0585 hold on 0586 0587 % before 0588 minAD = min(obj.cals.ADmn(:)); 0589 maxAD = max(obj.cals.ADmn(:)); 0590 0591 subplot(2,1,1) 0592 hold on 0593 title(sprintf('Before: Mean cal pulse in %s: (%d good cals of %d found)', ... 0594 obj.cals.params.cal_file, ... 0595 length(find(obj.cals.goodcaltrials == 1)), ... 0596 length(obj.cals.goodcaltrials) )); 0597 plot(obj.cals.ADmn'); 0598 plot([obj.cals.params.lo_curs_pts; obj.cals.params.lo_curs_pts], ... 0599 [minAD, maxAD], 'k-', 'LineWidth',1.5); 0600 plot([obj.cals.params.hi_curs_pts; obj.cals.params.hi_curs_pts], ... 0601 [minAD, maxAD], 'k-', 'LineWidth',1.5); 0602 plot(obj.cals.stepfunc, 'r-+', 'LineWidth',2); 0603 ylabel('AD units'); 0604 0605 % show the channel labels 0606 legend(obj.digheader.chndes, 'Location', 'northeastoutside'); 0607 before_pos = get(gca, 'Position'); 0608 hold off 0609 0610 % after 0611 subplot(2,1,2) 0612 hold on 0613 title(sprintf('Calibration file %s: After', obj.cals.params.cal_file)) 0614 plot(obj.cals.uVmn'); 0615 plot([obj.cals.params.lo_curs_pts; obj.cals.params.lo_curs_pts], ... 0616 [-1, 1].*obj.cals.params.cal_uV/2.0, 'k-', 'LineWidth',1.5); 0617 plot([obj.cals.params.hi_curs_pts; obj.cals.params.hi_curs_pts], ... 0618 [-1, 1].*obj.cals.params.cal_uV/2.0, 'k-', 'LineWidth',1.5); 0619 plot(obj.cals.stepfunc, 'r-+', 'LineWidth',2); 0620 ylabel('uV'); 0621 xlabel('Samples'); 0622 after_pos = get(gca, 'Position'); 0623 0624 % set right edge of after plot to match right edge of before 0625 after_pos(3) = before_pos(3); 0626 set(gca, 'Position', after_pos); 0627 hold off; 0628 0629 0630 end 0631 0632 0633 end % public methods 0634 0635 % These are utilities that no users need to call 0636 methods (Access='private') 0637 0638 function digheader = loaddigheader(obj) 0639 0640 syscmd = sprintf('raw2asci %s -header', obj.crwf); 0641 [status, rval] = system(syscmd); 0642 if (status ~= 0) 0643 fprintf('mklinux error\n %s', rval); 0644 return 0645 end 0646 0647 digheader = struct(); % initialize, e.g., 0648 0649 % evtno = 0x17A5 0650 % epleng = 0 0651 % nchans = 32 0652 % sums = 0 0653 % tpfuncs = 0 0654 % pp10uv = 0 0655 % verpos = 0 0656 % odelay = 8 0657 % totevnt = 0 0658 % ctickt = 400 0659 % dh_clkt = 0 0660 % ccoder = 0 0661 % presam = 0 0662 % trfuncs = 0 0663 % totrr = 0 0664 % totrej = 0 0665 % sbcode = 0 0666 % cprecis = 1 0667 % ovf_errs = 0 0668 % decfact = 0 0669 % dh_flag = 0 0670 % dh_item = 0 0671 % rfcnts[0] = 0 0672 % rfcnts[1] = 0 0673 % rfcnts[2] = 0 0674 % rfcnts[3] = 0 0675 % rfcnts[4] = 0 0676 % rfcnts[5] = 0 0677 % rfcnts[6] = 0 0678 % rfcnts[7] = 0 0679 % rftypes0 = 0680 % rftypes1 = 0681 % rftypes2 = 0682 % rftypes3 = 0683 % rftypes4 = 0684 % rftypes5 = 0685 % rftypes6 = 0686 % rftypes7 = 0687 % chndes0 = lle 0688 % chndes1 = lhz 0689 % chndes2 = MiPf 0690 % chndes3 = LLPf 0691 % chndes4 = RLPf 0692 % chndes5 = LMPf 0693 % chndes6 = RMPf 0694 % chndes7 = LDFr 0695 % chndes8 = RDFr 0696 % chndes9 = LLFr 0697 % chndes10 = RLFr 0698 % chndes11 = LMFr 0699 % chndes12 = RMFr 0700 % chndes13 = LMCe 0701 % chndes14 = RMCe 0702 % chndes15 = MiCe 0703 % chndes16 = MiPa 0704 % chndes17 = LDCe 0705 % chndes18 = RDCe 0706 % chndes19 = LDPa 0707 % chndes20 = RDPa 0708 % chndes21 = LMOc 0709 % chndes22 = RMOc 0710 % chndes23 = LLTe 0711 % chndes24 = RLTe 0712 % chndes25 = LLOc 0713 % chndes26 = RLOc 0714 % chndes27 = MiOc 0715 % chndes28 = A2 0716 % chndes29 = HEOG 0717 % chndes30 = rle 0718 % chndes31 = rhz 0719 % subdes = Sub#=ac08 List=C1 Yes=RH 05-15-13 0720 % sbcdes = 0721 % condes = 0722 % expdes = anocong 0723 % pftypes0 = 0724 % pftypes1 = 0725 % pftypes2 = 0726 % flags = 0x0000 0727 % nrawrecs = 0 0728 % idxofflow = 0 0729 % idxoffhi = 0 0730 0731 while (length(rval) > 0) 0732 0733 % MATLAB string parsing bullshit ... char(10) is linux newline 0734 % for (c = [1:length(rval)]) 0735 % if (strcmp(rval(c), char(10))) 0736 % fprintf('%d: %s\n', c, rval(c)) 0737 % error('found newline') 0738 % end 0739 % end 0740 0741 [l, rval] = strtok(rval, char(10)); 0742 0743 % ignore blanks 0744 if (strcmp(l,'') == 1 ) 0745 continue 0746 end 0747 0748 % ignore HEADER: 0749 if (strcmp(l, 'HEADER:')==1) 0750 continue 0751 end 0752 0753 % handle comment lines 0754 if ( strcmp(l(1),'*') == 1 ) 0755 if ~ isfield(digheader,'comment') 0756 digheader.comment = sprintf('%s\n', l); 0757 else 0758 digheader.comment = sprintf('%s%s\n', digheader.comment, l); 0759 end 0760 continue 0761 end 0762 0763 % parse and set the rest of the field = value lines 0764 [fieldname, l] = strtok(l); % fieldname 0765 [dummy, l ] = strtok(l); % skip the = 0766 0767 fieldname = strtrim(fieldname); % tidy up 0768 fieldval = strtrim(l); % string field value 0769 0770 % translate the square brackets and 0-base to 1-base ... ugh 0771 if ( strfind(fieldname, 'rfc')) 0772 for (rfc = [0:7]) 0773 fieldname = strrep(fieldname, ... 0774 sprintf('[%d]',rfc ), ... 0775 sprintf('(%d)',rfc+1 )); 0776 end 0777 end 0778 0779 % convert the 0-base channel descriptions back into 1-base array indices 0780 % chndes0, chndes1, ... -> chndes(1), chndes(2) ... 0781 if ( strfind(fieldname, 'chndes') ) 0782 index = str2num(fieldname(7:length(fieldname))); % 0-base string rep 0783 fieldname = regexprep(fieldname, ... 0784 '([\d]*)', ... 0785 sprintf('{%d}', index+1)); 0786 end 0787 0788 % build the cmdstring and run it ... 0789 0790 % these are string values (or 0xnumeric ...) 0791 if ( length(strfind(fieldname, 'des')) == 1 | ... 0792 length(regexp(fieldval, '^0x')) == 1 | ... 0793 length(fieldval) == 0 ) 0794 % strings for descriptions 0795 cmdstr = sprintf('digheader.%s=''%s'';', fieldname, fieldval); 0796 else 0797 % numeric values for everythings else 0798 cmdstr = sprintf('digheader.%s=%d;', fieldname, str2num(fieldval)); 0799 end 0800 0801 % disp(cmdstr); % debugging 0802 eval(cmdstr); 0803 end 0804 digheader; 0805 end 0806 0807 function [ceeg, eventtrack, cevents] = loadcrw(obj, crwf) 0808 % load continuous data and convert to single precision 0809 [d, e, t] = erpio2('get_cdata', crwf, 'msec', 0,0,0); 0810 ceeg = single(d); % raw AD samples 0811 eventtrack = e; % events (zero-filled) 0812 0813 % update build event array w/ cont. data events 0814 cevents = [t(e~=0); e(e~=0)]'; 0815 0816 end 0817 0818 % ------------------------------------------------------------ 0819 % utility functions 0820 % ------------------------------------------------------------ 0821 function samp = ms2samp(obj,ms) 0822 % convert ms to samples at sampling rate for this object 0823 samp = round(ms / (1000.0/obj.digheader.srate)); 0824 end 0825 0826 function ms = samp2ms(obj,samp) 0827 % convert samples to ms at sampling rate for this object 0828 ms = (samp* (1000/obj.digheader.srate)); 0829 end 0830 0831 end % private methods 0832 0833 end % classdef 0834 0835 0836 0837 0838