nt_dss1

[todss,pwr0,pwr1]=nt_dss1(x,w,keep1,keep2) - evoked-biased DSS denoising

  todss: denoising matrix
  pwr0: power per component (raw)
  pwr1: power per component (averaged)
  todssX: denoising matrix - crossvalidated
  dX: distance of each trial from mean of others - crossvalidated

  x: data to denoise (time * channels * trials)
  w: weight
  keep1: (in DSS0) number of PCs to retain (default: all)
  keep2: (in DSS0) ignore PCs smaller than keep2 (default: 10.^-12)

  The data mean is NOT removed prior to processing.

 NoiseTools

0001 function [todss,pwr0,pwr1,todssX,dX,d]=nt_dss1(x,w,keep1,keep2)
0002 %[todss,pwr0,pwr1]=nt_dss1(x,w,keep1,keep2) - evoked-biased DSS denoising
0003 %
0004 %  todss: denoising matrix
0005 %  pwr0: power per component (raw)
0006 %  pwr1: power per component (averaged)
0007 %  todssX: denoising matrix - crossvalidated
0008 %  dX: distance of each trial from mean of others - crossvalidated
0009 %
0010 %  x: data to denoise (time * channels * trials)
0011 %  w: weight
0012 %  keep1: (in DSS0) number of PCs to retain (default: all)
0013 %  keep2: (in DSS0) ignore PCs smaller than keep2 (default: 10.^-12)
0014 %
0015 %  The data mean is NOT removed prior to processing.
0016 %
0017 % NoiseTools
0018 
0019 if nargin<4; keep2=10.^-12; end
0020 if nargin<3; keep1=[]; end
0021 if nargin<2; w=[]; end
0022 if nargin<1; error('!'); end
0023 
0024 if ndims(x)<3; error('x should be 3D'); end
0025 if ~isa(x,'double'); warning('x is not double precision'); end
0026 
0027 x=x(:,:,:); % collapse higher dims
0028 
0029 [m,n,o]=size(x);
0030 
0031 if isempty(w)   % average over trials (--> bias function for DSS)
0032     [c0,nc0]=nt_cov(x);
0033     c0=c0/nc0;
0034     [c1,nc1]=nt_cov(mean(x,3)); 
0035     c1=c1/nc1;
0036 else
0037     % weighted average over trials (--> bias function for DSS)
0038     xx=nt_wmean(x,w,3);
0039     ww=min(w,[],2);
0040     % covariance of raw and biased data
0041     [c0,nc0]=nt_cov(x,[],w);
0042     c0=c0/nc0;
0043     [c1,nc1]=nt_cov(xx,[],ww); 
0044     c1=c1/nc1;
0045 end
0046 
0047 % derive DSS matrix
0048 [todss,pwr0,pwr1]=nt_dss0(c0,c1,keep1,keep2);
0049 
0050 % crossvalidated
0051 if nargout>3
0052     if nargin>1; error('!'); end % not implemented with weights
0053     [nsamples,nchans,ntrials]=size(x);
0054     % covariance per trial
0055     CC=zeros(nchans,nchans,ntrials);
0056     for iTrial=1:ntrials
0057         CC(:,:,iTrial)=nt_cov(x(:,:,iTrial));
0058     end
0059     % leave-one-out DSS
0060     todssX=zeros(nchans,nchans,ntrials);
0061     dX=zeros(nchans,ntrials);
0062     for iTrial=1:ntrials
0063         others=setdiff(1:ntrials,iTrial);
0064         C0=sum(CC(:,:,others),3)/nsamples/(ntrials-1);
0065         C1=nt_cov(mean(x(:,:,others),3))/nsamples;
0066         todssX(:,:,iTrial)=nt_dss0(C0,C1);
0067         z=nt_mmat(x,todssX(:,:,iTrial));
0068         tmp=mean(z(:,:,others),3);
0069         dX(:,iTrial)=mean( (nt_normcol(z(:,:,iTrial))-nt_normcol(tmp)).^2 );
0070     end
0071     d=zeros(nchans,ntrials);
0072     z=nt_mmat(x,todss);
0073     d=squeeze(mean(bsxfun(@minus,nt_normcol(z),nt_normcol(mean(z,3))).^2));
0074 end
0075 
0076 % if no output arguments, just plot
0077 if nargout==0
0078     figure(100); clf; 
0079     subplot 221; 
0080     plot(pwr1./pwr0,'.-');
0081     xlabel('component'); ylabel('score');
0082     z=nt_mmat(x,todss(:,1:3));
0083     for iComp=1:3
0084         subplot(2,2,iComp+1);
0085         nt_bsplot(z(:,iComp,:));
0086         title(iComp);
0087         xlabel('sample');
0088     end
0089     
0090     clear todss pwr0 pwr1
0091 end
0092