[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);
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