# Why is the sample spectrum considered inconsistent?

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?

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

• Ok but if the variance will go to zero as you say, then the standard periodogram is consistent. But literature says it is not. Commented Apr 28 at 15:59
• You can easily confirm for yourself that the standard DFT with the equation I gave does go to zero. Create a bandlimited time domain waveform of noise which will have finite energy and a consistent variance, and try for yourself. If you don't scale it as I have done, it won't go to zero. I have another post I'll try to locate that details that specifically. Commented Apr 28 at 16:03
• Here's the post: dsp.stackexchange.com/a/93479/21048 showing the variance of each DFT bin depending on how we scale. Note how with 1/N scaling the variance in each bin is the total variance of the time domain signal divided by the total number of samples ($\sigma^2/N$). So given a total power of white noise in the time domain waveform, the power in each bin goes down at $1/N$, this is consistent as the resolution of each bin is also decreasing by $1/N$. Going to zero as $N$ goes to infinity. I updated the answer to try to make that clearer. Commented Apr 28 at 16:12
• Exactly and that has nothing to do with a direct DFT squared vs Welch. There is a distinction between creating a power spectrum or a power spectral density, which you can create with either approach. The total variance in the time domain is the total variance in the frequency domain, that is Parseval's Theorem and both appraoches are consistent with this (as long as you scale the way I described, along with additional scaling for windowing if that is done). Commented Apr 28 at 18:29
• @DanBoschen I can’t give a detailed response at the moment. I’ll update my answer later tonight with a description and tag you. Commented Apr 28 at 22:33

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

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

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.

• This is very interesting and insightful! You are making me review some other related posts where I didn't see any bias (see dsp.stackexchange.com/a/87731/21048) but I think your important point is the case of "band limited" which I didn't detail there (I will add that and reference/upvote this if a bias does appear). I am also seeing details with 1/N vs 1 / sqrt(N) that I need to confirm that may also affect this answer. Commented Apr 29 at 11:12
• @DanBoschen if you’re looking at LFM cases, if you use a decent window with relatively low averaging, you will likely not notice it as much. Bias is definitely more pronounced in line spectra cases, furthering your point on not using Welch’s for that type of signal. I used fairly “unlikely” parameters to emphasize my points in both figures, so practically speaking this bias will be less noticeable. Commented Apr 29 at 13:49
• ok I think that is a very good point to include within your answer, for those reading and applying this to practical applications. Did you see any flaws in my answer (since it was selected, I want to be sure it is indeed correct). I do want to try my prior experiments with white noise for the bandlimited case to see first hand the bias in case it is something I would be concerned about - thanks for this. Commented Apr 29 at 13:55
• @DanBoschen okay, I will add some comments about that. The only thing that I would change on a technicality is saying Welch’s has lower error in your last paragraph. Bias is pretty hard to quantify for Welch’s. However, practically speaking most implementations will likely have much less variance with minimal added bias, so this is a technicality. Commented Apr 29 at 14:12
• I see- it is likely confusion with what I refer to as “error”, what I mean is lower noise in the estimate- I will fix that terminology. See my post here where I demonstrated that specifically how “noisy” the DFT is vs the Welch Method— the difference is very significant which is the main point I wanted to make dsp.stackexchange.com/questions/87723/… Commented Apr 29 at 16:04