Ornstein-Uhlenbeck power spectrum (source code)

= Ornstein-Uhlenbeck power spectrum
{c}

Let $X$ be the centered stationary solution of $dX_t=-AX_t\,dt+B\,dW_t$, where $A$ is a real $d\times d$ <matrix> with eigenvalues of positive real part and $W$ is an $m$-dimensional <Brownian motion>. Use the two-sided angular-frequency convention
$$
S(\omega)=\int_{\mathbb R}e^{-i\omega t}\mathbb E[X_tX_0^T]dt,
\qquad \Sigma=\frac1{2\pi}\int_{\mathbb R}S(\omega)d\omega.
$$
The <Ornstein-Uhlenbeck power spectrum> is
$$
S(\omega)=(A+i\omega I)^{-1}BB^T(A^T-i\omega I)^{-1}.
$$
The stationary linear filter has transfer <matrix> $H(\omega)=(A+i\omega I)^{-1}$, so the output spectral <matrix> is $HBB^TH^*$, explaining the adjoint and sign of the second factor. This <matrix> is Hermitian positive semidefinite for every frequency.

For scalar $dX=-\lambda X\,dt+\sigma\,dW$ with $\lambda>0$, the spectrum is $\sigma^2/(\lambda^2+\omega^2)$, proportional to a <Cauchy distribution> <probability density function>. Its inverse-frequency <integral> gives <variance> $\sigma^2/(2\lambda)$, checking the normalization. With ordinary frequency $f$ and inverse transform $\int e^{2\pi ift}S_f(f)df$, use $S_f(f)=S(2\pi f)$. A one-sided positive-frequency spectrum doubles this density away from zero. These conventions must be matched when comparing a theoretical spectrum with measured spectral power.