cont_data_calibrate

 cont_data_calibrate() - Extract calibration scaling factor based on
                         calibration pulses in continuous EEG data matrix.
                         Cal pulses amplitude relative to a baseline is
                         measured for each pulse and then averaged using
                         a trimmed mean (or regular mean if requested).  
                         From this mean, a scaling factor is derived.
                         Multiplying the continuous data by this scaling 
                         factor will convert it to units of microvolts.
                         

 Usage:
 >>[scale_fact, used_cal_ids]=cont_data_calibrate(crwfile,logfile,bdffile,use_chanlocs,params);

 Inputs:
  crwfile       - Name of Kutaslab crw file
  logfile       - Name of Kutaslab log file
  bdf_file      - Name of Kutaslab bdf file
  use_chanlocs  - EEGLAB chanlocs structure indicating electrode coordinates 
                  and labels. This should be the EEG.chanlocs field of
                  the EEG variable containing the data you want to scale
                  to microvolts
  params        - Variable with the following fields:
      cal_npts:     # of time points to include in cal window (e.g., 10)
      low_cursor:   center of pre-pulse baseline window in msec (e.g., -50)
      hi_cursor:    center of post-pulse baseline window in msec (e.g.,
                    50)
      cal_amp:      cal pulse amplitude (in uV)
      cal_polarity: [1 or -1] cal pulse polarity. 1=+ pulses, -1=- pulses.
      trm_pcnt:     percentage (e.g., 5) of most extreme cal pulses to 
                    discard (used when computing trimmed mean of cal
                    pulses. If trm_pcnt=0, the standard mean is computed.
      cal_plot_opt: If non-zero, a figure with calibration diagnostic
                    plots will be produced
      cal_plot_avg: If non-zero and cal_plot_op is non-zero, the average 
                    of used cal pulses will be computed and displayed. 
                    THIS TAKES A LONG TIME.

 Outputs:
   scale_fact   - A vector of values by which to scale each channel of
                  EEG.data to convert it to microvolts. To convert the 
                  data in EEG.data to microvolts, do the following:
                  >> EEG.data=EEG.data.*repmat(scale_fact,1,n_tpts);
   used_cal_ids - A 2D (channel x cal pulse) matrix that indicates which
                  cal pulses were included in the estimation of the average
                  cal pulse amplitude. An value of "1" means that first cal
                  pulse detected in cont_data was used, etc...

 Notes:
 -This function assumes that any event that falls in the calibration pulse
 bin (as determined by the function import_blf.m) is a cal pulse.  The
 calibration pulse bin is the first bin with a condition code of 0.

 Author:
 David Groppe
 Kutaslab, 6/2011

0001 % cont_data_calibrate() - Extract calibration scaling factor based on
0002 %                         calibration pulses in continuous EEG data matrix.
0003 %                         Cal pulses amplitude relative to a baseline is
0004 %                         measured for each pulse and then averaged using
0005 %                         a trimmed mean (or regular mean if requested).
0006 %                         From this mean, a scaling factor is derived.
0007 %                         Multiplying the continuous data by this scaling
0008 %                         factor will convert it to units of microvolts.
0009 %
0010 %
0011 % Usage:
0012 % >>[scale_fact, used_cal_ids]=cont_data_calibrate(crwfile,logfile,bdffile,use_chanlocs,params);
0013 %
0014 % Inputs:
0015 %  crwfile       - Name of Kutaslab crw file
0016 %  logfile       - Name of Kutaslab log file
0017 %  bdf_file      - Name of Kutaslab bdf file
0018 %  use_chanlocs  - EEGLAB chanlocs structure indicating electrode coordinates
0019 %                  and labels. This should be the EEG.chanlocs field of
0020 %                  the EEG variable containing the data you want to scale
0021 %                  to microvolts
0022 %  params        - Variable with the following fields:
0023 %      cal_npts:     # of time points to include in cal window (e.g., 10)
0024 %      low_cursor:   center of pre-pulse baseline window in msec (e.g., -50)
0025 %      hi_cursor:    center of post-pulse baseline window in msec (e.g.,
0026 %                    50)
0027 %      cal_amp:      cal pulse amplitude (in uV)
0028 %      cal_polarity: [1 or -1] cal pulse polarity. 1=+ pulses, -1=- pulses.
0029 %      trm_pcnt:     percentage (e.g., 5) of most extreme cal pulses to
0030 %                    discard (used when computing trimmed mean of cal
0031 %                    pulses. If trm_pcnt=0, the standard mean is computed.
0032 %      cal_plot_opt: If non-zero, a figure with calibration diagnostic
0033 %                    plots will be produced
0034 %      cal_plot_avg: If non-zero and cal_plot_op is non-zero, the average
0035 %                    of used cal pulses will be computed and displayed.
0036 %                    THIS TAKES A LONG TIME.
0037 %
0038 % Outputs:
0039 %   scale_fact   - A vector of values by which to scale each channel of
0040 %                  EEG.data to convert it to microvolts. To convert the
0041 %                  data in EEG.data to microvolts, do the following:
0042 %                  >> EEG.data=EEG.data.*repmat(scale_fact,1,n_tpts);
0043 %   used_cal_ids - A 2D (channel x cal pulse) matrix that indicates which
0044 %                  cal pulses were included in the estimation of the average
0045 %                  cal pulse amplitude. An value of "1" means that first cal
0046 %                  pulse detected in cont_data was used, etc...
0047 %
0048 % Notes:
0049 % -This function assumes that any event that falls in the calibration pulse
0050 % bin (as determined by the function import_blf.m) is a cal pulse.  The
0051 % calibration pulse bin is the first bin with a condition code of 0.
0052 %
0053 % Author:
0054 % David Groppe
0055 % Kutaslab, 6/2011
0056 
0057 function [scale_fact, used_cal_ids]=cont_data_calibrate(crwfile,logfile,bdffile,use_chanlocs,params)
0058 
0059 global VERBLEVEL
0060 
0061 %% Import Log File Info
0062 VerbReport(' ',1,VERBLEVEL);
0063 % If ascii version of logfile not specified or if new log file created with artifact flags,
0064 % create temporary ascii file
0065 
0066 % Misc info that needs to be read from crw file
0067 srate=erpio2('get_srate',crwfile);
0068 
0069 VerbReport(sprintf('Creating temporary ascii version of logfile %s',logfile),1,VERBLEVEL);
0070 %remove temp.asci if it exists
0071 asci_cmnd='rm -f temp.asci';
0072 unix(asci_cmnd);
0073 asci_cmnd=sprintf('log2asci %s temp.asci -d %d',logfile,srate);
0074 [s, w]=unix(asci_cmnd);
0075 if s~=0,
0076     disp('Could not create ascii version of logfile using log2asci.');
0077     error('According to log2asci: %s',w);
0078 end
0079 asciilogfile='temp.asci';
0080 tempascii=1;
0081 VerbReport(sprintf('Getting log information from %s', asciilogfile), 1, VERBLEVEL);
0082 
0083 %import log information for events caught by blf file
0084 log_info=import_entire_log(asciilogfile);
0085 
0086 %If a temp ascii version of logfile was created, delete it
0087 if  tempascii,
0088     VerbReport(sprintf('Removing temporary ascii version of logfile.'),1,VERBLEVEL);
0089     asci_cmnd='rm -f temp.asci';
0090     [s, w]=unix(asci_cmnd);
0091     if s~=0,
0092         disp('Could not remove ascii version of logfile: temp.asci.');
0093         error('According to rm: %s',w);
0094     end
0095 end
0096 
0097 % Parse log file up by condition code
0098 [cc_bounds, cc_bounds_codes]=parse_log_by_cc(log_info);
0099 
0100 %% Figure out which electrodes to use
0101 % Get channel names from crw file
0102 n_chans=erpio2('get_nchans',crwfile);
0103 errstring=erpio2('get_errstr'); % check for error
0104 if ~isempty(errstring),
0105     error(errstring);
0106 end
0107 
0108 crw_chanlabels=cell(1,n_chans);
0109 for c=1:n_chans,
0110     chan_label=erpio2('get_chandes',crwfile,c-1); %have to subtract 1 because channel indexing starts at 0
0111     errstring=erpio2('get_errstr'); % check for error
0112     if ~isempty(errstring),
0113         error(errstring);
0114     end
0115     crw_chanlabels{c}=chan_label;
0116 end
0117 
0118 % get channel names from EEGLAB chanlocs variable
0119 n_chanloc_chans=length(use_chanlocs);
0120 chanloc_labels=cell(1,n_chanloc_chans);
0121 crwfile2chanloc_ids=zeros(1,n_chans); %vector used to map channels in crw file to channels in EEG variable
0122 for a=1:n_chanloc_chans,
0123     chanloc_labels{a}=use_chanlocs(a).labels;
0124     [in_crw_file, id]=ismember(chanloc_labels{a},crw_chanlabels);
0125     if ~in_crw_file,
0126         error('The cal pulse crw file %s does not have a channel for EEG variable channel %s.', ...
0127             crwfile,chanloc_labels{a});
0128     end
0129     crwfile2chanloc_ids(id)=a;
0130 end
0131 
0132 
0133 %% Load data corresponding to a condition code of 0 from crw File
0134 cont_data=[];
0135 eventcodes=[];
0136 timestamps=[];
0137 n_cc_bounds=size(cc_bounds,1);
0138 for a=1:n_cc_bounds,
0139     if ~cc_bounds_codes(a)
0140         % Only import data corresponding to condition codes of 0
0141         [cont_dataTEMP, eventcodesTEMP, timestampsTEMP]=erpio2('get_cdata', ...
0142             crwfile,'event',cc_bounds(a,1)-1,cc_bounds(a,2)-1,0); %note you have to subtract 1 because logfile event indices
0143         %start at 0
0144         
0145         cont_data=[cont_data cont_dataTEMP];
0146         eventcodes=[eventcodes eventcodesTEMP];
0147         timestamps=[timestamps timestampsTEMP];
0148     end
0149 end
0150 clear cont_dataTEMP eventcodesTEMP timestampsTEMP
0151 
0152 %% Map cont_data to the order and number of channels in EEG variable
0153 % Select only unignored channels
0154 use_chans=find(crwfile2chanloc_ids>0);
0155 cont_data=cont_data(use_chans,:);
0156 % Now make sure order of channels matches order of chans in EEG variable
0157 cont_data=cont_data(crwfile2chanloc_ids(use_chans),:);
0158 %Note, you could speed this up by changing the order of channels and
0159 %selecting channels after cal pulses have been extracted
0160 
0161 
0162 %% Import blf/bdf file information
0163 [blf_info, bin_LongDesc, blffile, blf_evnum, rtfile, cal_bin, cond_LongDesc, bin_CondCode]= ...
0164     import_blf(bdffile,logfile,srate,[]);
0165 
0166 
0167 %% Extract calibration pulses
0168 % cal_bin identifies calibration pulses
0169 [n_chan, n_tpts]=size(cont_data);
0170 bndry_times_in_ms=1000*log_info(log_info(:,2)==-16384,5); %times at which the dig machine was paused
0171 cal_blf_ids=find(blf_info(:,cal_bin+1));
0172 cal_log_ids=blf_evnum(cal_blf_ids);
0173 %cc0_ids=find(log_info(:,3)==0); %all events with a condition code of 0
0174 %n_cc0_ids=length(cc0_ids);
0175 n_cal_ids=length(cal_log_ids);
0176 diff=zeros(n_chan,n_cal_ids);
0177 pulse_log_ids=zeros(1,n_cal_ids);
0178 pulse_ct=0;
0179 if VERBLEVEL>1,
0180     fprintf('Searching data for cal pulses (i.e., items in Bin %d). Cal pulses found: ',cal_bin);
0181 end
0182 
0183 back_step=floor((params.cal_npts-1)/2);
0184 fwd_step=ceil((params.cal_npts-1)/2);
0185 for a=cal_log_ids,
0186     if (VERBLEVEL>1) && ~mod(pulse_ct,50),
0187         fprintf('%d ',pulse_ct);
0188     end
0189     
0190     % Did this event fall in a bin
0191     log_id=a+1; %have to add 1, since log item #'s start at 0
0192     
0193     % Time index of event
0194     event_time=1000*log_info(log_id,5); %need to convert log_info times into msec
0195     pre_cntr_id=find_tpt(event_time+params.low_cursor,timestamps); %center of pre-pulse baseline window
0196     post_cntr_id=find_tpt(event_time+params.hi_cursor,timestamps); %center of post-pulse baseline window
0197     
0198     % Make sure no recording boundaries occur between beginning of baseline
0199     % and end of cal pulse
0200     
0201     % Extract data
0202     if (pre_cntr_id-back_step)<1,
0203         watchit(sprintf('Cal pulse #%d is too close to boundary to be used.',find(cal_log_ids==a)));
0204     elseif (post_cntr_id+fwd_step)>n_tpts,
0205         watchit(sprintf('Cal pulse #%d is too close to boundary to be used.',find(cal_log_ids==a)));
0206     else
0207         min_time_needed=timestamps(pre_cntr_id-back_step);
0208         max_time_needed=timestamps(post_cntr_id+fwd_step);
0209         if sum( (bndry_times_in_ms>=min_time_needed).*(bndry_times_in_ms<=max_time_needed))
0210             watchit(sprintf('Cal pulse #%d is too close to boundary to be used.',find(cal_log_ids==a)));
0211         else
0212             %Use cal pulse
0213             mn_prepulse=mean(cont_data(:,(pre_cntr_id-back_step):(pre_cntr_id+fwd_step)),2);
0214             mn_postpulse=mean(cont_data(:,(post_cntr_id-back_step):(post_cntr_id+fwd_step)),2);
0215             
0216             % Pre/post difference
0217             pulse_ct=pulse_ct+1;
0218             diff(:,pulse_ct)=mn_postpulse-mn_prepulse;
0219             pulse_log_ids(pulse_ct)=a+1;
0220         end
0221     end
0222 end
0223 
0224 if VERBLEVEL>1
0225     fprintf('\n');
0226 end
0227 
0228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0229 % THROW AN ERROR AND DIE IF NO CALS
0230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0231 if pulse_ct==0,
0232     error('No calibration pulses found. Cal pulses are events with a condition code of 0 that fall in a bin according to the specified bdf/blf file.');
0233 end
0234 
0235 
0236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0237 % CONVERTING TO uV ...
0238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0239 VerbReport(sprintf(['cont_data_calibrate: converting raw A/D to calibrated' ...
0240     ' uV with the following parameters']),1,VERBLEVEL);
0241 VerbReport(sprintf('# of time points in mean windows: %d',params.cal_npts),1,VERBLEVEL);
0242 VerbReport(sprintf('pre-pulse time window center (in msec): %d',params.low_cursor),1,VERBLEVEL);
0243 VerbReport(sprintf('post-pulse time window center (in msec): %d',params.hi_cursor),1,VERBLEVEL);
0244 VerbReport(sprintf('cal pulse amplitude: %d uV',params.cal_amp),1,VERBLEVEL);
0245 if params.cal_polarity==1,
0246     VerbReport('cal pulse polarity: positive',1, ...
0247         VERBLEVEL);
0248 elseif params.cal_polarity==-1,
0249     VerbReport('cal pulse polarity: negative',1, ...
0250         VERBLEVEL);
0251 else
0252     error('cal_polarity must be 1 or -1');
0253 end
0254 VerbReport(sprintf('Percent of extreme pulses to trim: %d',params.trm_pcnt),VERBLEVEL);
0255 
0256 
0257 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0258 % COMPUTE TRIMMED MEANS
0259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0260 
0261 VerbReport(sprintf('Averaging cal pulses'),2,VERBLEVEL);
0262 % FIGURE OUT WHICH PULSES TO TRIM FOR EACH ELECTRODE
0263 [diff srt_id]=sort(diff(:,1:pulse_ct)');
0264 cut_pls=round((params.trm_pcnt*pulse_ct)/200); %number of max and min pulses
0265 %to remove (approx. params.trm_pcnt of pulses)
0266 nUse=length(cut_pls+1:pulse_ct-cut_pls);
0267 true_trm_pcnt=100*(1-nUse/pulse_ct)/2;
0268 VerbReport(sprintf('Removing max %.2f percent and min %.2f percent of pulses', ...
0269     true_trm_pcnt,true_trm_pcnt),2,VERBLEVEL);
0270 
0271 % FIND THE TRIMMED MEAN PULSE AMP
0272 useamps=diff(cut_pls+1:pulse_ct-cut_pls,:);
0273 scale_fact=params.cal_polarity*mean(useamps)';
0274 scale_fact=params.cal_amp./scale_fact;
0275 
0276 % Keep track of which cals were actually used
0277 used_cal_ids=zeros(n_chan,nUse);
0278 for dg=1:n_chan,
0279     used_cal_ids(dg,:)=srt_id(cut_pls+1:pulse_ct-cut_pls,dg)';
0280 end
0281 
0282 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0283 % OPTIONAL DIAGNOSTIC PLOTS
0284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0285 if params.cal_plot_opt
0286     if params.cal_plot_avg,
0287         n_subplot=5;
0288     else
0289         n_subplot=4;
0290     end
0291     
0292     %make the following optional inputs?
0293     show_prepulse=100; % Time before cal pulse to begin plotting in diagnostic plots (in msec)
0294     show_postpulse=700; % Time after cal pulse to begin plotting in diagnostic plots (in msec)
0295     
0296     % PLOT FIRST AND LAST CAL PULSES TO MAKE SURE THERE WAS NO INITIAL OR
0297     % FINAL CRAZINESS (e.g., amplifier locked, outlier value)
0298     
0299     % Collect waveforms of first pulse for plotting diagnostics
0300     show_pulse_ct=1;
0301     show_prepulse_tpt=[]; %set to empty to potentially flag failure
0302     tme=-show_prepulse:1000/srate:show_postpulse;
0303     n_show_tpts=length(tme);
0304     while show_pulse_ct<=pulse_ct,
0305         a=pulse_log_ids(show_pulse_ct);
0306         event_time=1000*log_info(a,5); %need to convert log_info times into msec
0307         show_prepulse_tpt=find_tpt(event_time-show_prepulse,timestamps);
0308         show_postpulse_tpt=find_tpt(event_time+show_postpulse,timestamps);
0309         if length(show_prepulse_tpt:show_postpulse_tpt)~=n_show_tpts
0310             watchit(sprintf('Used pulse #%d is too close to boundary to plot.',show_pulse_ct));
0311             show_prepulse_tpt=[]; %set to empty to potentially flag failure
0312         else
0313             break;
0314         end
0315         show_pulse_ct=show_pulse_ct+1;
0316     end
0317     if isempty(show_prepulse_tpt)
0318         error('Could not plot a cal pulse.  They are all too close to log file boundaries.');
0319     end
0320     firstpulse=cont_data(:,show_prepulse_tpt:show_postpulse_tpt);
0321     
0322     % Collect waveforms of last pulse for plotting diagnostics
0323     show_pulse_ct=pulse_ct;
0324     show_prepulse_tpt=[];  %set to empty to potentially flag failure
0325     while show_pulse_ct>=1,
0326         a=pulse_log_ids(show_pulse_ct);
0327         event_time=1000*log_info(a,5); %need to convert log_info times into msec
0328         show_prepulse_tpt=find_tpt(event_time-show_prepulse,timestamps);
0329         show_postpulse_tpt=find_tpt(event_time+show_postpulse,timestamps);
0330         if length(show_prepulse_tpt:show_postpulse_tpt)~=n_show_tpts
0331             watchit(sprintf('Used pulse #%d is too close to boundary to plot.',show_pulse_ct));
0332             show_prepulse_tpt=[]; %set to empty to potentially flag failure
0333         else
0334             break;
0335         end
0336         show_pulse_ct=show_pulse_ct-1;
0337     end
0338     if isempty(show_prepulse_tpt)
0339         error('Could not plot a cal pulse.  They are all too close to log file boundaries.');
0340     end
0341     lastpulse=cont_data(:,show_prepulse_tpt:show_postpulse_tpt);
0342     
0343     figure;
0344     subplot(n_subplot,1,1);
0345     plot(tme,firstpulse');
0346     h=title('First Cal Pulse (Uncalibrated)');
0347     set(h,'fontsize',14,'fontweight','bold');
0348     hold on;
0349     v=axis;
0350     axis([tme(1) tme(end) v(3:4)]);
0351     h=plot([0 0],v(3:4),'k');
0352     set(h,'linewidth',3);
0353     h=xlabel('Time (ms)');
0354     set(gca,'ygrid','on');
0355     set(h,'fontsize',14);
0356     
0357     subplot(n_subplot,1,2);
0358     plot(tme,lastpulse');
0359     h=title('Last Cal Pulse (Uncalibrated)');
0360     set(h,'fontsize',14,'fontweight','bold');
0361     hold on;
0362     v=axis;
0363     axis([tme(1) tme(end) v(3:4)]);
0364     h=plot([0 0],v(3:4),'k');
0365     set(h,'linewidth',3);
0366     set(gca,'ygrid','on');
0367     h=xlabel('Time (ms)');
0368     set(h,'fontsize',14);
0369     
0370     % PLOT BOXPLOT OF ALL CAL PULSES ACROSS ALL CHANNELS
0371     subplot(n_subplot,1,3);
0372     boxplot(diff.*repmat(scale_fact',pulse_ct,1));
0373     set(gca,'ygrid','on','xtick',1:n_chan,'xticklabel',1:n_chan);
0374     hh=xlabel({' ','Electrode Number'});
0375     ylabel('Microvolts');
0376     h=title('Boxplots of All Pulses (Calibrated)');
0377     set(h,'fontsize',14,'fontweight','bold');
0378     
0379     % PLOT BOXPLOT OF TRIMMED CAL PULSES ACROSS ALL CHANNELS
0380     subplot(n_subplot,1,4);
0381     boxplot(useamps.*repmat(scale_fact',nUse,1));
0382     set(gca,'ygrid','on','xtick',1:n_chan,'xticklabel',1:n_chan);
0383     hh2=xlabel({' ','Electrode Number'}); %The extra lines are to make the xlabel
0384     %appear below the x-tick labels.  I tried set 'pos' to a lower value
0385     %but it wasn't working for some reason.
0386     ylabel('Microvolts');
0387     h=title(['Boxplots of ' num2str(200*cut_pls/pulse_ct,3) ' Percent Trimmed Pulses (Calibrated)']);
0388     set(h,'fontsize',14,'fontweight','bold');
0389     
0390     if params.cal_plot_avg,
0391         % Compute mean of used cal pulses
0392         fprintf('Computing average used cal pulse (i.e., excluding trimmed pulses).\n');
0393         cal_erp=zeros(n_chan,n_show_tpts);
0394         cal_erp_ct=zeros(n_chan,1);
0395         for c=1:n_chan
0396             fprintf('Channel %d of %d.\n',c,n_chan);
0397             for b=1:nUse,
0398                 a=pulse_log_ids(used_cal_ids(c,b)); %log_info index of event
0399                 
0400                 % Time index of event
0401                 event_time=1000*log_info(a,5); %need to convert log_info times into msec
0402                 start_tpt=find_tpt(event_time-show_prepulse,timestamps); %center of pre-pulse baseline window
0403                 stop_tpt=find_tpt(event_time+show_postpulse,timestamps); %center of post-pulse baseline window
0404                 if length(start_tpt:stop_tpt)~=n_show_tpts,
0405                     watchit(sprintf('Pulse %d, Channel %d cannot be included in average since plotting epoch duration extends beyond recording boundaries.', ...
0406                         b,c));
0407                 else
0408                     cal_erp(c,:)=cal_erp(c,:)+cont_data(c,start_tpt:stop_tpt);
0409                     cal_erp_ct(c)=cal_erp_ct(c)+1;
0410                 end
0411             end
0412         end
0413         %convert sums to averages and scale to microvolts
0414         cal_erp=cal_erp./repmat(cal_erp_ct./scale_fact,1,n_show_tpts);
0415         %remove pre-pulse baseline
0416         center_tpt=find_tpt(params.low_cursor,tme);
0417         bsln_start_tpt=center_tpt-floor((params.cal_npts-1)/2);
0418         bsln_stop_tpt=center_tpt+ceil((params.cal_npts-1)/2);
0419         bsln_mn=mean(cal_erp(:,bsln_start_tpt:bsln_stop_tpt),2);
0420         cal_erp=cal_erp-repmat(bsln_mn,1,n_show_tpts);
0421         
0422         % PLOT AVERAGE CAL PULSE OF USED CAL PULSES
0423         subplot(n_subplot,1,5);
0424         plot(tme,cal_erp');
0425         h=title('Average of Used Cal Pulse (Baselined, Trimmed, Calibrated)');
0426         set(h,'fontsize',14,'fontweight','bold');
0427         hold on;
0428         v=axis;
0429         axis([tme(1) tme(end) v(3:4)]);
0430         h=plot([0 0],v(3:4),'k');
0431         set(h,'linewidth',3);
0432         set(gca,'ygrid','on');
0433         h=xlabel('Time (ms)');
0434         set(h,'fontsize',14);
0435     end
0436 end