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?

• You have to start out with a discrete-time white Gaussian signal. Sampling a continuous-time white process is mathematically ill-defined, because the autocorrelation function of that process is described by a Dirac delta distribution. Since the autocorrelation of the sampled process is a sampled version of the autocorrelation of the original continuous process, you would need to consider a sampled version of the Dirac delta distribution, which is not defined. – Matt L. Jun 16 '15 at 20:08
• @MattL. "[The] autocorrelation of the sampled process is a sampled version of the autocorrelation of the original continuous process...". This isn't obvious to me, actually. Explaining that would be a useful self Q&A. – DanielSank Jun 16 '15 at 21:47
• Pay attention that the answers will hold for any unitary transformation of White Gaussian Noise. – Royi Aug 19 '19 at 19:42
• @Royi I disagree with your edit. Can you provide a link indicating that the upper case you've used in the title is consistent with a site policy? – DanielSank Aug 19 '19 at 19:44
• Restored your style. The major thing in the edit was adding relevant tags. – Royi Aug 19 '19 at 19:46

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

• All this is fine and dandy but when dealing with multiple random variables, and especially Gaussian random variables, the covariances are of some importance as is the question of which of the random variables are independent. Could you address this issue in your answer? (Marginally Gaussian random variables need not be jointly Gaussian too; are your $2N$ random variables jointly Gaussian? are they independent? – Dilip Sarwate Jun 16 '15 at 19:01
• @DilipSarwate this is a good question. Unfortunately I don't know the answer (yet). I am going through what you might call "self study" of stochastic signal processing and I haven't yet understood why the values of physical processes at different times are frequently modeled as jointly Gaussian (or even what that really means). I suspect it has to do with the differential equations governing the underlying process, but again I just don't know yet. If you care to make a self Q&A it would be really useful. Otherwise I will eventually ask the relevant questions on this site. – DanielSank Jun 16 '15 at 21:49
• @DilipSarwate I noticed you used the assumption of a Gaussian process in your answer to this other question. You even noted that a "Gaussian process" is not the same thing as just saying that the distribution of the process at a fixed $t$ is Gaussian distributed. This suggests that Gaussian processes are common in Nature/engineering. Is this true? If so, can you give me a hint as to where I can learn why? – DanielSank Jun 16 '15 at 21:58
• @DanielSank According to the central limit theorem, the combinaison of a very large number of independent random variables will always produce a normal distribution, no matter the original distribution of individual random variables. Since the normal distribution is very well studied, it is often assumed that the process being observed fits the central limit theorem. This is not always the case, (such as photons on a CCD, for example), but it tend to be a safe approximation for many macroscopic physics problems. – PhilMacKay Dec 19 '17 at 14:37
• @anishtain4 Here's a single (long!) line of Python that simulates thr process: import numpy as np; np.std(np.real(np.sum(np.random.normal(0, 1, (10000, 10000)) * np.exp(1.0j * 2 * np.pi * np.linspace(0, 1, 10000) * 50), axis=1))). When I do this, I get the output 70, which is equal to $\sqrt{10,000/2}$ as it should be. Perhaps you can compare your simulation to that line. – DanielSank Nov 21 '18 at 18:14

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.