nt_detrend

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
 
  y: detrended data
  w: updated weights
  r: basis matrix used

  x: raw data
  order: order of polynomial or number of sin/cosine pairs
  w0: weights
  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
  thresh: threshold for outliers [default: 3 sd]
  niter: number of iterations [default: 3]
  wsize: window size for local detrending [default: all]

 This NEW (circa Oct 2019) version of detrend allows detrending to be
 applied to smaller overlapping windows, which are then stitched together
 using overlap-add. This is not described in the paper.

0001 function [y,w,r]=nt_detrend(x,order,w0,basis,thresh,niter,wsize)
0002 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
0003 %
0004 %  y: detrended data
0005 %  w: updated weights
0006 %  r: basis matrix used
0007 %
0008 %  x: raw data
0009 %  order: order of polynomial or number of sin/cosine pairs
0010 %  w0: weights
0011 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0012 %  thresh: threshold for outliers [default: 3 sd]
0013 %  niter: number of iterations [default: 3]
0014 %  wsize: window size for local detrending [default: all]
0015 %
0016 % This NEW (circa Oct 2019) version of detrend allows detrending to be
0017 % applied to smaller overlapping windows, which are then stitched together
0018 % using overlap-add. This is not described in the paper.
0019 
0020 nt_greetings;
0021 
0022 %% arguments
0023 if nargin<2||isempty(order); error('!'); end
0024 if nargin<3; w0=[]; end
0025 if nargin<4||isempty(basis); basis='polynomials'; end
0026 if nargin<5||isempty(thresh); thresh=3; end
0027 if nargin<6||isempty(niter); niter=3; end
0028 if nargin<7; wsize=[]; end
0029 
0030 if iscell(x)
0031     if ~isempty(w0); error('not implemented'); end
0032     y={}; w={}; r={};
0033     for iTrial=1:numel(x);
0034         [y{iTrial},w{iTrial},r{iTrial}]=nt_detrend(x{iTrial},order,w0,basis,thresh,niter,wsize);
0035     end
0036     return
0037 end
0038 
0039 w=w0;
0040 if ~isempty(w)
0041     w=w(:);
0042     if numel(w)<size(x,1)
0043         % assume that w contains indices, set them to 1
0044         idx=w;
0045         w=zeros(size(x,1),1);
0046         w(idx)=1;
0047     end
0048 end
0049 
0050 if isempty(wsize) || ~wsize
0051     % standard detrending (trend fit to entire data)
0052     [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter);
0053 else
0054     wsize=2*floor(wsize/2);
0055     
0056     % do some sanity checks because many parameters
0057     if numel(order)>1; error('!'); end
0058     if ~isempty(w) && ~(size(w,1)==size(x,1)) ; disp(size(w)); error('!'); end
0059     if ~(isempty(basis) || isstring(basis) || ~(isnumeric(basis)&&size(basis,1)==size(x,1))); error('!'); end
0060     if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0061     if numel(thresh)>1; error('!'); end
0062     if numel(niter)>1; error('!'); end
0063 
0064     dims=size(x); nchans=dims(2);
0065     x=x(:,:); % concatenates dims >= 2
0066     w=w(:,:);
0067     if isempty(w); w=ones(size(x)); end
0068     if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0069 
0070 
0071     % (1) divide into windows, (2) detrend each, (3) stitch together, (4)
0072     % estimate w
0073 
0074     for iIter=1:niter
0075 
0076         y=zeros(size(x));
0077         trend=zeros(size(x));
0078         a=zeros(size(x,1),1);
0079 
0080     %     figure(1); clf
0081 
0082         offset=0;
0083         while true
0084             start=offset+1;
0085             stop=min(size(x,1),offset+wsize);
0086 
0087             % if not enough valid samples grow window:
0088             MAXGROW=1000; % to avoid endless calculation (necessary?)
0089             counter=MAXGROW;
0090             while any( sum(w(start:stop,:)) < wsize )... % smaller than wsize?
0091                     || any(mean(w(start:stop,:))  < 1/4) % window has less than quarter ones?
0092                 if counter <= 0 ; break; end 
0093                 start=max(1,start-wsize/2);
0094                 stop=min(size(x,1),stop+wsize/2);
0095                 counter=counter-1;
0096             end
0097             if rem(stop-start+1,2)==1; stop=stop-1; end
0098             wsize2=stop-start+1;
0099 
0100             % detrend this window
0101             NITER=1;
0102             if ~sum(w(start:stop)); error('!'); end
0103             yy=nt_detrend_helper(x(start:stop,:),order,w(start:stop,:),basis,thresh,NITER);
0104 
0105             % triangular weighting
0106             if start==1
0107                 b=[ones(1,wsize2/2)*wsize2/2, wsize2/2:-1:1]';
0108             elseif stop==size(x,1)
0109                 b=[1:wsize2/2, ones(1,wsize2/2)*wsize2/2]';
0110             else
0111                 b=[1:wsize2/2, wsize2/2:-1:1]';
0112             end
0113 
0114             % overlap-add to output
0115             y(start:stop,:)=y(start:stop,:)+bsxfun(@times,yy,b);
0116             trend(start:stop,:)=trend(start:stop,:)+bsxfun(@times,x(start:stop,:)-yy,b);
0117 
0118             a(start:stop,1)=a(start:stop)+b;
0119 
0120             offset=offset+wsize/2;
0121             if offset>size(x,1)-wsize/2; break; end
0122         end
0123         
0124         if stop<size(x,1); y(end,:)=y(end-1,:); a(end,:)=a(end-1,:); end; % last sample can be missed
0125         
0126         y=bsxfun(@times,y,1./a); y(find(isnan(y)))=0;
0127         trend=bsxfun(@times,trend,1./a);  trend(find(isnan(trend)))=0;
0128 
0129         % find outliers
0130         d=x-trend; 
0131 
0132 
0133         if ~isempty(w); d=bsxfun(@times,d,w); end
0134         ww=ones(size(x));
0135         ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0136         clear d
0137 
0138         % update weights
0139         if isempty(w); 
0140             w=ww;
0141         else
0142             w=min(w,ww);
0143         end
0144         clear ww;
0145 
0146     end % for iIter...
0147     
0148     w=[];r=[]; % not informative
0149     
0150 end % if isempty(wsize)...
0151 
0152 if ~nargout
0153     % don't return, just plot
0154     subplot 411; plot(x); title('raw'); xlim([1,size(x,1)])
0155     subplot 412; plot(y); title('detrended'); xlim([1,size(x,1)])
0156     subplot 413; plot(x-y); title('trend'); xlim([1,size(x,1)])
0157     subplot 414; nt_imagescc(w'); title('weight');
0158     clear y w r
0159 end
0160 end
0161 
0162 %% Original version of detrend (no windows) is called by new version (windows)
0163 function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
0164 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robustly remove trend
0165 %
0166 %  y: detrended data
0167 %  w: updated weights
0168 %  r: basis matrix used
0169 %
0170 %  x: raw data
0171 %  order: order of polynomial or number of sin/cosine pairs
0172 %  w: weights
0173 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0174 %  thresh: threshold for outliers [default: 3 sd]
0175 %  niter: number of iterations [default: 3]
0176 %
0177 % Noise tools
0178 % See nt_regw().
0179 %
0180 % The data are fit to the basis using weighted least squares. The weight is
0181 % updated by setting samples for which the residual is greater than 'thresh'
0182 % times its std to zero, and the fit is repeated at most 'niter'-1 times.
0183 %
0184 % The choice of order (and basis) determines what complexity of the trend
0185 % that can be removed.  It may be useful to first detrend with a low order
0186 % to avoid fitting outliers, and then increase the order.
0187 %
0188 % Examples:
0189 % Fit linear trend, ignoring samples > 3*sd from it, and remove:
0190 %   y=nt_detrend(x,1);
0191 % Fit/remove polynomial order=5 with initial weighting w, threshold = 4*sd:
0192 %   y=nt_detrend(x,5,w,[],4);
0193 % Fit/remove linear then 3rd order polynomial:
0194 %   [y,w]=nt_detrend(x,1);
0195 %   [yy,ww]=nt_detrend(y,3);
0196 %
0197 nt_greetings;
0198 
0199 % arguments
0200 if nargin<2; error('!'); end
0201 if nargin<3; w=[]; end
0202 if nargin<4||isempty(basis); basis='polynomials'; end
0203 if nargin<5||isempty(thresh); thresh=3; end
0204 if nargin<6||isempty(niter); niter=3; end
0205 
0206 if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0207 
0208 dims=size(x);
0209 x=x(:,:); % concatenates dims >= 2
0210 w=w(:,:);
0211 
0212 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0213 
0214 %% regressors
0215 if isnumeric(basis)
0216     r=basis;
0217 else
0218     switch basis
0219         case 'polynomials'
0220             r=zeros(size(x,1),numel(order));
0221             lin=linspace(-1,1,size(x,1));
0222             for k=1:order
0223                 r(:,k)=lin.^k;
0224             end
0225         case 'sinusoids'
0226             r=zeros(size(x,1),numel(order)*2);
0227             lin=linspace(-1,1,size(x,1));
0228             for k=1:order
0229                 r(:,2*k-1)=sin(2*pi*k*lin/2);
0230                 r(:,2*k)=cos(2*pi*k*lin/2);
0231             end
0232         otherwise
0233             error('!');
0234     end
0235 end
0236 
0237 
0238 % remove trends
0239 % The tricky bit is to ensure that weighted means are removed before
0240 % calculating the regression (see nt_regw).
0241 
0242 for iIter=1:niter
0243     
0244     % weighted regression on basis
0245     [~,y]=nt_regw(x,r,w);
0246     
0247     % find outliers
0248     d=x-y; 
0249     if ~isempty(w); d=bsxfun(@times,d,w); end
0250     ww=ones(size(x));
0251     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0252      
0253     % update weights
0254     if isempty(w); 
0255         w=ww;
0256     else
0257         w=min(w,ww);
0258     end
0259     clear ww;    
0260 end
0261 y=x-y;
0262 y=reshape(y,dims);
0263 w=reshape(w,dims);
0264 
0265 end
0266 
0267 %% test code
0268 function test_me
0269 if 0
0270     % basic
0271     x=(1:100)'; x=x+ randn(size(x));
0272     WSIZE=30;
0273     y=nt_detrend2(x,1,[],[],[],[],WSIZE);
0274     figure(1); clf; plot([x,y]);
0275 end
0276 if 0
0277     % basic
0278     x=(1:100)'; x=x+ randn(size(x));
0279     x(40:50)=0;
0280     WSIZE=30;
0281     [yy,ww]=nt_detrend2(x,1,[],[],2,[],WSIZE);
0282     [y,w]=nt_detrend(x,1);
0283     figure(1); clf; subplot 211; 
0284     plot([x,y,yy]);
0285     subplot 212; plot([w,ww],'.-');
0286 end
0287 if 0
0288     % detrend biased random walk
0289     x=cumsum(randn(1000,1)+0.1);
0290     WSIZE=200;
0291     [y1,w1]=nt_detrend(x,1,[]);
0292     [y2,w2]=nt_detrend(x,1,[],[],[],[],WSIZE);
0293     figure(1); clf; 
0294     plot([x,y1,y2]); legend('before', 'global', 'local');
0295 end
0296 if 0
0297     % weights
0298     x=linspace(0,100,1000)';
0299     x=x+3*randn(size(x));
0300     x(1:100,:)=100;
0301     w=ones(size(x)); w(1:100,:)=0;
0302     y=nt_detrend2(x,3,[],[],[],[],WSIZE);
0303     yy=nt_detrend2(x,3,w,[],[],[],WSIZE);
0304     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0305 end
0306 if 0
0307     [p,x]=nt_read_data('/users/adc/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_43.mat'); x=x'; x=x(:,1:128); %x=x(1:10000,:);
0308     %[p,x]=nt_read_data('/data/meg/arzounian/ADC_DA_140521_p20/ADC_DA_140521_p20_01_calib'); x=x'; x=x(1:10000,:);
0309     
0310     x=nt_demean(x);
0311     figure(1);
0312     nt_detrend(x,3);
0313     figure(2);
0314     WSIZE=1000;
0315     nt_detrend2(x(:,:),3,[],[],[],[],WSIZE);
0316 end
0317 end
0318 
0319