python - 使用 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]

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

我怎样才能找到和 E[x(t)]一样的能量或 E[x_tf]使用 Pxx ?我试图计算:

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', ...)返回 power spectral density (PSD) 的估计值离散信号 x[n]。对 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 的值(value)和 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

我们现在想了解如何 S_xx_welch , 以上通过 scipy.signal.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在值(value)上应该都非常接近:) 正如上一节末尾所讨论的，我们确实预计 E_welch 之间会有一些细微的差异。与 E 相比和 E_fft ，但这归因于从 Welch 方法得出的值减少了随机误差(即更准确)。