I remember a result about the cross-correlation function of a stationary signal and that same signal, filtered through a linear time-invariant filter. There might have been an assumption that the signal was Gaussian. What is the cross-correlation function?
$$E\left\{ x(t) y(t) \right\} = \,\,\,?$$
where $y(t)$ is the output of an LTI filter with the input $x(t)$.
If the linear system (does not have to be time-invariant) is bounded-input bounded-output (BIBO) stable (meaning that if $|x(t)| \leq M$ for all $t$ for some finite positive number $M$, then there is another finite positive number $M^\prime$ such that $|y(t)| \leq M^\prime$ for all $t$), then when the input is a finite-variance random process (a.k.a. second-order proecess), one can interchange the order of the expectation operation and the integration/sum for the convolution. That is, if $$Y(t) = \int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds$$ where $h(t)$ is the impulse response of the LTI system, then $$\begin{align*} E[Y(t)] &= E\left[\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right ]\\ &= \int_{-\infty}^{\infty} E\left[h(s)X(t-s)\,\mathrm ds\right ]\\ &= \int_{-\infty}^{\infty} h(s)E\left[X(t-s)\right ]\,\mathrm ds\\ \end{align*} $$ and so if the input is a WSS process with $E[X(t)] = \mu_X$ for all $t$, $$E[Y(t)] = \mu_X\int_{-\infty}^{\infty} h(s)\,\mathrm ds = \mu_Y$$ is also a constant. In fact, $\{Y(t)\}$ is also a WSS process. If we define the crosscorrelation function $R_{X,Y}(\tau)$ as $R_{X,Y}(\tau) = E[X(t-\tau)Y(t)]$, 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*R_X.$$ For completeness, $$R_Y = \tilde{h}*h*R_X = (\tilde{h}*h)*R_X = R_h*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).$$
