# How can we compute statistics of the DFT of a random signal?

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'm not a theoretical statistics person, so maybe you could lay out where you think your argument in the independent case goes awry? – Bill Barth Dec 24 '13 at 15:40
• It would be important to state what exactly you mean by "in the case that $x_n$ are correlated". I ask this because this can be interpreted as saying that "the choice of $x_n$ depends on $x_{n-1}$", implying a causal correlation, or "the $n$-dimensional point $x$ is chosen from a distribution such as $N(0,\Sigma)$ with a covariance matrix $\Sigma$". I suspect that the answer to your question depends on your definition of correlation. – Wolfgang Bangerth Dec 26 '13 at 14:22
• What I mean is basically that the distribution for $x_n$ depends in some way on the value of $x_{n-1}$. In my case where the noise is filtered Gaussian white noise, I think what will happens is that the mean of the distribution for $x_n$ depends on $x_{n-1}$ and the width of the distribution for $x_n$ depends on the time between points. – DanielSank Dec 26 '13 at 20:58
• Note that for real inputs $x_k$ you get $X_k=\bar X_{N-k}$. For everything else, the Fourier transform is just a linear transformation, so the usual rules for the linear transformation of Gaussian vectors apply. – Lutz Lehmann Dec 29 '13 at 19:14
• @LutzL, that sounds like a useful observation. Can you comment on or direct me to information on how to deal with correlation? – DanielSank Dec 29 '13 at 23:29

I asked this question on the Signal Processing exchange as well. Eventually I found the answer myself and posted it there. Here I copy/paste my answer as posted there.

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