使用 scipy.signal.welch 找不到正确的能量

Question

对于给定的离散时间信号x(t)，间距为dt（等于1/fs，fs 是采样率），能量为：

E[x(t)] = sum(abs(x)**2.0)/fs

然后我做一个x(t)的DFT：

x_tf = np.fft.fftshift( np.fft.fft( x ) ) / ( fs * ( 2.0 * np.pi ) ** 0.5 )

并再次计算能量：

E[x_tf] = sum( abs( x_tf ) ** 2.0 ) * fs * 2 * np.pi / N

（这里的因子fs*2*np.pi/N = 脉动间距dk，fftfreq 的文档提供了有关频域间距的更多详细信息），我有同样的能量：

E[x(t)] = E[x_tf]

但是...当我使用scipy.signal.welch 计算x(t) 的功率谱密度时，我找不到合适的能量。 scipy.signal.welch 返回频率向量f 和能量Pxx（或每个频率的能量，取决于我们在scipy.signal.welch 的参数中输入的scaling）。

如何使用Pxx 找到与E[x(t)] 或E[x_tf] 相同的能量？我试图计算：

E_psd = sum(Pxx_den) / nperseg

其中nperseg 是 Welch 算法每个段的长度，fs 和 np.sqrt(2*np.pi) 等因素被取消，并用 nperseg 重新调整 E[x(t)]，但没有任何成功（订单数量级小于E[x(t)])

我使用以下代码生成信号：

#Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.

fs = 10e3   #sampling rate, dt = 1/fs
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

我做了以下得到功率谱密度：

f, Pxx_den = signal.welch(x, fs )

Answer 1

这种明显差异的解决方法在于仔细理解和应用

• 连续与离散傅里叶变换，以及
• 给定信号的能量、功率和功率谱密度

我也一直在为这个确切的问题而苦恼，所以我将在下面的讨论中尽可能明确。

离散傅里叶变换 (DFT)

满足某些可积性条件的连续信号 x(t) 具有傅里叶变换 X(f)。然而，当使用 离散 信号 x[n] 时，通常使用离散时间傅里叶变换 (DTFT)。我将 DTFT 表示为 X_{dt}(f)，其中dt 等于相邻样本之间的时间间隔。回答您的问题的关键要求您认识到 DTFT 不 等于相应的傅立叶变换！其实这两者是有关系的

X_{dt}(f) = (1 / dt) * X(f)

此外，离散傅里叶变换 (DFT) 只是 DTFT 的 离散 样本。当然，DFT 是 Python 在使用 np.fft.fft(...) 时返回的结果。因此，您计算的 DFT 不 等于傅立叶变换！

功率谱密度

scipy.signal.welch(..., scaling='density', ...) 返回离散信号 x[n] 的power spectral density (PSD) 的估计值。对 PSD 的完整讨论有点超出本文的范围，但对于简单的周期性信号（例如您的示例中的信号），PSD S_{xx}(f) 给出为

S_{xx} = |X(f)|^2 / T

其中|X(f)|是信号的傅立叶变换，T 是信号的总持续时间（以时间计）（如果您的信号 x(t) 是一个随机过程，我们必须对系统的许多实现进行整体平均。 ..)。信号中的总功率只是 S_{xx} 在系统频率带宽上的积分。使用上面的代码，我们可以编写

import scipy.signal

# Estimate PSD `S_xx_welch` at discrete frequencies `f_welch`
f_welch, S_xx_welch = scipy.signal.welch(x, fs=fs)

# Integrate PSD over spectral bandwidth
# to obtain signal power `P_welch`
df_welch = f_welch[1] - f_welch[0]
P_welch = np.sum(S_xx_welch) * df_welch

要联系您的np.fft.fft(...) 计算（返回 DFT），我们必须使用上一节中的信息，即

X[k] = X_{dt}(f_k) = (1 / dt) * X(f_k)

