blociobi.m

剑桥大学用于线性和双线性动力学系统辨识的工具箱

function Z = blociobi(U,Y,i,l,m)
%BLOCIOBI   Assembles Y_i (U_i) for use in the Bilinear Identification Toolbox.
%
%   BLOCIOBI(U,Y,i,l,m) returns Y_{i-1|0} (U_{i-1|0} if U==Y and l==m), where U and Y are
%   block Hankel matrices generated by input and output data, i>0 is an
%   integer, l is the number of outputs and m is the number of inputs of
%   the system.
%
%   See also BILINID, BLOCPLBI, KRONS, BLOCHANK, KHATRI, KRON.
%
% CUED System Identification Toolbox.
% Cambridge University Engineering Department.
% Copyright (C) 1998-2002. All Rights Reserved.

% Version 1.00, Date: 01/06/2002
% Created by H. Chen and E.C. Kerrigan.

[RowY,ColY]=size(Y);
[RowU,ColU]=size(U);

if ColY ~= ColU 
  error('The number of columns of Y and U have to be the same.')
end	
if rem(RowU,i) ~=0 | rem(RowY,i) ~=0 ,
  error('The input/output matrices have to be block matrices of appropriate size.')
end	  

di=0;
for j=1:i,
  di=di+l*(m+1)^(j-1);
end  

% Old code

%y=zeros(di,ColY);
%if i == 1,
%  C=Y(1:l,:);
%else
%  C=Y(1:l,:);
%  for j=1:i-1,
%	A=Y(j*l+1:(j+1)*l,:);
%	B=khatri(U(j*m+1:(j+1)*m,:),C);
%	C=[A; C;B ];
%  end  
%  y(1:di,:)=C;
%end

% More efficient code - avoids unneccesary dynamic memory allocation. See
% the Matlab User's Guide on efficient memory management.

Z = zeros(di,ColY); % preallocate just enough memory
Z(1:l,:) = Y(1:l,:);

if i > 1
  A = zeros(di-l*(m+1)^(i-1),ColY); % preallocate just enough memory
  p = l; % set counter
  for j=1:i-1
	A(1:p,:) = Z(1:p,:); % avoids dynamic memory allocation
	Z(1:l,:) = Y(j*l+1:(j+1)*l,:);
	q = m*p; % number of rows in khatri(U(j*m+1:(j+1)*m,:),Z(1:p,:))
	Z(l+1:l+p,:) = A(1:p,:);
	Z(l+p+1:l+p+q,:) = khatri(U(j*m+1:(j+1)*m,:),A(1:p,:));
	p = p+l+q; % update counter
  end  
end

% *** last line of blociobi.m  ***