Home > NoiseTools > nt_detrend_robust_old.m

# nt_detrend_robust_old

## PURPOSE

[y,w]=nt_detrend_robust(x,order,w,basis,thresh,percentage) - 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,percentage) - remove polynomial or sinusoidal trend

      y: detrended data

      x: raw data
      order: order of polynomial
      w: weight
      basis: 'polynomials' [default] or 'sinusoids'
      thresh: threshold for outliers [default: 2 sd]
      niter: number of iterations [default: 2]

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_detrend_robust [y,w]=nt_detrend_robust(x,order,w,basis,thresh,niter) - remove polynomial or sinusoidal trend
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1

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,percentage) - remove polynomial or sinusoidal trend
    0003 %
    0004 %  y: detrended data
    0005 %
    0006 %  x: raw data
    0007 %  order: order of polynomial
    0008 %  w: weight
    0009 %  basis: 'polynomials' [default] or 'sinusoids'
    0010 %  thresh: threshold for outliers [default: 2 sd]
    0011 %  niter: number of iterations [default: 2]
    0012
    0013
    0014
    0015 %% arguments
    0016 if nargin<2; error('!'); end
    0017 if nargin<3; w=[]; end
    0018 if nargin<4||isempty(basis); basis='polynomials'; end
    0019 if nargin<5||isempty(thresh); thresh=1.5; end
    0020 if nargin<6||isempty(niter); niter=4; end
    0021
    0022 dims=size(x);
    0023 x=x(:,:); % concatenates dims >= 2
    0024
    0025 %% regressors
    0026 switch basis
    0027     case 'polynomials'
    0028         r=zeros(size(x,1),numel(order));
    0029         lin=linspace(-1,1,size(x,1));
    0030         for k=1:order
    0031             r(:,k)=lin.^k;
    0032         end
    0033     case 'sinusoids'
    0034         r=zeros(size(x,1),numel(order)*2);
    0035         lin=linspace(-1,1,size(x,1));
    0036         for k=1:order
    0037             r(:,2*k-1)=sin(2*pi*k*lin/2);
    0038             r(:,2*k)=cos(2*pi*k*lin/2);
    0039         end
    0040     otherwise
    0041         error('!');
    0042 end
    0043 r=nt_demean(r);
    0044 r=nt_normcol(r);
    0045
    0046 %% give weights same geometry as data
    0047 if size(w)==size(x)
    0048     ; % OK
    0049 elseif isempty(w)
    0050     w=ones(size(x));
    0051 else
    0052     w=repmat(w,1,size(x,2)); % to each channel its own weight
    0053 end
    0054
    0055 figure(1); clf;
    0056 figure(2); clf;
    0057 %% remove trends
    0058 for iChan=1:size(x,2) % detrend each channel individually
    0059     xx=x(:,iChan);
    0060     ww=w(:,iChan);
    0061     for iIter=1:niter
    0062
    0063
    0064         www=ww;
    0065         for iRegressor = 1:size(r,2) % remove regressors one by one
    0066             % weight data, remove mean
    0067             xxx=xx.*www;
    0068             avg=sum(xxx)/sum(www);
    0069             xxx=xxx-avg;
    0070             xxx=xxx.*www;
    0071             % weight basis, remove mean, PCA, normalize
    0072             rr=r(:,iRegressor);
    0073             rr=nt_demean(rr,www);
    0074             rr=nt_normcol(rr,www);
    0075             a=xxx' * rr / sum(www); % cross product with basis
    0076             y=rr*a'; % project
    0077
    0078             y=y+avg; % restore mean
    0079             d=xx-y;
    0080             xx=d;
    0081             www=min(www,(abs(d)<=thresh*std(d)));
    0082             figure(1)
    0083             subplot 311; plot([xx,y]);
    0084             subplot 312; plot(d); hold on
    0085             subplot 313; plot(www); ylim([-0.1 1.1]); pause
    0086             disp(mean(d.^2.*www))
    0087         end
    0088         ww=min(ww,www);
    0089         xx=xx-y;
    0090         %figure(2); plot([xx]); hold on
    0091     end
    0092     x(:,iChan)=xx;
    0093     w(:,iChan)=ww;
    0094 end
    0095
    0096 x=reshape (x,dims);
    0097
    0098 %% test code
    0099 if 0
    0100     x=(1:100)';
    0101     nt_detrend_robust(x,2,[],[],[],[],[]);
    0102
    0103 end