因此，要从 FFT 计算中计算功率谱密度（或总功率），我们需要认识到

S_{xx} = |X[k]|^2 * (dt ^ 2) / T

# Compute DFT
Xk = np.fft.fft(x)

# Compute corresponding frequencies
dt = time[1] - time[0]
f_fft = np.fft.fftfreq(len(x), d=dt)

# Estimate PSD `S_xx_fft` at discrete frequencies `f_fft`
T = time[-1] - time[0]
S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T

# Integrate PSD over spectral bandwidth to obtain signal power `P_fft`
df_fft = f_fft[1] - f_fft[0]
P_fft = np.sum(S_xx_fft) * df_fft

P_welch 和 P_fft 的值应该非常接近，并且接近信号中的预期功率，可以计算为

# Power in sinusoidal signal is simply squared RMS, and
# the RMS of a sinusoid is the amplitude divided by sqrt(2).
# Thus, the sinusoidal contribution to expected power is
P_exp = (amp / np.sqrt(2)) ** 2 

# For white noise, as is considered in this example,
# the noise is simply the noise PSD (a constant)
# times the system bandwidth. This was already
# computed in the problem statement and is given
# as `noise_power`. Simply add to `P_exp` to get
# total expected signal power.
P_exp += noise_power

注意：P_welch 和 P_fft 不会完全相等，甚至可能在数值精度范围内也不相等。这归因于与功率谱密度的估计相关的随机误差这一事实。为了减少此类错误，Welch 的方法将您的信号分成几个段（其大小由 nperseg 关键字控制），计算每个段的 PSD，并对 PSD 进行平均以获得对信号的更好估计PSD（平均分段越多，产生的随机误差越小）。 FFT 方法实际上等同于仅对一个大段进行计算和平均。因此，我们预计P_welch 和P_fft 之间存在一些差异，但我们应该预计P_welch 更准确。

信号能量

正如你所说，信号能量可以从 Parseval 定理的离散版本中获得

# Energy obtained via "integrating" over time
E = np.sum(x ** 2)

# Energy obtained via "integrating" DFT components over frequency.
# The fact that `E` = `E_fft` is the statement of 
# the discrete version of Parseval's theorem.
N = len(x)
E_fft = np.sum(np.abs(Xk) ** 2) / N

我们现在想了解上面通过scipy.signal.welch(...) 计算得出的S_xx_welch 与信号中的总能量E 之间的关系。从上面，S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T。重新排列这个表达式中的术语，我们看到np.abs(Xk) ** 2 = (T / (dt ** 2)) * S_xx_fft。此外，

从上面我们知道np.sum(S_xx_fft) = P_fft / df_fft 和P_fft 和P_welch 大致相等。进一步，P_welch = np.sum(S_xx_welch) / df_welch 让我们得到

np.sum(S_xx_fft) = (df_welch / df_fft) * np.sum(S_xx_welch)

此外，S_xx_fft = ((np.abs(Xk) * dt) ** 2) / T。将S_xx_fft 代入上述等式并重新排列各项，我们得到

np.sum(np.abs(Xk) ** 2) = (T / (dt ** 2)) * (df_welch / df_fft) * np.sum(S_xx_welch)

上式中的左侧 (LHS) 现在看起来应该非常接近根据 DFT 分量计算得出的信号总能量的表达式。现在，请注意T / dt = N，其中N 是信号中的采样点数。除以N，我们现在有一个LHS，根据定义，等于上面计算的E_fft。因此，我们可以通过 Welch 的 PSD 得到信号中的总能量

# Signal energy from Welch's PSD
E_welch = (1. / dt) * (df_welch / df_fft) * np.sum(S_xx_welch)

E、E_fft 和 E_welch 的值都应该非常接近 :) 正如上一节末尾所讨论的，我们确实预计 E_welch 与 E 和E_fft，但这是因为从 Welch 方法得出的值减少了随机误差（即更准确）。