1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function svm = svplot(Sl,w)
//svplot singular-value sigma-plot.
// svm = svplot(sl,w) computes for the linear dynamical system
// sl, the singular values of its transfer function matrix:
// -1
// g(jw) = c(jw*i-a) b+d
//
// or
// -1
// g(exp(jw)) = c(exp(jw)*i-a) b+d
//
// evaluated over the frequency range specified by w.
// sl is a sylin list (see syslin) representing the system
// [a,b,c,d] in state-space form.
// the i-th column of the output matrix svm contains the singular
// values of g(exp(jw)) for the i-th frequency value.
// svm = svplot(sl) is equivalent to
// svm = svplot(sl,logspace(-3,3)) (continuous) or
// svm = svplot(sl,logspace(-3,pi)) (discrete).
//!
[a,b,c,d]=abcd(Sl);
// Reduce a to Hessenberg form
[q,a] = hess(a); b = q'*b; c = c*q;
// Compute the singular values of the frequency response
select Sl.dt
case []
warning(msprintf(gettext("%s: Input argument #%d is assumed continuous time.\n"),"svplot",1));
if argn(2) == 1
w = logspace(-3,3);
end
nf = max(size(w)); nsv = min(size(d)); j = sqrt(-1);
svm(nsv,nf) = 0;
for i = 1:nf
svm(:,i) = svd(c*((j*w(i)*eye()-a)\b)+d);
end
case "c"
if argn(2) == 1
w = logspace(-3,3);
end
nf = max(size(w)); nsv = min(size(d)); j = sqrt(-1);
svm(nsv,nf) = 0;
for i = 1:nf
svm(:,i) = svd(c*((j*w(i)*eye()-a)\b)+d);
end
case "d"
if argn(2) == 1
w = logspace(-3,%pi);
end
nf = max(size(w)); nsv = min(size(d)); j = sqrt(-1);
svm(nsv,nf) = 0;
for i = 1:nf
svm(:,i) = svd(c*((exp(j*w(i))*eye()-a)\b)+d);
end
else T=Sl("dt");
if argn(2) == 1
w = logspace(-3,%pi);
end
nf = max(size(w)); nsv = min(size(d)); j = sqrt(-1);
svm(nsv,nf) = 0;
for i = 1:nf
svm(:,i) = svd(c*((exp(j*w(i)*T)*eye()-a)\b)+d);
end
end
endfunction
|
