What happens to the 0 and Nyquist frequency when using Welch's method?

Question

Welch's method is a way to get better power spectral density (PSD) estimations than simple, naive periodograms. It has two main components:

• Cutting the input into many segments and average their individual PSD to reduce variance.
• Windowing each segment to reduce leakage from one frequency bin to the other.

My question relates to the first point: when we cut the input signal into many components, we essentially estimate the PSD of a shorter signal, which means that the frequency bins are wider. How do we then interpret the value at the Nyquist frequency? (and at the 0 frequency as well).

In the example below, I start with a perfectly flat one-sided PSD and compute a corresponding signal assuming a random phase for each bin (except for 0 and Nyquist because those are always real values for real signals). I then estimate the PSD of this signal using Welch's method with a lot of short segments. As you can see on the plot, the estimate is good except for the 0 and Nyquist frequency where it is half the expected valued. Why is that? How do I interpret these values? I can see why segmenting the full signal would yield a different DC component for each segment and therefor alter the average DC component, but what is happening to the Nyquist frequency?

This is not implementation-specific, it seems to be a fundamental side-effect of segmenting the input signal, as I have tried other windows as well as my own custom implementation of the method. I've included in the plot an average of 512 PSD across 512 independent signals sharing the same original PSD (in blue). This clearly doesn't drop to half the average at 0 and the Nyquist frequency.

corresponding python code:

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
N = 8192
freq = np.fft.rfftfreq(N)
def get_signal(spectrum: np.ndarray) -> np.ndarray:
    spectrum = init_spectrum.copy() + 0j
    spectrum[1:-1] /= 2  # 1-sided PSD -> pseudo-2-sided PSD
    spectrum = np.sqrt(spectrum)  # Power -> amplitude
    spectrum[1:-1] = np.exp(np.random.rand(N // 2 - 1)  2j * np.pi)  # random phase
    return np.sqrt(N) * np.fft.irfft(spectrum)
np.random.seed(1)
init_spectrum = np.ones(N // 2 + 1)
spec_full = []
for _ in range(512):
    full_signal = get_signal(init_spectrum)
    f_full, spec = signal.welch(full_signal, nperseg=N, detrend=False)
    spec_full.append(spec)
spec_full = np.mean(spec_full, axis=0)
f_welch, spectrum = signal.welch(full_signal, nperseg=N // 64, noverlap=N // 128, detrend=False)
plt.plot(f_full, spec_full, label="windowed, no segmentation")
plt.plot(f_welch, spectrum, label="Welch's method")
plt.plot(freq, init_spectrum, label="initial PSD")
plt.xlabel("frequency")
plt.ylabel("spectral intensity")
plt.legend()
plt.show()

Answer 1

This is all related to how you account for negative frequencies.

For real signals ($x[n] \in \mathbb{R}$) the spectrum has complex conjugate symmetry, i.e.

$$X[-k] = X^*[k]$$

Most people don't bother plotting the spectrum at the negative frequencies since it's a mirror image of the positives ones. The energy is still there and need to be accounted for.

scipy.signal.welch() gives you two options for that:

return_onesided: bool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. Defaults to True, but for complex data, a two-sided spectrum is always returned.

If you choose the default return_onesided=True, the function will add the energy at the negative frequency to the bin of the corresponding positive frequency. That's why you see 6 dB less gain at DC and Nyquist: they don't have a negative equivalent.

If you choose return_onesided=False you will get a two-sided PSD with (more or less) equal values everwhere.

Answer 2

I agree with Hilmar's answer, and also wanted to point out an issue in your code that leads to the issue you are describing. And also showing that is it not due to segmenting the data.

You start with a real signal that has a flat spectrum. Because it is real, the FT at negative frequencies is the conjugate of the FT at positive frequencies, as Hilmar pointed out.

If one considers the initial PSD as the PSD resulting from a signal with flat spectrum, in which the power of the negative frequencies have been added to the positive ones, your init_spectrum should be:

init_spectrum = np.ones(N // 2 + 1)
init_spectrum[0] = 0.5
init_spectrum[-1] = 0.5

Then all the 3 spectra align:

Here is the full code (note I also changed the argument in get_signal to init_spectrum):

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
N = 8192
freq = np.fft.rfftfreq(N)
def get_signal(init_spectrum: np.ndarray) -> np.ndarray:
    spectrum = init_spectrum.copy() + 0j
    spectrum[1:-1] /= 2  # 1-sided PSD -> pseudo-2-sided PSD
    spectrum = np.sqrt(spectrum)  # Power -> amplitude
    spectrum[1:-1] = np.exp(np.random.rand(N // 2 - 1)  2j * np.pi)  # random phase
    return np.sqrt(N) * np.fft.irfft(spectrum)
np.random.seed(1)
init_spectrum = np.ones(N // 2 + 1)
init_spectrum[0] = 0.5
init_spectrum[-1] = 0.5
spec_full = []
for _ in range(512):
    full_signal = get_signal(init_spectrum)
    f_full, spec = signal.welch(full_signal, nperseg=N, detrend=False)
    spec_full.append(spec)
spec_full = np.mean(spec_full, axis=0)
f_welch, spectrum = signal.welch(full_signal, nperseg=N // 64, noverlap=N // 128, detrend=False)
plt.plot(f_full, spec_full, label="windowed, no segmentation", lw=4)
plt.plot(f_welch, spectrum, label="Welch's method")
plt.plot(freq, init_spectrum, '--', label="initial PSD", c='r')
plt.xlabel("frequency")
plt.ylabel("spectral intensity")
plt.xlim(-0.01, 0.05)
plt.legend()
plt.show()