nt_subspace_prune

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

  Y: denoised data
  
  X: data to denoise (nsamples X nchans X ntrials)
  npass: number of passes [default: 1]
  thresh: threshold power ratio between segment & all [default: 10]

0001 function [Y,scores]=nt_subspace_prune(X,npass,thresh)
0002 %[Y]=nt_subspace_prune(X,npass,thresh) - local cleaning matrices
0003 %
0004 %  Y: denoised data
0005 %
0006 %  X: data to denoise (nsamples X nchans X ntrials)
0007 %  npass: number of passes [default: 1]
0008 %  thresh: threshold power ratio between segment & all [default: 10]
0009 
0010 if nargin<2||isempty(npass); npass=1; end
0011 if nargin<3||isempty(thresh); thresh=1; end
0012 
0013 if isnumeric(X)
0014     if ndims(X)<3; error('!'); end
0015     tmp={};
0016     for iTrial=1:size(X,3); 
0017         tmp{iTrial}=X(:,:,iTrial); 
0018     end
0019     X=tmp;
0020     [Y,scores]=nt_subspace_prune(X,npass,thresh);
0021     tmp=zeros(size(Y{1},1),size(Y{2},2),numel(Y));
0022     for iTrial=1:numel(X); 
0023         tmp(:,:,iTrial)=Y{iTrial}; 
0024     end
0025     Y=tmp;
0026     return
0027 end
0028 
0029 ntrials=numel(X);
0030 nchans=size(X{1},2);
0031 
0032 scores=zeros(ntrials,1,npass);
0033 original_power=nt_wpwr(X);
0034 for iPass=1:npass
0035     
0036     X=nt_demean2(X);
0037     
0038     % covariance matrices
0039     [C0,tw]=nt_cov(X); 
0040     C0=C0/tw;
0041     CC1=zeros(nchans,nchans,ntrials);
0042     for iTrial=1:numel(X)
0043         a=X{iTrial};
0044         %disp(size(a))
0045         CC1(:,:,iTrial)=nt_cov(a)/size(a,1);
0046     end
0047     
0048     % denoising matrices
0049     nochange=1; % to abort when no effect`
0050     MM=zeros(nchans,nchans,ntrials);
0051     for iTrial=1:ntrials
0052         C1=CC1(:,:,iTrial);
0053         [todss,pwr0,pwr1]=nt_dss0(C0,C1,[],0);
0054         fromdss=pinv(todss);
0055         
0056         if 0
0057             figure(9); plot(pwr1./pwr0, '.-'); drawnow; pause(1)
0058         end
0059         
0060         if pwr1(1)./pwr0(1) < thresh
0061             MM(:,:,iTrial)=eye(nchans);
0062         else
0063             nochange=0;
0064             MM(:,:,iTrial)=todss(:,2:end)*fromdss(2:end,:);
0065         end
0066         scores(iTrial,1,iPass)=pwr1(1)./pwr0(1);
0067     end
0068     if nochange; break; end
0069 
0070     %figure(10); clf; plot(scores(:,:)); pause
0071     
0072     % denoise
0073     Y=X;
0074     for iTrial=1:ntrials
0075         tmp=X{iTrial};
0076         X{iTrial}=X{iTrial}*MM(:,:,iTrial);
0077         
0078         if 0
0079             figure(10); clf
0080             subplot 121; plot(tmp)
0081             subplot 122; plot(X{iTrial});
0082         end
0083     end
0084 
0085     %disp(nt_wpwr(X)/original_power);
0086 end
0087     
0088 Y=X;