Spectral factorization of linear systems
H = G'*S*G
H = H'. In this factorization,
Sis a symmetric matrix and
Gis a square, stable, and minimum-phase system with unit (identity) feedthrough.
G'is the conjugate of
G, which has transfer function G(–s)T in continuous time, and G(1/z)T in discrete time.
Consider the following system.
G0 = ss(zpk([-1 -5 1+2i 1-2i],[-100 1+2i 1-2i -10],1e3)); H = G0'*G0;
G0 has a mix of stable and unstable dynamics.
H is a self-conjugate system whose dynamics consist of the poles and zeros of
G0 and their reflections across the imaginary axis. Use spectral factorization to separate the stable poles and zeros into
G and the unstable poles and zeros into
[G,S] = spectralfact(H);
G is stable and minimum phase, by checking that all its poles and zeros fall in the left half-plane (
Re(s) < 0).
p = pole(G)
p = 4×1 complex 102 × -0.0100 + 0.0200i -0.0100 - 0.0200i -0.1000 + 0.0000i -1.0000 + 0.0000i
z = zero(G)
z = 4×1 complex -1.0000 + 2.0000i -1.0000 - 2.0000i -1.0000 + 0.0000i -5.0000 + 0.0000i
G also has unit feedthrough.
ans = 1
H is SISO,
S is a scalar. If
H were MIMO, the dimensions of
S would match the I/O dimensions of
S = 1000000
H = G'*S*G by comparing the original system to the difference between the original and factored systems.
sigmaplot throws a warning because the difference is very small.
Hf = G'*S*G; sigmaplot(H,H-Hf)
Warning: The frequency response has poor relative accuracy. This may be because the response is nearly zero or infinite at all frequencies, or because the state-space realization is ill conditioned. Use the "prescale" command to investigate further.
Suppose that you have the following 2-output, 2-input state-space model,
A = [-1.1 0.37; 0.37 -0.95]; B = [0.72 0.71; 0 -0.20]; C = [0.12 1.40 1.49 1.41]; D = [0.67 0.7172; -1.2 0]; F = ss(A,B,C,D);
Suppose further that you have a symmetric 2-by-2 matrix,
R = [0.65 0.61 0.61 -3.42];
Compute the spectral factorization of the system given by
H = F'*R*F, without explicitly computing
[G,S] = spectralfact(F,R);
G is a minimum-phase system with identity feedthrough.
ans = 2×2 1 0 0 1
F is has two inputs and two outputs, both
S are 2-by-2 matrices.
G'*S*G = F'*R*F by comparing the original factorization to the difference between the two factorizations. The singular values of the difference are far below those of the original system.
Ff = F'*R*F; Gf = G'*S*G; sigmaplot(Ff,Ff-Gf)
Consider the following discrete-time system.
F = zpk(-1.76,[-1+i -1-i],-4,0.002);
F has poles and zeros outside the unit circle. Use
spectralfact to compute a system
G with stable poles and zeros, such that
G'*G = F'*F.
G = spectralfact(F,)
G = -3.52 z (z+0.5682) ------------------ (z^2 + z + 0.5) Sample time: 0.002 seconds Discrete-time zero/pole/gain model.
G has no poles or zeroes outside the unit circle.
G does have an additional zero at
z = 0, which is a reflection of the unstable zero at
z = Inf in
G'*G = F'*F by comparing the original factorization to the difference between the two factorizations. The singular values of the difference are far below those of the original factorization.
Ff = F'*F; Gf = G'*G; sigmaplot(Ff,Ff-Gf)
H— Self-conjugate LTI model
Self-conjugate LTI model, specified as a
zpk model. Self-conjugate means that is equal
to its conjugate,
H' is the transfer function H(–s)T in
continuous time and H(1/z)T in
H can be SISO or MIMO, provided it has
as many outputs as inputs.
H can be continuous
or discrete with the following restrictions:
In continuous time,
H must be
biproper with no poles or zeros at infinity or on the imaginary axis.
In discrete time,
H must have
no poles or zeros on the unit circle.
F factor of the factored form
= F'*R*F, specified as a
F cannot have
more inputs than outputs.
G— LTI factor
LTI factor, returned as a
G is a stable,
minimum-phase system that satisfies:
H = G'*S*G, if you use the syntax
G'*S*G = F'*R*F, if you use the
[G,S] = spectralfact(F,R).
G'*G = F'*F, if you use the syntax
S— Numeric factor
Numeric factor, returned as a symmetric matrix that satisfies:
spectralfact assumes that
self-conjugate. In some cases when
H is not self-conjugate,
do not satisfy
H = G'*S*G. Therefore, verify that
your input model is in fact self-conjugate before using
One way to verify
H is to compare
- H' on a singular value plot.
H is self-conjugate, the
- H' line on the plot lies far below the