vb_temporal_pca

 Temporal PCA
 [z,indxT,s] = vb_temporal_pca(x, Twin, smin)
 --- Input
 x : time series signal [NX x T]
 x(n,t) : n-th component signal at time t , t=1:T, n=1:NX
 Twin = [Nbck , Nfwd] : time shift length for backward and forward direction
  if Twin is scalar, Nbck = Nfwd = Twin is assumed
   smin : >= 1  : number of components for threshould PCA (muse be integer, default 40)
           < 1  : contribution rate for threshould PCA
 --- Output
  z : Time shifted signal after orthogonalization by PCA
  indxT : time index for z
  z(:,:) <-> x(:,indxT)
  indxT = (1+Nbck):(T-Nfwd))
  s : eigen value

 2007-8-21 Masa-aki Sato
 2015-4-20 Yusuke Fujiwara
 2015-6-12 Yusuke Takeda

0001 function    [z,indxT,s] = vb_temporal_pca(x, Twin, smin)
0002 % Temporal PCA
0003 % [z,indxT,s] = vb_temporal_pca(x, Twin, smin)
0004 % --- Input
0005 % x : time series signal [NX x T]
0006 % x(n,t) : n-th component signal at time t , t=1:T, n=1:NX
0007 % Twin = [Nbck , Nfwd] : time shift length for backward and forward direction
0008 %  if Twin is scalar, Nbck = Nfwd = Twin is assumed
0009 %   smin : >= 1  : number of components for threshould PCA (muse be integer, default 40)
0010 %           < 1  : contribution rate for threshould PCA
0011 % --- Output
0012 %  z : Time shifted signal after orthogonalization by PCA
0013 %  indxT : time index for z
0014 %  z(:,:) <-> x(:,indxT)
0015 %  indxT = (1+Nbck):(T-Nfwd))
0016 %  s : eigen value
0017 %
0018 % 2007-8-21 Masa-aki Sato
0019 % 2015-4-20 Yusuke Fujiwara
0020 % 2015-6-12 Yusuke Takeda
0021 
0022 [NX,T] = size(x);
0023 
0024 D  = sum(Twin) + 1;    % time shift length
0025 TD = T - sum(Twin);    % Effective time length for process
0026 
0027 if TD <= 0, error('Signal time length is too short'); end;
0028 
0029 % subtract bias component
0030 x = vb_repadd(x, - mean(x,2));
0031 
0032 % --- Time shifted embedding index
0033 %
0034 %  z(:,t) = [ x(:,t-Nbck);
0035 %                ...;
0036 %             x(:,t);
0037 %                ...;
0038 %             x(:,t+Nfwd) ]
0039 
0040 indxT = (1+Twin(1)):(T-Twin(2));
0041 indxD = -Twin(1):Twin(2);
0042 indx  = repmat( indxD(:),[1 TD]) + repmat(indxT,[D 1]);
0043 indx  = indx(:);
0044 
0045 X = reshape( x(:,indx), [NX*D TD]);
0046 
0047 % --- SVD ---
0048 % X = U * S * V , U:[N x N], S=diag(s):[N x N], V:[N x T]
0049 %     U' * U = V * V' = eye(N)
0050 % --- PCA ---
0051 %  X * X'  = U * S^2 * U'
0052 % (X * X') * U = U * S^2
0053 %  V = S^(-2) * U' * X
0054 % [U, D] = eig(A) : A * U = U * D , D = diag(d)
0055 
0056 %[U,S,V] = svd( X * X' );
0057 
0058 [U,S] = eig( X * X' );
0059 s  = diag(S);
0060 if smin<1
0061     ix = find( s/max(s) > smin );
0062 else
0063     Npcs = length(s);
0064     ix = (Npcs - smin + 1):Npcs;
0065 end
0066 
0067 % --- Orthogonalize time shifted signal
0068 z  = U(:,ix)' * X;
0069 
0070 % Normalization
0071 z  = vb_repmultiply( z, 1./sqrt(sum(z.^2,2)));
0072 
0073 err1 = 1 - sum(s(ix))/sum(s);
0074 err2 = max(max(abs( z*z' - eye(size(z,1)))));
0075 
0076 fprintf('Dim of noise space = %d, ',length(ix))
0077 fprintf('PCA err = %g, Ortho err = %g\n',err1,err2)
0078 
0079 return
0080