Home > NoiseTools > nt_cluster_jd.m

nt_cluster_jd

PURPOSE ^

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

SYNOPSIS ^

function [IDX,TODSS,SCORE]=nt_cluster_jd(x,dsr,smooth,flags,init,verbose)

DESCRIPTION ^

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

  IDX: cluster ownership (IDX{1}: low amp, IDX{2{: high amp)
  TODSS: DSS matrix (1st column --> discriminating component)
  SCORE: score (smaller means better contrast)

  x: data (time*channels)
  dsr: downsample ratio for cross product series
  smooth: further smooth cross-product series
  flags: 'norm', 'norm2': give each dsr-sized slice the same weight
  init: provide initial clustering
  verbose: display & plot (default=no)

 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_jd(x,dsr,smooth,flags,init,verbose)
0002 %[IDX,todss,SCORE]=nt_cluster_jd(x,dsr,smooth,flags,init,verbose) - cluster with joint diagonalization
0003 %
0004 %  IDX: cluster ownership (IDX{1}: low amp, IDX{2{: high amp)
0005 %  TODSS: DSS matrix (1st column --> discriminating component)
0006 %  SCORE: score (smaller means better contrast)
0007 %
0008 %  x: data (time*channels)
0009 %  dsr: downsample ratio for cross product series
0010 %  smooth: further smooth cross-product series
0011 %  flags: 'norm', 'norm2': give each dsr-sized slice the same weight
0012 %  init: provide initial clustering
0013 %  verbose: display & plot (default=no)
0014 %
0015 % See nt_bias_cluster, nt_cluster1D
0016 
0017 
0018 if nargin<2; error('!'); end
0019 if nargin<3 ||isempty(smooth); smooth=1; end
0020 if nargin<4 ||isempty(flags); flags=[]; end
0021 if nargin<5; init=[]; end
0022 if nargin<6||isempty(verbose); verbose=0; end
0023 
0024 if ~nargout; disp('entering nt_cluster_jd...'); end
0025 
0026 if ndims(x)>2 || size(x,2) ==1;
0027     error('should be 2D matrix');
0028 end
0029 
0030 
0031 %{
0032  Calculate the time series of cross products (terms of the covariance matrix).
0033  This time series has coarser temporal resolution than x by a factor dsr.
0034 %}
0035 [xx,ind]=nt_xprod(x,'lower',dsr);  
0036 
0037 % figure(2); clf;
0038 % subplot 211;
0039 % plot(xx)
0040 
0041 % option: give each slice the same weight (counters amplitude variations)
0042 if find(strcmp(flags,'norm'))
0043     xx=nt_normrow(xx);
0044 end
0045 if find(strcmp(flags,'norm2'))
0046     xx=norm2(xx,size(x,2),ind);
0047 end
0048 
0049 % subplot 212;
0050 % plot(xx);
0051 % pause;
0052 
0053 smooth
0054 size(xx)
0055 xx=nt_smooth(xx,smooth,[],1);
0056 
0057 %{
0058 Cluster each column the time series of cross products, 
0059 choose the column with best score (reduction in energy), 
0060 and use it's cluster index to initialize the first JD analysis.
0061 %}
0062 if isempty(init)
0063     [C,A,score]=nt_cluster1D(xx);
0064     [~,idx]=min(score); % select column with best score (tightest clusters)
0065     A=A(:,idx); 
0066         
0067     % upsample the cluster ownership index so we can apply it to x
0068     A=repmat(A',[dsr,1]);
0069     A=A(:);
0070     A(end:size(x,1))=A(end);
0071     IDX{1}=find(A==0);
0072 else
0073     IDX{1}=init;
0074 end
0075 
0076 % initial DSS to contrast clusters
0077 c0=nt_cov(x);
0078 c1=nt_cov(x(IDX{1},:));
0079 [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0080 todss2=todss(:,[1 end]); % keep only first and last components
0081 z=nt_mmat(x,todss2);
0082 
0083 PLOT_FIG2=0;
0084 if PLOT_FIG2
0085     figure(2);  clf; set(gcf, 'name','in nt_cluster_jd');
0086     A=zeros(size(x,1),1); A(IDX{1})=1;
0087     subplot 411; plot(x); title('data');
0088     subplot 412; plot(A,'.-'); title('initial cluster map');
0089     subplot 413; plot(z(:,1)); title('initial DSS1');
0090     subplot 414; plot(z(:,2)); title('initial DSS2');
0091     drawnow; pause;
0092 end
0093 
0094 % iterate until stable
0095 old_IDX=IDX{1};
0096 for k=1:10
0097 
0098     [zz,ind]=nt_xprod(z,[],dsr);
0099     zz=zz(:,1:2);       % keep only the squares
0100     zz=log2(zz+eps);    % log to make it sensitive to power ratio
0101     [C,A]=nt_cluster1D(zz);
0102     [~,idx]= max(diff(C));    % the best component is the one with greatest ratio
0103     A=A(:,idx);
0104     
0105     %plot(nt_smooth(A,smooth, [],1));
0106     
0107     A=double(nt_smooth(A,smooth, [],1)>=1/smooth); % extend ownership to include effect of smoothing
0108 
0109     % upsample the cluster ownership index so we can apply it to x
0110     A=repmat(A',[dsr,1]); % upsample
0111     A=A(:); 
0112     A(end:size(x,1))=A(end);
0113     IDX{1}=find(A==0); % 0: low values, 1: high values
0114     
0115     
0116     % DSS to contrast clusters
0117     c0=nt_cov(x)/size(x,1);
0118     c1=nt_cov(x(IDX{1},:))/size(x(IDX{1},:),1);
0119     [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0120     z=nt_mmat(x,todss(:,[1 end])); % keep first and last
0121 
0122     if ~nargout||verbose; 
0123         disp(['low amp cluster: ', num2str((100*numel(IDX{1})/size(x,1)), 2), ' % of samples, power ratio: ' num2str(pwr1(end)/pwr0(end), 3)]); 
0124     end
0125 
0126     if PLOT_FIG2
0127         figure(2);  
0128         subplot 515; plot(A,'.-'); title('final cluster map');
0129     end
0130     if all(size(old_IDX)==size(IDX{1})) && all(old_IDX==IDX{1}); break; end
0131     old_IDX=IDX{1};
0132 end 
0133 IDX{2}=setdiff(1:size(x,1), IDX{1});
0134 
0135 % score
0136 [TODSS,pwr0,pwr1]=nt_dss0(c0,c1);
0137 SCORE=pwr1(end)/pwr0(end);
0138 
0139 nargout
0140 
0141 if nargout==0||verbose;
0142     
0143     % no output, just plot
0144 
0145     z1=nt_mmat(x,TODSS(:,1));
0146     z2=nt_mmat(x,TODSS(:,end));
0147 
0148     figure(101); clf ;
0149     subplot 221;
0150     plot(pwr1./pwr0,'.-'); xlabel('component'); ylabel('score'); title('DSS cluster [low amp] vs all');
0151     subplot 222;
0152     wsize=min(1024,size(z2,1));
0153     %nt_spect_plot(z1/sqrt(mean(z1(:).^2)),wsize,[],[],1);
0154     hold on
0155     nt_spect_plot(z2/sqrt(mean(z2(:).^2)),wsize,[],[],1);
0156     nt_spect_plot(x/sqrt(mean(x(:).^2)),wsize,[],[],1);
0157     xlim([0 .5]);
0158     nt_linecolors([],[1 3 2]);
0159     legend('last','all'); legend boxoff
0160     hold off
0161 
0162     z=nt_mmat(x,todss); 
0163     z=nt_normcol(z);
0164     subplot 223; imagescc(nt_cov(z(IDX{1},:))); title('cluster [low amp]'); 
0165     subplot 224; imagescc(nt_cov(z)-nt_cov(z(IDX{1},:))); title('cluster [high amp]');
0166 
0167     
0168     figure(102); clf
0169     if 0
0170         subplot 311;
0171         plot(x); hold on
0172         xx=x; xx(IDX{1},:)=nan;
0173         plot(xx,'k');
0174         axis tight
0175         title('black: cluster [high amp]');
0176         subplot 312;
0177         plot(z1); axis tight
0178         title('first DSS component');
0179         subplot 313;
0180         plot(z2); axis tight
0181         title('last DSS component');
0182     else
0183         subplot 311;
0184         plot(x); hold on
0185         xx=x; xx(IDX{1},:)=nan;
0186         plot(xx,'k');
0187         axis tight
0188         title('black: cluster [high amp]');
0189         subplot 312;
0190         plot(z2); axis tight
0191         title('last DSS component');
0192         subplot 313;
0193         nt_sgram(z2,128,1); axis tight
0194         title('last DSS component');
0195     end
0196     
0197     if 1 
0198         figure(105); clf
0199         subplot 211;
0200         %sr=44100; ERBspect(z1,sr);
0201         nt_sgram(z1,1024,32,[],1);
0202         title('first');
0203         subplot 212;
0204         %ERBspect(z2,sr);
0205         nt_sgram(z2,1024,32,[],1);
0206         title('last');
0207     end
0208     if nargout==0; clear IDX SCORE TODSS; end
0209     
0210 end
0211 
0212 function y=norm2(x,n,ind)
0213 [I,J]=ind2sub([n,n],ind);
0214 for k=1:size(x,1)
0215     a=x(k,1:n);
0216     b=sqrt(a(I).*a(J));
0217     y(k,:)=x(k,:)./b;
0218 end
0219 
0220     
0221     
0222

Generated on Mon 28-Nov-2016 20:12:47 by m2html © 2005