Home > matlabmk > trm_raw2nrm.m

trm_raw2nrm

PURPOSE ^

trm_raw2nrm() - Normalize EEG via estimated trimmed mean amplitude of

SYNOPSIS ^

function [usedcal_ids, scale_factor]=trm_raw2nrm(firstpulse,lastpulse,pulseamps,meancals,tme,Results,prewind,postwind)

DESCRIPTION ^

 trm_raw2nrm() - Normalize EEG via estimated trimmed mean amplitude of
                 calibration pulses
              
 Usage:
  >> usedcal_ids=trm_raw2nrm(firstpulse,lastpulse,pulseamps,meancals,tme,Results,prewind,postwind)

 Required Global Variables:
   SUBJECT.DATA.data = 3D matrix of EEG data, Chan x Sample x Epoch, in raw A/D units
   VERBLEVEL         = matlabMK level of verbosity (i.e., tells functions
                        how much to report about what they're doing during
                        runtime)

 Inputs:
   firstpulse  = 2D matrix of the first cal pulse at all channels (for
                  plotting diagnostics)
   lastpulse   = 2D matrix of the last cal pulse at all channels (for
                  plotting diagnostics)
   pulseamps   = matrix of amplitudes of each cal pulse (in A/D units, NOT in uV)
                  relative to the prepulse baseline. Matrix dimensions are
                  number of channels X number of pulses. Used by trm_raw2nrm.m to  
                  normalize EEG.
   meancals    = 2D matrix of the mean of all pulses (for plotting
                  diangostics)
   tme         = vector of peri-stimulus times (in msec) corresponding to
                  each time point of data (for plotting diagnostics)
   Results     = the input parser "p.Results" variable from crw2set.  The fields of Results
                        used in this function are:
                           Results.cal_npts
                           Results.low_cursor
                           Results.hi_cursor
                           Results.cal_amp
                           Results.cal_polarity
                           Results.trm_pcnt
                           Results.cal_plot_opt
   prewind     = Beginning and end of pre-pulse baseline window (in msec)
                  used for measuring pulse amplitude (for plotting
                  diagnostics)
   postwind    = Beginning and end of post-pulse baseline window (in msec)
                  used for measuring pulse amplitude (for plotting
                  diagnostics)

 Outputs:
   usedcal_ids  = 2D matrix, Chan x Number of Used Pulses, containing the column
                   numbers of pulseamps that were used to calibrate each
                   channel.  If Results.trm_pcnt>0, some extreme pulses
                   are not used.
   scale_factor = vector containing the values by which to scale each
                   channel to make the data in units of microvolts. For example, 
                   if you divide the first channel's data by the first
                   element of scale factor, its data will now be in
                   microvolts.

 Explanation of Results parameters and example values:
 
  Results.cal_npts     =  samples in window centered on low_cursor, hi_cursor
  Results.low_cursor   =  in ms relative to cal pulse onset
  Results.hi_cursor    =  in ms relative to cal pulse onset
  Results.cal_polarity =  +1 or -1
  Results.cal_amp      =  how big the calibration pulse is (in uV)
  Resutls.trm_pcnt     =  the percentage of extreme pulse amplitudes to trim.  This 
                           prevents outlier pulses from affecting mean pulse amplitude 
                           estimate.  The max trm_pcnt/2 and the min trm_pcnt/2 
                            percent of values are discarded.
  Results.cal_plot_opt =  if 1, create figures to help verify that
                           calibration has worked, else, no plots

  EXAMPLE:
  cal_npts     =  10
  low_cursor   = -50
  hi_cursor    =  50
  cal_amp      =  10
  cal_polarity =   1
  trm_pcnt     =   5
  cal_plot_opt =   1

 Similar functionality to normerp except with trm_pcnt option.

 Note: This function is based on Tom Urbach's mkp_raw2nrm.m function.  The two key differences 
 are that this function allows the trimming of extreme cal pulses and it assumes the data should 
 be converted into microvolts (i.e., you cannot arbitrarily set the point per 10 uV option).  
 Tom notes that scaling the data to microvolts is OK for processing in Matlab because it
 uses double precision floating point. However, when writing the data to an ERP avg file, the data
 should not be written as microvolts since it will be low precision.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [usedcal_ids, scale_factor]=trm_raw2nrm(firstpulse,lastpulse,pulseamps,meancals,tme,Results,prewind,postwind)
