[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);
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