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

Question

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

但是，使用简单的正弦测试函数，当每个 window 使用的点数时，两种方法之间存在系统差异（scipy.signal.welch 的 nperseg 关键字； NFFT关键字为mlab.psd)为even，如下图为每window

64分的情况

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

当每个window使用的点数为奇数时，这两个函数有更好的一致性，如下图相同信号，但每 window

处理 65 个点

向信号添加其他特征（例如噪声）会稍微减弱这种影响，但它仍然存在，当每个 [=49= 使用偶数点时，两种方法之间的相对误差约为 10% ]. （我意识到上面每个 window 65 点计算出的 PSD 在最高频率处的相对误差相对较大。但是，这是在信号的奈奎斯特频率上，我不太关心这样的特征高频。当每个 window 的点数为偶数时，我更关心大部分信号带宽的大而系统的相对误差。

这种差异有什么原因吗？就准确性而言，一种方法是否优于另一种方法？我正在使用 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

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

生成的

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

这会导致为偶数长度生成不同的 window。从 this section of the Window function entry on Wikipedia 开始，与 Matplotlib 的 window_hanning 相比，这可能会给您带来轻微的性能优势，后者总是 returns 对称版本。

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

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

使用此 window 作为参数（在两个函数中使用 window=win）将给出以下图表，您可能会注意到两个 PSD 估计函数之间的一致性更接近：

或者，假设您不想要周期性 window 属性，您也可以使用：

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

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

Why do the power spectral density estimates from matplotlib.mlab.psd and scipy.signal.welch differ when the number of points per window is even?

signal-processing

discrete-mathematics

scipy

python-2.7

mlab