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:

- nt_cov [c,tw]=nt_cov(x,shifts,w) - time shift covariance
- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_lower_to_full b=nt_lower_to_full(a,ind) - transform lower diagonal to full covariance
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1
- nt_normrow y=nt_normcol(x) - normalize each row so its msq is 1
- nt_sgram [s,f,t]=nt_sgram(x,window,noverlap,nfft,sr,flags) - spectrogram
- nt_spect_plot nt_spect_plot - plot power spectrum
- nt_xprod [y,ind]=nt_xprod(x,flag,dsratio,order) - form all crossproducts

This function is called by:

## SUBFUNCTIONS

- function y=norm2(x,n,ind)

## 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