0001 function [x,w,ww]=nt_star3(x,thresh,closest,w)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015 26 Jan 2017
0016 Regress also on polynomial to handle trends.
0017
0018
0019
0020
0021 nt_greetings;
0022
0023 PCA_THRESH=10^-15;
0024 NSMOOTH=10;
0025 MINPROP=0.3;
0026 NITER=3;
0027 VERBOSE=1;
0028
0029 if nargin<1; error; end
0030 if nargin<2 || isempty(thresh); thresh=1; end
0031 if nargin<3; closest=[]; end
0032 if ~isempty(closest)&&size(closest,1)~=size(x,2);
0033 error('closest array should have as many rows as channels of x');
0034 end
0035 if nargin<4 || isempty(w); w=1; end
0036
0037 if nargout==0
0038 [y,w,ww]=nt_star3(x,thresh,closest,w);
0039 disp([mean(w(:)), mean(ww(:))])
0040 figure(1); clf;
0041 subplot 311; plot(x);
0042 subplot 312; plot(y);
0043 subplot 313; plot(w, '.-'); ylim([0 1.1]);
0044 clear x w ww
0045 return
0046 end
0047
0048
0049 [nsample,nchan,~]=size(x);
0050 x=nt_unfold(x);
0051
0052 p0=nt_wpwr(x);
0053 mn=mean(x);
0054 x=nt_demean(x);
0055 nn=sqrt(mean(x.^2));
0056 x=nt_normcol(x);
0057 p00=nt_wpwr(x);
0058
0059
0060 iNan=find(all(isnan(x)));
0061 iZero=find(all(x==0));
0062 x(:,iNan)=randn(size(x,1),numel(iNan));
0063 x(:,iZero)=randn(size(x,1),numel(iZero));
0064
0065 order=1;
0066 r=zeros(size(x,1),numel(order));
0067 lin=linspace(-1,1,size(x,1));
0068 for k=1:order
0069 r(:,k)=lin.^k;
0070 end
0071
0072 x=nt_demean(x);
0073 c0=nt_cov(x);
0074
0075
0076 Find time intervals where at least one channel is excentric --> w==0.
0077
0078
0079 iIter=NITER;
0080 while iIter>0
0081
0082
0083 w=ones(size(x,1),1);
0084 for iChan=1:nchan
0085
0086
0087 if ~isempty(closest);
0088 oChan=closest(iChan,:);
0089 else
0090 oChan=setdiff(1:nchan,iChan);
0091 end
0092 oChan(oChan>nchan)=[];
0093
0094
0095 [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan));
0096 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0097 topcs=topcs(:,idx);
0098
0099
0100 b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs);
0101 y=x(:,oChan)*(topcs*b');
0102 dx=abs(y-x(:,iChan));
0103 dx=dx+eps;
0104
0105 d=dx/mean(dx(find(w)));
0106 if NSMOOTH>0;
0107 d=filtfilt(ones(NSMOOTH,1)/NSMOOTH,1,d);
0108 end
0109
0110 d=d/thresh;
0111 w=min(w,(d<1));
0112
0113 end
0114
0115 prop=mean(w);
0116 if VERBOSE>0; disp(['proportion artifact free: ', num2str(prop)]); end
0117
0118 if iIter==NITER && prop<MINPROP
0119 thresh=thresh*1.1;
0120 if VERBOSE>0; disp(['Warning: nt_star increasing threshold to ', num2str(thresh)]); end
0121 w=ones(size(w));
0122 else
0123 iIter=iIter-1;
0124 end
0125
0126 x=nt_demean(x,w);
0127 c0=nt_cov(x,[],w);
0128 end
0129
0130
0131 We now know which part contains channel-specific artifacts (w==0 for artifact part),
0132 and we have an estimate of the covariance matrix of the artifact-free part.
0133
0134
0135
0136 Find which channel is most excentric at each time point.
0137 Here we use an excentricity measure based on the absolute value of the signal,
0138 rather than the difference between signal and projection.
0139
0140
0141 xx=abs(x);
0142 xx=bsxfun(@times,xx, 1 ./ sqrt(mean(xx(find(w),:).^2)));
0143 if NSMOOTH>0;
0144 xx=filtfilt(ones(NSMOOTH,1)/NSMOOTH,1,xx);
0145 end
0146 [~,rank]=sort(xx','descend');
0147 rank=rank';
0148 rank(find(w),:)=0;
0149
0150 depth=1;
0151 depth=min(depth,nchan-1);
0152 ww=ones(size(x));
0153 for iDepth=1:depth
0154
0155
0156 Fix each channel by projecting on other channels.
0157
0158
0159 iFixed=nchan;
0160 nFixed=0;
0161 for iChan=1:nchan
0162
0163 bad_samples=find(iChan==rank(:,iDepth));
0164 if iDepth ~=1;
0165 bad_samples(find(xx(bad_samples,iChan)<thresh)) =[];
0166 end
0167
0168 nFixed=nFixed+numel(bad_samples);
0169 if isempty(bad_samples);
0170 iFixed=iFixed-1;
0171 continue;
0172 end
0173 ww(bad_samples,iChan)=0;
0174
0175
0176 if ~isempty(closest);
0177 oChan=closest(iChan,:);
0178 else
0179 oChan=setdiff(1:nchan,iChan);
0180 end
0181 oChan(oChan>nchan)=[];
0182
0183
0184 [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan));
0185 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0186 topcs=topcs(:,idx);
0187
0188
0189 b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs);
0190 y=x(bad_samples,oChan)*(topcs*b');
0191
0192 x(bad_samples,iChan)=y(:);
0193
0194 end
0195
0196 if VERBOSE>0;
0197 disp(['depth: ', num2str(iDepth), ', n fixed channels: ',num2str(iFixed),...
0198 ', n fixed samples: ', num2str(nFixed), ', ratio: ',num2str(nt_wpwr(x)/p00)]);
0199 end
0200 end
0201
0202 x=nt_demean(x);
0203 x=bsxfun(@times,x,nn);
0204 x=bsxfun(@plus,x,mn);
0205
0206 x=nt_fold(x,nsample);
0207 w=nt_fold(w,nsample);
0208 ww=nt_fold(ww,nsample);
0209
0210
0211
0212
0213 x(:,iNan)=nan;
0214 x(:,iZero)=0;
0215
0216 if VERBOSE>0; disp(['power ratio: ', num2str(nt_wpwr(x)/p0)]); end