python-2.7 - 当每个窗口的点数为偶数时，为什么 matplotlib.mlab.psd 和 scipy.signal.welch 的功率谱密度估计不同？

Question

matplotlib.mlab.psd(...)和 scipy.signal.welch(...)两者都实现 Welch 的平均周期图方法来估计信号的功率谱密度 (PSD)。假设将等效参数传递给每个函数，则返回的 PSD 应该是等效的。

然而，使用简单的正弦测试函数，当每个窗口使用的点数( nperseg 的 scipy.signal.welch 关键字； NFFT 的 mlab.psd 关键字)是偶数时，两种方法之间存在系统差异，如下面显示了每个窗口 64 个点的情况



上图显示了通过两种方法计算的 PSD，而下图显示了它们的相对误差(即两个 PSD 的绝对差除以其绝对平均值)。较大的相对误差表明两种方法之间存在较大的分歧。

当每个窗口使用的点数为奇数时，这两个函数的一致性要好得多，如下所示，对于相同的信号，但每个窗口处理了 65 个点



向信号添加其他特征(例如噪声)会在一定程度上减少这种影响，但它仍然存在，当每个窗口使用偶数个点时，两种方法之间的相对误差约为 10%。 (我意识到用上面每个窗口 65 个点计算的 PSD 在最高频率处的相对误差相对较大。但是，这是在信号的奈奎斯特频率处，我不太关心如此高频率的特征。我当每个窗口的点数为偶数时，更关心信号带宽的大部分上的大的和系统的相对误差)。

这种差异有原因吗？就准确性而言，一种方法是否优于另一种方法？我正在使用 scipy 版本 0.16.0 和 matplotlib 版本 1.4.3。

用于生成上述图的代码包含在下面。对于纯正弦信号，设置 A_noise = 0 ;对于嘈杂的信号，设置 A_noise到一个有限值。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch
from matplotlib import mlab

# Sampling information
Fs = 200.
t0 = 0
tf = 10
t = np.arange(t0, tf, 1. / Fs)

# Pure sinusoidal signal
f0 = 27.
y = np.cos(2 * np.pi * f0 * t)

# Add in some white noise
A_noise = 1e-3
y += (A_noise * np.random.randn(len(y)))

nperseg = 64    # even
# nperseg = 65  # odd

if nperseg > len(y):
    raise ValueError('nperseg > len(y); preventing unintentional zero padding')

# Compute PSD with `scipy.signal.welch`
f_welch, S_welch = welch(
    y, fs=Fs, nperseg=nperseg, noverlap=(nperseg // 2),
    detrend=None, scaling='density', window='hanning')

# Compute PSD with `matplotlib.mlab.psd`, using parameters that are
# *equivalent* to those used in `scipy.signal.welch` above
S_mlab, f_mlab = mlab.psd(
    y, Fs=Fs, NFFT=nperseg, noverlap=(nperseg // 2),
    detrend=None, scale_by_freq=True, window=mlab.window_hanning)

fig, axes = plt.subplots(2, 1, sharex=True)

# Plot PSD computed via both methods
axes[0].loglog(f_welch, S_welch, '-s')
axes[0].loglog(f_mlab, S_mlab, '-^')
axes[0].set_xlabel('f')
axes[0].set_ylabel('PSD')
axes[0].legend(['scipy.signal.welch', 'mlab.psd'], loc='upper left')

# Plot relative error
delta = np.abs(S_welch - S_mlab)
avg = 0.5 * np.abs(S_welch + S_mlab)
relative_error = delta / avg
axes[1].loglog(f_welch, relative_error, 'k')
axes[1].set_xlabel('f')
axes[1].set_ylabel('relative error')

plt.suptitle('nperseg = %i' % nperseg)
plt.show()

Answer 1

虽然参数可能看起来是等效的，但对于均匀大小的窗口，窗口参数可能略有不同。更具体地说，除非提供了特定的窗口向量，否则 scipy 使用的窗口 welch函数是用

win = get_window(window, nperseg)

它使用默认参数 fftbins=True ，并根据 scipy documentation :

If True, create a “periodic” window ready to use with ifftshift and be multiplied by the result of an fft (SEE ALSO fftfreq).

这导致对于偶数长度生成不同的窗口。来自 this section of the Window function entry on Wikipedia ，与 Matplotlib 的 window_hanning 相比，这可能会给您带来轻微的性能优势。它总是返回对称版本。

要使用相同的窗口，您可以为两个 PSD 估计函数明确指定窗口向量。例如，您可以使用以下方法计算此窗口:

win = scipy.signal.get_window('hanning',nperseg)

使用此窗口作为参数(在两个函数中都带有 window=win)将给出下图，您可能会注意到两个 PSD 估计函数之间的一致性更接近:



或者，假设您不想要周期窗口属性，您也可以使用:

win = mlab.window_hanning(np.ones(nperseg)) # or
win = scipy.signal.get_window('hanning',nperseg,fftbins=False)