Home > NoiseTools > nt_cov_cluster.m

nt_cov_cluster

PURPOSE ^

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

SYNOPSIS ^

function [A,C]=nt_cov_cluster(x,dsr,flags)

DESCRIPTION ^

[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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [A,C]=nt_cov_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 
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 % time series of cross-products
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 % give each slice the same weight (counters amplitude variations)
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 % cluster the time series (2 clusters)
0039 NCLUSTERS=2;
0040 [C,A]=vl_kmeans(xx',NCLUSTERS,'algorithm', 'elkan','initialization','plusplus',...
0041     'numrepetitions', 100);
0042 
0043 % make sure the first cluster is biggest
0044 if numel(find(A==1))<numel(A)/2;
0045     C=fliplr(C);
0046     A=3-A;
0047 end
0048 
0049 % upsample the cluster ownership index
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 % full covariance matrices from lower diagonal values
0059 c0=squeeze(nt_lower_to_full(C(:,1)',ind));   
0060 c1=squeeze(nt_lower_to_full(C(:,2)',ind));   
0061 end
0062 
0063 % DSS to find components maximally different between clusters
0064 [todss,pwr0,pwr1]=nt_dss0(c0+c1,c1);
0065 
0066 
0067 if nargout==0;
0068     % no output, just plot
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

Generated on Fri 23-May-2014 19:34:17 by m2html © 2005