当每个窗口的点数均匀时,为什么matplotlib.mlab.psd和scipy.signal.welch的功率谱密度估计值会有所不同？

Question

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

但是,使用简单的正弦测试函数,当每个窗口使用的点数(nperseg关键字for scipy.signal.welch; NFFT关键字for 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

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

win = get_window(window, nperseg)

它使用默认参数fftbins=True,并根据scipy文档:

如果为True,则创建一个"周期性"窗口,准备与ifftshift一起使用,并乘以fft的结果(另请参见fftfreq).

这导致了均匀长度的不同生成窗口.从维基百科上的Window函数条目的这一部分来看,这可以比Matplotlib提供轻微的性能优势,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)