An adapted and enhanced version of
an answer of mine on dsp.SE
Suppose that
$$\begin{align*}
Y(t) &= \int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds \tag{1}
\end{align*}
$$
If we define the crosscorrelation function $R_{X,Y}(\tau)$ as
$$R_{X,Y}(\tau) = E[X(t-\tau)Y(t)],\tag{2}$$ then
$$\begin{align*}
R_{X,Y}(\tau) &= E\left[X(t-\tau)\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right]\\
&= \int_{-\infty}^{\infty} h(s)E[X(t-\tau)X(t-s)]\,\mathrm ds\\
&= \int_{-\infty}^{\infty} h(s)R_X(\tau-s)\,\mathrm ds.
\end{align*}$$
In short,
$R_{X,Y} = h\star R_X.$
Next, consider
$$\begin{align}
R_Y(\tau) &= E[Y(t-\tau)Y(t)]\\
&= E\left[\int_{-\infty}^{\infty} h(s)X(t-\tau-s)\,\mathrm ds \,Y(t)\right]
&{\scriptstyle{\text{substituting from} ~ (1)}}\\
&= \int_{-\infty}^{\infty} h(s) E[X(t-\tau-s)Y(t)]\,\mathrm ds\\
&= \int_{-\infty}^{\infty} h(s) R_{X,Y}(\tau+s)\,\mathrm ds\\
&= \int_{-\infty}^{\infty} \tilde{h}(-s) R_{X,Y}(\tau+s)\,\mathrm ds
&{\scriptstyle{\tilde{h}(t) = h(-t)\ \forall \, t ~\text{is the time-reversed impulse response}}}\\
&= \int_{-\infty}^{\infty} \tilde{h}(\lambda) R_{X,Y}(\tau-\lambda)\,\mathrm d\lambda
&{\scriptstyle{\text{substitute}~ \lambda = -s}}
\end{align}$$
that is, $R_Y = \tilde{h}\star R_{X,Y}$, and it follows that
$$R_Y = \tilde{h}\star h \star R_X = (\tilde{h}\star h)\star R_X = R_h\star R_X$$
where $\tilde{h}(t) = h(-t)$ for all $t$ is the time-reversed impulse response
and $R_h = \tilde{h}*h$ is the autocorrelation function of the
deterministic signal $h(t)$. Translated to the frequency domain, this gives the
power spectral density relationship
$$S_Y(f) = |H(f)|^2 S_X(f).$$
