Home > NoiseTools > nt_bias_cluster.m

# nt_bias_cluster

## PURPOSE

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

## SYNOPSIS

function [c0,c1,A,todss,pwr0,pwr1]=nt_bias_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

      See: nt_cluster1D, nt_cluster_jd.

## 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_imagescc nt_imagescc - plot image with symmetric scaling
- nt_linecolors nt_colorlines(h,permutation) - apply different colors to lines of plot
- 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,normrow_flag) - form all crossproducts

This function is called by:

## SUBFUNCTIONS

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

## SOURCE CODE

    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