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