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