dsigma.m

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

function [svout,w] = dsigma(Z1,Z2,Z3,Z4,Z5,Z6,Z7)
%function [svout,w] = dsigma(a,b,c,d,Ts,w,invflag)
%DSIGMA Singular value frequency response of discrete-time linear systems.
%	DSIGMA(A,B,C,D,Ts) (or optional SIGMA(SS_,Ts) in RCT) produces a 
%       singular value plot of matrix:   -1
%                    G(w) = C(exp(jwT)I-A) B + D 
%	as a function of frequency.  The singular values are an extension
%	of Bode magnitude response for MIMO system.  The frequency range 
%       and number of points are chosen automatically. For square systems, 
%       DSIGMA(A,B,C,D,Ts,'inv') produces the singular values of the inverse 
%       matrix     -1                    -1      -1
%	          G (w) = [ C(exp(jwT)I-A) B + D ]
%	DSIGMA(A,B,C,D,Ts,W) or DSIGMA(A,B,C,D,Ts,W,'inv') use the user-
%	supplied frequency vector W which must contain the frequencies, in
%	radians/sec, at which the singular value response is to be 
%	evaluated. Aliasing will occur at frequencies greater than the 
%	Nyquist frequency (pi/Ts).  When invoked with left hand arguments,
%           [SV,W] = DSIGMA(A,B,C,D,Ts,...)
%	or  [SV,W] = SIGMA(SS_,Ts,...)    (for Robust Control Toolbox user)
%	returns the frequency vector W and the matrix SV with MIN(NU,NY)
%	columns and length(W) rows, where NU is the number of inputs and
%	NY is the number of outputs.  No plot is drawn on the screen. 
%	The singular values are returned in descending order.
%
%	See also: LOGSPACE, SEMILOGX, DNICHOLS, DNYQUIST and DBODE.

%	Clay M. Thompson  7-10-90
%	Revised A.Grace 2-12-91, 6-21-92
%	Revised W.Wang 7/24/92
%	Copyright (c) 1986-93 by the MathWorks, Inc.

if nargin==0, eval('dexresp(''dsigma'',1)'), return, end

if nargin <= 4  %if not enough parameters, check rct data structure
  cmd=[];
  %The following if..end is add to match up RCT format change
  if exist('mkargs') == 2, %If RCT installed
    inargs='(a,b,c,d,Ts,w,invflag)';
    eval('cmd=mkargs(inargs,nargin,''ss'');')
    eval(cmd);
  end;
else                       % If RCT not installed
  % assign a=Z1;...,w=Z6;invflag=Z7;
  a=Z1; b=Z2; c=Z3; d=Z4;
  if nargin>6,
    invflag=Z7; w=Z6; Ts=Z5;
  elseif nargin>5,
    w=Z6; Ts=Z5;
  elseif nargin>4,
    Ts=Z5;
  end
end;

if nargin==7,		% Trap call to RCT function
	if ~isstr(invflag),
		eval('svout = dsigma2(a,b,c,d,Ts,w,invflag);')
		return
	end
end

error(nargchk(5,7,nargin));
error(abcdchk(a,b,c,d));
if ~(length(d) | (length(b) &  length(c)))
	return;
end

% Determine status of invflag
if nargin==5, 
	invflag = [];
	w = [];
elseif (nargin==6)
	if (isstr(w)),
		invflag = w;
		w = [];
		[ny,nu] = size(d);
		if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
	else
		invflag = [];
	end

else
	[ny,nu] = size(d);
	if (ny~=nu), error('The state space system must be square when using ''inv''.'); end
end

% Generate frequency range if one is not specified.

% If frequency vector supplied then use Auto-selection algorithm
% Fifth argument determines precision of the plot.
if ~length(w)  
	w=dfrqint(a,b,c,d,Ts,30);
end

[nx,na] = size(a);
[no,ns] = size(c);
nw = max(size(w));

% Balance A
[t,a] = balance(a);
b = t \ b;
c = c * t;

% Reduce A to Hessenberg form:
[p,a] = hess(a);

% Apply similarity transformations from Hessenberg
% reduction to B and C:
b = p' * b;
c = c * p;

s = exp(sqrt(-1)*w*Ts);
I=eye(length(a));
if nx > 0,
	for i=1:length(s)
		if isempty(invflag),
		    sv(:,i)=svd(c*((s(i)*I-a)\b) + d);
		else
		    sv(:,i)=svd(inv(c*((s(i)*I-a)\b) + d));
		end
	end
else
	for i=1:length(s)
		if isempty(invflag),
		    sv(:,i)=svd(d);
		else
		    sv(:,i)=svd(inv(d));
		end
	end
end

% If no left hand arguments then plot graph.
if nargout==0
	semilogx(w, 20*log10(sv), [w(1), w(length(w))], [0,0], 'w:')
	xlabel('Frequency (rad/sec)')
	ylabel('Singular Values dB')
	return % Suppress output
end
svout =sv;