h2ox.m

A Matlab toolbox for exact linear time-invariant system identification is presented. The emphasis is

% H2OX - From an impulse response to an extended 
%        observability matrix and a state sequence.
%
% [O,X] = h2ox(h,lmax,n,tol)
%
% H - first T>2*LMAX samples of the impulse response
%     for MIMO (M-inputs, P-outputs) system, H is PxMxT
%     for SISO system, H can be Tx1 vector
% LMAX - upper bound on the system lag
% 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.
% O - extended observability matrix
% X - state sequence

function [O,X] = h2ox(h,lmax,n,tol)

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

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

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

% Rank revealing factorization
sqrt_s = sqrt(1./s(1:n))';
O = sqrt_s(ones(size(U,1),1),:) .* U(:,1:n);
if nargout > 1
  X = (sqrt_s(ones(size(V,1),1),:) .* V(:,1:n))';
end