Home > matlabmk > robust_locdetrend.m

robust_locdetrend

PURPOSE ^

robust_locdetrend() - Detrends the data in an EEGLAB EEG variable using

SYNOPSIS ^

function EEG=robust_locdetrend(EEG,robust,wind_duration,step_size,make_plots)

DESCRIPTION ^

 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/2011

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % 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

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