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).$$
