How to compute the statistics of the DFT of correlated noise?

I would like to know how to compute the statistics of the discrete Fourier transform of a noise signal. To illustrate what I mean, I will first explain in detail a computation I have managed to do myself.

Suppose we have a discrete time series of values $x_n$ with $n$ from 0 to $N-1$. Each $x_n$ is a random variable, uncorrelated with the others, and Gaussian distributed with width $\sigma$. If I define the discrete Fourier transform

$$X_k = \frac{1}{N}\sum_{n=0}^{N-1} x_n e^{-2 \pi i n k / N}$$

then I find that $X_k$ is a complex random variable with real and imaginary parts Gaussian distributed with width $\sigma/\sqrt{2 N}$. I did the computation by using the fact that the distribution of a sum is the convolution of the distributions, etc.

Now I want to know how to do this computation in the case that $x_n$ are correlated. How does one approach this problem? I can make the assumption that the process is Markovian.

I had originally asked this on the Computation Science site, but I think here is a better fit.

Answer

After many weeks I give the answer to my own question.

There is a limit in which we can solve this problem in a reasonably simple way. Suppose we sum enough points in our DFT that the central limit theorem guarantees that the distribution of the sum's real and imaginary parts are Gaussian distributed. Then we only need to compute the variance. If we specialize the case of the variances of the real part of the DFT we can write

$$\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}\langle x_n x_m \rangle \cos(2\pi n k / N)\cos(2 \pi m l / N)$$ The thing to note is that $\langle x_n x_m \rangle$ is just the correlation function of the time domain samples, which we can denote $\rho(n-m)$. Putting this in, we get

$$\langle (\textrm{Re} X_k) (\textrm{Re} X_l) \rangle = \frac{1}{N^2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1} \rho(n-m) \cos(2\pi n k / N)\cos(2 \pi m l / N)$$

This form in terms of the correlation function is useful because the correlation function is frequently known; By the Wiener-Khinchin theorem the correlation function is the Fourier transform of the spectral density.

This sum can be computed numerically, or even analytically for some particular forms of $\rho$. To compute covariance with imaginary parts of $X$ just put $\sin$ instead of $\cos$.

I found the idea for this in J. Schoukens and J. Renneboog, IEEE Transactions on Instrumentation and Measurement, Vol. IM-35, No. 3, September (1986).