nt_sca

[M,y,score,proportion]=nt_sca(x,ncomp) - shared component analysis

 x: data (time X channels) (can be cell array)
 ncomp: keep only ncomp components (faster)

 M: SCA transform matrix
 z: transformed data
 score: sharedness score, per component
 prop: proportion of power accounted for, per component

0001 function [M,y,score,proportion]=nt_sca(x,ncomp)
0002 %[M,y,score,proportion]=nt_sca(x,ncomp) - shared component analysis
0003 %
0004 % x: data (time X channels) (can be cell array)
0005 % ncomp: keep only ncomp components (faster)
0006 %
0007 % M: SCA transform matrix
0008 % z: transformed data
0009 % score: sharedness score, per component
0010 % prop: proportion of power accounted for, per component
0011 
0012 
0013 if nargin<1; error('!'); end
0014 if nargin<2; ncomp=[]; end
0015 
0016 if iscell(x)
0017     xx=[];
0018     for iTrial=1:numel(x)
0019         xx=[xx;x{iTrial}];
0020     end
0021     x=xx; clear xx;
0022 end
0023 x=nt_demean(x);
0024 
0025 if isempty(ncomp); ncomp=size(x,2); end
0026 
0027 
0028 
0029 THRESH=10^-12;
0030 
0031 %{
0032 Todo: 
0033 - allow xvalidation
0034 - operate on covariance matrices
0035 %}
0036 
0037 T=eye(size(x,2)); % current transform
0038 M=eye(size(x,2)); % result
0039 C0=nt_cov(x); % initial covariance
0040 
0041 score=[];
0042 for iComp=1:ncomp
0043     C=T'*C0*T; % current covariance
0044     N=diag(1./sqrt(diag(C))); % normalizing matrix
0045     N(find(isnan(N)))=0;
0046     C=N'*C*N; % normalize current covariance
0047     [topcs,ev]=nt_pcarot(C); % PCA
0048     frompcs=pinv(topcs);
0049     M(:,iComp)=T*N*topcs(:,1); % keep first PC
0050     T=T*N*(topcs(:,2:end)*frompcs(2:end,:)); % project out first PC
0051     score(iComp)=ev(1);
0052 end
0053 
0054 if ncomp<size(x,2) 
0055     % fill rest of transform matrix with leftover (unprocessed) dimensions
0056     C=T'*C0*T; % current covariance
0057     N=diag(1./sqrt(diag(C))); % normalizing matrix
0058     N(find(isnan(N)))=0;
0059     C=N'*C*N; % normalize current covariance
0060     topcs=nt_pcarot(C); % PCA
0061     T=T*topcs;
0062     M(:,ncomp+1:end)=T(:,1:(size(x,2)-ncomp));
0063 end
0064 
0065 prop=diag(M'*C*M);
0066 
0067 %figure(1); clf; plot(prop); pause
0068 
0069 if nargout>1; y=x*M; end
0070 
0071 % test code
0072 if 0
0073     x=randn(1000,10);
0074     [M,y]=nt_sca(x);
0075 end
0076 
0077 if 0
0078     % data are 11 chan:
0079     % 10 chan share same source (sine),
0080     % 1 chan is different source (noise) with higher variance
0081     x=randn(1000,10);
0082     s=sin(2*pi*(1:1000)'/1000);
0083     x=bsxfun(@plus,x,s); % add same to all
0084     x=[x,10*randn(1000,1)]; % extra channel with large variance
0085     %[y,M,score]=nt_sca_old(x);
0086     MM=nt_sca(x); y=x*MM;
0087     yy=nt_pca(x);
0088     figure(1); clf;  plot(nt_normcol(s)'*nt_normcol(y)/size(s,1)); 
0089     hold on; plot(nt_normcol(s)'*nt_normcol(yy)/size(s,1));    legend('sca','pca');
0090 end
0091 
0092 
0093 if 0
0094     % two shared sources
0095     x=randn(1000,10);
0096     s=sin(2*pi*(1:1000)'/1000);
0097     s2=sin(2*pi*2*(1:1000)'/1000);
0098     x=x+s*rand(1,10); % add same to all
0099     x=x+s2*rand(1,10); % add same to all
0100     x=[x,10*randn(1000,3)]; % extra channel with large variance
0101     %[y,M,score]=nt_sca_old(x);
0102     MM=nt_sca(x); yyy=x*MM;
0103     yy=nt_pca(x);
0104     figure(1); clf; 
0105     subplot 121; %bar(abs(nt_normcol(s)'*nt_normcol(y)/size(s,1)));
0106     hold on; bar(abs(nt_normcol(s)'*nt_normcol(yyy)/size(s,1))); 
0107     hold on; bar(abs(nt_normcol(s)'*nt_normcol(yy)/size(s,1)));    legend('sca','pca'); title('source 1');
0108     subplot 122; %bar(abs(nt_normcol(s2)'*nt_normcol(y)/size(s,1)));
0109     hold on;bar(abs(nt_normcol(s2)'*nt_normcol(yyy)/size(s,1))); 
0110     hold on; bar(abs(nt_normcol(s2)'*nt_normcol(yy)/size(s,1)));    legend('sca','pca'); title('source 2');
0111 end
0112     
0113