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) && ~(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
0100 if stop<size(x,1); y(end,:)=y(end-1,:); a(end,:)=a(end-1,:); end;
0101
0102 y=bsxfun(@times,y,1./a); y(find(isnan(y)))=0;
0103 trend=bsxfun(@times,trend,1./a); trend(find(isnan(trend)))=0;
0104
0105
0106 d=x-trend;
0107
0108
0109 if ~isempty(w); d=bsxfun(@times,d,w); end
0110 ww=ones(size(x));
0111 ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0112 clear d
0113
0114
0115 if isempty(w);
0116 w=ww;
0117 else
0118 w=min(w,ww);
0119 end
0120 clear ww;
0121
0122 end
0123 end
0124
0125 if ~nargout
0126
0127 subplot 411; plot(x); title('raw'); xlim([1,size(x,1)])
0128 subplot 412; plot(y); title('detrended'); xlim([1,size(x,1)])
0129 subplot 413; plot(x-y); title('trend'); xlim([1,size(x,1)])
0130 subplot 414; nt_imagescc(w'); title('weight');
0131 clear y w r
0132 end
0133 end
0134
0135
0136 function [y,w,r]=nt_detrend_helper(x,order,w,basis,thresh,niter)
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
0168
0169
0170 nt_greetings;
0171
0172
0173 if nargin<2; error('!'); end
0174 if nargin<3; w=[]; end
0175 if nargin<4||isempty(basis); basis='polynomials'; end
0176 if nargin<5||isempty(thresh); thresh=3; end
0177 if nargin<6||isempty(niter); niter=3; end
0178
0179 if thresh==0; error('thresh=0 is not what you want...'); end
0180
0181 dims=size(x);
0182 x=x(:,:);
0183 w=w(:,:);
0184
0185 if size(w,2)==1; w=repmat(w,1,size(x,2)); end
0186
0187
0188 if isnumeric(basis)
0189 r=basis;
0190 else
0191 switch basis
0192 case 'polynomials'
0193 r=zeros(size(x,1),numel(order));
0194 lin=linspace(-1,1,size(x,1));
0195 for k=1:order
0196 r(:,k)=lin.^k;
0197 end
0198 case 'sinusoids'
0199 r=zeros(size(x,1),numel(order)*2);
0200 lin=linspace(-1,1,size(x,1));
0201 for k=1:order
0202 r(:,2*k-1)=sin(2*pi*k*lin/2);
0203 r(:,2*k)=cos(2*pi*k*lin/2);
0204 end
0205 otherwise
0206 error('!');
0207 end
0208 end
0209
0210
0211
0212
0213
0214
0215 for iIter=1:niter
0216
0217
0218 [~,y]=nt_regw(x,r,w);
0219
0220
0221 d=x-y;
0222 if ~isempty(w); d=bsxfun(@times,d,w); end
0223 ww=ones(size(x));
0224 ww(find(abs(d)>thresh*repmat(std(d),size(x,1),1))) = 0;
0225
0226
0227 if isempty(w);
0228 w=ww;
0229 else
0230 w=min(w,ww);
0231 end
0232 clear ww;
0233 end
0234 y=x-y;
0235 y=reshape(y,dims);
0236 w=reshape(w,dims);
0237
0238 end
0239
0240
0241 function test_me
0242 if 0
0243
0244 x=(1:100)'; x=x+ randn(size(x));
0245 WSIZE=30;
0246 y=nt_detrend2(x,1,[],[],[],[],WSIZE);
0247 figure(1); clf; plot([x,y]);
0248 end
0249 if 0
0250
0251 x=(1:100)'; x=x+ randn(size(x));
0252 x(40:50)=0;
0253 WSIZE=30;
0254 [yy,ww]=nt_detrend2(x,1,[],[],2,[],WSIZE);
0255 [y,w]=nt_detrend(x,1);
0256 figure(1); clf; subplot 211;
0257 plot([x,y,yy]);
0258 subplot 212; plot([w,ww],'.-');
0259 end
0260 if 0
0261
0262 x=cumsum(randn(1000,1)+0.1);
0263 WSIZE=200;
0264 [y1,w1]=nt_detrend(x,1,[]);
0265 [y2,w2]=nt_detrend2(x,1,[],[],[],[],WSIZE);
0266 figure(1); clf;
0267 plot([x,y1,y2]); legend('before', 'after');
0268 end
0269 if 0
0270
0271 x=linspace(0,100,1000)';
0272 x=x+3*randn(size(x));
0273 x(1:100,:)=100;
0274 w=ones(size(x)); w(1:100,:)=0;
0275 y=nt_detrend2(x,3,[],[],[],[],WSIZE);
0276 yy=nt_detrend2(x,3,w,[],[],[],WSIZE);
0277 figure(1); clf; plot([x,y,yy]); legend('before', 'unweighted','weighted');
0278 end
0279 if 0
0280 [p,x]=nt_read_data('/data/meg/theoldmanandthesea/eeg/mc/MC_aespa_speech_43.mat'); x=x'; x=x(:,1:128);
0281
0282
0283 x=nt_demean(x);
0284 figure(1);
0285 nt_detrend(x,3);
0286 figure(2);
0287 WSIZE=1000;
0288 nt_detrend2(x(:,:),3,[],[],[],[],WSIZE);
0289 end
0290 end
0291
0292