uy2y0.m

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

% UY2Y0 - From data to sequential free responses.
%
% [y0, delta] = uy2y0(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 block
% DELTA - optional desired length of the impulse response 
%   If skipped the response the stop condition is 
%   |y0(delta)| > EPS.
% Y0 - matrix of J, DELTA samples long sequential free resp.
%   the i-th response is Y(:,i) 

function [y0,delta] = uy2y0(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
fu = zeros(L*m,j); 
fu(1:lmax*m,:) = blkhank(u,lmax,j);
fy = zeros(L*p,j);
fy(1:lmax*p,:) = blkhank(y,lmax,j);

% Iteration
maxi = ceil(delta/l);
y0   = zeros(delta*p,j);
for i = 1:maxi
  % Solve the system
  g = pinvUY1 * [fu;fy(1:lmax*p,:)];
  fy(lmax*p+1:end,:) = Y2 * g;
  % Store the result
  y0((i-1)*l*p+1:i*l*p,:) = fy(lmax*p+1:end,:);
  % Check exit condition (if delta is not given)
  if (find_delta & norm(fy(lmax*p+1:end,1)) < EPS)
    delta = i;
    break;
  end
  % Shift for the next iteration
  fu(1:m*l,:) = []; 
  fu = [fu; zeros(m*l,j)];
  fy(1:p*l,:) = []; 
  fy = [fy; zeros(p*l,j)]; 
end 

y0 = y0(1:delta*p,:);