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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
|
// Copyright (C) 2018 - IIT Bombay - FOSSEE
//
// 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-en.txt
// Author:[insert name]
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
function [multp, indx] = mpoles (p, tol, reorder)
//Calling sequence:
// [multp, idxp] = mpoles (p)
// [multp, idxp] = mpoles (p, tol)
// [multp, idxp] = mpoles (p, tol, reorder)
// Identify unique poles in p and their associated multiplicity.
//
// The output is ordered from largest pole to smallest pole.
//
// If the relative difference of two poles is less than tol then they are
// considered to be multiples. The default value for tol is 0.001.
//
// If the optional parameter reorder is zero, poles are not sorted.
//
// The output multp is a vector specifying the multiplicity of the poles.
// multp(n) refers to the multiplicity of the Nth pole
// p(idxp(n)).
//
//
//
// test case:
// p = [2 3 1 1 2];
// [m, n] = mpoles (p)
// n =
// 2.
// 5.
// 1.
// 4.
// 3.
// m =
// 1.
// 1.
// 2.
// 1.
// 2.
//
[nargout,nargin]=argn();
if (nargin < 1 | nargin > 3)
error("wrong number of input arguments");
end
if (nargin < 2 | isempty (tol))
tol = 0.001;
end
if (nargin < 3 | isempty (reorder))
reorder = %t;
end
Np = length (p);
// Force the poles to be a column vector.
p = p(:);
// Sort the poles according to their magnitidues, largest first.
if (reorder)
// Sort with smallest magnitude first.
[p, ordr] = gsort (p,"r","i");
// Reverse order, largest maginitude first.
n = Np:-1:1;
p = p(n);
ordr = ordr(n);
else
ordr = 1:Np;
end
// Find pole multiplicty by comparing the relative differnce in the
// poles.
multp = zeros (Np, 1);
indx = [];
n = find (multp == 0, 1);
while (n)
dp = abs (p-p(n));
if (p(n) == 0.0)
if (or (abs (p) > 0 & (abs(p)<%inf)))
p0 = mean (abs (p(abs (p) > 0 & abs(p)<%inf)));
else
p0 = 1;
end
else
p0 = abs (p(n));
end
k = find (dp < tol * p0)';
// Poles can only be members of one multiplicity group.
// if (length(indx))
// mk=members(k,indx);
// k = (~ bool2s(mk~=0));
// end
m = 1:length (k);
multp(k) = m';
indx = [indx; k];
n = find (multp == 0, 1);
end
multp = multp(indx);
indx = ordr(indx);
endfunction
|
