How to calculate the PSD of a stochastic process

I'll reply myself with a full derivation of the PSD of the harmonic oscillator. It's only half of the answer to my original question; the other half can be found here, after I asked the same question on mathoverflow

The full equation takes the expression $$ m \ddot{x} + \gamma \dot{x} + \delta x = \sigma \eta(t) $$ Performing a change of variables, $x = e^{\frac{-\gamma}{2m}t}x_1$, we get $$ m \ddot{x}_1 + \left(\delta - \frac{\gamma^2}{4m}\right)x_1 = \sigma e^{\frac{\gamma}{2m}t}\eta(t) $$ thus eliminating $\dot{x}$. Setting $a = \frac{\delta}{m} - \frac{\gamma^2}{4m^2}$, $b = \frac{\sigma}{m}$ and rewriting the equation as a first order linear system with $$ X = \begin{pmatrix} x_2\\ v_2 \end{pmatrix} $$ where $v_2 = \dot{x_2}$, we get in Ito's notation \begin{equation} X = \begin{pmatrix}\label{eq:complete} 0 & 1 \\ -a & 0 \end{pmatrix} \cdot X \mathrm{d}t + \begin{pmatrix} 0\\ b e^{\frac{\gamma}{2 m}t} \end{pmatrix} \cdot \mathrm{d} W_t \end{equation} The solution of a linear homogeneous SDE is $$ X_t = e^{\int_0^t A(t)\mathrm{d}t}\cdot X_0 + e^{\int_0^t A(t)\mathrm{d}t}\cdot \int_0^t e^{-\int A(s)\mathrm{d}s}\sigma(s)\mathrm{d}W_s $$ where $A(t)$ is the (generally vector) coefficient of $X$. For this SDE a fundamental matrix solution of the associated homogeneous noise-free system is $$ \Phi(t) = \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \\ -\sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} $$ The determinant of this matrix is 1, so its inverse matrix will be $$ \Phi^{-1}(t) = e^{-\int A(\tau)\mathrm{d}\tau} = \det \Phi(t)^{-1}\cdot \begin{pmatrix} \cos \sqrt{a}t & -\sin \sqrt{a}t/\sqrt{a} \\ \sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} = \begin{pmatrix} \cos \sqrt{a}t & -\sin \sqrt{a}t/\sqrt{a} \\ \sqrt{a} \sin\sqrt{a}t & \cos \sqrt{a}t \end{pmatrix} $$ and hence we can solve the complete system. We are interested in the first component of $X$, the position (we will only calculate the PSD of $x$, although this can be done considering the full matrix system)

\begin{equation*} \begin{split} x_1(t) = & \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \begin{pmatrix} x_1(0)\\ v_1(0) \end{pmatrix} + \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \int_0^{t} b e^{\frac{\gamma}{2m}r}\cdot \begin{pmatrix} -\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r \end{pmatrix} \mathrm{d}W_r \end{split} \end{equation*} Finally, $x_1(t) = e^{\frac{\gamma t}{2m}}x(t)$, so \begin{equation} \begin{split} x(t) = & e^{-\frac{\gamma t}{2m}}\begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \begin{pmatrix} x(0)\\ v(0) + \frac{\gamma}{2m}x(0) \end{pmatrix} +\\ & e^{-\frac{\gamma t}{2m}} \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \int_0^{t} b e^{\frac{\gamma}{2m}r}\cdot \begin{pmatrix} -\sin \sqrt{a}r/\sqrt{a}\\ \cos \sqrt{a}r \end{pmatrix} \mathrm{d} W_r \end{split} \end{equation}

We see that, after a transient time, only the term depending on $\mathrm{d} W_r$ remains, so the first moment of the process is zero. Now, applying Ito's isometry to calculate the covariance we get $$ \mathbb{E} \left[ \int_0^{T} f(u)\,\mathrm{d} W_u \int_0^S f(v) \,\mathrm{d} W_v \right] = $$ $$ b^2 e^{-\frac{\gamma (t + s)}{2m}} \begin{pmatrix} \cos \sqrt{a}t & \sin \sqrt{a}t/\sqrt{a} \end{pmatrix} \cdot \mathbb{E} \left[ \int_0^{\min(t,s)} e^{\frac{\gamma}{m}u} \begin{pmatrix} \frac{\sin^2 \sqrt{a}u}{a} & -\frac{\sin \sqrt{a}u \cos \sqrt{a}u}{\sqrt{a}}\\ -\frac{\sin \sqrt{a}u \cos \sqrt{a}u}{\sqrt{a}} & \cos^2 \sqrt{a}u \end{pmatrix} \, \mathrm{d} u \right] \cdot \begin{pmatrix} \cos \sqrt{a}s \\ \sin \sqrt{a}s/\sqrt{a} \end{pmatrix} $$

This is a quite uninteresting calculation. After simplification one gets $$ R(\tau) = \frac{b^2 m^2 e^{\frac{-\Gamma |\tau|}{2}}\left( 2 \sqrt{a}m \cos(\sqrt{a}|\tau|) + \gamma \sin(\sqrt{a}|\tau|) \right)}{\sqrt{a}(\gamma^3 + 4a\gamma m^2)} $$ where we have defined the normalized damping constant $\Gamma = \frac{\gamma}{m}$.

From the expression of the autocorrelation we see that $R(t,\tau) = R(\tau)$: therefore, the process is wide-sense stationary and the conditions to apply the Wiener-Khinchin theorem are satisfied. The Fourier transform of this autocorrelation function is the power spectral density $$ S(f) = \frac{16b^2}{\Gamma^4 + 8(a + 4\pi^2f^2)\Gamma^2 + 16(a-4\pi^2f^2)^2} $$ which, after replacing the variables and some rearranging takes the simpler expression $$ S(\omega) = \frac{\sigma^2/m^2}{(\omega_0^2 - \omega^2)^2 + \Gamma^2 \omega^2} $$ where we have replaced the unitary ordinary frequency Fourier transform (in terms of $f$) by the non-unitary angular frequency Fourier transform.