0001 function [IDX,TODSS,SCORE]=nt_cluster_jd(x,dsr,smooth,flags,init,verbose)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
0038
0039
0040
0041
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
0050
0051
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);
0065 A=A(:,idx);
0066
0067
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
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]);
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
0095 old_IDX=IDX{1};
0096 for k=1:10
0097
0098 [zz,ind]=nt_xprod(z,[],dsr);
0099 zz=zz(:,1:2);
0100 zz=log2(zz+eps);
0101 [C,A]=nt_cluster1D(zz);
0102 [~,idx]= max(diff(C));
0103 A=A(:,idx);
0104
0105
0106
0107 A=double(nt_smooth(A,smooth, [],1)>=1/smooth);
0108
0109
0110 A=repmat(A',[dsr,1]);
0111 A=A(:);
0112 A(end:size(x,1))=A(end);
0113 IDX{1}=find(A==0);
0114
0115
0116
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]));
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
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
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
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
0201 nt_sgram(z1,1024,32,[],1);
0202 title('first');
0203 subplot 212;
0204
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