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

Generated on Wed 12-Oct-2016 15:09:44 by m2html © 2005