robust_locdetrend() - Detrends the data in an EEGLAB EEG variable using
local linear regression and a linear regression.
Regression can be performed using ordinary least
squares regression (which is fast) or a method
robust to outliers ("least angles regression"
as implementd by the MATLAB function fit.m).
The way local linear regression works is that a line is
fit to all the data in a window, the window is then
moved a step, and the process is repeated. The
lines from all windows are then averaged together
via a weighted average to produce the estimated
slow drift that is subtracted from the recording.
Required Input:
EEG - An EEGLAB EEG variable containing epoched data.
Optional Inputs:
robust - ['yes' or 'no'] If 'yes', least angles regression (a
regression algorithm robust to outliers is used).
Otherwise, ordinary least squares regression is used.
{default: 'yes'}
wind_duration - The duration of the regression window in MILLISECONDS.
For local linear regression, this should be less than
the length of the epoch. The smaller this is, the
higher the frequencies that will be affected by the
filter. For strictly linear regression (i.e., NOT local
linear regression), leave this empty or make it equal
to the epoch length. {default: 0}
step_size - The amount to increment the regression window at each
step of the local linear regression process in
MILLISECONDS. The smaller this is, the higher the
frequencies that will be affected by the filter. As a
rule of thumb, this should not exceed wind_duration/2.
For strictly linear regression (i.e., NOT local linear
regression), leave this empty or make it 0. {default:
0)
make_plots - ['yes' or 'no'] If 'yes', a plot will be made for each
epoch and channel to illustrate the trend estimate
and effect of detrending. The plot for an epoch will be
held for two seconds and then it will be erased and
replaced by the plot for the next epoch. {default:
'no'}
Outputs:
EEG - Same as input EEG variable but with the following
fields modified:
EEG.data
EEG.history
EEG.saved
EEG.icaact
EEG.icspectra
EEG.icfreqs
Examples:
-To perform robust linear detrending
>>EEG=robust_locdetrend(EEG);
-To perform ordinary least squares linear detrending
>>EEG=robust_locdetrend(EEG,'no');
-To perform robust local linear detrending with a one second window and
a half second step between windows
>>EEG=robust_locdetrend(EEG,'yes',1000,500);
Notes:
When the length of the data is such that some data points at the end of
the epoch don't fall in a regression window. The last regression line is
simply extended to the end of the data.
Author:
David Groppe
Kutaslab, 3/20110001 % robust_locdetrend() - Detrends the data in an EEGLAB EEG variable using 0002 % local linear regression and a linear regression. 0003 % Regression can be performed using ordinary least 0004 % squares regression (which is fast) or a method 0005 % robust to outliers ("least angles regression" 0006 % as implementd by the MATLAB function fit.m). 0007 % The way local linear regression works is that a line is 0008 % fit to all the data in a window, the window is then 0009 % moved a step, and the process is repeated. The 0010 % lines from all windows are then averaged together 0011 % via a weighted average to produce the estimated 0012 % slow drift that is subtracted from the recording. 0013 % 0014 % Required Input: 0015 % EEG - An EEGLAB EEG variable containing epoched data. 0016 % 0017 % Optional Inputs: 0018 % robust - ['yes' or 'no'] If 'yes', least angles regression (a 0019 % regression algorithm robust to outliers is used). 0020 % Otherwise, ordinary least squares regression is used. 0021 % {default: 'yes'} 0022 % wind_duration - The duration of the regression window in MILLISECONDS. 0023 % For local linear regression, this should be less than 0024 % the length of the epoch. The smaller this is, the 0025 % higher the frequencies that will be affected by the 0026 % filter. For strictly linear regression (i.e., NOT local 0027 % linear regression), leave this empty or make it equal 0028 % to the epoch length. {default: 0} 0029 % step_size - The amount to increment the regression window at each 0030 % step of the local linear regression process in 0031 % MILLISECONDS. The smaller this is, the higher the 0032 % frequencies that will be affected by the filter. As a 0033 % rule of thumb, this should not exceed wind_duration/2. 0034 % For strictly linear regression (i.e., NOT local linear 0035 % regression), leave this empty or make it 0. {default: 0036 % 0) 0037 % make_plots - ['yes' or 'no'] If 'yes', a plot will be made for each 0038 % epoch and channel to illustrate the trend estimate 0039 % and effect of detrending. The plot for an epoch will be 0040 % held for two seconds and then it will be erased and 0041 % replaced by the plot for the next epoch. {default: 0042 % 'no'} 0043 % 0044 % Outputs: 0045 % EEG - Same as input EEG variable but with the following 0046 % fields modified: 0047 % EEG.data 0048 % EEG.history 0049 % EEG.saved 0050 % EEG.icaact 0051 % EEG.icspectra 0052 % EEG.icfreqs 0053 % 0054 % Examples: 0055 % -To perform robust linear detrending 0056 % >>EEG=robust_locdetrend(EEG); 0057 % 0058 % -To perform ordinary least squares linear detrending 0059 % >>EEG=robust_locdetrend(EEG,'no'); 0060 % 0061 % -To perform robust local linear detrending with a one second window and 0062 % a half second step between windows 0063 % >>EEG=robust_locdetrend(EEG,'yes',1000,500); 0064 % 0065 % Notes: 0066 % When the length of the data is such that some data points at the end of 0067 % the epoch don't fall in a regression window. The last regression line is 0068 % simply extended to the end of the data. 0069 % 0070 % Author: 0071 % David Groppe 0072 % Kutaslab, 3/2011 0073 0074 %%%%%%%%%%%%%%%% REVISION LOG %%%%%%%%%%%%%%%%% 0075 % 4/22/2011-'robust' and 'make_plots' options added 0076 % 0077 0078 %%%%%%%%%%%%%%%% Future Work %%%%%%%%%%%%%%%%% 0079 % Make work on continuous data 0080 % Add verblevel? 0081 0082 0083 function EEG=robust_locdetrend(EEG,robust,wind_duration,step_size,make_plots) 0084 0085 0086 %% Default Parameters 0087 if nargin<2, 0088 robust=1; %1=use robust regression 0089 elseif strcmpi(robust,'yes'), 0090 robust=1; 0091 elseif strcmpi(robust,'no'), 0092 robust=0; 0093 else 0094 error('Argument ''robust'' should be ''yes'' or ''no''.'); 0095 end 0096 if (nargin<3) 0097 wind_duration=0; 0098 end 0099 if (nargin<4) 0100 step_size=0; 0101 end 0102 if (nargin<5) 0103 make_plots=0; %set to 1 to illustrate trend estimation for each epoch and channel of data 0104 elseif strcmpi(make_plots,'yes'), 0105 make_plots=1; 0106 elseif strcmpi(make_plots,'no'), 0107 make_plots=0; 0108 else 0109 error('Argument ''make_plots'' should be ''yes'' or ''no''.'); 0110 end 0111 0112 %% Check Parameters 0113 %conver wind_duration and step_size to seconds 0114 wind_duration=wind_duration/1000; 0115 step_size=step_size/1000; 0116 if wind_duration<0, 0117 error('wind_duration needs to be greater than or equal to 0.'); 0118 end 0119 if (step_size<0) 0120 error('step_size needs to be greater than or equal 0.'); 0121 elseif (step_size>wind_duration), 0122 error('step_size needs to be smaller than or equal to wind_duration.'); 0123 end 0124 0125 %% Update history of EEG variable 0126 if ischar(EEG.history), 0127 temp=EEG.history; 0128 EEG.history=[]; 0129 EEG.history{1}=temp; 0130 end 0131 EEG.history{length(EEG.history)+1}=['EEG=robust_locdetrend(EEG,' num2str(wind_duration) ... 0132 ',' num2str(step_size) ');']; 0133 0134 if ~wind_duration 0135 %doing strictly linear detrending (not local linear detrending) 0136 wind_duration=length(EEG.times)/EEG.srate; 0137 step_size=wind_duration; 0138 end 0139 0140 % Are data continuous or epoched? 0141 s=size(EEG.data); 0142 if length(s)==3, 0143 continuous_data=0; 0144 nt=s(2); %number of time points per epoch 0145 Fs=EEG.srate; 0146 data_duration=nt./Fs; 0147 if (wind_duration>data_duration) 0148 error('wind_duration needs to be equal to or smaller than the duration of your epochs, %g.', ... 0149 data_duration); 0150 end 0151 else 0152 continuous_data=1; 0153 end 0154 0155 if make_plots 0156 h1=figure; 0157 end 0158 0159 %turn off robust regression warnings 0160 warning('off','curvefit:fit:iterationLimitReached'); 0161 warning('off','curvefit:fit:nonDoubleYData'); 0162 if continuous_data, 0163 error('This function can''t handle continuous data yet.'); 0164 else 0165 n=wind_duration*Fs; % number of time points in window 0166 dn=step_size*Fs; % number of time points to increment window at each step 0167 xwt=((1:n)-n/2)/(n/2); 0168 wt=(1-abs(xwt).^3).^3; %weighting function for combining regression lines 0169 nwin=1+floor((nt-n)/dn); % number of time windows per epoch 0170 if ~nwin 0171 error('No time windows span these data.'); 0172 end 0173 0174 if robust, 0175 fo_ = fitoptions('method','LinearLeastSquares','Robust','LAR'); 0176 ft_ = fittype('poly1'); 0177 end 0178 0179 for channel=1:s(1), 0180 fprintf('Robust detrending channel: %s\n',EEG.chanlocs(channel).labels); 0181 for epoch=1:s(3), 0182 y=EEG.data(channel,:,epoch)'; 0183 0184 y_line=zeros(nt,1); % trend estimate 0185 norm=y_line; %sum of weights at each time point (we'll divide by this to average the linear fits in each window 0186 0187 yfit=zeros(nwin,n); %linear fit for each moving window 0188 0189 if make_plots 0190 figure(h1); clf; 0191 subplot(1,2,1); 0192 %plot raw data in Subplot 1 0193 plot([EEG.times(1) EEG.times(end)],[0 0],'k'); hold on; 0194 plot(EEG.times,y,'b'); 0195 end 0196 0197 for j=1:nwin 0198 use_time_pt_ids=dn*(j-1)+1:dn*(j-1)+n; %data within moving time window 0199 tseg=y(use_time_pt_ids); 0200 0201 if robust, 0202 cf_=fit([1:n]',tseg,ft_,fo_); %robust regression 0203 a=cf_.p1; %slope 0204 b=cf_.p2; %intercept 0205 else 0206 %OLS regression parameters 0207 y1=mean(tseg); %the mean of the segment of data in the window 0208 y2=mean((1:n)'.*tseg)*2/(n+1); 0209 a=(y2-y1)*6/(n-1); %a=slope, 0210 b=y1-a*(n+1)/2; %b=intercept 0211 end 0212 0213 yfit(j,:)=(1:n)*a+b; %regression line 0214 y_line((j-1)*dn+(1:n))=y_line((j-1)*dn+(1:n))+(yfit(j,:).*wt)'; %weighted sum of regression lines 0215 norm((j-1)*dn+(1:n))=norm((j-1)*dn+(1:n))+wt'; %sum of the weights for combining regression lines 0216 0217 if make_plots 0218 figure(h1); 0219 subplot(1,2,1); 0220 hh=plot(EEG.times(use_time_pt_ids),yfit(j,:),'r'); 0221 set(hh,'linewidth',2); 0222 end 0223 0224 end 0225 mask=find(norm>0); %data points at the end might have a weight of 0, ignore them when normalizing 0226 y_line(mask)=y_line(mask)./norm(mask); 0227 indx=(nwin-1)*dn+n-1; %index of time points at the end that have a weight of 0 0228 npts=length(y)-indx+1; 0229 y_line(indx:end)=(n+1:n+npts)'*a+b; %flesh out final points with last estimated slope & intercept 0230 EEG.data(channel,:,epoch)=[y-y_line]'; %detrended data 0231 0232 if make_plots 0233 figure(h1); 0234 subplot(1,2,1); 0235 hh=plot(EEG.times,y_line,'g--'); 0236 set(hh,'linewidth',2); 0237 v=axis(); 0238 axis([EEG.times(1) EEG.times(end) v(3:4)]); 0239 [LEGH, OBJH]=legend('EEG','Local Linear Fits','Weighted Mean of Fits','location','best'); 0240 set(OBJH(4),'color','b','linewidth',2); 0241 set(OBJH(6),'color','r') 0242 set(OBJH(8),'color','g','linestyle','--'); 0243 hh=xlabel('Time (ms)'); 0244 hh=ylabel('\muV'); 0245 title('Raw Data'); 0246 0247 %plot detrended data in Subplot 2 0248 subplot(1,2,2); 0249 plot([EEG.times(1) EEG.times(end)],[0 0],'k'); hold on; 0250 hh=plot(EEG.times,EEG.data(channel,:,epoch),'b'); %detrended data 0251 %v=axis(); 0252 axis([EEG.times(1) EEG.times(end) v(3:4)-mean(v(3:4))]); 0253 hh=xlabel('Time (ms)'); 0254 hh=ylabel('\muV'); 0255 title('Local Linear Detrended Data'); 0256 0257 pause(2); 0258 end 0259 end 0260 end 0261 end 0262 %turn robust regression warnings back on 0263 warning('on','curvefit:fit:iterationLimitReached'); 0264 warning('on','curvefit:fit:nonDoubleYData'); 0265 0266 EEG.saved='no'; 0267 0268 % Recompute ICA activations and remove IC spectral info if they exist 0269 if ~isempty(EEG.icaact), 0270 EEG=recomp_act(EEG); 0271 0272 if isfield(EEG,'icspectra'), 0273 EEG.icspectra=[]; 0274 end 0275 0276 if isfield(EEG,'icfreqs'), 0277 EEG.icfreqs=[]; 0278 end 0279 end 0280 0281 function EEG=recomp_act(EEG) 0282 %function EEG=recomp_act(EEG) 0283 % 0284 % Recompute ICA activations 0285 0286 s=size(EEG.data); 0287 if length(s)==3, 0288 %epoched data 0289 EEG.icaact=reshape(EEG.data,s(1),s(2)*s(3)); 0290 EEG.icaact=EEG.icaweights*EEG.icasphere*EEG.icaact; 0291 EEG.icaact=reshape(EEG.icaact,s(1),s(2),s(3)); 0292 else 0293 %continuous data 0294 EEG.icaact=EEG.icaweights*EEG.icasphere*EEG.data; 0295 end 0296 0297 0298