0001 function [y,w,r]=ntdetrend(x,order,w,basis,thresh,niter,wsize)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 nt_greetings;
0021
0022
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
0032 [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter);
0033 else
0034 wsize=2*floor(wsize/2);
0035
0036 if numel(order)>1; error('!'); end
0037 if ~(isempty(w) || (isnumeric(w) && size(w,1)==size(x,1))) ; disp(size(w)); error('!'); end
0038 if ~(isempty(basis) || isstring(basis) || ~(isnumeric(basis)&&size(basis,1)==size(x,1))); error('!'); end
0039 if thresh==0; error('thresh=0 is not what you want...'); end
0040 if numel(thresh)>1; error('!'); end
0041 if numel(niter)>1; error('!'); end
0042
0043 dims=size(x); nchans=dims(2);
0044 x=x(:,:);
0045 w=w(:,:);
0046 if isempty(w); w=ones(size(x)); end
0047 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0048
0049
0050
0051
0052
0053 for iIter=1:niter
0054
0055 y=zeros(size(x));
0056 trend=zeros(size(x));
0057 a=zeros(size(x,1),1);
0058
0059
0060
0061 offset=0;
0062 while true
0063 start=offset+1;
0064 stop=min(size(x,1),offset+wsize);
0065
0066
0067 counter=5;
0068 while any (sum(min(w(start:stop),2))) <wsize
0069 if counter <= 0 ; break; end
0070 start=max(1,start-wsize/2);
0071 stop=min(size(x,1),stop+wsize/2);
0072 counter=counter-1;
0073 end
0074 if rem(stop-start+1,2)==1; stop=stop-1; end
0075 wsize2=stop-start+1;
0076
0077
0078 NITER=1;
0079 yy=nt_detrend_helper(x(start:stop,:),order,w(start:stop,:),basis,thresh,NITER);
0080
0081
0082 if start==1
0083 b=[ones(1,wsize2/2)*wsize2/2, wsize2/2:-1:1]';
0084 elseif stop==size(x,1)
0085 b=[1:wsize2/2, ones(1,wsize2/2)*wsize2/2]';
0086 else
0087 b=[1:wsize2/2, wsize2/2:-1:1]';
0088 end
0089
0090
0091 y(start:stop,:)=y(start:stop,:)+bsxfun(@times,yy,b);
0092 trend(start:stop,:)=trend(start:stop,:)+bsxfun(@times,x(start:stop,:)-yy,b);
0093
0094 a(start:stop,1)=a(start:stop)+b;
0095
0096 offset=offset+wsize/2;
0097 if offset>size(x,1)-wsize/2; break; end
0098 end
0099 y=bsxfun(@times,y,1./a); y(find(isnan(y)))=0;
0100 trend=bsxfun(@times,trend,1./a); trend(find(isnan(trend)))=0;
0101
0102
0103 d=x-trend;
0104
0105
0106 if ~isempty(w); d=bsxfun(@times,d,w); end
0107 ww=ones(size(x));
0108 ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0109 clear d
0110
0111
0112 if isempty(w);
0113 w=ww;
0114 else
0115 w=min(w,ww);
0116 end
0117 clear ww;
0118
0119 end
0120 end
0121
0122 if ~nargout
0123
0124 subplot 411; plot(x); title('raw'); xlim([1,size(x,1)])
0125 subplot 412; plot(y); title('detrended'); xlim([1,size(x,1)])
0126 subplot 413; plot(x-y); title('trend'); xlim([1,size(x,1)])
0127 subplot 414; nt_imagescc(w'); title('weight');
0128 clear y w r
0129 end
0130 end
0131
0132
0133 function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167 nt_greetings;
0168
0169
0170 if nargin<2; error('!'); end
0171 if nargin<3; w=[]; end
0172 if nargin<4||isempty(basis); basis='polynomials'; end
0173 if nargin<5||isempty(thresh); thresh=3; end
0174 if nargin<6||isempty(niter); niter=3; end
0175
0176 if thresh==0; error('thresh=0 is not what you want...'); end
0177
0178 dims=size(x);
0179 x=x(:,:);
0180 w=w(:,:);
0181
0182 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0183
0184
0185 if isnumeric(basis)
0186 r=basis;
0187 else
0188 switch basis
0189 case 'polynomials'
0190 r=zeros(size(x,1),numel(order));
0191 lin=linspace(-1,1,size(x,1));
0192 for k=1:order
0193 r(:,k)=lin.^k;
0194 end
0195 case 'sinusoids'
0196 r=zeros(size(x,1),numel(order)*2);
0197 lin=linspace(-1,1,size(x,1));
0198 for k=1:order
0199 r(:,2*k-1)=sin(2*pi*k*lin/2);
0200 r(:,2*k)=cos(2*pi*k*lin/2);
0201 end
0202 otherwise
0203 error('!');
0204 end
0205 end
0206
0207
0208
0209
0210
0211
0212 for iIter=1:niter
0213
0214
0215 [~,y]=nt_regw(x,r,w);
0216
0217
0218 d=x-y;
0219 if ~isempty(w); d=bsxfun(@times,d,w); end
0220 ww=ones(size(x));
0221 ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0222
0223
0224 if isempty(w);
0225 w=ww;
0226 else
0227 w=min(w,ww);
0228 end
0229 clear ww;
0230 end
0231 y=x-y;
0232 y=reshape(y,dims);
0233 w=reshape(w,dims);
0234
0235 end
0236
0237
0238 function test_me
0239 if 0
0240
0241 x=(1:100)'; x=x+ randn(size(x));
0242 WSIZE=30;
0243 y=nt_detrend2(x,1,[],[],[],[],WSIZE);
0244 figure(1); clf; plot([x,y]);
0245 end
0246 if 0
0247
0248 x=(1:100)'; x=x+ randn(size(x));
0249 x(40:50)=0;
0250 WSIZE=30;
0251 [yy,ww]=nt_detrend2(x,1,[],[],2,[],WSIZE);
0252 [y,w]=nt_detrend(x,1);
0253 figure(1); clf; subplot 211;
0254 plot([x,y,yy]);
0255 subplot 212; plot([w,ww],'.-');
0256 end
0257 if 0
0258
0259 x=cumsum(randn(1000,1)+0.1);
0260 WSIZE=200;
0261 [y1,w1]=nt_detrend(x,1,[]);
0262 [y2,w2]=nt_detrend2(x,1,[],[],[],[],WSIZE);
0263 figure(1); clf;
0264 plot([x,y1,y2]); legend('before', 'after');
0265 end
0266 if 0
0267
0268 x=linspace(0,100,1000)';
0269 x=x+3*randn(size(x));
0270 x(1:100,:)=100;
0271 w=ones(size(x)); w(1:100,:)=0;
0272 y=nt_detrend2(x,3,[],[],[],[],WSIZE);
0273 yy=nt_detrend2(x,3,w,[],[],[],WSIZE);
0274 figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0275 end
0276 if 0
0277 [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_43.mat'); x=x'; x=x(:,1:128);
0278
0279
0280 x=nt_demean(x);
0281 figure(1);
0282 nt_detrend(x,3);
0283 figure(2);
0284 WSIZE=1000;
0285 nt_detrend2(x(:,:),3,[],[],[],[],WSIZE);
0286 end
0287 end
0288
0289