0001 function [x,w]=ntdetrend_robust(x,order,w,basis,thresh,niter)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 if nargin<2; error('!'); end
0017 if nargin<3; w=[]; end
0018 if nargin<4||isempty(basis); basis='polynomials'; end
0019 if nargin<5||isempty(thresh); thresh=1.5; end
0020 if nargin<6||isempty(niter); niter=4; end
0021
0022 dims=size(x);
0023 x=x(:,:);
0024
0025
0026 switch basis
0027 case 'polynomials'
0028 r=zeros(size(x,1),numel(order));
0029 lin=linspace(-1,1,size(x,1));
0030 for k=1:order
0031 r(:,k)=lin.^k;
0032 end
0033 case 'sinusoids'
0034 r=zeros(size(x,1),numel(order)*2);
0035 lin=linspace(-1,1,size(x,1));
0036 for k=1:order
0037 r(:,2*k-1)=sin(2*pi*k*lin/2);
0038 r(:,2*k)=cos(2*pi*k*lin/2);
0039 end
0040 otherwise
0041 error('!');
0042 end
0043 r=nt_demean(r);
0044 r=nt_normcol(r);
0045
0046
0047 if size(w)==size(x)
0048 ;
0049 elseif isempty(w)
0050 w=ones(size(x));
0051 else
0052 w=repmat(w,1,size(x,2));
0053 end
0054
0055 figure(1); clf;
0056 figure(2); clf;
0057
0058 for iChan=1:size(x,2)
0059 xx=x(:,iChan);
0060 ww=w(:,iChan);
0061 for iIter=1:niter
0062
0063
0064 www=ww;
0065 for iRegressor = 1:size(r,2)
0066
0067 xxx=xx.*www;
0068 avg=sum(xxx)/sum(www);
0069 xxx=xxx-avg;
0070 xxx=xxx.*www;
0071
0072 rr=r(:,iRegressor);
0073 rr=nt_demean(rr,www);
0074 rr=nt_normcol(rr,www);
0075 a=xxx' * rr / sum(www);
0076 y=rr*a';
0077
0078 y=y+avg;
0079 d=xx-y;
0080 xx=d;
0081 www=min(www,(abs(d)<=thresh*std(d)));
0082 figure(1)
0083 subplot 311; plot([xx,y]);
0084 subplot 312; plot(d); hold on
0085 subplot 313; plot(www); ylim([-0.1 1.1]); pause
0086 disp(mean(d.^2.*www))
0087 end
0088 ww=min(ww,www);
0089 xx=xx-y;
0090
0091 end
0092 x(:,iChan)=xx;
0093 w(:,iChan)=ww;
0094 end
0095
0096 x=reshape (x,dims);
0097
0098
0099 if 0
0100 x=(1:100)';
0101 nt_detrend_robust(x,2,[],[],[],[],[]);
0102
0103 end