Home > NoiseTools > nt_detrend.m

nt_detrend

PURPOSE ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robust remove polynomial or other trend

SYNOPSIS ^

function [y,w,r]=ntdetrend(x,order,w,basis,thresh,niter)

DESCRIPTION ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robust remove polynomial or other trend
 
  y: detrended data
  w: updated weights
  r: basis matrix used

  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(x,1); 
 Fit/remove polynomial trend with initial weighting w, threshold = 4*sd
   y=nt_detrend(x,5,w,3);
 Fit/remove linear then 3rd order polynomial:
   [y,w]=nt_detrend(x,1);
   [yy,ww]=nt_detrend(y,3);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 27-Feb-2017 15:36:07 by m2html © 2005