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:

- nt_detrend_robust [y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_regw [b,z]=nt_regw(y,x,w) - weighted regression

This function is called by:

- nt_detrend_robust [y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend
- nt_inpaint y=nt_inpaint(x,w) - interpolate missing samples of matrix
- nt_star2 [y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression

## 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