0001 function [x,w]=nt_sns2(x,thresh)
0002
0003
0004
0005
0006
0007
0008
0009
0010 nt_greetings;
0011 if nargin<1; error; end
0012 if nargin<2 || isempty(thresh); thresh=1; end
0013 PCA_THRESH=10^-15;
0014 [nsample,nchan,~]=size(x);
0015
0016 x=nt_unfold(x);
0017 mn=mean(x);
0018 x=nt_demean(x);
0019 nn=sqrt(mean(x.^2));
0020 x=nt_normcol(x);
0021
0022
0023 For each channel, find sections for which it does not fit the
0024 subspace spanned by other sensors. The calculation is repeated
0025 and the projection matrix is refined at each step.
0026
0027
0028 w=ones(size(x));
0029
0030 NITER=3;
0031 for k=1:NITER
0032
0033
0034 ww=min(w,[],2);
0035 x=nt_demean(x,ww);
0036 [c0,tw]=nt_cov2(x,ww);
0037 c0=c0./tw; c0(isnan(c0))=0;
0038
0039 for iChan=1:nchan
0040
0041
0042 oChan=setdiff(1:nchan,iChan);
0043
0044 [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan));
0045 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0046 topcs=topcs(:,idx);
0047 b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs);
0048 y(:,iChan)=(x(:,oChan)*topcs)*b';
0049
0050
0051 for iChan2=oChan
0052 oChan2=setdiff(oChan,iChan2);
0053 iBad=find(w(:,iChan2)==0);
0054 if ~isempty(iBad)
0055 [topcs,eigenvalues]=nt_pcarot(c0(oChan2,oChan2));
0056 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0057 topcs=topcs(:,idx);
0058 b=c0(iChan,oChan2)*topcs/(topcs'*c0(oChan2,oChan2)*topcs);
0059 y(iBad,iChan)=(x(iBad,oChan2)*topcs)*b';
0060 end
0061 end
0062
0063
0064
0065 dd=y(:,iChan)-x(:,iChan);
0066
0067 d=mahal(dd,dd)/thresh;
0068
0069 w(:,iChan)=(d<1);
0070
0071 end
0072
0073 disp(mean(w(:)))
0074 end
0075
0076 x=y;
0077
0078
0079 To do:
0080 Record all the DC shifts introduced when x mean is removed, so as to
0081 restore accurately.
0082
0083
0084 x=nt_demean(x);
0085 x=bsxfun(@times,x,nn);
0086 x=bsxfun(@plus,x,mn);
0087
0088 x=nt_fold(x,nsample);
0089 w=nt_fold(w,nsample);
0090
0091
0092