Home > NoiseTools > nt_detrend.m

nt_detrend

PURPOSE ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend

SYNOPSIS ^

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

DESCRIPTION ^

[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - 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: 3]
  wsize: window size for local detrending [default: all]

 This NEW (circa Oct 2019) version of detrend allows detrending to be
 applied to smaller overlapping windows, which are then stitched together
 using overlap-add.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [y,w,r]=ntdetrend(x,order,w,basis,thresh,niter,wsize)
0002 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - 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: 3]
0014 %  wsize: window size for local detrending [default: all]
0015 %
0016 % This NEW (circa Oct 2019) version of detrend allows detrending to be
0017 % applied to smaller overlapping windows, which are then stitched together
0018 % using overlap-add.
0019 
0020 nt_greetings;
0021 
0022 %% arguments
0023 if nargin<2; error('!'); end
0024 if nargin<3; w=[]; end
0025 if nargin<4||isempty(basis); basis='polynomials'; end
0026 if nargin<5||isempty(thresh); thresh=3; end
0027 if nargin<6||isempty(niter); niter=3; end
0028 if nargin<7; wsize=[]; end
0029 
0030 if isempty(wsize) || ~wsize;
0031     % standard detrending (trend fit to entire data)
0032     [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter);
0033 else
0034     wsize=2*floor(wsize/2);
0035     % do some sanity checks because many parameters
0036     if numel(order)>1; error('!'); end
0037     if ~isempty(w) && ~(size(w,1)==size(x,1)) ; disp(size(w)); error('!'); end
0038     if ~(isempty(basis) || isstring(basis) || ~(isnumeric(basis)&&size(basis,1)==size(x,1))); error('!'); end
0039     if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0040     if numel(thresh)>1; error('!'); end
0041     if numel(niter)>1; error('!'); end
0042 
0043     dims=size(x); nchans=dims(2);
0044     x=x(:,:); % concatenates dims >= 2
0045     w=w(:,:);
0046     if isempty(w); w=ones(size(x)); end
0047     if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0048 
0049 
0050     % (1) divide into windows, (2) detrend each, (3) stitch together, (4)
0051     % estimate w
0052 
0053     for iIter=1:niter
0054 
0055         y=zeros(size(x));
0056         trend=zeros(size(x));
0057         a=zeros(size(x,1),1);
0058 
0059     %     figure(1); clf
0060 
0061         offset=0;
0062         while true
0063             start=offset+1;
0064             stop=min(size(x,1),offset+wsize);
0065 
0066             % if not enough valid samples grow window:
0067             counter=5;
0068             while any (sum(min(w(start:stop),2))) <wsize
0069                 if counter <= 0 ; break; end 
0070                 start=max(1,start-wsize/2);
0071                 stop=min(size(x,1),stop+wsize/2);
0072                 counter=counter-1;
0073             end
0074             if rem(stop-start+1,2)==1; stop=stop-1; end
0075             wsize2=stop-start+1;
0076 
0077             % detrend this window
0078             NITER=1;
0079             yy=nt_detrend_helper(x(start:stop,:),order,w(start:stop,:),basis,thresh,NITER);
0080 
0081             % triangular weighting
0082             if start==1
0083                 b=[ones(1,wsize2/2)*wsize2/2, wsize2/2:-1:1]';
0084             elseif stop==size(x,1)
0085                 b=[1:wsize2/2, ones(1,wsize2/2)*wsize2/2]';
0086             else
0087                 b=[1:wsize2/2, wsize2/2:-1:1]';
0088             end
0089 
0090             % overlap-add to output
0091             y(start:stop,:)=y(start:stop,:)+bsxfun(@times,yy,b);
0092             trend(start:stop,:)=trend(start:stop,:)+bsxfun(@times,x(start:stop,:)-yy,b);
0093 
0094             a(start:stop,1)=a(start:stop)+b;
0095 
0096             offset=offset+wsize/2;
0097             if offset>size(x,1)-wsize/2; break; end
0098         end
0099         
0100         if stop<size(x,1); y(end,:)=y(end-1,:); a(end,:)=a(end-1,:); end; % last sample can be missed
0101         
0102         y=bsxfun(@times,y,1./a); y(find(isnan(y)))=0;
0103         trend=bsxfun(@times,trend,1./a);  trend(find(isnan(trend)))=0;
0104 
0105         % find outliers
0106         d=x-trend; 
0107 
0108 
0109         if ~isempty(w); d=bsxfun(@times,d,w); end
0110         ww=ones(size(x));
0111         ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0112         clear d
0113 
0114         % update weights
0115         if isempty(w); 
0116             w=ww;
0117         else
0118             w=min(w,ww);
0119         end
0120         clear ww;
0121 
0122     end % for iIter...
0123 end % if isempty(wsize)...
0124 
0125 if ~nargout
0126     % don't return, just plot
0127     subplot 411; plot(x); title('raw'); xlim([1,size(x,1)])
0128     subplot 412; plot(y); title('detrended'); xlim([1,size(x,1)])
0129     subplot 413; plot(x-y); title('trend'); xlim([1,size(x,1)])
0130     subplot 414; nt_imagescc(w'); title('weight');
0131     clear y w r
0132 end
0133 end
0134 
0135 %% Original version of detrend (no windows) is called by new version (windows)
0136 function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
0137 %[y,w,r]=nt_detrend(x,order,w,basis,thresh,niter) - robustly remove trend
0138 %
0139 %  y: detrended data
0140 %  w: updated weights
0141 %  r: basis matrix used
0142 %
0143 %  x: raw data
0144 %  order: order of polynomial or number of sin/cosine pairs
0145 %  w: weights
0146 %  basis: 'polynomials' [default] or 'sinusoids', or user-provided matrix
0147 %  thresh: threshold for outliers [default: 3 sd]
0148 %  niter: number of iterations [default: 3]
0149 %
0150 % Noise tools
0151 % See nt_regw().
0152 %
0153 % The data are fit to the basis using weighted least squares. The weight is
0154 % updated by setting samples for which the residual is greater than 'thresh'
0155 % times its std to zero, and the fit is repeated at most 'niter'-1 times.
0156 %
0157 % The choice of order (and basis) determines what complexity of the trend
0158 % that can be removed.  It may be useful to first detrend with a low order
0159 % to avoid fitting outliers, and then increase the order.
0160 %
0161 % Examples:
0162 % Fit linear trend, ignoring samples > 3*sd from it, and remove:
0163 %   y=nt_detrend(x,1);
0164 % Fit/remove polynomial order=5 with initial weighting w, threshold = 4*sd:
0165 %   y=nt_detrend(x,5,w,[],4);
0166 % Fit/remove linear then 3rd order polynomial:
0167 %   [y,w]=nt_detrend(x,1);
0168 %   [yy,ww]=nt_detrend(y,3);
0169 %
0170 nt_greetings;
0171 
0172 % arguments
0173 if nargin<2; error('!'); end
0174 if nargin<3; w=[]; end
0175 if nargin<4||isempty(basis); basis='polynomials'; end
0176 if nargin<5||isempty(thresh); thresh=3; end
0177 if nargin<6||isempty(niter); niter=3; end
0178 
0179 if thresh==0; error('thresh=0 is not what you want...'); end % common mistake
0180 
0181 dims=size(x);
0182 x=x(:,:); % concatenates dims >= 2
0183 w=w(:,:);
0184 
0185 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0186 
0187 %% regressors
0188 if isnumeric(basis)
0189     r=basis;
0190 else
0191     switch basis
0192         case 'polynomials'
0193             r=zeros(size(x,1),numel(order));
0194             lin=linspace(-1,1,size(x,1));
0195             for k=1:order
0196                 r(:,k)=lin.^k;
0197             end
0198         case 'sinusoids'
0199             r=zeros(size(x,1),numel(order)*2);
0200             lin=linspace(-1,1,size(x,1));
0201             for k=1:order
0202                 r(:,2*k-1)=sin(2*pi*k*lin/2);
0203                 r(:,2*k)=cos(2*pi*k*lin/2);
0204             end
0205         otherwise
0206             error('!');
0207     end
0208 end
0209 
0210 
0211 % remove trends
0212 % The tricky bit is to ensure that weighted means are removed before
0213 % calculating the regression (see nt_regw).
0214 
0215 for iIter=1:niter
0216     
0217     % weighted regression on basis
0218     [~,y]=nt_regw(x,r,w);
0219     
0220     % find outliers
0221     d=x-y; 
0222     if ~isempty(w); d=bsxfun(@times,d,w); end
0223     ww=ones(size(x));
0224     ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0225      
0226     % update weights
0227     if isempty(w); 
0228         w=ww;
0229     else
0230         w=min(w,ww);
0231     end
0232     clear ww;    
0233 end
0234 y=x-y;
0235 y=reshape(y,dims);
0236 w=reshape(w,dims);
0237 
0238 end
0239 
0240 %% test code
0241 function test_me
0242 if 0
0243     % basic
0244     x=(1:100)'; x=x+ randn(size(x));
0245     WSIZE=30;
0246     y=nt_detrend2(x,1,[],[],[],[],WSIZE);
0247     figure(1); clf; plot([x,y]);
0248 end
0249 if 0
0250     % basic
0251     x=(1:100)'; x=x+ randn(size(x));
0252     x(40:50)=0;
0253     WSIZE=30;
0254     [yy,ww]=nt_detrend2(x,1,[],[],2,[],WSIZE);
0255     [y,w]=nt_detrend(x,1);
0256     figure(1); clf; subplot 211; 
0257     plot([x,y,yy]);
0258     subplot 212; plot([w,ww],'.-');
0259 end
0260 if 0
0261     % detrend biased random walk
0262     x=cumsum(randn(1000,1)+0.1);
0263     WSIZE=200;
0264     [y1,w1]=nt_detrend(x,1,[]);
0265     [y2,w2]=nt_detrend2(x,1,[],[],[],[],WSIZE);
0266     figure(1); clf; 
0267     plot([x,y1,y2]); legend('before', 'after');
0268 end
0269 if 0
0270     % weights
0271     x=linspace(0,100,1000)';
0272     x=x+3*randn(size(x));
0273     x(1:100,:)=100;
0274     w=ones(size(x)); w(1:100,:)=0;
0275     y=nt_detrend2(x,3,[],[],[],[],WSIZE);
0276     yy=nt_detrend2(x,3,w,[],[],[],WSIZE);
0277     figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0278 end
0279 if 0
0280     [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_43.mat'); x=x'; x=x(:,1:128); %x=x(1:10000,:);
0281     %[p,x]=nt_read_data('/data/meg/arzounian/ADC_DA_140521_p20/ADC_DA_140521_p20_01_calib'); x=x'; x=x(1:10000,:);
0282     
0283     x=nt_demean(x);
0284     figure(1);
0285     nt_detrend(x,3);
0286     figure(2);
0287     WSIZE=1000;
0288     nt_detrend2(x(:,:),3,[],[],[],[],WSIZE);
0289 end
0290 end
0291 
0292

Generated on Thu 07-May-2020 08:42:54 by m2html © 2005