0001 function [c0,c1,A,todss,pwr0,pwr1]=nt_bias_cluster(x,dsr,flags)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 SMOOTH=1;
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
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
0034 xx=filter(ones(SMOOTH,1),1,xx);
0035
0036
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
0045 NCLUSTERS=2;
0046 [C,A]=vl_kmeans(xx',NCLUSTERS,'algorithm', 'elkan','initialization','plusplus',...
0047 'numrepetitions', 100);
0048
0049
0050 if numel(find(A==1))<numel(A)/2;
0051 C=fliplr(C);
0052 A=3-A;
0053 end
0054
0055
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
0065 c0=squeeze(nt_lower_to_full(C(:,1)',ind));
0066 c1=squeeze(nt_lower_to_full(C(:,2)',ind));
0067 end
0068
0069
0070 [todss,pwr0,pwr1]=nt_dss0(c0+c1,c1);
0071
0072
0073 if nargout==0;
0074
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; imagescc(c0); title('cluster A');
0117 subplot 122; 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