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:

- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_detrend [y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
- nt_greetings nt_greetings - display message the first time the toolbox is used
- nt_imagescc nt_imagescc - plot image with symmetric scaling
- nt_read_data [p,data]=nt_read_data(fname,flag) - read data from file
- nt_regw [b,z]=nt_regw(y,x,w) - weighted regression

This function is called by:

- nt_deboing y=nt_deboing(x,events) - fit, remove ringing associated with events
- nt_detrend [y,w,r]=nt_detrend(x,order,w,basis,thresh,niter,wsize) - robustly remove trend
- nt_inpaint function y=nt_inpaint(x,w) - weighted interpolation based on correlation structure
- nt_outliers [w,y]=nt_outliers(x,w,thresh,niter) - detect outliers based on weighted correlation structure
- nt_star2 [y,w,ww]=nt_star2(x,thresh,closest,w) - sensor noise suppression

## SUBFUNCTIONS

- function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
- function test_me

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