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     %disp(iIter);
0084     %nt_whoss;
0085     
0086     % weighted regression on basis
0087     [~,y]=nt_regw(x,r,w);
0088     
0089     % find outliers
0090     d=x-y; 
0091     if ~isempty(w); d=bsxfun(@times,d,w); end
0092     ww=ones(size(x));
0093     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0094     clear d
0095     
0096     % update weights
0097     if isempty(w); 
0098         w=ww;
0099     else
0100         w=min(w,ww);
0101     end
0102     clear ww;
0103     
0104 end
0105 
0106 y=x-y;
0107 y=reshape(y,dims);
0108 w=reshape(w,dims);
0109 
0110 if ~nargout
0111     % don't return, just plot
0112     figure(1); clf;
0113     subplot 411; plot(x); title('raw');
0114     subplot 412; plot(y); title('detrended');
0115     subplot 413; plot(x-y); title('trend');
0116     subplot 414; nt_imagescc(w'); title('weight');
0117     clear y w r
0118 end
0119 
0120 
0121 
0122 
0123 
0124 %% test code
0125 if 0
0126     % basic
0127     x=(1:100)'; x=x+ randn(size(x));
0128     y=nt_detrend(x,1);
0129     figure(1); clf; plot([x,y]);
0130 end
0131 if 0
0132     % detrend biased random walk
0133     x=cumsum(randn(1000,1)+0.1);
0134     y=nt_detrend(x,3,[]);
0135     figure(1); clf; plot([x,y]); legend('before', 'after');
0136 end
0137 if 0
0138     % weights
0139     x=linspace(0,100,1000)';
0140     x=x+3*randn(size(x));
0141     x(1:100,:)=100;
0142     w=ones(size(x)); w(1:100,:)=0;
0143     y=nt_detrend(x,3,[],[],100);
0144     yy=nt_detrend(x,3,w);
0145     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0146 end
0147 if 0
0148     [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,:);
0149     %[p,x]=nt_read_data('/data/meg/arzounian/ADC_DA_140521_p20/ADC_DA_140521_p20_01_calib'); x=x'; x=x(1:10000,:);
0150     
0151     x=nt_demean(x);
0152     nt_detrend(x,10);   
0153 end
0154

Generated on Tue 09-Oct-2018 10:58:04 by m2html © 2005