soltritoep.m

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

function [T,H] = soltritoep(SS,PP,k)
%SOLTRITOEP   Solves for a lower triangular block Toeplitz matrix.
%
%   [T,H] = soltritoep(S,P,k) 
%
%   The matrix
% 
%        [ H(0)    0    ...  0    0  ]
%        [ H(1)   H(0)  ...  0    0  ]
%    T = [  :                     :  ]
%        [H(k-2) H(k-3) ... H(0)  0  ]
%        [H(k-1) H(k-2) ... H(1) H(0)]
%
%   is a least squares lower triangular block Toeplitz matrix solution to 
%   S = T*P. The matrix H = [H(0);H(1);...;H(k-1)].
%
%   The number of columns in S and P have to be the same and the number of
%   rows of S and P have to be divisible by k.
%
%   See also BLOCTOEP, TOEPLITZ.
%
% 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.

% See Section 6 of N.L.C. Chui and J.M. Maciejowski. "Subspace
% identification - a Markov parameter approach". Technical Report
% CUED/F-INFENG/TR.337, December 1998, University of Cambridge, UK.

% SS = [S1;S2;...;Sk]
% PP = [P1;P2;...;Pk]

p = size(SS,1)/k; % number of rows in each Si
m = size(PP,1)/k; % number of rows in each Pi
n = size(SS,2);   % number of columns in SS and PP

if ~(fix(p) == p & fix(m) == m)
  error('The number of rows of S are P are not divisible by k.')
end

if ~(n == size(PP,2))
  error('The number of columns of S and P differ.')
end

% Rearrange SS and PP such that
% S = [S1 S2 ...  Sk ]  and
% P = [P1 P2 ...  Pk ]  is an upper triangular block Toeplitz matrix
%     [0  P1 ... Pk-1] 
%     [:           : ]
%     [0  0  ...  P1 ]

S = zeros(p,n);
for i = 1:k
  S(:,(i-1)*n+1:i*n) = SS((i-1)*p+1:i*p,:); 
end
clear SS
P = bloctoep([PP(1:m,:);zeros(m*k-m,n)],PP,k);
clear PP

% Solve for HH in S = HH*P
HH = S/P; % HH = [H0 H2 ... Hk-1]

% Reshape to get H = [H0;H2;...;Hk-1] and T
H = zeros(k*p,m);
for i = 1:k
  H((i-1)*p+1:i*p,:) = HH(:,(i-1)*m+1:i*m);
end
clear HH

T = bloctoep(H,[H(1:p,:); zeros(p*k-p,m)],k);

% *** last line of soltritoep.m ***