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

Question

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 1

Math tools

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) \, . $$

Calculation

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}$.

Answer 2

I would like to give another take on @DanielSank's answer. We first suppose that $v_{n} \sim \mathcal{CN}(0, \sigma^{2})$ and is i.i.d. Its Discrete Fourier Transform is then:

$$ V_{k} = \frac{1}{N} \sum_{n=0}^{N-1} v_{n} e^{-j 2 \pi \frac{n}{N} k}$$.

We want to calculate the distribution of $V_{k}$ To start, we note that since $v_{n}$ is white Gaussian noise, it is circularly symmetric, so the real and imaginary parts of its Fourier Transform will distributed the same. Therefore, we only need to calculate the distribution of the real part and then combine it with the imaginary part.

So we separate $V_{k}$ into its real and imaginary parts. We have:

$$ V_{k} = \frac{1}{N} \sum_{n=0}^{N-1} v_{n} e^{-j 2 \pi \frac{n}{N} k} $$ $$ V_{k} = \frac{1}{N} \sum_{n=0}^{N-1} \big( R\{v_{n}\} + j I\{v_{n} \} \big) \cdot \big( cos(2 \pi \frac{n}{N} k) + j sin(2 \pi \frac{n}{N} k) \big)$$ $$ V_{k} = R\{V_{k}\}_{1} + R\{V_{k}\}_{2} + jI\{V_{k}\}_{1} + jI\{V_{k}\}_{2} $$ $$ V_{k} = R\{V_{k}\} + jI\{V_{k}\}$$

Where:

$$ R\{V_{k}\} = R\{V_{k}\}_{1} + R\{V_{k}\}_{2} $$ $$ I\{V_{k}\} = I\{V_{k}\}_{1} + I\{V_{k}\}_{2} $$

And:

$$ R\{V_{k}\}_{1} = \frac{1}{N} \sum_{n=0}^{N-1} R\{ v_{n} \} cos(2 \pi \frac{n}{N} k) $$

$$ R\{V_{k}\}_{2} = - \frac{1}{N} \sum_{n=0}^{N-1} I\{ v_{n} \} sin(2 \pi \frac{n}{N} k) $$

$$ I \{V_{k}\}_{1} = \frac{1}{N} \sum_{n=0}^{N-1} R\{ v_{n} \} sin(2 \pi \frac{n}{N} k) $$

$$ I\{V_{k}\}_{2} = \frac{1}{N} \sum_{n=0}^{N-1} I\{ v_{n} \} cos(2 \pi \frac{n}{N} k) $$

Now we work on deriving the distribution of $R\{V_{k}\}_{1}$ and $R\{V_{k}\}_{2}$. As in @DanielSank's answer, we define:

$$ x_{n,k} = \frac{1}{N} cos(2 \pi \frac{n}{N} k) R\{v_{n}\} = \frac{1}{N} c_{n,k} R\{v_{n}\}$$

Thus we can write: $$ R\{ V_{k} \}_{1} = \sum_{n=0}^{N-1} x_{n,k} $$

This allows us the easily apply the following facts about linear combinations of Gaussian random variables. Namely, we know that:

1. When $x \sim \mathcal{CN}(0, \sigma^{2})$ then $R\{ x \} \sim \mathcal{N}(0, \frac{1}{2} \sigma^{2})$
2. When $x \sim \mathcal{N}(\mu, \sigma^{2})$ then $cx \sim \mathcal{N}(c \mu, c^{2} \sigma^{2})$

Together, these imply that $x_{n,k} \sim \mathcal{N}(0, \frac{c^{2}_{n,k}}{2N^{2}} \sigma^{2})$. Now we work on the sum. We know that:

1. When $x_{n} \sim \mathcal{N}(\mu_{n}, \sigma^{2}_{n})$ then $y = \sum_{n=0}^{N-1} x_{n} \sim \mathcal{N}(\sum_{n=0}^{N-1} \mu_{n}, \sum_{n=0}^{N-1} \sigma^{2}_{n}) $
2. $\sum_{n=0}^{N-1} c^{2}_{n,k} = \frac{N}{2}$

These imply that:

$$ R\{V_{k}\}_{1} \sim \mathcal{N}(0, \sum_{n=0}^{N-1} \frac{c_{n,k}^{2}}{2N^{2}} \sigma^{2}) = \mathcal{N}(0 , \frac{\frac{N}{2}}{2N^{2}} \sigma^{2} = \mathcal{N}(0, \frac{\sigma^{2}}{4N})$$

So we have shown that:

$$ R\{V_{k}\}_{1} \sim \mathcal{N}(0, \frac{\sigma^{2}}{4N}) $$

Now we apply the same argument to $R\{V_{k}\}_{2}$. Abusing our notation, we rewrite:

$$ x_{n,k} = \frac{1}{N} sin(2 \pi \frac{n}{N} k) I\{v_{n}\} = \frac{1}{N} s_{n,k} I\{v_{n}\}$$

Repeating the same argument, and noting that the Gaussian is a symmetric distribution (so we can ignore the sign difference), gives us:

$$ R\{V_{k}\}_{2} \sim \mathcal{N}(0, \frac{\sigma^{2}}{4N}) $$

Since $\sum_{n=0}^{N-1} s^{2}_{n,k} = \frac{N}{2}$ as well. So therefore since $R\{V_{k}\} = R\{V_{k}\}_{1} + R\{V_{k}\}_{2}$, we get:

$$ R\{V_{k}\} \sim \mathcal{N}(0, \frac{\sigma^{2}}{4N} + \frac{\sigma^{2}}{4N}) = \mathcal{N}(0, \frac{\sigma^{2}}{2N}) $$

So we have shown that:

$$ R\{V_{k}\} \sim \mathcal{N}(0, \frac{\sigma^{2}}{2N}) $$

By circular symmetry, we also know then that:

$$ I\{V_{k}\} \sim \mathcal{N}(0, \frac{\sigma^{2}}{2N}) $$

So since $V_{k} = R\{V_{k}\} + jI\{V_{k}\}$, we finally arrive at:

$$ V_{k} \sim \mathcal{CN}(0, \frac{\sigma^{2}}{N}) $$

Therefore taking the DFT divides the variance by the length of the DFT window -- assuming the window is rectangular of course -- which is the same result as in @DanielSank's answer.