Description of nt

0001 function [Y,excentricity,removed,cov1,cov2]=nt_LSP(X,npass,thresh,tol,guard)
0002 %[Y,scores,removed]=nt_LSP(X,npass,thresh,tol,guard) - local subspace pruning
0003 %
0004 %  Y: denoised data
0005 %  scores: record of excentricity scores for each trial and each pass
0006 %  removed: components removed
0007 %
0008 %  X: data to denoise (nsamples X nchans X ntrials matrix or array of 2D matrices)
0009 %  npass: number of passes [default: 10]
0010 %  thresh: threshold excentricity score [default: 10]
0011 %  tol: tolerance factor to speed up calculation [default: 0.5]
0012 %  guard: don't modify channels with correlation below this limit [default: 0.1]
0013 %
0014 % For each trial, JD is used to contrast it with all other trials.  If the
0015 % power ratio ('score') of the first component is above threshold, that
0016 % component is discarded from that trial.
0017 %
0018 % NoiseTools.
0019 
0020 if nargin<2||isempty(npass); npass=10; end
0021 if nargin<3||isempty(thresh); thresh=10; end
0022 if nargin<4||isempty(tol); tol=0.5; end
0023 if nargin<5||isempty(guard); guard=0.1; end
0024 
0025 if isnumeric(X)
0026     % transfer 3D matrix to array of 2D
0027     if ndims(X)<3; error('!'); end
0028     tmp={};
0029     for iTrial=1:size(X,3); 
0030         tmp{iTrial}=X(:,:,iTrial); 
0031     end
0032     X=tmp;
0033     
0034     % process
0035     [Y,excentricity,removed,cov1,cov2]=nt_lsp(X,npass,thresh,tol);
0036     
0037     % transfer back to 3D matrix
0038     tmp=zeros(size(Y{1},1),size(Y{1},2),numel(Y));
0039     tmp2=zeros(size(Y{1},1),size(removed{1},2),numel(Y));
0040     for iTrial=1:numel(X) 
0041         tmp(:,:,iTrial)=Y{iTrial}; 
0042         tmp2(:,:,iTrial)=removed{iTrial}; 
0043     end
0044     Y=tmp;
0045     removed=tmp2;
0046     return
0047 end
0048 
0049 ntrials=numel(X);
0050 nchans=size(X{1},2);
0051 [C00,tw]=nt_cov(X);
0052 C00=C00/tw;
0053 
0054 % create matrix array to save removed component/trials
0055 removed={};
0056 for iTrial=1:ntrials
0057     removed{iTrial}=zeros(size(X{iTrial},1),npass);
0058 end
0059 
0060 excentricity=[]; D=[]; score_fig=figure(10);
0061 cov1=zeros(size(X{1},2),size(X{1},2),npass+1);
0062 cov2=zeros(size(X{1},2),size(X{1},2),npass+1);
0063 cov1(:,:,1)=nt_cov(X); 
0064 %disp(size(X{1}));
0065 %disp(size(nt_trial2mat(X)));
0066 cov2(:,:,1)=nt_cov(mean(nt_trial2mat(X),3));
0067 for iPass=1:npass
0068     
0069     X=nt_demean2(X);
0070     
0071 %     for k=1:ntrials
0072 %         X{k}=X{k}/sqrt(mean(X{k}(:).^2));
0073 %     end
0074     
0075     % covariance of full data
0076     [C0,tw]=nt_cov(X); 
0077     %[C0,tw]=nt_cov_smr(X);
0078     C0=C0/tw;
0079     
0080     
0081     % mix with original estimate
0082     alpha=0.01;
0083     C0=alpha*C00+(1-alpha)*C0;
0084     
0085     % find most excentric trial
0086     iBest=0; best_score=0; 
0087     CC1=zeros(nchans,nchans,ntrials);
0088     for iTrial=1:numel(X)
0089         a=X{iTrial};
0090         % covariance of this trial
0091         C1=nt_cov(a)/size(a,1);       
0092         % contrast this trial with rest
0093         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0094         % is this trial the most excentric?
0095         if pwr1(1)/pwr0(1)>best_score
0096             iBest=iTrial;
0097             best_score=pwr1(1)/pwr0(1);
0098         end
0099         excentricity(iPass,iTrial)=pwr1(1)/pwr0(1);
0100 %         if pwr1(1)<pwr0(1)
0101 %             figure(1); clf; plot([pwr1', pwr0']); title([pwr0(1), pwr1(1)]); drawnow; pause
0102 %         end
0103     end
0104     
0105     % remove most excentric component of most excentric trials
0106     if best_score>thresh
0107         
0108         % find other trials for which this component is large
0109         a=X(iBest);
0110         C1=nt_cov(a)/(size(a,1)*size(a,3));
0111         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0112         z=nt_mmat(X,todss(:,1));
0113         for k=1:ntrials
0114             p(k)=mean(z{k}.^2);
0115         end
0116         iRemove=find(p/max(p)>tol & p>thresh);
0117         if isempty(iRemove); break; end
0118         
0119         disp(numel(iRemove))
0120         
0121         % update DSS to fit all trials to be removed
0122         a=X(iRemove);
0123         C1=nt_cov(a)/(size(a,1)*size(a,3));
0124         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0125         fromdss=pinv(todss);
0126 
0127         
0128         D=todss(:,2:end)*fromdss(2:end,:); % denoising matrix
0129         X0=X;
0130         X(iRemove)=nt_mmat(X(iRemove),D);
0131         
0132         if 0
0133             a=nt_unfold(nt_cell2mat(X(iRemove)));
0134             b=nt_unfold(nt_cell2mat(X0(iRemove)));
0135             r=diag(corr(b,a-b));
0136             thresh_r=guard; % don't modify channels with correlation below this
0137             mask=abs(r)<thresh_r;
0138 
0139             for k=iRemove
0140                 X{k}(:,mask)=X0{k}(:,mask); % revert masked channels
0141             end
0142 
0143             for k=iRemove
0144                 removed{k}(:,iPass)=nt_mmat(X{k},todss(:,1));
0145             end
0146         end
0147     else
0148         break;
0149     end
0150 
0151     if ~isreal(excentricity); return; end
0152     
0153     if nt_verbose 
0154         set(0,'currentfigure',score_fig); clf; 
0155         imagesc(excentricity); 
0156         h=colorbar; set(get(h,'label'),'string','excentricity score');
0157         xlabel('trial'); ylabel('pass'); drawnow
0158     end
0159     cov1(:,:,iPass+1)=nt_cov(nt_cell2mat(X));
0160     cov2(:,:,iPass+1)=nt_cov(mean(nt_cell2mat(X),3));
0161 end
0162 
0163 for iTrial=1:ntrials
0164     removed{iTrial}=removed{iTrial}(:,1:iPass);
0165 end
0166     
0167 Y=X;
nt_lsp

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
