Autocorrelation of convolution integral Work out the autocorrelation
$r_Y(\tau) = E[Y(t)Y(t+\tau)]$
with 
$Y(t) = \int_{-\infty}^{\infty} h(t-u) x(u)$
and
$X$ a WSS, ergodic process
I always get: $h(t)* h(t+\tau) * r_X(\tau)$  (with $*$ convolution)
My approach
(after multiplication)
$\int h(t-u)  ( \int h(t+\tau -u')r_X(u'-u) \text{d}u' ) \text{d}u$
The inner integral is the definition of a convolution hence
$= \int h(t-u)  ( h(t+\tau)*r_x(t+\tau-u) ) \text{d}u$
This gives again a convolution integral with time variable t. so,
$= h(t)* h(t+\tau) * r_X(\tau) $
Correct solution:
$  r_X(\tau) * h(\tau) * h(-\tau) $
Question
What am I doing wrong?
 A: 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).$$
