# 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.

• I would just look at each $X_k$ as the linear combination of random variables. The exponentials are, relative to a given $X_k$, just constants. – Jim Clay Dec 24 '13 at 14:36
• @Jim Clay, indeed this is just a linear combination of random variables. In the case that they are independent I just used the fact that the distribution of a linear combination of random variables is the Convolution of the distributions of the elements being summed. In the case that they're correlated I'm not sure how to proceed. – DanielSank Dec 24 '13 at 17:55
• Since it is Markovian you can make the random variables conditionally independent, but that would quickly become a mess for large numbers of points. I would try the conditionally independent route with a four point DFT and see how it works out. – Jim Clay Dec 24 '13 at 18:09
• I'm trying to work up an analytic formula. – DanielSank Dec 24 '13 at 18:43

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).