Home > NoiseTools > nt_detrend_robust.m

nt_detrend_robust

PURPOSE ^

[y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend

SYNOPSIS ^

function [x,w]=ntdetrend_robust(x,order,w,basis,thresh,niter)

DESCRIPTION ^

[y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend
 
  y: detrended data
  w: updated weights

  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: 5]

 Noise tools
 See nt_regw().

 The data are fit to the basis using weighted least squares. The weight is
 updated by setting samples for which the residual is greater than 'thresh' 
 times its std to zero, and the fit is repeated at most 'niter'-1 times.

 The choice of order (and basis) determines what complexity of the trend
 that can be removed.  It may be useful to first detrend with a low order
 to avoid fitting outliers, and then increase the order.

 Examples:
 Fit linear trend, ignoring samples > 3*sd from it, and remove:
   y=nt_detrend_robust(x,1); 
 Fit/remove polynomial trend with initial weighting w, threshold = 4*sd
   y=nt_detrend_robust(x,5,w,3);
 Fit/remove linear then 3rd order polynomial:
   [y,w]=nt_detrend_robust(x,1);
   [yy,ww]=nt_detrend_robust(y,3);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [x,w]=ntdetrend_robust(x,order,w,basis,thresh,niter)
0002 %[y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend
0003 %
0004 %  y: detrended data
0005 %  w: updated weights
0006 %
0007 %  x: raw data
0008 %  order: order of polynomial or number of sin/cosine pairs
0009 %  w: weights
0010 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0011 %  thresh: threshold for outliers [default: 3 sd]
0012 %  niter: number of iterations [default: 5]
0013 %
0014 % Noise tools
0015 % See nt_regw().
0016 %
0017 % The data are fit to the basis using weighted least squares. The weight is
0018 % updated by setting samples for which the residual is greater than 'thresh'
0019 % times its std to zero, and the fit is repeated at most 'niter'-1 times.
0020 %
0021 % The choice of order (and basis) determines what complexity of the trend
0022 % that can be removed.  It may be useful to first detrend with a low order
0023 % to avoid fitting outliers, and then increase the order.
0024 %
0025 % Examples:
0026 % Fit linear trend, ignoring samples > 3*sd from it, and remove:
0027 %   y=nt_detrend_robust(x,1);
0028 % Fit/remove polynomial trend with initial weighting w, threshold = 4*sd
0029 %   y=nt_detrend_robust(x,5,w,3);
0030 % Fit/remove linear then 3rd order polynomial:
0031 %   [y,w]=nt_detrend_robust(x,1);
0032 %   [yy,ww]=nt_detrend_robust(y,3);
0033 %
0034 nt_greetings;
0035 
0036 %% arguments
0037 if nargin<2; error('!'); end
0038 if nargin<3; w=[]; end
0039 if nargin<4||isempty(basis); basis='polynomials'; end
0040 if nargin<5||isempty(thresh); thresh=3; end
0041 if nargin<6||isempty(niter); niter=4; end
0042 
0043 dims=size(x);
0044 x=x(:,:); % concatenates dims >= 2
0045 
0046 if size(w,2)==1; w=repmat(w,1,size(w,2)); end
0047 
0048 %% regressors
0049 if isnumeric(basis)
0050     r=basis;
0051 else
0052     switch basis
0053         case 'polynomials'
0054             r=zeros(size(x,1),numel(order));
0055             lin=linspace(-1,1,size(x,1));
0056             for k=1:order
0057                 r(:,k)=lin.^k;
0058             end
0059         case 'sinusoids'
0060             r=zeros(size(x,1),numel(order)*2);
0061             lin=linspace(-1,1,size(x,1));
0062             for k=1:order
0063                 r(:,2*k-1)=sin(2*pi*k*lin/2);
0064                 r(:,2*k)=cos(2*pi*k*lin/2);
0065             end
0066         otherwise
0067             error('!');
0068     end
0069 end
0070 %r=nt_normcol(nt_demean(r));
0071 
0072 %% remove trends
0073 
0074 % The tricky bit is to ensure that weighted means are removed before
0075 % calculating the regression (see nt_regw).
0076 
0077 for iIter=1:niter
0078     % regress on basis
0079     [~,y]=nt_regw(x,r,w);
0080     % residual
0081     d=x-y;
0082     if ~isempty(w); d=bsxfun(@times,d,w); end
0083     %figure(1); clf; plot(d); pause
0084     % find outliers
0085     ww=ones(size(x));
0086     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0087     % update weights
0088     if isempty(w); 
0089         w=ww;
0090     else
0091         w=min(w,ww);
0092     end
0093     
0094 end
0095 
0096 x=x-y;
0097 x=reshape(x,dims);
0098 
0099 
0100 %% test code
0101 if 0
0102     % basic
0103     x=(1:100)' + randn(size(x));
0104     y=nt_detrend_robust(x,1);
0105     figure(1); clf; plot([x,y]);
0106 end
0107 if 0
0108     % detrend biased random walk
0109     x=cumsum(randn(1000,1)+0.1);
0110     y=nt_detrend_robust(x,3,[]);
0111     figure(1); clf; plot([x,y]); legend('before', 'after');
0112 end
0113 if 0
0114     % weights
0115     x=cumsum(rand(100,1));
0116     x(1:10,:)=100;
0117     w=ones(size(x)); w(1:10,:)=0;
0118     y=nt_detrend_robust(x,3,[]);
0119     yy=nt_detrend_robust(x,3,w);
0120     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0121 end
0122

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005