Why is the sample spectrum considered inconsistent?

Question

In multiple sources of the literature, e.g. here or here, it is stated that the sample spectrum (i.e. calculating the spectrum of an input signal by simply squaring the magnitude of the DFT) is inconsistent, i.e. its variance does not go to zero as the record length N goes to infinity. This is why we use modified periodogram methods like Welch.

However, in other places, e.g. this answer on dsp se, it is shown that the variance of the dft does indeed decrease as N increase.

There seem to be a contradiction here. Which is true?

Answer 1

It depends on the scaling used in the DFT but this doesn't mean we can't use the DFT vs Welch for computing a power spectrum. More on the benefits of Welch further below, but first I will cover the scaling: If the signal has finite energy, the variance of each DFT bin will go to zero as $N$ goes to infinity when using a DFT scaled by the total number of samples $N$, or specifically:

$$X[k] = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j2\pi n k/N}$$

(And makes sense since the resolution bandwidth of each DFT bin also goes to zero as $N$ goes to infinity).

If we scale by $1\sqrt{N}$, the variance in the frequency domain will equal the variance in the time domain, which is consistent with Parseval's Theorem, regardless of $N$: Given a discrete time waveform $x[n]$, the variance of $x[n]$ will equal the variance of $X[k]$ which is specifically:

$$\sigma_k^2 = \sum_{n=0}^{N-1}|X[k]|^2$$

Where here $X[k]$ is given as:

$$X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x[n]e^{-j2\pi n k/N}$$

This keeps the DFT unitary and is consistent with Parseval's but does not result in a accurate power spectrum (see details on this at the bottom of my answer to DSP.SE #87723. showing how $1/N$ is the correct scaling, consistent for both tones and total power.)

Given a power spectrum, and the resolution bandwidth that was used, we can scale the result to be a power spectral density by dividing by the resolution bandwidth. This can be done with a DFT, and is typically adjusted as such in the standard implementations of the Welch algorithm found in MATLAB, Octave and Python's scipy.signal.

Regardless of the above points which show how a DFT can be properly used for creating a power spectrum, this is not why we used a modified periodogram like Welch. We can indeed use a squared 1/N DFT to compute the power spectrum as long as we deal with the limitations mentioned in the OP's link such as aliasing, spectral leakage, and lack of stationarity in the longer term, and we can convert it to a power spectral density with proper scaling for resolution bandwidth. The primary motivation for using the Welch spectral estimation technique is that the resulting noise variation in the estimate itself of the PSD will be less, at the expense of resolution bandwidth. Also given the Welch method uses computations of the DFT with shorter time blocks and with windowing, it addresses the stationarity issues (shorter time blocks) and spectral leakage (windowing) we may have with a single very long DFT (however we can also compute the spectrum with a shorter DFT and with windowing to address those). I demonstrate this in more detail in my answer to DSP.SE #79255. Ultimately it is the lower noise result in the resulting spectrum that is the primary advantage of the Welch algorithm for computing the power spectral density of noise and waveforms that have their power distributed over many frequencies (as opposed to single tones and waveforms with their spectral content less than the resolution bandwidth of the measurement, in which cases the Welch method is not useful over using the DFT directly).

Answer 2

The Stoica and Moses book gives a thorough derivation of this, so I’m not going to repeat it here. For filtered white noise, the output periodogram is

\begin{equation} \hat{\phi}_{y}(\omega) =|H(\omega)|^{2} \hat{\phi}_{e}(\omega) + \mathcal{O}(\frac{1}{\sqrt{N}}) \end{equation}

where $\hat{\phi}_{e}(\omega)$ is the noise PSD, and $\mathcal{O}\frac{1}{\sqrt{N}}$ is an unkown constant $\in \: [-\frac{1}{\sqrt{N}},\frac{1}{\sqrt{N}}]$ depending on the draw of the white noise, or something like that.

