sfact

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)

See also

• gtild — tilde operation
• fspecg — stable factorization of continuous time dynamical systems

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