Home > NoiseTools > nt_detrend.m

nt_detrend

PURPOSE ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robustly remove 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) - robustly remove 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 order=5 with initial weighting w, threshold = 4*sd:
   y=nt_detrend(x,5,w,[],4);
 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) - robustly remove 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 order=5 with initial weighting w, threshold = 4*sd:
0030 %   y=nt_detrend(x,5,w,[],4);
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 if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0045 
0046 dims=size(x);
0047 x=x(:,:); % concatenates dims >= 2
0048 w=w(:,:);
0049 
0050 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0051 
0052 %% regressors
0053 if isnumeric(basis)
0054     r=basis;
0055 else
0056     switch basis
0057         case 'polynomials'
0058             r=zeros(size(x,1),numel(order));
0059             lin=linspace(-1,1,size(x,1));
0060             for k=1:order
0061                 r(:,k)=lin.^k;
0062             end
0063         case 'sinusoids'
0064             r=zeros(size(x,1),numel(order)*2);
0065             lin=linspace(-1,1,size(x,1));
0066             for k=1:order
0067                 r(:,2*k-1)=sin(2*pi*k*lin/2);
0068                 r(:,2*k)=cos(2*pi*k*lin/2);
0069             end
0070         otherwise
0071             error('!');
0072     end
0073 end
0074 %r=nt_normcol(nt_demean(r));
0075 
0076 %% remove trends
0077 
0078 % The tricky bit is to ensure that weighted means are removed before
0079 % calculating the regression (see nt_regw).
0080 
0081 for iIter=1:niter
0082     
0083     % weighted regression on basis
0084     [~,y]=nt_regw(x,r,w);
0085     
0086     % find outliers
0087     d=x-y;
0088     if ~isempty(w); d=bsxfun(@times,d,w); end
0089     ww=ones(size(x));
0090     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0091     
0092     % update weights
0093     if isempty(w); 
0094         w=ww;
0095     else
0096         w=min(w,ww);
0097     end
0098     
0099 end
0100 
0101 y=x-y;
0102 y=reshape(y,dims);
0103 w=reshape(w,dims);
0104 
0105 if ~nargout
0106     % don't return, just plot
0107     figure(1); clf;
0108     subplot 411; plot(x); title('raw');
0109     subplot 412; plot(y); title('detrended');
0110     subplot 413; plot(x-y); title('trend');
0111     subplot 414; nt_imagescc(w'); title('weight');
0112     clear y w r
0113 end
0114 
0115 
0116 
0117 
0118 
0119 %% test code
0120 if 0
0121     % basic
0122     x=(1:100)'; x=x+ randn(size(x));
0123     y=nt_detrend(x,1);
0124     figure(1); clf; plot([x,y]);
0125 end
0126 if 0
0127     % detrend biased random walk
0128     x=cumsum(randn(1000,1)+0.1);
0129     y=nt_detrend(x,3,[]);
0130     figure(1); clf; plot([x,y]); legend('before', 'after');
0131 end
0132 if 0
0133     % weights
0134     x=linspace(0,100,1000)';
0135     x=x+3*randn(size(x));
0136     x(1:100,:)=100;
0137     w=ones(size(x)); w(1:100,:)=0;
0138     y=nt_detrend(x,3,[],[],100);
0139     yy=nt_detrend(x,3,w);
0140     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0141 end
0142 if 0
0143     [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,:);
0144     %[p,x]=nt_read_data('/data/meg/arzounian/ADC_DA_140521_p20/ADC_DA_140521_p20_01_calib'); x=x'; x=x(1:10000,:);
0145     
0146     x=nt_demean(x);
0147     nt_detrend(x,10);   
0148 end
0149

Generated on Thu 30-Nov-2017 17:26:18 by m2html © 2005