Home > NoiseTools > nt_detrend.m

nt_detrend

PURPOSE ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend

SYNOPSIS ^

function [y,w,r]=ntdetrend(x,order,w,basis,thresh,niter,wsize)

DESCRIPTION ^

[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
  w: 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Tue 18-Feb-2020 11:23:12 by m2html © 2005