Covariance between real and imaginary parts of Fourier transform of a stationary time series Since Fourier transform of a random stationary process in time (in the case of existence) is not necessarily real, my question is what is the relation between the covariance of real and imaginary parts of the transformed random variable.

More specifically, according to Wiener-Khinchin theorem for the time function $x(t)$ with some conditions the Fourier transform $X(\omega)$ exists. 

The question is can one recover  $$\mathbb{E}\left[\Re\{X(\omega)\},\Re\{X(\omega')\}\right]-\mathbb{E}\left[\Re\{X(\omega)\}\right]\mathbb{E}\left[\Re\{X(\omega')\}\right]$$ 
and $$\mathbb{E}\left[\Im\{X(\omega)\},\Im\{X(\omega')\}\right]-\mathbb{E}\left[\Im\{X(\omega)\}\right]\mathbb{E}\left[\Im\{X(\omega')\}\right]$$ 

from $$\mathbb{E}\left[\bar{x}(t),x(s)\right]-\mathbb{E}\left[\bar{x}(t)\right]\mathbb{E}\left[x(s)\right]\;?$$

I would be happy to be guided to a reference as well.

This is the most general case but the simpler case where the input signal is real is the main concern for me. I would be happy to be guided to a reference as well.
Note that in general the Fourier transform of a stationary process $x(t)$ does not exist. The Wiener-Khinchin theorem only states that under certain conditions the power spectral density of $x(t)$ exists, and it can be computed as the Fourier transform of the autocorrelation function of $x(t)$.

Having said that, if for some reason one assumes that the Fourier transform of $x(t)$ exists, then you can do the math and see if you get a useful result. So, let's see. I assume that $x(t)$ is real-valued. This is not at all necessary, but just simplifies things a bit.

Since $\text{Re}\{X(\omega)\}=\frac12[X(\omega)+X^*(\omega)]$ and $\text{Im}\{X(\omega)\}=\frac{1}{2j}[X(\omega)-X^*(\omega)]$, we can compute the desired expectations from $E\{X(\omega)X(\omega')\}$ and $E\{X^*(\omega)X(\omega')\}$. For the first of these expectations we have

$$E\{X(\omega)X(\omega')\}=E\left\{\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt
\int_{-\infty}^{\infty}x(t')e^{-j\omega' t'}dt'  \right\}=\\
=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}E\{x(t)x(t')\}e^{-j\omega t}e^{-j\omega't'}dtdt'$$

With the substution $\tau=t'-t$ and with the autocorrelation function $R_x(\tau)=E\{x(t)x(t+\tau)\}$ (assuming wide-sense stationarity of $x(t)$) we obtain

$$E\{X(\omega)X(\omega')\}=\int_{-\infty}^{\infty}R_x(\tau)e^{-j\omega'\tau}d\tau\int_{-\infty}^{\infty}e^{-j(\omega+\omega')t}dt$$

Interpreting the second integral as a distribution, and with the power spectrum $S_x(\omega)=\mathcal{F}\{R_x(\tau)\}$ we finally get

$$E\{X(\omega)X(\omega')\}=2\pi\delta(\omega+\omega')S_x(\omega')=
2\pi\delta(\omega+\omega')S_x(-\omega)=2\pi\delta(\omega+\omega')S_x(\omega)\tag{1}$$

because $S_x(\omega)$ is an even function (we assumed $x(t)$ to be real-valued). The other expectation $E\{X^*(\omega)X(\omega')\}$ can be derived in a completely analogous manner. The result is

$$E\{X^*(\omega)X(\omega')\}=2\pi\delta(\omega-\omega')S_x(\omega)\tag{2}$$

From (1) and (2) we can obtain the desired expectations

$$E\{\text{Re}\{X(\omega)\}\text{Re}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')+\delta(\omega+\omega')]\\
E\{\text{Im}\{X(\omega)\}\text{Im}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')-\delta(\omega+\omega')]$$

It remains to compute the quantities $E\{\text{Re}\{X(\omega)\}$ and $E\{\text{Im}\{X(\omega)\}$. We can easily derive them from the expectations $E\{X(\omega)\}$ and $E\{X^*(\omega)\}$:

$$E\{X(\omega)\}=E\left\{ \int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt \right\}=
\int_{-\infty}^{\infty}E\{x(t)\}e^{-j\omega t}dt=
\mu_x\int_{-\infty}^{\infty}e^{-j\omega t}dt=2\pi\mu_x\delta(\omega)$$

with $\mu_x=E\{x(t)\}$, where we have again assumed wide-sense stationarity of $x(t)$. Obviously, we get the same result for $E\{X^*(\omega)\}$:

$$E\{X^*(\omega)\}=2\pi\mu_x\delta(\omega)$$

This results in

$$E\{\text{Re}\{X(\omega)\}\}=2\pi\mu_x\delta(\omega)\quad\text{and}\quad
E\{\text{Im}\{X(\omega)\}\}=0$$

So we have

$$E\{\text{Re}\{X(\omega)\}\}\cdot E\{\text{Re}\{X(\omega')\}\}=4\pi^2\mu_x^2\delta(\omega)\delta(\omega')\\
E\{\text{Im}\{X(\omega)\}\}\cdot E\{\text{Im}\{X(\omega')\}\}=0\tag{3}$$

Combining (3) with (1) and (2) you obtain the desired result. From the result you can see that the covariances of the real parts of $X(\omega)$ and $X(\omega')$ vanish everywhere except for $\omega=\pm\omega'$. The same is true for the covariances of the imaginary parts. The difference is an additional term for $\omega=\omega'=0$ for the covariance of the real parts.