uy2hy0.m

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

% UY2HY0 - From data to impulse and sequential free responses.
%
% [h, y0, delta] = uy2hy0(u, y, lmax, l, delta, j)
%
% U,Y - input (TxM matrix) and output (TxP matrix)
% LMAX  - upper bound for the system lag
% L     - optional number of samples in a computed block
% DELTA - optional desired length of the responses 
%   If skipped the response is computed till 
%   ||H(delta)||_F > EPS
% J - optional number of free responses
% H  - first 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 J, DELTA samples long sequential free resp.
%   the i-th response is Y(:,i) 

function [h,y0,delta] = uy2h(u,y,lmax,l,delta,j)

% Constants
EPS  = 1e-5; % tolerance for the convergence of the response
IMAX = 100;  % maximum response length
[T,m] = size(u);
p = size(y,2);

% Check inputs
if (nargin < 3)
  error('Not enough input arguments');
end
if (nargin < 4 | isempty(l))
  l = lmax;
end
L = lmax + l; % # of rows in the Hankel matrix
if (nargin < 5 | isempty(delta))
  find_delta = 1;
  delta = IMAX; 
else
  find_delta = 0;
end
if (nargin < 6 | isempty(j))
  j = T - L + 1;
end

% SVD of the Hankel matrix of the data
r   = triu(qr([blkhank(u,L);blkhank(y,L)]'))';
in  = 1:L*m+lmax*p;
r11 = r(in,in);
r21 = r(in(end)+1:end,in);
pinvUY1 = pinv(r11);
Y2      = r21;

% Initialization for H
fu_h = zeros(L*m,m); 
fu_h(lmax*m+1:(lmax+1)*m,:) = eye(m);
fy_h = zeros(L*p,m);

% Initialization for Y0
fu_y0 = zeros(L*m,j); 
fu_y0(1:lmax*m,:) = blkhank(u,lmax,j);
fy_y0 = zeros(L*p,j);
fy_y0(1:lmax*p,:) = blkhank(y,lmax,j);

% Iteration
imax = ceil(delta/l);
h    = zeros(p*l*imax,m);
y0   = zeros(delta*p,j);
for i = 1:imax  
  % Solve the system for H
  g = pinvUY1 * [fu_h; fy_h(1:p*lmax,:)];
  fy_h(p*lmax+1:end,:)  = Y2 * g;
  % Solve the system for Y0
  g = pinvUY1 * [fu_y0; fy_y0(1:p*lmax,:)];
  fy_y0(p*lmax+1:end,:) = Y2 * g;
  % Store the results
  h((i-1)*p*l+1:i*p*l,:)  = fy_h(p*lmax+1:end,:);  
  y0((i-1)*l*p+1:i*l*p,:) = fy_y0(p*lmax+1:end,:);
  % Check exit condition (if delta is not given)
  if (find_delta & i*l*min(m,p) > ...
      2*lmax & norm(fy_h(p*lmax+1:end,:),'fro') < EPS)
    delta = i*l;
    break;
  end
  % Shift FU_H, FY_H for the next iteration
  fu_h(1:m*l,:) = []; 
  fu_h = [fu_h; zeros(m*l,m)];
  fy_h(1:p*l,:) = []; 
  fy_h = [fy_h; zeros(p*l,m)];
  % Shift FU_Y0, FY_Y0 for the next iteration
  fu_y0(1:m*l,:) = []; 
  fu_y0 = [fu_y0; zeros(m*l,j)];
  fy_y0(1:p*l,:) = []; 
  fy_y0 = [fy_y0; zeros(p*l,j)];  
end 

if ( m > 1 | p > 1 ) % MIMO
  % Construct a 3d output array
  h_ = h;
  h = zeros(p,m,delta);
  for i = 1:delta
    h(:,:,i) = h_((i-1)*p+1:i*p,:);
  end
else % SISO
  h = h(1:delta,:);
end
  
y0 = y0(1:delta*p,:);