nt_split_jd

[idx,score_vector,todss]=nt_split_dss(x,thresh,depth) - segmentation based on joint diagonalization

  idx: index at which to split
  score_vector: from nt_split
  todss: DSS matrix

  x: data
  thresh: truncation threshold for PCA
  depth: recursion depth

0001 function [idx,score_vector,todss]=nt_split_jd(x,thresh,depth);
0002 %[idx,score_vector,todss]=nt_split_dss(x,thresh,depth) - segmentation based on joint diagonalization
0003 %
0004 %  idx: index at which to split
0005 %  score_vector: from nt_split
0006 %  todss: DSS matrix
0007 %
0008 %  x: data
0009 %  thresh: truncation threshold for PCA
0010 %  depth: recursion depth
0011 
0012 if nargin<3||isempty(depth); depth=1; end
0013 if nargin<2||isempty(thresh); thresh=0; end
0014 if isempty(x); error('!'); end
0015 if ndims(x)>2; error('x should be 2D'); end
0016 
0017 [m,n]=size(x);
0018 
0019 % initial PCA to remove below-threshold dimensions
0020 topcs=nt_pca0(x);
0021 z=nt_mmat(x,topcs);
0022 keep=find(mean(z.^2)/mean(z(:,1).^2)>thresh);
0023 z=z(:,keep);
0024 
0025 if m<=2; warning('m==',num2str(m)); end
0026 
0027 idx=ceil(m/2); % initial split into two arbitrary intervals
0028 % iterate until stable
0029 for k=1:10
0030     c1=nt_cov(z(1:idx,:));
0031     c0=c1+nt_cov(z(idx+1:end,:));
0032     todss=nt_dss0(c0,c1);
0033     zz=nt_mmat(z,todss(:,[1,end]));
0034     zz=nt_normcol(zz);
0035     old_idx=idx;
0036     [idx,score_vector]=nt_split(nt_normcol(nt_demean(zz.^2)));
0037     %figure(1); clf; subplot 211; plot(zz(:,1));subplot 212; plot(zz(:,2)); idx, pause
0038     if idx==old_idx; break; end
0039     disp(num2str([idx,old_idx]));
0040 end
0041 todss=topcs(:,keep)*todss; % to return
0042 
0043 if depth>1
0044     [a]=nt_split_jd(x(1:idx,:),thresh, depth-1);
0045     [b]=nt_split_jd(x(idx+1:end,:), thresh,depth-1);
0046     idx=[a,idx,idx+b];
0047     idx=unique(idx);
0048 end
0049 
0050 % disp(['depth, ndims: ' num2str([depth, size(z,2)])])
0051 % disp(['idx: ' num2str([idx])])
0052 % nt_split(nt_normcol(nt_demean(zz.^2)));
0053 
0054 disp(['nt_split_jd nargout: ', num2str(nargout)])
0055 
0056 if nargout==0;
0057     disp(['split at ', num2str(idx)]);
0058     disp(['(%: ', num2str(100*idx/m, '  %.01f'), ')'])
0059     nd=zeros(1,size(x,1));
0060     figure(201); clf
0061 
0062     subplot 312
0063     plot(score_vector);  xlim([1 size(x,1)]); ylim([0 1]); drawnow
0064     nt_mark(idx);
0065     ylabel('score')
0066 
0067     subplot 313
0068     colors='brgcmyk';
0069     hold on
0070     idx2=unique([idx,size(x,1)]); % add on the last sample
0071     old_idx=0;
0072     for iInterval=1:numel(idx2)
0073         z=nt_pca(x(old_idx+1:idx2(iInterval),:));
0074         dim= numel(find( mean(z.^2)/mean(z(:,1).^2) > thresh));
0075         nd(old_idx+1:idx2(iInterval))=dim;
0076         old_idx=idx2(iInterval);
0077     end
0078     z=nt_pca(x);
0079     nd0=numel(find( mean(z.^2)/mean(z(:,1).^2)>thresh));drawnow;
0080     hold on
0081     plot(nd0*ones(size(x,1),1), 'k--')
0082     plot(nd,'b'); 
0083     xlim([1 size(x,1)]); ylim([0 max([nd0,nd])+1]); drawnow
0084     nt_mark(idx);
0085     ylabel('ndims');
0086     
0087     subplot 311
0088     plot(x);  xlim([1 size(x,1)]); drawnow
0089     nt_mark(idx);
0090     if numel(idx)>1; disp(['smallest interval: ', num2str(min(diff(idx)))]); end
0091 end
0092