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
0029 y=zeros(size(x));
0030 d=ones(size(x,1),nchan);
0031
0032 NITER=2;
0033 for k=1:NITER
0034
0035
0036 dd=max(d,[],2);
0037 c0=nt_cov(x,[],dd<=1);
0038
0039 d0=d;
0040 for iChan=1:nchan
0041
0042 oChan=setdiff(1:nchan,iChan);
0043
0044 [topcs,eigenvalues]=nt_pcarot(c0(oChan,oChan));
0045 PCA_THRESH=10^-6;
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
0051 y0=y;
0052
0053
0054 competitors=zeros(size(d0(:,oChan)));
0055 iCompetitors=d0(:,oChan)> 1 ...
0056 & d0(:,oChan) > repmat(d0(:,iChan),1,nchan-1);
0057 competitors(find(iCompetitors))=d0(find(iCompetitors));
0058 idx=find(~min(competitors,[],2));
0059 [~,worst]=max(competitors,[],2);
0060 worst(idx)=0;
0061 for k=1:numel(oChan);
0062
0063 iReplace=find(k==worst);
0064 if iReplace
0065 ooChan=setdiff(1:nchan,[iChan,oChan(k)]);
0066 [topcs,eigenvalues]=nt_pcarot(c0(ooChan,ooChan));
0067 PCA_THRESH=10^-6;
0068 idx=find(eigenvalues/max(eigenvalues) > PCA_THRESH);
0069 topcs=topcs(:,idx);
0070 b=c0(iChan,ooChan)*topcs/(topcs'*c0(ooChan,ooChan)*topcs);
0071 y(iReplace,iChan)=(x(iReplace,ooChan)*topcs)*b';
0072 end
0073 end
0074
0075
0076 ddd=y(:,iChan)-x(:,iChan);
0077
0078 d(:,iChan)=mahal(ddd,ddd)/thresh;
0079
0080 figure(1); clf
0081 FOCUS=861600+(1:500); subplot 311; plot(x(FOCUS,:)); title('raw');
0082 subplot 312; plot([x(FOCUS,iChan), y0(FOCUS,iChan), y(FOCUS,iChan)]); title (num2str(iChan)); legend('raw','proj','fixed');
0083 subplot 313; plot(d(FOCUS,:)); pause
0084
0085 end
0086
0087 end
0088
0089 dd=max(d,[],2);
0090
0091 c0=nt_cov(x,[],dd<1);
0092
0093
0094
0095 For each channel, find the part for which it dominates the artifact.
0096 Replace that part based on projection on the other channels.
0097
0098
0099 idx=find(dd>1);
0100 xx=x(idx,:);
0101 xx=abs(nt_normcol(xx));
0102
0103 for iChan=1:nchan
0104
0105
0106
0107
0108
0109
0110
0111 oChan=setdiff(1:nchan,iChan);
0112 b=c0(iChan,oChan)/c0(oChan,oChan);
0113 idx2=find(xx(:,iChan)>max(xx(:,oChan),[],2));
0114
0115 x(idx(idx2),iChan)=x(idx(idx2),oChan)*b';
0116
0117 idx2=find(d(:,iChan)>1 & d(:,iChan)>max(d(:,oChan),[],2));
0118
0119
0120 end
0121
0122
0123 To do:
0124 - avoid problems with rank-deficient data
0125 - exclude borderline channels from projection
0126
0127
0128 x=bsxfun(@times,x,nn);
0129 x=bsxfun(@plus,x,mn);
0130
0131 x=nt_fold(x,nsample);
0132 d=nt_fold(d,nsample);
0133
0134
0135