使用 scipy.signal.welch 找不到合适的能量
Can't find the right energy using scipy.signal.welch
对于给定的离散时间信号x(t)，间距为dt（等于1/fs，fs是采样率），能量为：
E[x(t)] = sum(abs(x)**2.0)/fs
然后我做一个 DFT x(t):
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 returns 频率向量 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) 等因素被抵消，并重新调整 E[x(t)] 与 nperseg，但没有任何成功（数量级小于 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 )
解决这种明显差异的方法在于仔细理解和应用
- 连续与离散傅立叶变换，以及
- 给定信号的能量、功率和功率谱密度
我也一直在为这个确切的问题而苦苦挣扎，所以我会在下面的讨论中尽可能明确。
离散傅立叶变换 (DFT)
满足一定可积性条件的连续信号x(t)具有傅立叶变换X(f)。然而，当处理 discrete 信号 x[n] 时，通常通常使用离散时间傅里叶变换 (DTFT)。我将 DTFT 表示为 X_{dt}(f)，其中 dt 等于相邻样本之间的时间间隔。回答你这个问题的关键是你要认识到DTFT不等于对应的傅立叶变换！事实上，两者是相关的
X_{dt}(f) = (1 / dt) * X(f)
此外，离散傅里叶变换 (DFT) 只是 DTFT 的 离散 样本。当然，DFT 就是使用 np.fft.fft(...) 时的 Python return。因此，您计算的 DFT 不 等于傅里叶变换！
功率谱密度
scipy.signal.welch(..., scaling='density', ...) return 是对离散信号 x[n] 的 power spectral density (PSD) 的估计。 PSD 的完整讨论有点超出了这个 post 的范围，但是对于一个简单的周期信号（例如您的示例中的信号），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(...) 计算（return 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 更准确。
信号能量
正如您所说，信号能量可以从帕塞瓦尔定理的离散版本中获得
# 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 的值应该都非常接近 :) 正如在上一节末尾讨论的那样，我们确实预计 [= =46=] 与 E 和 E_fft 相比，但这归因于这样一个事实，即从韦尔奇方法得出的值减少了随机误差（即更准确）。
对于给定的离散时间信号x(t)，间距为dt（等于1/fs，fs是采样率），能量为：
E[x(t)] = sum(abs(x)**2.0)/fs
然后我做一个 DFT x(t):
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 returns 频率向量 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) 等因素被抵消，并重新调整 E[x(t)] 与 nperseg，但没有任何成功（数量级小于 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 )
解决这种明显差异的方法在于仔细理解和应用
- 连续与离散傅立叶变换，以及
- 给定信号的能量、功率和功率谱密度
我也一直在为这个确切的问题而苦苦挣扎，所以我会在下面的讨论中尽可能明确。
离散傅立叶变换 (DFT)
满足一定可积性条件的连续信号x(t)具有傅立叶变换X(f)。然而，当处理 discrete 信号 x[n] 时，通常通常使用离散时间傅里叶变换 (DTFT)。我将 DTFT 表示为 X_{dt}(f)，其中 dt 等于相邻样本之间的时间间隔。回答你这个问题的关键是你要认识到DTFT不等于对应的傅立叶变换！事实上，两者是相关的
X_{dt}(f) = (1 / dt) * X(f)
此外，离散傅里叶变换 (DFT) 只是 DTFT 的 离散 样本。当然，DFT 就是使用 np.fft.fft(...) 时的 Python return。因此，您计算的 DFT 不 等于傅里叶变换！
功率谱密度
scipy.signal.welch(..., scaling='density', ...) return 是对离散信号 x[n] 的 power spectral density (PSD) 的估计。 PSD 的完整讨论有点超出了这个 post 的范围，但是对于一个简单的周期信号（例如您的示例中的信号），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(...) 计算（return 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 更准确。
信号能量
正如您所说，信号能量可以从帕塞瓦尔定理的离散版本中获得
# 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。因此，我们可以通过
# Signal energy from Welch's PSD
E_welch = (1. / dt) * (df_welch / df_fft) * np.sum(S_xx_welch)
E、E_fft 和 E_welch 的值应该都非常接近 :) 正如在上一节末尾讨论的那样，我们确实预计 [= =46=] 与 E 和 E_fft 相比，但这归因于这样一个事实，即从韦尔奇方法得出的值减少了随机误差（即更准确）。