hy02xbal.m

来自「A Matlab toolbox for exact linear time-i」· M 代码 · 共 50 行

% HY02XBAL - From impulse response and sequential free 
%            responses to a balanced state sequence.
%
% XBAL = hy02xbal(h,y0,delta,n,tol)
%
% H - first T>2*DELTA samples of the impulse response
%     for MIMO (M-inputs, P-outputs) system, H is PxMxT
%     for SISO system, H can be Tx1 vector
% Y0 - matrix of DELTA samples long sequential free resp. 
%      (Y0(:,i) is such a response)
% DELTA - optional finite time balancing parameter
%     DELTA>LMAX+1, LMAX - upper bound on the system lag 
%     default: DELTA = floor(T/2) (the maximum)
% N    - optional order of the system
% TOL - tolerance for order selection, default 1e-7
%   If TOL < 0, the SVs are plotted and an input is asked.
% XBAL - finite time DELTA balanced state sequence

function xbal = hy02xbal(h,y0,delta,n,tol)

% Remove H(0) from H
if length(size(h)) == 2
  T = size(h,1);
  h = h(2:end);
else
  T = size(h,3);
  h = h(:,:,2:end);
end

if nargin < 3 | isempty(delta)
  delta = floor(T/2);
end

% SVD of the Hankel matrix of the Markov parameters
[U,S,V] = svd(blkhank(h,delta),0);
s = diag(S);

% If not given, determine the order
if (nargin < 5 | isempty(tol))
  tol = [];
end
if (nargin < 4 | isempty(n))
  n = order(s,tol);
end

% Compute balanced state sequence
sqrt_s = sqrt(1./s(1:n));
xbal = (sqrt_s(:,ones(1,size(U,1))) .* U(:,1:n)') ...
       * y0(1:size(U,1),:);