# 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)```