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: This function is called by:

SUBFUNCTIONS ^

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

Generated on Mon 30-Jan-2017 18:59:11 by m2html © 2005