Home > NoiseTools > nt_cluster_jd.m

nt_cluster_jd

PURPOSE ^

[A,todss]=nt_cluster_jd(x,dsr,flags) - cluster with joint diagonalization

SYNOPSIS ^

function [A,todss]=nt_cluster_jd(x,dsr,flags)

DESCRIPTION ^

[A,todss]=nt_cluster_jd(x,dsr,flags) - cluster with joint diagonalization

  A: map of cluster ownership
  todss: 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: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [A,todss]=nt_cluster_jd(x,dsr,flags)
0002 %[A,todss]=nt_cluster_jd(x,dsr,flags) - cluster with joint diagonalization
0003 %
0004 %  A: map of cluster ownership
0005 %  todss: result of DSS
0006 %
0007 %  x: data (time*channels)
0008 %  dsr: downsample ratio for cross product series
0009 %  flags: 'norm': give each dsr-sized slice the same weight
0010 
0011 
0012 if nargin<3; flags=[]; end
0013 if nargin<2; error('!'); end
0014 
0015 if ~exist('vl_kmeans');
0016     disp('vl_kmeans() not found, download from http://www.vlfeat.org');
0017 end
0018 if ndims(x)>2; 
0019     error('x should be time*channels');
0020 end
0021 
0022 disp('entering nt_cluster_jd...');
0023 
0024 % initial clustering
0025 [C0,C1,A,todss]=nt_bias_cluster(x,dsr,flags);
0026 todss2=todss(:,[1 end]); % keep only first and last components
0027 z=nt_mmat(x,todss2);
0028 
0029 figure(2);clf
0030 subplot 511; plot(x); title('data');
0031 subplot 512; plot(A,'.-'); title('initial cluster map');
0032 subplot 513; plot(z(:,1)); title('DSS1');
0033 subplot 514; plot(z(:,2)); title('DSS2');
0034 drawnow;
0035 
0036 % iterate until stable
0037 old_A=A;
0038 for k=1:10
0039     [~,~,A]=nt_bias_cluster(nt_normcol(z),dsr,flags); %cluster first & last components
0040     disp(['cluster sizes: ', num2str([sum(A==1), sum(A==2)])]);
0041     c1=nt_cov(x(find(A==1),:));
0042     c0=c1+nt_cov(x(find(A==2),:));
0043     [todss,pwr0,pwr1]=nt_dss0(c0,c1);
0044     todss2=todss(:,[1 end]); % keep only first and last
0045     z=nt_mmat(x,todss2);
0046     %figure(3); clf; plot(double(A)); pause
0047     subplot 515; plot(A,'.-'); title('final cluster map');
0048     if old_A==A; break; end
0049     old_A=A;
0050 end 
0051 
0052 if nargout==0;
0053     % no output, just plot
0054     disp(['cluster1: ',num2str(100*numel(find(A==1))/numel(A)), '%']);
0055 
0056     z1=nt_mmat(x,todss(:,1));
0057     z2=nt_mmat(x,todss(:,end));
0058     z=nt_mmat(x,todss); 
0059     z=nt_normcol(z);
0060     e1=mean(z(find(A==1),:).^2);
0061     e2=mean(z(find(A==2),:).^2);
0062 
0063     figure(100); clf
0064     plot(x); hold on
0065     xx=x;
0066     xx(find(A==2),:)=nan;
0067     plot(xx,'k');
0068     axis tight
0069     title('black: cluster2');
0070     
0071     figure(101); clf
0072     subplot 121;
0073     plot(pwr1./pwr0,'.-'); xlabel('component'); ylabel('score'); title('DSS cluster A vs B');
0074     subplot 122;
0075     nt_spect_plot(z1,1024,[],[],1);
0076     hold on
0077     nt_spect_plot(z2,1024,[],[],1);
0078     xlim([0 .5])
0079     nt_colorlines([],[1 3]);
0080     legend('first','last'); legend boxoff
0081     hold off
0082 
0083     
0084     figure(102); clf
0085     subplot 211;
0086     plot(z1); axis tight
0087     title('first DSS component')
0088     subplot 212;
0089     plot(z2); axis tight
0090     title('last DSS component');
0091     
0092     figure(103); clf
0093     plot([e1',e2'], '.-'); legend('cluster A', 'cluster B'); title ('power per component');
0094     
0095     figure(104);
0096     subplot 121; imagescc(nt_cov(z(find(A==1),:))); title('cluster A'); 
0097     subplot 122; imagescc(nt_cov(z(find(A==2),:))); title('cluster B');
0098     
0099     if 0 
0100         figure(105); clf
0101         subplot 211;
0102         nt_sgram(z1,1024,32,[],1);
0103         title('first');
0104         subplot 212;
0105         nt_sgram(z2,1024,32,[],1);
0106         title('last');
0107     end
0108     clear c0 c1 A todss pwr0 pwr1
0109     
0110 end
0111 
0112 function y=norm2(x,n,ind)
0113 [I,J]=ind2sub([n,n],ind);
0114 for k=1:size(x,1)
0115     a=x(k,1:n);
0116     b=sqrt(a(I).*a(J));
0117     y(k,:)=x(k,:)./b;
0118 end
0119 
0120     
0121     
0122

Generated on Mon 10-Nov-2014 14:40:42 by m2html © 2005