nt_bias_cluster

[c0,c1,A,todss,pwr0,pwr1]=nt_bias_cluster(x,dsr,flags) - cluster covariance

  c0,c1: covariance matrices of clusters
  A: map of cluster ownership
  todss,pwr0,pwr1: result of DSS

  x: data (time*channels)
  dsr: downsample ratio for cross product series
  flags: 'norm': give each dsr-sized slice the same weight

  See: nt_cluster1D, nt_cluster_jd.

0001 function [c0,c1,A,todss,pwr0,pwr1]=nt_bias_cluster(x,dsr,flags)
0002 %[c0,c1,A,todss,pwr0,pwr1]=nt_bias_cluster(x,dsr,flags) - cluster covariance
0003 %
0004 %  c0,c1: covariance matrices of clusters
0005 %  A: map of cluster ownership
0006 %  todss,pwr0,pwr1: result of DSS
0007 %
0008 %  x: data (time*channels)
0009 %  dsr: downsample ratio for cross product series
0010 %  flags: 'norm': give each dsr-sized slice the same weight
0011 %
0012 %  See: nt_cluster1D, nt_cluster_jd.
0013 
0014 SMOOTH=1;%2; % smooth the time series of cross products
0015 
0016 if nargin<3 ||isempty(flags); flags=[]; end
0017 if nargin<2; error('!'); end
0018 
0019 if ~exist('vl_kmeans');
0020     disp('vl_kmeans() not found, download from http://www.vlfeat.org');
0021 end
0022 if ndims(x)>2; 
0023     error('x should be time*channels');
0024 end
0025 
0026 % time series of cross-products
0027 if find(strcmp(flags,'nodiag'))
0028     [xx,ind]=nt_xprod(x,'nodiag',dsr);
0029 else
0030     [xx,ind]=nt_xprod(x,'lower',dsr);
0031 end
0032 
0033 % smooth
0034 xx=filter(ones(SMOOTH,1),1,xx); 
0035 
0036 % give each slice the same weight (counters amplitude variations)
0037 if find(strcmp(flags,'norm'))
0038     xx=nt_normrow(xx);
0039 end
0040 if find(strcmp(flags,'norm2'))
0041     xx=norm2(xx,size(x,2),ind);
0042 end
0043 
0044 % cluster the time series (2 clusters)
0045 NCLUSTERS=2;
0046 [C,A]=vl_kmeans(xx',NCLUSTERS,'algorithm', 'elkan','initialization','plusplus',...
0047     'numrepetitions', 100);
0048 
0049 % make sure the first cluster is biggest
0050 if numel(find(A==1))<numel(A)/2;
0051     C=fliplr(C);
0052     A=3-A;
0053 end
0054 
0055 % upsample the cluster ownership index
0056 A=repmat(A,[dsr,1]);
0057 A=A(:);
0058 A(end:size(x,1))=A(end);
0059 
0060 if 1
0061 c0=nt_cov(x(find(A==1),:));
0062 c1=nt_cov(x(find(A==2),:));
0063 else
0064 % full covariance matrices from lower diagonal values
0065 c0=squeeze(nt_lower_to_full(C(:,1)',ind));   
0066 c1=squeeze(nt_lower_to_full(C(:,2)',ind));   
0067 end
0068 
0069 % DSS to find components maximally different between clusters
0070 [todss,pwr0,pwr1]=nt_dss0(c0+c1,c1);
0071 
0072 
0073 if nargout==0;
0074     % no output, just plot
0075     disp(['cluster1: ',num2str(100*numel(find(A==1))/numel(A)), '%']);
0076 
0077     z1=nt_mmat(x,todss(:,1));
0078     z2=nt_mmat(x,todss(:,end));
0079     z=nt_mmat(x,todss); 
0080     z=nt_normcol(z);
0081     e1=mean(z(find(A==1),:).^2);
0082     e2=mean(z(find(A==2),:).^2);
0083 
0084     figure(100); clf
0085     plot(x); hold on
0086     x(find(A==2),:)=nan;
0087     plot(x,'k');
0088     axis tight
0089     title('black: cluster2');
0090     
0091     figure(101); clf
0092     subplot 121;
0093     plot(pwr1./pwr0,'.-'); xlabel('component'); ylabel('score'); title('DSS cluster A vs B');
0094     subplot 122;
0095     nt_spect_plot(z1,1024,[],[],1);
0096     hold on
0097     nt_spect_plot(z2,1024,[],[],1);
0098     xlim([0 .5])
0099     nt_linecolors([],[1 3]);
0100     legend('first','last'); legend boxoff
0101     hold off
0102 
0103     
0104     figure(102); clf
0105     subplot 211;
0106     plot(z1); axis tight
0107     title('first DSS component')
0108     subplot 212;
0109     plot(z2); axis tight
0110     title('last DSS component');
0111     
0112     figure(103); clf
0113     plot([e1',e2'], '.-'); legend('cluster A', 'cluster B'); title ('power per component');
0114     
0115     figure(104);
0116     subplot 121; nt_imagescc(c0); title('cluster A'); 
0117     subplot 122; nt_imagescc(c1); title('cluster B');
0118     
0119     if 0 
0120         figure(105); clf
0121         subplot 211;
0122         nt_sgram(z1,1024,32,[],1);
0123         title('first');
0124         subplot 212;
0125         nt_sgram(z2,1024,32,[],1);
0126         title('last');
0127     end
0128     clear c0 c1 A todss pwr0 pwr1
0129     
0130 end
0131 
0132 function y=norm2(x,n,ind)
0133 [I,J]=ind2sub([n,n],ind);
0134 for k=1:size(x,1)
0135     a=x(k,1:n);
0136     b=sqrt(a(I).*a(J));
0137     y(k,:)=x(k,:)./b;
0138 end
0139 
0140     
0141     
0142