expected value of two LTI output signals multiplied (cross correlation) I have an input signal x (assumed to be iid Gaussian with $\mu=0$, $\sigma^2$) which is fed into two linear systems:


$y_1 = h_1 * x$
$y_2 = h_2 * x$


Now I would like to calculate $\mathbb{E}[y_1 y_2]$. What is the proper way to do this - preferably in the frequency domain?

Background: I want to calculate something like $\operatorname{var}(y_1 + y_2)$. This expands to $\operatorname{var}(y_1)+\operatorname{var}(y_2)+2\mathbb{E}[y_1 y_2]$. Now I am very familar with the first two components - the variance. I know that the variance of a stochastic signal is given by the autocorrelation sequence at position 0. With the Wiener-Kinchin theorem this translates to:

$$
\operatorname{var}(y_1) = r_{yy}(0) = \int_0^{\infty} \Phi_{yy}(f) df
$$

and then:

$$
\cdots = \sigma_x^2 \int_0^{\infty} |H_1(f)|^2 df
$$

Now for my problem - with $\mathbb{E}[y_1 y_2]$ I arrive at $\mathbb{E}[x^2 \cdots]$ (giving a variance) but the filter is not square magnitude.
Assumimg that the two linear systems are BIBO-stable, the random processes $\{Y_1(t)\}$ and $\{Y_2(t)\}$ are zero-mean WSS Gaussian processes with autocorrelation functions and power spectral densities given by 
\begin{align}
R_{Y_1} &amp;= \sigma^2 (h_1 \star \tilde{h}_1)\\
R_{Y_2} &amp;= \sigma^2 (h_2 \star \tilde{h}_2)\\
S_{Y_1} &amp;= \sigma^2 |H_1|^2\\
S_{Y_2} &amp;= \sigma^2 |H_2|^2
\end{align}
In fact, the processes are also jointly Gaussian and jointly WSS processes with cross-correlation function $$R_{Y_1, Y_2}(\tau) = E[Y_1(t), Y_2(t+\tau)] = \sigma^2 (h_1 \star \tilde{h}_2)$$ and cross-power spectral density 
$$S_{Y_1,Y_2}(f) = \sigma^2 H_1(f)H_2^*(f).$$
The OP wants to find $E[Y_1(t)Y_2(t)]$ which is given by
\begin{align}E[Y_1(t)Y_2(t)] &amp;= R_{Y_1,Y_2}(0)\\
&amp;= \sigma^2 h_1\star \tilde{h}_2\big|_{\tau=0}\tag{1}\\
&amp;= \sigma^2\int_{-\infty}^{\infty} H_1(f)H_2^*(f) \,\mathrm df\tag{2}
\end{align}
since the OP prefers the frequency-domain calculation. Personally, given only $h_1$ and $h_2$ (and not $H_1$ and $H_2$), I would say that it is easier/cheaper to just grind out $(1)$ rather than use the frequency-domain calculation $(2)$ but the OP might have specific reasons for opting for the frequency-domain calculation. In particular, $\operatorname{var}(Y_1+Y_2)$ which is what the OP seems to really want to find is just
$$\operatorname{var}(Y_1+Y_2) = \sigma^2\begin{cases}\displaystyle\int_{-\infty}^\infty |h_1(t)+h_2(t)|^2 \,\mathrm dt,\\
\displaystyle\sum_{n=-\infty}^\infty |h_1[n]+h_2[n]|^2\end{cases}$$
which seems easier than the frequency-domain version of the same calculation, but ymmv.