0002 % trm_raw2nrm() - Normalize EEG via estimated trimmed mean amplitude of
0003 %                 calibration pulses
0004 %
0005 % Usage:
0006 %  >> usedcal_ids=trm_raw2nrm(firstpulse,lastpulse,pulseamps,meancals,tme,Results,prewind,postwind)
0007 %
0008 % Required Global Variables:
0009 %   SUBJECT.DATA.data = 3D matrix of EEG data, Chan x Sample x Epoch, in raw A/D units
0010 %   VERBLEVEL         = matlabMK level of verbosity (i.e., tells functions
0011 %                        how much to report about what they're doing during
0012 %                        runtime)
0013 %
0014 % Inputs:
0015 %   firstpulse  = 2D matrix of the first cal pulse at all channels (for
0016 %                  plotting diagnostics)
0017 %   lastpulse   = 2D matrix of the last cal pulse at all channels (for
0018 %                  plotting diagnostics)
0019 %   pulseamps   = matrix of amplitudes of each cal pulse (in A/D units, NOT in uV)
0020 %                  relative to the prepulse baseline. Matrix dimensions are
0021 %                  number of channels X number of pulses. Used by trm_raw2nrm.m to
0022 %                  normalize EEG.
0023 %   meancals    = 2D matrix of the mean of all pulses (for plotting
0024 %                  diangostics)
0025 %   tme         = vector of peri-stimulus times (in msec) corresponding to
0026 %                  each time point of data (for plotting diagnostics)
0027 %   Results     = the input parser "p.Results" variable from crw2set.  The fields of Results
0028 %                        used in this function are:
0029 %                           Results.cal_npts
0030 %                           Results.low_cursor
0031 %                           Results.hi_cursor
0032 %                           Results.cal_amp
0033 %                           Results.cal_polarity
0034 %                           Results.trm_pcnt
0035 %                           Results.cal_plot_opt
0036 %   prewind     = Beginning and end of pre-pulse baseline window (in msec)
0037 %                  used for measuring pulse amplitude (for plotting
0038 %                  diagnostics)
0039 %   postwind    = Beginning and end of post-pulse baseline window (in msec)
0040 %                  used for measuring pulse amplitude (for plotting
0041 %                  diagnostics)
0042 %
0043 % Outputs:
0044 %   usedcal_ids  = 2D matrix, Chan x Number of Used Pulses, containing the column
0045 %                   numbers of pulseamps that were used to calibrate each
0046 %                   channel.  If Results.trm_pcnt>0, some extreme pulses
0047 %                   are not used.
0048 %   scale_factor = vector containing the values by which to scale each
0049 %                   channel to make the data in units of microvolts. For example,
0050 %                   if you divide the first channel's data by the first
0051 %                   element of scale factor, its data will now be in
0052 %                   microvolts.
0053 %
0054 % Explanation of Results parameters and example values:
0055 %
0056 %  Results.cal_npts     =  samples in window centered on low_cursor, hi_cursor
0057 %  Results.low_cursor   =  in ms relative to cal pulse onset
0058 %  Results.hi_cursor    =  in ms relative to cal pulse onset
0059 %  Results.cal_polarity =  +1 or -1
0060 %  Results.cal_amp      =  how big the calibration pulse is (in uV)
0061 %  Resutls.trm_pcnt     =  the percentage of extreme pulse amplitudes to trim.  This
0062 %                           prevents outlier pulses from affecting mean pulse amplitude
0063 %                           estimate.  The max trm_pcnt/2 and the min trm_pcnt/2
0064 %                            percent of values are discarded.
0065 %  Results.cal_plot_opt =  if 1, create figures to help verify that
0066 %                           calibration has worked, else, no plots
0067 %
0068 %  EXAMPLE:
0069 %  cal_npts     =  10
0070 %  low_cursor   = -50
0071 %  hi_cursor    =  50
0072 %  cal_amp      =  10
0073 %  cal_polarity =   1
0074 %  trm_pcnt     =   5
0075 %  cal_plot_opt =   1
0076 %
0077 % Similar functionality to normerp except with trm_pcnt option.
0078 %
0079 % Note: This function is based on Tom Urbach's mkp_raw2nrm.m function.  The two key differences
0080 % are that this function allows the trimming of extreme cal pulses and it assumes the data should
0081 % be converted into microvolts (i.e., you cannot arbitrarily set the point per 10 uV option).
0082 % Tom notes that scaling the data to microvolts is OK for processing in Matlab because it
0083 % uses double precision floating point. However, when writing the data to an ERP avg file, the data
0084 % should not be written as microvolts since it will be low precision.
0085 %
0086 
0087 % REVISION HISTORY
0088 %
0089 % 2009/8/17 - Fixed plot of normalized pulses to accomodate amps besides 10 uV.  Pulse windows are
0090 % now labeled with a legend instead of colored text. -DG
0091 %
0092 % 2009/8/21 - Modified to be more memory efficient. Now function only requires cal pulse amplitudes,
0093 % not the entire pulse. -DG
0094 %
0095 % 2009/10/07 - X-tick labels added to boxplots -DG
0096 %
0097 
0098 global VERBLEVEL;
0099 global SUBJECT;
0100 
0101 
0102 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0103 % THROW AN ERROR AND DIE FOR EMPTY CAL BIN
0104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0105 if isempty(pulseamps)
0106     error('trm_raw2nrm: Cals not found');
0107 end;
0108 
0109 
0110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0111 % CONVERTING TO uV ...
0112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0113 VerbReport(sprintf(['trm_raw2nrm: converting raw A/D to calibrated' ...
0114     ' uV with following parameters']),1,VERBLEVEL);
0115 VerbReport(sprintf('# of time points in mean windows: %d',Results.cal_npts),1,VERBLEVEL);
0116 VerbReport(sprintf('pre-pulse time window center (in msec): %d',Results.low_cursor),1,VERBLEVEL);
0117 VerbReport(sprintf('post-pulse time window center (in msec): %d',Results.hi_cursor),1,VERBLEVEL);
0118 VerbReport(sprintf('cal pulse amplitude: %d uV',Results.cal_amp),1,VERBLEVEL);
0119 if Results.cal_polarity==1,
0120     VerbReport('cal pulse polarity: positive',1, ...
0121         VERBLEVEL);
0122 elseif Results.cal_polarity==-1,
0123     VerbReport('cal pulse polarity: negative',1, ...
0124         VERBLEVEL);
0125 else
0126    error('cal_polarity must be 1 or -1'); 
0127 end
0128 VerbReport(sprintf('Percent of extreme pulses to trim: %d',Results.trm_pcnt),VERBLEVEL);
0129 %unit_type = 'uV'; %?? might be important someday for writing from Matlab to .avg files
0130 
0131 
0132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0133 % BEGIN PROCESSING IN EARNEST ...
0134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0135   
0136 VerbReport(sprintf('Averaging cal data'),2,VERBLEVEL);
0137 % FIGURE OUT WHICH PULSES TO TRIM FOR EACH ELECTRODE
0138 nCals=size(pulseamps,2);
0139 [nChan, nSamp, nEpochs]=size(SUBJECT.DATA.data);
0140 [pulseamps srt_id]=sort(pulseamps');
0141 cut_pls=round((Results.trm_pcnt*nCals)/200);  %number  of max and min pulses
0142 %to remove (approx. Results.trm_pcnt of pulses)
0143 nUse=length(cut_pls+1:nCals-cut_pls);
0144 true_trm_pcnt=100*(1-nUse/nCals)/2;
0145 VerbReport(sprintf('Removing max %.2f percent and min %.2f percent of pulses', ...
0146     true_trm_pcnt,true_trm_pcnt),2,VERBLEVEL);
0147 
0148 % FIND THE TRIMMED MEAN PULSE AMP
0149 useamps=pulseamps(cut_pls+1:nCals-cut_pls,:);
0150 trimmed_mn=Results.cal_polarity*mean(useamps)';
0151 
0152 % Keep track of which cals were actually used
0153 usedcal_ids=zeros(nChan,nUse);
0154 for dg=1:nChan,
0155     usedcal_ids(dg,:)=srt_id(cut_pls+1:nCals-cut_pls,dg)';
0156 end
0157   
0158 
0159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0160 % OPTIONALLY PLOT BEFORE ...
0161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0162 if Results.cal_plot_opt
0163     nrmfig = figure;
0164   
0165     cp = subplot(2,2,1);
0166     hold on
0167     title('Mean of all Loaded Pulses (No Trimming)');
0168     plot(tme,meancals');
0169     v=axis;
0170     axis([tme(1) tme(end) v(3:4)]);
0171     plot([1 1]*prewind(1),[min(min(meancals)) max(max(meancals))]*1.2,'k--');
0172     plot([1 1]*prewind(2),[min(min(meancals)) max(max(meancals))]*1.2,'k--');
0173     plot([1 1]*postwind(1),[min(min(meancals)) max(max(meancals))]*1.2,'b--');
0174     plot([1 1]*postwind(2),[min(min(meancals)) max(max(meancals))]*1.2,'b--');
0175     set(get(cp, 'YLabel'), 'String', 'A/D Units');
0176     set(get(cp, 'xLabel'), 'String', 'Time (msec)');
0177     set(gca,'ygrid','on');
0178     hold off
0179     
0180     cp = subplot(2,2,3);
0181     hold on
0182     plot(tme,mean(SUBJECT.DATA.data,3)');
0183     v=axis;
0184     axis([tme(1) tme(end) v(3:4)]);
0185     title('Mean Data Across all Epochs Before Norming');
0186     set(get(cp, 'xLabel'), 'String', 'Time (msec)');
0187     set(get(cp, 'YLabel'), 'String', 'A/D Units');
0188     set(gca,'ygrid','on');
0189     hold off
0190 end
0191 
0192 
0193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0194 % CONVERT TO MICROVOLTS, NORMERP TYPE A/D UNITS OR DIE ...
0195 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0196 
0197 %?? I think this line might be important when writing ERPs from Matlab to
0198 % .avg
0199 SUBJECT.DATA.pp10uv = 80;  % = 10uV * factor of 8 HANSEN-STYLE A/D ADJUSTMENT
0200 
0201 % APPLY CALIBRATION VALUE TO CONVERT SAMPLES TO uV ALL CHANNELS
0202 scale_factor=trimmed_mn/Results.cal_amp;
0203 cal_vals_ad_array = repmat(scale_factor,1,nSamp);
0204 
0205 % TRANSFORM EEG ...
0206 for e = 1: nEpochs
0207     %VerbReport(sprintf('Normalizing epoch %4d/%-5d', e, nEpochs),2,VERBLEVEL);
0208     SUBJECT.DATA.data(:,:,e) = squeeze(SUBJECT.DATA.data(:,:,e))./ ...
0209         cal_vals_ad_array;  
0210 end
0211   
0212 % TRANSFORM CALS ...
0213 meancals=meancals./cal_vals_ad_array;
0214 firstpulse=firstpulse./cal_vals_ad_array;
0215 lastpulse=lastpulse./cal_vals_ad_array;
0216 pulseamps=(pulseamps')./repmat(scale_factor,1,nCals);
0217 useamps=(useamps')./repmat(scale_factor,1,nUse);
0218 
0219 
0220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0221 % OPTIONALLY PLOT AFTER ...
0222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0223 if Results.cal_plot_opt    
0224     %%%%%%% Plotting post-calibration pulses and ERPs
0225     figure(nrmfig)
0226 
0227     subplot(2,2,2)
0228     hold on
0229     title('Mean of all Loaded Pulses After Norming (No Trimming)');
0230 
0231     % PLOT NORMALIZED CAL PULSES
0232     % plot time windows first so that they show up in legend
0233     plot([1 1]*prewind(1),[min(min(meancals)) max(max(meancals))]*1.2,'k--');
0234     plot([1 1]*postwind(1),[min(min(meancals)) max(max(meancals))]*1.2,'b--');
0235     plot(tme,meancals');
0236     v=axis;
0237     axis([tme(1) tme(end) v(3:4)]);
0238     xlabel('Time (msec)');
0239     ylabel('Microvolts');
0240     plot([1 1]*prewind(1),[min(min(meancals)) max(max(meancals))]*1.2,'k--');
0241     plot([1 1]*prewind(2),[min(min(meancals)) max(max(meancals))]*1.2,'k--');
0242     plot([1 1]*postwind(1),[min(min(meancals)) max(max(meancals))]*1.2,'b--');
0243     plot([1 1]*postwind(2),[min(min(meancals)) max(max(meancals))]*1.2,'b--');
0244     legend('Pre-Pulse Time Window','Post-Pulse Time Window','location','east');
0245     set(gca,'ygrid','on');
0246     hold off
0247 
0248     % PLOT NORMALIZED ERPS
0249     subplot(2,2,4)
0250     hold on
0251     plot(tme,squeeze(mean(SUBJECT.DATA.data,3))');
0252     title('Mean Data Across all Epochs After Norming');
0253     v=axis;
0254     axis([tme(1) tme(end) v(3:4)]);
0255     xlabel('Time (msec)');
0256     ylabel('Microvolts');
0257     set(gca,'ygrid','on');
0258     hold off
0259     
0260     h=textsc('Note: Each electrode is graphed as a separate line in each subplot', ...
0261          'title');
0262     set(h,'fontsize',14,'fontweight','bold');
0263     
0264     
0265     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0266     % NEW FIGURE
0267     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0268     % plot first and last cal pulses to make sure there was no initial
0269     % craziness (e.g., amplifier locked, outlier value)
0270     figure;
0271     subplot(4,1,1);
0272     plot(tme,firstpulse');
0273     h=title('First Cal Pulse');
0274     set(h,'fontsize',14,'fontweight','bold');
0275     hold on;
0276     v=axis;
0277     axis([tme(1) tme(end) v(3:4)]);
0278     h=plot([0 0],v(3:4),'k');
0279     set(h,'linewidth',3);
0280     h=xlabel('Time (ms)');
0281     set(gca,'ygrid','on');
0282     set(h,'fontsize',14);
0283     
0284     subplot(4,1,2);
0285     plot(tme,lastpulse');
0286     h=title('Last Cal Pulse');
0287     set(h,'fontsize',14,'fontweight','bold');
0288     hold on;
0289     v=axis;
0290     axis([tme(1) tme(end) v(3:4)]);
0291     h=plot([0 0],v(3:4),'k');
0292     set(h,'linewidth',3);
0293     set(gca,'ygrid','on');
0294     h=xlabel('Time (ms)');
0295     set(h,'fontsize',14);
0296 
0297     % PLOT BOXPLOT OF ALL CAL PULSES ACROSS ALL CHANNELS
0298     subplot(4,1,3);
0299     boxplot(pulseamps');
0300     set(gca,'ygrid','on','xtick',1:nChan,'xticklabel',1:nChan);
0301     hh=xlabel({' ','Electrode Number'});
0302     ylabel('Microvolts'); 
0303     h=title(['Boxplots of Mean Pulse Amplitude from ' int2str(postwind(1)) ...
0304          ' to ' int2str(postwind(2)) ' ms Post-Pulse']);
0305     set(h,'fontsize',14,'fontweight','bold');
0306     
0307     
0308     % PLOT BOXPLOT OF TRIMMED CAL PULSES ACROSS ALL CHANNELS
0309     subplot(4,1,4);
0310     boxplot(useamps');
0311     set(gca,'ygrid','on','xtick',1:nChan,'xticklabel',1:nChan);
0312     hh2=xlabel({' ','Electrode Number'}); %The extra lines are to make the xlabel
0313     %appear below the x-tick labels.  I tried set 'pos' to a lower value
0314     %but it wasn't working for some reason.
0315     ylabel('Microvolts'); 
0316     h=title(['Boxplots of ' num2str(200*cut_pls/nCals,3) ' Percent Trimmed Pulses']);
0317     set(h,'fontsize',14,'fontweight','bold');
0318   
0319 end    
0320 VerbReport(' ',1,VERBLEVEL);
0321 
0322   
0323

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005