[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 %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