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
125
126
127
128
|
// 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
// Original Source : https://octave.sourceforge.io/
// Modifieded by: Abinash Singh Under FOSSEE Internship
// Last Modified on : 3 Feb 2024
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
function [multp, idxp] = mpoles (p, tol, reorder)
if (nargin < 1)
error("mpoles: Invalid number of input arguments");
end
if ~( type(p)== 1)
error ("mpoles: P must be a single or double floating point vector");
end
if (nargin < 2 || isempty (tol))
tol = 0.001;
elseif (~(isscalar (tol) && isreal (tol) && tol > 0))
error ("mpoles: TOL must be a real scalar greater than 0");
end
if (nargin < 3 || isempty (reorder))
reorder = %t;
elseif (~(isscalar (reorder) && isreal (reorder)))
error ("mpoles: REORDER must be a numeric or logical scalar");
end
Np = length (p);
p = p(:); // force poles to be a column vector
if (reorder)
//// sort with largest magnitude first
[_, order] = gsort (abs(p));
p = p(order);
else
order = (1:Np).';
end
//// Create vector of tolerances for use in algorithm.
vtol = zeros (Np, 1, typeof(p));
p_nz = (p ~= 0); // non-zero poles
vtol(~p_nz) = tol; // use absolute tolerance for zero poles
//// Find pole multiplicity by comparing relative difference of poles.
multp = zeros (Np, 1, typeof(p));
idxp = [];
n = find (multp == 0, 1);
while (n)
dp = abs (p - p(n));
vtol(p_nz) = tol * abs (p(n));
k = find (dp < vtol);
//// Poles can only be members of one multiplicity group.
if (length (idxp))
k = k(~ismember (k, idxp));
end
m = 1:length (k);
multp(k) = m;
// disp("k")
// disp(k)
// disp("idxp")
// disp(idxp)
idxp = [idxp; k(:)];
n = find (multp == 0, 1);
end
multp = multp(idxp);
idxp = order(idxp);
endfunction
function y = ismember(a, b) // FIXME : do this in a more efficient .
y = zeros(size(a,1),size(a,2));
for i = 1:length(a)
for j = 1:length(b)
if a(i) == b(j)
y(i) = 1;
break;
end
end
end
endfunction
/*
// test
[mp, ip] = mpoles ([0 0], 0.01);
assert_checkequal (mp, [1; 2]);
// test
[mp, ip] = mpoles ([-1e4, -0.1, 0]);
assert_checkequal (mp, [1; 1; 1]);
assert_checkequal (ip, [1; 2; 3]);
// Test single inputs
// test
[mp, ip] = mpoles ([-1e4, -0.1, 0]);
assert_checkequal (mp,[1; 1; 1]);
assert_checkequal (ip, [1; 2; 3]);
// Test relative tolerance criteria
// test
[mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1);
assert_checkequal (mp, [1; 1; 1]);
[mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1 + %eps);
assert_checkequal (mp, [1; 1; 2]);
// Test absolute tolerance criteria with a zero pole
// test
[mp, ip] = mpoles ([0, -0.1, 0.3], .1);
assert_checkequal (mp, [1; 1; 1]);
[mp, ip] = mpoles ([0, -0.1, 0.3], .1 + %eps);
assert_checkequal (mp, [1; 1; 2]);
//// Test input validation
error <Invalid call> mpoles ()
error <P must be a single or double floating point vector> mpoles (uint8 (1))
error <TOL must be a real scalar greater than 0> mpoles (1, [1, 2])
error <TOL must be a real scalar greater than 0> mpoles (1, 1i)
error <TOL must be a real scalar greater than 0> mpoles (1, 0)
error <REORDER must be a numeric or logical scalar> mpoles (1, 1, [1, 2])
error <REORDER must be a numeric or logical scalar> mpoles (1, 1, {1})
*/
|
