Home > NoiseTools > nt_cluster_jd_old.m

# nt_cluster_jd_old

## PURPOSE

[IDX,todss,SCORE]=nt_cluster_jd2(x,dsr,flags,init) - cluster with joint diagonalization

## SYNOPSIS

function [IDX,TODSS,SCORE]=nt_cluster_jd2(x,dsr,flags,init)

## DESCRIPTION

    [IDX,todss,SCORE]=nt_cluster_jd2(x,dsr,flags,init) - cluster with joint diagonalization

      IDX: indices of cluster ownership
      TODSS: DSS matrices to emphasize both clusters
      SCORE: scores for both clusters (smaller means better contrast)

      x: data (time*channels)
      dsr: downsample ratio for cross product series
      flags: 'norm', 'norm2': give each dsr-sized slice the same weight
      init: provide initial clustering

     See nt_bias_cluster, nt_cluster1D

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_cluster1D [C,A,score]=nt_cluster1D_b(x) - cluster 1D data into 2 clusters
- 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_linecolors nt_colorlines(h,permutation) - apply different colors to lines of plot
- 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 [IDX,TODSS,SCORE]=nt_cluster_jd2(x,dsr,flags,init)
    0002 %[IDX,todss,SCORE]=nt_cluster_jd2(x,dsr,flags,init) - cluster with joint diagonalization
    0003 %
    0004 %  IDX: indices of cluster ownership
    0005 %  TODSS: DSS matrices to emphasize both clusters
    0006 %  SCORE: scores for both clusters (smaller means better contrast)
    0007 %
    0008 %  x: data (time*channels)
    0009 %  dsr: downsample ratio for cross product series
    0010 %  flags: 'norm', 'norm2': give each dsr-sized slice the same weight
    0011 %  init: provide initial clustering
    0012 %
    0013 % See nt_bias_cluster, nt_cluster1D
    0014
    0015
    0016 if nargin<2; error('!'); end
    0017 if nargin<3 ||isempty(flags); flags=[]; end
    0018 if nargin<4; init=[]; end
    0019
    0020 if ~nargout; disp('entering nt_cluster_jd...'); end
    0021
    0022 % time series of cross products (terms of covariance matrix)
    0023 [xx,ind]=nt_xprod(x,'lower',dsr);
    0024
    0025 % option: give each slice the same weight (counters amplitude variations)
    0026 if find(strcmp(flags,'norm'))
    0027     xx=nt_normrow(xx);
    0028 end
    0029 if find(strcmp(flags,'norm2'))
    0030     xx=norm2(xx,size(x,2),ind);
    0031 end
    0032
    0033 % initial clustering of the time series of cross products
    0034 if isempty(init)
    0035     [C,A,score]=nt_cluster1D(xx);
    0036     [~,idx]=min(score);
    0037     A=A(:,idx);
    0038     % upsample the cluster ownership index
    0039     A=repmat(A',[dsr,1]);
    0040     A=A(:);
    0041     A(end:size(x,1))=A(end);
    0042     IDX{1}=find(A==1);
    0043     IDX{2}=find(A==2);
    0044 else
    0045     IDX{1}=setdiff(1:size(x,1),init);
    0046     IDX{2}=init;
    0047 end
    0048
    0049 % initial DSS to contrast clusters
    0050 c0=nt_cov(x);
    0051 c1=nt_cov(x(IDX{2},:));
    0052 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
    0053 todss2=todss(:,[1 end]); % keep only first and last components
    0054 z=nt_mmat(x,todss2);
    0055
    0056 PLOT_FIG2=0;
    0057 if PLOT_FIG2
    0058     figure(2);  clf; set(gcf, 'name','in nt_cluster_jd');
    0059     A=zeros(size(x,1),1); A(IDX{1})=1;
    0060     subplot 511; plot(x); title('data');
    0061     subplot 512; plot(A,'.-'); title('initial cluster map');
    0062     subplot 513; plot(z(:,1)); title('initial DSS1');
    0063     subplot 514; plot(z(:,2)); title('initial DSS2');
    0064     drawnow; %pause(1);
    0065 end
    0066
    0067 % iterate until stable
    0068 old_IDX=IDX;
    0069 for k=1:10
    0070
    0071     [zz,ind]=nt_xprod(z,[],dsr);
    0072
    0073     % option: give each slice the same weight (counters amplitude variations)
    0074     if find(strcmp(flags,'norm'))
    0075         zz=nt_normrow(zz);
    0076     end
    0077     if find(strcmp(flags,'norm2'))
    0078         zz=norm2(zz,size(z,2),ind);
    0079     end
    0080
    0081     if 0 % expand power terms to improve distribution of small values
    0082         EXPONENT=0.5;
    0083         zz(1,1)=(zz(1,1)).^EXPONENT;
    0084         zz(2,2)=(zz(2,2)).^EXPONENT;
    0085     end
    0086
    0087     % 1D clustering of each of the 3 time series
    0088     [C,A,score]=nt_cluster1D(zz);
    0089     [~,idx]=min(score); % choose best
    0090     A=A(:,idx);
    0091     A=repmat(A',[dsr,1]);
    0092     A=A(:);
    0093     A(end:size(x,1))=A(end);
    0094     IDX{1}=find(A==1);
    0095     IDX{2}=find(A==2);
    0096
    0097     if ~nargout; disp(['cluster sizes: ', num2str([numel(IDX{1}), numel(IDX{2})])]); end
    0098
    0099     % DSS to contrast clusters
    0100     c1=nt_cov(x(IDX{1},:));
    0101     c0=c1+nt_cov(x(IDX{2},:));
    0102     [todss,pwr0,pwr1]=nt_dss0(c0,c1);
    0103     z=nt_mmat(x,todss(:,[1 end])); % keep first and last
    0104
    0105     if PLOT_FIG2
    0106         figure(2);
    0107         subplot 515; plot(A,'.-'); title('final cluster map');
    0108     end
    0109     if all(size(old_IDX{1})==size(IDX{1})) && all(old_IDX{1}==IDX{1}); break; end
    0110     old_IDX=IDX;
    0111 end
    0112
    0113 % scores for each cluster
    0114 c1=nt_cov(x(IDX{1},:));
    0115 c0=c1+nt_cov(x(IDX{2},:));
    0116 [TODSS{1},pwr0,pwr1]=nt_dss0(c0,c1);
    0117 SCORE(1)=pwr1(end)/pwr0(end);
    0118
    0119 c1=nt_cov(x(IDX{2},:));
    0120 c0=c1+nt_cov(x(IDX{1},:));
    0121 [TODSS{2},pwr0,pwr1]=nt_dss0(c0,c1);
    0122 SCORE(2)=pwr1(end)/pwr0(end);
    0123
    0124 SCORE=SCORE(:);
    0125
    0126 if nargout==0;
    0127
    0128     % no output, just plot
    0129     disp(['cluster1: ',num2str(100*numel(find(A==1))/numel(A)), '%']);
    0130
    0131     figure(100); clf
    0132     subplot 311;
    0133     plot(x); hold on
    0134     xx=x;
    0135     xx(IDX{1},:)=nan;
    0136     plot(xx,'k');
    0137     axis tight
    0138     title('black: cluster B');
    0139
    0140     z1=nt_mmat(x,todss(:,1));
    0141     z2=nt_mmat(x,todss(:,end));
    0142
    0143     figure(101); clf ;
    0144     subplot 221;
    0145     plot(pwr1./pwr0,'.-'); xlabel('component'); ylabel('score'); title('DSS cluster A vs all');
    0146     subplot 222;
    0147     wsize=min(1024,size(z1,1));
    0148     nt_spect_plot(z1,wsize,[],[],1);
    0149     hold on
    0150     wsize=min(1024,size(z2,1));
    0151     nt_spect_plot(z2,wsize,[],[],1);
    0152     xlim([0 .5])
    0153     nt_linecolors([],[1 3]);
    0154     legend('first','last'); legend boxoff
    0155     hold off
    0156
    0157     z=nt_mmat(x,todss);
    0158     z=nt_normcol(z);
    0159     subplot 223; imagescc(nt_cov(z(IDX{1},:))); title('cluster A');
    0160     subplot 224; imagescc(nt_cov(z(IDX{2},:))); title('cluster B');
    0161
    0162
    0163     figure(100);
    0164     subplot 312;
    0165     plot(z1); axis tight
    0166     title('first DSS component')
    0167     subplot 313;
    0168     plot(z2); axis tight
    0169     title('last DSS component');
    0170
    0171 %     figure(103); clf
    0172 %     e1=mean(z(find(A==1),:).^2);
    0173 %     e2=mean(z(find(A==2),:).^2);
    0174 %     plot([e1',e2'], '.-'); legend('cluster A', 'cluster B'); title ('power per component');
    0175
    0176
    0177     if 0
    0178         figure(105); clf
    0179         subplot 211;
    0180         nt_sgram(z1,1024,32,[],1);
    0181         title('first');
    0182         subplot 212;
    0183         nt_sgram(z2,1024,32,[],1);
    0184         title('last');
    0185     end
    0186     clear IDX SCORE TODSS
    0187
    0188 end
    0189
    0190 function y=norm2(x,n,ind)
    0191 [I,J]=ind2sub([n,n],ind);
    0192 for k=1:size(x,1)
    0193     a=x(k,1:n);
    0194     b=sqrt(a(I).*a(J));
    0195     y(k,:)=x(k,:)./b;
    0196 end
    0197
    0198
    0199
    0200