blocplbi.m

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

function Y = blocplbi(U,i,m)
%BLOCPLBI   Assembles U+ for use in the Bilinear Identification Toolbox.
%
%   BLOCPLBI(U,i,m) returns U^{+}_{i-1|0}, where U is a block Hankel
%   matrix, i>0 is an integer and m is the number of inputs of the system.
%
%   See also BILINID, BLOCIOBI, 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.

% Old code

%[RowU,ColU]=size(U);
%if rem(RowU,i) ~= 0 | rem(RowU,m) ~= 0 ,
%  error('The input/output  matrix  has to be a N block matrix')
%end	  
%d=(m+1)^(i)-1;
%y=zeros(d,ColU);
%if i == 1,
%  C=U(1:m,:);
%else
%  C=U(1:m,:);
%  for j=1:i-1,
%	A=U(j*m+1:(j+1)*m,:);
%	B=khatri(U(j*m+1:(j+1)*m,:),C);
%	C=[C; A;B ];
%  end  
%  y(1:d,:)=C;
%end


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

[RowU,ColU]=size(U);
if rem(RowU,i) ~= 0 | rem(RowU,m) ~= 0
  error('The number of rows of U need to be divisible by i and m.')
end	  

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

Y = zeros(ei,ColU); % preallocate just enough memory
Y(1:m,:) = U(1:m,:);

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

% *** last line of blocplbi.m  ***