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.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