5
$\begingroup$

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.

Blue: Average of 512 PSD estimate with no segmentation is noisy but clearly doesn't dip at x=0 or x=0.5. Orange: PSD estimate with Welch's method is a bit noisy around y=1 except for x=0 and x=0.5 where the estimate is ~0.5. Green: initial flat spectrum at y=1.

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

$\endgroup$

2 Answers 2

1
$\begingroup$

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:

enter image description here

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()
$\endgroup$
2
  • $\begingroup$ The wording of your answer is what made it click for me: I was assuming that a flat spectrum meant "flat in the 1-sided" picture, but that doesn't really make much sense. $\endgroup$
    – dedebenui
    Commented Feb 21 at 8:12
  • $\begingroup$ Very glad to hear! $\endgroup$
    – Chachni
    Commented Feb 21 at 16:48
4
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Interesting there's no correction factor for this. Or at least an option to do the correction factor. Matlab does the same thing too. $\endgroup$
    – Baddioes
    Commented Dec 18, 2023 at 16:21
  • $\begingroup$ Thank you for the answer. I figured that one-sided vs two-sided was playing a role, which is why I did take care of this in my example, and why the blue line shows the correct spectrum. The drop at 0 and Nyquist in the orange line seem more related to the fact that the data is segmented then averaged. $\endgroup$
    – dedebenui
    Commented Dec 20, 2023 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.