What this means is that, for filtered white noise, which can describe a large variety of signals, the variance of the periodogram approaches the variance of the white noise plus a constant proportional to $\frac{1}{\sqrt{N}}$ (with scaling of $\frac{1}{\sqrt{N}}$ for the FFT). If using $\frac{1}{N}$ in the FFT computation, you get an $N^{2}$ in the periodogram denominator, which will produced a biased estimate of the periodogram. Haven’t worked out the math in a while, but if I remember correctly it will get rid of the constant term due to the noise input asymptotically, making it consistent, but this results in a scale factor bias.

It’s important to note that Bartlett’s, Welch’s, etc, while consistent, are biased, not asymptotically unbiased, like the periodogram. To achieve this, you need MVM or a modern parametric technique.

EDIT #1 Discussion of statistical properties of the periodogram and its derivatives

There's a couple ways to dissect this. Often times the periodogram is thought of as unbiased, but it's not, it's asymptotically unbiased, and this is actually a very important point. For white noise, the periodogram is virtually unbiased, except for cases with very small $N$. However, for most other cases, the periodogram is biased, even for large $N$, especially for signals with high dynamic range, because of the sidelobes of the Fejér kernel. This is easily verified by generating an AR PSD as they have a closed form solution, and the exact PSD can be known from the AR parameters. A picture is shown below.

As you can see, the periodogram follows the true PSD with relatively zero bias until the high frequency content where it starts to deviate from the true PSD. To be fair, other PSD methods including mismatched parametric models also suffer this bias. In line spectral cases I believe Eigenvector produces an unbiased result. Regardless, it is important to note that in cases of high dynamic range, many PSD methods, particularly most filterbank methods, will suffer bias.

Notice also that Welch's method suffers a severe amount of bias (a decent amount of averaging was applied to emphasize). At the peaks, Welch's method underestimates the true PSD value, whereas elsewhere, Welch's method overestimates the true PSD value. This is because bias = resolution for spectral estimators. The less resolution you have (compared to full resolution), the more bias you will have 1. Because by definition Welch's method (or any periodogram derivative method) cannot be full resolution, it will always have inherent bias. This is also not just the case for line spectra in noise (even though the above picture was of an AR process, ie filtered white noise, which covers a wide variety of cases, other interesting things happen when using other examples). Take for example a linear FM chirp. Below is a picture of this case.

Because of the bias, Welch's slightly underestimates the bandwidth of the chirp, although it has better first sidelobe response. Playing around with the values of the window, segment size, and overlap size can get you a combination of 3 dB mainlobe width and PSL that gives the appearance of being unbiased and consistent, but this is trial and error and may not always work, and mathematically isn't consistent since resolution = bias.

With respect to my point in the comments about the error, the easiest way (IMO) to quantify the error, and the method used in 1 for a lot of their discussions, is MSE. The MSE of a PSD is $\sigma^{2} + (\text{bias})^{2}$. When empirically comparing single draws of different PSDs, this is a pretty hard thing to quantify, especially for a wide variety of different cases. For example, you could potentially make the argument that for the AR process the periodogram has less MSE, whereas in the LFM case Welch's may have less MSE. It's really hard to say, but my point was primarily to say that if you changed error to variance you would be correct across the board.

At the end of the day, without getting into far more complicated techniques, the only way I can think of to guarantee an unbiased and consistent PSD estimator is to have a matched parametric model. Hopefully this discussion clears up some things on what I was trying to get at!

EDIT 2 Practical Implications The pictures shown in this post don't necessarily reflect a likely selection of parameters for Welch's method. As such, they are meant to be extreme cases to demonstrate the bias that Welch's introduces. In practical scenarios with more reasonable parameter selection, it is unlikely that this amount of bias will be seen, although some bias will be present.

1 Stoica, Petre, and Randolph L. Moses. Spectral analysis of signals. Vol. 452. Upper Saddle River, NJ: Pearson Prentice Hall, 2005.