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.
