What are the statistics of the discrete Fourier transform of white Gaussian noise?

Consider a white Gaussian noise signal $ x \left( t \right) $.
If we sample this signal and compute the discrete Fourier transform, what are the statistics of the resulting Fourier amplitudes?

Answer

We can do the calculation using some basic elements of probability theory and Fourier analysis. There are three elements (we denote the probability density of a random variable $X$ at value $x$ as $P_X(x)$):

1. 
   

   Given a random variable $X$ with distribution $P_X(x)$, the distribution of the scaled variable $Y = aX$ is $P_Y(y) = (1/a)P_X(y/a)$.

   
   


2. 
   

   The probability distribution of a sum of two random variables is equal to the convolution of the probability distributions of the summands. In other words, if $Z = X + Y$ then $P_Z(z) = (P_X \otimes P_Y)(z)$ where $\otimes$ indicates convolution.

   
   


3. 
   

   The Fourier transform of the convolution of two functions is equal to the product of the Fourier transforms of those two functions. In other words:

   
   

$$\int dx \, (f \otimes g)(x) e^{-i k x} = \left( \int dx \, f(x) e^{-ikx} \right) \left( \int dx \, g(x) e^{-ikx} \right) \, . $$

Denote the random process as $x(t)$. Discrete sampling produces a sequence of values $x_n$ which we assume to be statistically uncorrelated. We also assume that for each $n$ $x_n$ is Gaussian distributed with standard deviation $\sigma$. We denote the Gaussian function with standard deviation $\sigma$ by the symbol $G_\sigma$ so we would say that $P_{x_n}(x) = G_{\sigma}(x)$.

The discrete Fourier transform amplitudes are defined as $$X_k \equiv \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N} \, .$$ Focusing for now on just the real part we have $$\Re X_k = \sum_{n=0}^{N-1} x_n \cos(2 \pi n k /N) \, .$$ This is just a sum, so by rule #2 the probability distribution of $\Re X_k$ is equal to the multiple convolution of the probability distributions of the terms being summed. We rewrite the sum as $$\Re X_k = \sum_{n=0}^{N-1} y_n$$ where $$y_n \equiv x_n \cos(2\pi n k / N) \, .$$ The cosine factor is a deterministic scale factor. We know that the distribution of $x_n$ is $G_\sigma$ so we can use rule #1 from above to write the distribution of $y_n$: $$P_{y_n}(y) = \frac{1}{\cos(2 \pi n k / N)}G_\sigma \left( \frac{y}{\cos(2 \pi n k / N)} \right) = G_{\sigma c_{n,k}}(y)$$ where for brevity of notation we've defined $c_{n,k} \equiv \cos(2\pi n k / N)$.

Therefore, the distribution of $\Re X_k$ is the multiple convolution over the functions $G_{\sigma c_{n,k}}$: $$P_{\Re X_k}(x) = \left( G_{\sigma c_{0,k}} \otimes G_{\sigma c_{1,k}} \otimes \cdots \right)(x) \, .$$

It's not obvious how to do the multiple convolution, but using rule #3 it's easy. Denoting the Fourier transform of a function by $\mathcal{F}$ we have $$\mathcal{F}(P_{\Re X_k}) = \prod_{n=0}^{N-1} \mathcal{F}(G_{\sigma c_{n,k}}) \, .$$

The Fourier transform of a Gaussian with width $\sigma$ is another Gaussian with width $1/\sigma$, so we get \begin{align} \mathcal{F}(P_{\Re X_k})(\nu) &= \prod_{n=0}^{N-1} G_{1/\sigma c_{n,k}}\\ &= \prod_{n=0}^{N-1} \sqrt{\frac{\sigma^2 c_{n,k}^2}{2\pi}} \exp \left[ \frac{-\nu^2}{2 (1 / \sigma^2 c_{n,k}^2)}\right] \\ &= \left( \frac{\sigma^2}{2\pi} \right)^{N/2} \left(\prod_{n=0}^{N-1}c_{n,k} \right) \exp \left[ -\frac{\nu^2}{2} \sigma^2 \sum_{n=0}^{N-1} \cos(2\pi nk/N)^2 \right] \, . \end{align} All of the stuff preceding the exponential are independent of $\nu$ and are therefore normalization factors, so we ignore them. The sum is just $N/2$ so we get \begin{align} \mathcal{F}(P_{\Re X_k}) &\propto \exp \left[ - \frac{\nu^2}{2} \sigma^2 \frac{N}{2}\right]\\ &= G_{\sqrt{2 / \sigma^2 N}} \end{align} and therefore $$P_{\Re X_k} = G_{\sigma \sqrt{N/2}} \, .$$

We have therefore computed the probability distribution of the real part of the Fourier coefficient $X_k$. It is Gaussian distributed with standard deviation $\sigma \sqrt{N/2}$. Note that the distribution is independent of the frequency index $k$, which makes sense for uncorrelated noise. By symmetry the imaginary part should be distributed exactly the same.

Intuitively we expect adding more integration should reduce the width of the resulting noise distribution. However, we found that the standard deviation of the distribution of $X_k$ grows as $\sqrt{N}$. This is just due to our choice of normalization of the discrete Fourier transform. If we had instead normalized it like this $$X_k = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N}$$ then we would have found $$P_{\Re X_k} = G_{\sigma / \sqrt{2N}}$$ which agrees with the intuition that the noise distribution gets smaller as we add more data. With this normalization, a coherent signal would demodulate to a fixed amplitude phasor, so we recover the usual relation that the ratio of the signal to noise amplitudes scales as $\sqrt{N}$.