[Scilab-users] Simplify with threshold
Federico Miyara
fmiyara at fceia.unr.edu.ar
Fri Aug 28 18:12:06 CEST 2020
Sorry, the previous version yielded a normalized version of the
simplified rational. This one fixes that.
Federico Miyara
--
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200828/95f7f5eb/attachment.htm>
-------------- next part --------------
function [Hsimp, p, z] = simplify(H, th)
// Simplify rational function when poles and zeros are
// found whose values differ less than a given numerical
// threshold
//
// Usage:
// [Hsimp, p, z] = simplify(H, th)
// where
// H: Rational function
// th: Threshold; if omitted its default value is 1e-6
// Hsimp: Version of H with coincident poles and zeros
// removed
// p: Poles after simplification
// z: Zeros after simplification
//
// NOTE: The simplification is imperfect if poles and zeros
// are not identical, but the rational function computed
// at specific values of the variable will differ very
// little if the threshold is carefully selected. The
// threshold should be small, otherwise poles and zeros
// which differ considerably might pe cancelled out.
//
// ----------------------------
// Author: Federico Miyara
// Date: 2020-08-27
// Input argument error handling
if argn(2)<1 | argn(2)>2
error(msprintf(gettext(..
"%s: Wrong number of input arguments: %d or %d expected.\n"),..
"simplify",1,2));
end
if argn(2)<2
th = 0.000001;
end
if type(th)<>1 | (type(th)==1 & ~isreal(th))
error(msprintf(gettext(..
"%s: Argument #%d: Decimal number expected.\n"), ..
"simplify", 2))
end
if type(H)<>16
error(msprintf(gettext(..
"%s: Argument #%d: Rational fraction expected.\n"), ..
"simplify", 1))
end
if argn(1)>3
error(msprintf(gettext(..
"%s: Wrong number of output arguments: %d to %d expected.\n"),..
"simplify",0,3));
end
// Zeros and poles of H
z = roots(H.num);
p = roots(H.den);
// Compute the differences between all poles and all
// zeros. Find the pole and zero whose difference is
// minimum. If that difference is smaller than the
// threshold, remove that pole anz zero. Iterate
// until no pole and zero differing less than the
// threshold are left.
terminate = 0;
while terminate==0
// Number of zeros and poles
nz = length(z);
np = length(p);
// Initialize matrix of differences
delta = zeros(np,nz);
// Compute matrix of differences
for k=1:nz
for n=1:np
delta(n,k) = abs(z(k) - p(n));
end
end
// Find the minimum differenc
[mini, I] = min(delta);
// If the minimum is smaller than the threshold
// remove the corresponding pole and zero.
// Otherwise, terminate the iteration.
if mini<th
p(I(1)) = [];
z(I(2)) = [];
else
terminate = 1;
end
end
// Recreate simplified numerator and denominator
Nsimp = poly(z, varn(H.num));
Dsimp = poly(p, varn(H.den));
// Simplified version of rational (normalized)
Hsimp = Nsimp/Dsimp;
// Find constant to denormalize
xo = max(abs([z; p])) + 1;
K = horner(H,xo)/horner(Hsimp,xo);
// Final simplified version of rational
Hsimp = K*Hsimp;
endfunction
More information about the users
mailing list