0001 function [x,d]=nt_sns3(x,thresh)
0002
0003
0004
0005
0006
0007
0008
0009
0010 nt_greetings;
0011
0012 if nargin<1; error; end
0013 if nargin<2 || isempty(thresh); thresh=1; end
0014
0015 [nsample,nchan,~]=size(x);
0016 x=nt_unfold(x);
0017
0018 mn=mean(x);
0019 x=nt_demean(x);
0020 nn=sqrt(mean(x.^2));
0021 x=nt_normcol(x);
0022
0023
0024 For each channel, find sections for which it does not fit the
0025 subspace spanned by other sensors. The calculation is repeated
0026 and the projection matrix is refined at each step.
0027
0028 PCA_THRESH=10^-15;
0029
0030 d=zeros(size(x));
0031
0032
0033
0034 x=nt_demean(x);
0035 c0=nt_cov(x);
0036
0037 for iChan=1:nchan
0038
0039 oChan=setdiff(1:nchan,iChan);
0040
0041
0042 [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan));
0043 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0044 topcs=topcs(:,idx);
0045
0046
0047 b=c0(iChan,oChan)*topcs/(topcs'*c0(oChan,oChan)*topcs);
0048 y(:,iChan)=(x(:,oChan)*topcs)*b';
0049 dd=y(:,iChan)-x(:,iChan);
0050
0051 d(:,iChan)=mahal(dd,dd)/thresh;
0052
0053 end
0054
0055
0056
0057 Find which channel is most excentric at each moment.
0058
0059 w=max(d,[],2)<1;
0060 [~,iWorstChan]=sort(d','descend');
0061 iWorstChan=iWorstChan';
0062 iWorstChan(find(w),:)=nan;
0063
0064 x=nt_demean(x,w);
0065 c0=nt_cov(x(find(w),:));
0066 for iChan=1:nchan
0067
0068 idx=find(iWorstChan(:,1)==iChan);
0069
0070 ww=zeros(size(w));
0071 ww(idx)=1;
0072
0073 if numel(idx)>0
0074 c1=nt_cov(x(idx,:));
0075 todss=nt_dss0(c0,c1);
0076 z=nt_mmat(x,todss);
0077 fromdss=pinv(todss);
0078 x(idx,iChan)=nt_mmat(z(idx,2:end),fromdss(2:end,iChan));
0079 end
0080
0081
0082 end
0083
0084
0085 x=bsxfun(@times,x,nn);
0086 x=bsxfun(@plus,x,mn);
0087
0088 x=nt_fold(x,nsample);
0089 d=nt_fold(d,nsample);
0090