sfact
discrete time spectral factorization
Syntax
F = sfact(P)
Arguments
- P
Square real polynomial matrix
Description
If P
is a real polynomial matrix verifying
P(z) / z^m = P(1/z)' * z^m
where n=2*m
is the maximum degree (necessarily even) of
P(z)
elements and P(z)/z^n
is positive definite on the unit circle then
sfact(P)
yields a polynomial matrix F
such that
P(z) = F(z) * F(1/z)' * z^n
By this factorization the roots of det(P) are splitted in two sets strictly included in the exterior or the interior
of the unit circle or its exterior, the (antistable) roots of det(F)
or the (stable) roots of det(F(1/z)'*z^n)
respectively.
In the scalar case the first condition simplifies, as P
has to be palindromic (or self-reciprocal)
i.e. the coefficients of the monomials of degrees k
and n-k
must be equal.
In the contect of continous time (see the third example below), the factorization can be applied to a real polynomial
matrix P(s)
verifying
P(s) = P(-s)
and P(s)
positive definite on the imaginary axis by using the
continuous to discrete time bilinear transformation.
The spectral factorization has applications in system theory, optimal control, network theory and control system synthesis. The details of the algorithm can be found in the referenced paper.
Examples
// single polynomial example z = %z; p = (z -1/2) * (2 - z) // p is palindromic w = sfact(p); w * horner(w, 1/z) * z
// polynomial matrix example z = %z; // construct P verifying the conditions // gtild(F1,'d') is horner(F1, 1/z)'*z^3 F1 = [z-1/2, z+1/2, z^2+2; 1, z, -z; z^3+2*z, z, 1/2-z]; P = F1*gtild(F1,'d'); F = sfact(P) roots(det(P)) // roots of det(P) roots(det(gtild(F,'d'))) // stable roots roots(det(F)) // antistable roots clean(P-F*gtild(F,'d'))
// Example of continuous time use s = %s; z = %z; p = -3*(s+(1+%i))*(s+(1-%i))*(s+0.5)*(s-0.5)*(s-(1+%i))*(s-(1-%i)); p = real(p); // p(s) = polynomial in s^2 , looks for stable f such that p=f(s)*f(-s) w = horner(p,(1-z)/(1+z)); // bilinear transform w=p((1-z)/(1+z)) wn = w.num; // take the numerator fn = sfact(wn); // factor f = horner(fn,(1-s)/(1+s)).num; // back transform f = f/sqrt(horner(f*horner(f,-s),0)); f = f*sqrt(horner(p,0)); // normalization roots(f) // f is stable clean(f*horner(f,-s)-p) // f(s)*f(-s) is p(s)
References
J. Ježek, V. Kučera, Efficient Algorithm for Spectral Factorization, IFAC Proceedings Volumes, Volume 17, Issue 2, 1984, Pages 257-262, ISSN 1474-6670, https://doi.org/10.1016/S1474-6670(17)60979-0
Report an issue | ||
<< rowcompr | Polynomials | simp >> |