signextb.m

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

function [sys,Bind,par,T] = signextb(sys0) 
%SIGNEXTB   Corrects the signs and extracts the elements of the B matrix.
%
%   [SYS,Bind,par,T] = SIGNEXTB(SYS0) corrects the signs of the B elements in
%   the balanced system SYS0 and outputs the new balanced system in SYS. The
%   diagonal matrix T == inv(T) is the state transformation z = Tx used to
%   convert SYS0 to SYS. Bind(j) contains the index of the first non-zero
%   element in the j'th row of SYS.B. This positive element and the
%   remaining elements in the row are contained in par. 
%
%   See also FILLB.
%
% 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.

sys = sys0; % Copy all fields
B = sys0.b;
[n,m] = size(B);
T = ones(n,1);
Bind = T;
par = [];

for j = 1:n 
  if all(B(j,:) == 0)
	error('One of the rows of B (K if ALG=''mp'') == 0. This is not allowed.')
  end
  % Find the first non-zero element in the row
  k = min(find(B(j,:)));
  Bind(j) = k;
  % If the first non-zero element is negative then change the
  % corresponding sign of the element in T
  if sign(B(j,k)) == -1
	T(j) = -1;
  end
end

T = diag(T); % Note that inv(T) == T

if any(any(T == -1)) % If any of the signs of B needs to change
  sys = ss2ss(sys0,T);
%  sys.a = T  * sys0.a * T; 
%  sys.b = T  * B;
%  sys.c = sys0.c * T;
end

% Extract data from B row-wise and put into par
for j = 1:n
  k = Bind(j);
  par = [par; log(sys.b(j,k))]; % Guarantees that the first non-zero element is
                                % positive. Extracted using exp(.) - See FILLB
  for i = k+1:m
	par = [par; sys.b(j,i)];    % Any sign allowed
  end
end