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