nt_subspace_prune6

[Y]=nt_subspace_prune6(X,npass,thresh) - local cleaning matrices

  Y: denoised data
  
  X: data to denoise (nsamples X nchans X ntrials matrix or array of 2D matrices)
  npass: number of passes [default: 10]
  thresh: threshold power ratio between segment & all [default: 10]
  tol: tolerance factor to speed up calculation [default: 0.5]

0001 function [Y,scores,removed]=nt_subspace_prune6(X,npass,thresh,tol)
0002 %[Y]=nt_subspace_prune6(X,npass,thresh) - local cleaning matrices
0003 %
0004 %  Y: denoised data
0005 %
0006 %  X: data to denoise (nsamples X nchans X ntrials matrix or array of 2D matrices)
0007 %  npass: number of passes [default: 10]
0008 %  thresh: threshold power ratio between segment & all [default: 10]
0009 %  tol: tolerance factor to speed up calculation [default: 0.5]
0010 
0011 if nargin<2||isempty(npass); npass=10; end
0012 if nargin<3||isempty(thresh); thresh=10; end
0013 if nargin<4||isempty(tol); tol=0.5; end
0014 
0015 if isnumeric(X)
0016     % transfer 3D matrix to array of 2D
0017     if ndims(X)<3; error('!'); end
0018     tmp={};
0019     for iTrial=1:size(X,3); 
0020         tmp{iTrial}=X(:,:,iTrial); 
0021     end
0022     X=tmp;
0023     
0024     % process
0025     [Y,scores,removed]=nt_subspace_prune6(X,npass,thresh,tol);
0026         disp(nt_wpwr(Y)/nt_wpwr(X));
0027     
0028     % transfer back to 3D matrix
0029     tmp=zeros(size(Y{1},1),size(Y{1},2),numel(Y));
0030     tmp2=zeros(size(Y{1},1),size(removed{1},2),numel(Y));
0031     for iTrial=1:numel(X) 
0032         tmp(:,:,iTrial)=Y{iTrial}; 
0033         tmp2(:,:,iTrial)=removed{iTrial}; 
0034     end
0035     Y=tmp;
0036     removed=tmp2;
0037     return
0038 end
0039 
0040 ntrials=numel(X);
0041 nchans=size(X{1},2);
0042 [C00,tw]=nt_cov(X);
0043 C00=C00/tw;
0044 
0045 % matrix array to save removed component/trials
0046 removed={};
0047 for iTrial=1:ntrials
0048     removed{iTrial}=zeros(size(X{iTrial},1),npass);
0049 end
0050 
0051 original_power=nt_wpwr(X);
0052 scores=[]; D=[];
0053 for iPass=1:npass
0054     
0055     X=nt_demean2(X);
0056     
0057 %     for k=1:ntrials
0058 %         X{k}=X{k}/sqrt(mean(X{k}(:).^2));
0059 %     end
0060     
0061     % covariance of full data
0062     [C0,tw]=nt_cov(X); 
0063     C0=C0/tw;
0064     
0065     % mix with original estimate
0066     alpha=0.01;
0067     C0=alpha*C00+(1-alpha)*C0;
0068     
0069     % find most excentric trial
0070     iBest=0; best_score=0; 
0071     CC1=zeros(nchans,nchans,ntrials);
0072     for iTrial=1:numel(X)
0073         a=X{iTrial};
0074         % covariance of this trial
0075         C1=nt_cov(a)/size(a,1);       
0076         % contrast this trial with rest
0077         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0078         % is this trial the most excentric?
0079         if pwr1(1)/pwr0(1)>best_score
0080             iBest=iTrial;
0081             best_score=pwr1(1)/pwr0(1);
0082         end
0083         scores(iPass,iTrial)=pwr1(1)/pwr0(1);
0084         if pwr1(1)<pwr0(1);
0085             figure(1); clf; plot([pwr1', pwr0']); pause
0086         end
0087     end
0088     
0089     % remove most excentric component of most excentric trials
0090     if best_score>thresh
0091         
0092         % find other trials for which this component is large
0093         a=X(iBest);
0094         C1=nt_cov(a)/(size(a,1)*size(a,3));
0095         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0096         z=nt_mmat(X,todss(:,1));
0097         for k=1:ntrials
0098             p(k)=mean(z{k}.^2);
0099         end
0100         p=p/mean(p);
0101         iRemove=find(p>1/tol);
0102         
0103         %disp(numel(iRemove))
0104         
0105         % update DSS to fit all trials to be removed
0106         a=X(iRemove);
0107         C1=nt_cov(a)/(size(a,1)*size(a,3));
0108         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0109         fromdss=pinv(todss);
0110         
0111         D=todss(:,2:end)*fromdss(2:end,:);
0112         X0=X;
0113         X(iRemove)=nt_mmat(X(iRemove),D);
0114                 
0115         for k=iRemove
0116             removed{k}(:,iPass)=nt_mmat(X{k},todss(:,1));
0117         end
0118     else
0119         break;
0120     end
0121 
0122     if ~isreal(scores); return; end
0123     
0124     figure(10); clf; nt_imagescc(scores); colorbar;
0125     %disp(nt_wpwr(X)/original_power);
0126 end
0127 
0128 for iTrial=1:ntrials
0129     removed{iTrial}=removed{iTrial}(:,1:iPass);
0130 end
0131     
0132 Y=X;