Solution (source code)

= Solution

Split the Fourier integral at zero and use the two stationary covariance branches:
$$
s(\omega)=\int_0^\infty e^{-(A+i\omega I)\tau}\sigma\,d\tau
+\int_0^\infty e^{-(A-i\omega I)\tau}\!{}^T\sigma\,d\tau.
$$
Equivalently, Fourier transforming the <Multivariate Ornstein-Uhlenbeck process> equation gives
$$
(A+i\omega I)x(\omega)=b\Lambda(\omega).
$$
Unit white-noise covariance then yields the <Ornstein-Uhlenbeck power spectrum>
$$
s(\omega)=(A+i\omega I)^{-1}bb^T(A^T-i\omega I)^{-1},
$$
so
$$
\boxed{(A+i\omega I)s(\omega)(A^T-i\omega I)=bb^T}.
$$

Solved by gpt-5.6-sol high.