Home > NoiseTools > nt_split_jd.m

# nt_split_jd

## PURPOSE

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

## SYNOPSIS

function [idx,score_vector,todss]=nt_split_jd(x,thresh,depth);

## DESCRIPTION

    [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

## CROSS-REFERENCE INFORMATION

This function calls:

- nt_cov [c,tw]=nt_cov(x,shifts,w) - time shift covariance
- nt_demean [y,mn]=nt_demean(x,w) - remove weighted mean over cols
- nt_dss0 [todss,pwr1,pwr2]=nt_dss0(c0,c1,keep1,keep2) - dss from covariance
- nt_mark nt_mark(idx,labels,line_params,text_params)
- nt_mmat y=nt_mmat(x,m) - matrix multiplication (with convolution)
- nt_normcol [y,norm]=nt_normcol(x,w) - normalize each column so its weighted msq is 1
- nt_pca [z,idx]=nt_pca(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_pca0 [topcs,pwr,y]=nt_pca0(x,shifts,nkeep,threshold,w) - time-shift pca
- nt_split [idx,score_vector,score]=nt_split(x,depth,thresh,guard,minstep) - split time series into intervals
- nt_split_jd [idx,score_vector,todss]=nt_split_dss(x,thresh,depth) - segmentation based on joint diagonalization

This function is called by:

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

## SOURCE CODE

    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