tzero2.m

数字通信第四版原书的例程

function z = tzero(a,b,c,d)
%TZERO2	Transmission zeros.
%	TZERO2(A,B,C,D) returns the transmission zeros of the state-space
%	system:
%		.
%		x = Ax + Bu
%		y = Cx + Du
%
%	Care must be exercised to disregard any large zeros that may
%	actually be at infinity.  It is best to use FORMAT SHORT E so
%	large zeros don't swamp out the display of smaller finite zeros.
%
%	This M-file finds the transmission zeros using an algorithm based
%	on QZ techniques.  It is not as numerically reliable as the 
%	algorithm in TZERO.
%
%	See also: TZERO.

% For more information on this algorithm, see:
%   [1] A.J Laub and B.C. Moore, Calculation of Transmission Zeros Using
%   QZ Techniques, Automatica, 14, 1978, p557
%
% For a better algorithm, not implemented here, see:
%   [2] A. Emami-Naeini and P. Van Dooren, Computation of Zeros of Linear 
%   Multivariable Systems, Automatica, Vol. 14, No. 4, pp. 415-430, 1982.
%

%	J.N. Little 4-21-85
%	Copyright (c) 1986-93 by the MathWorks, Inc.

error(nargchk(4,4,nargin));
error(abcdchk(a,b,c,d));

[n,m] = size(b);
[r,n] = size(c);
aa = [a b;c d];

if m == r
	[mcols, nrows] = size(aa);
	bb = zeros(mcols, nrows);
	bb(1:n,1:n) = eye(n);
	z = eig(aa,bb);
	z = z(~isnan(z) & finite(z)); % Punch out NaN's and Inf's
else
	nrm = norm(aa,1);
	if m > r
		aa1 = [aa;nrm*(rand(m-r,n+m)-.5)];
		aa2 = [aa;nrm*(rand(m-r,n+m)-.5)];
	else
		aa1 = [aa nrm*(rand(n+r,r-m)-.5)];
		aa2 = [aa nrm*(rand(n+r,r-m)-.5)];
	end
	[mcols, nrows] = size(aa1);
	bb = zeros(mcols, nrows);
	bb(1:n,1:n) = eye(n);
	z1 = eig(aa1,bb);
	z2 = eig(aa2,bb);
	z1 = z1(~isnan(z1) & finite(z1));   % Punch out NaN's and Inf's
	z2 = z2(~isnan(z2) & finite(z2));
	nz = length(z1);
	for i=1:nz
		if any(abs(z1(i)-z2) < nrm*sqrt(eps))
			z = [z;z1(i)];
		end
	end
end