uy2h.m

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

% UY2H - From data to impulse response. (Recursive algorithm)
%
% [h, delta] = uy2h(u, y, lmax, l, delta)
%
% 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 impulse response 
%   If skipped the response is computed till 
%   ||H(delta)||_F > EPS.
% 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

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

% Constants
EPS  = 1e-5;
IMAX = 100;
[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

% Needed in the iterative computation
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
fu = zeros(L*m,m); 
fu(lmax*m+1:(lmax+1)*m,:) = eye(m);
fy = zeros(L*p,m);

% Iteration
imax = ceil(delta/l);
h = zeros(p*l*imax,m);
for i = 1:imax  
  % Solve the system
  g = pinvUY1 * [fu;fy(1:lmax*p,:)];
  fy(lmax*p+1:end,:) = Y2 * g;
  % Store the result
  h((i-1)*p*l+1:i*p*l,:) = fy(lmax*p+1:end,:);  
  % Check exit condition if delta is not given
  if find_delta 
    if ( i*l*min(m,p) > 2*lmax & ...
         norm(fy(lmax*p+1:end,:),'fro') < EPS)
      delta = i*l;
      break;
    end
  end
  % Shift for the next iteration
  fu(1:m*l,:) = []; 
  fu = [fu; zeros(m*l,m)];
  fy(1:p*l,:) = []; 
  fy = [fy; zeros(p*l,m)];
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