scipy.signal.welch 的缩放比例：频谱与密度

Question

我需要估计某些信号的功率谱密度，并为此目的使用 scipy.signal 提供的 welch 算法。即使经过数小时的研究我也找不到，与默认的“密度”设置相比，当我将 kwarg“缩放”设置为“频谱”时，输出到底有什么不同。

我认为'频谱'可能是'密度'乘以韦尔奇返回的两个相邻频率之间的频率差，这确实会将功率密度转换为功率。这仅在 1.5 倍以下是正确的，我真的很难理解它的来源。有什么想法吗？

import numpy as np
from scipy.signal import welch
dt = 1/4
tt = np.arange(0,1000,dt)
w = 1/5
data = 2 * np.sin(2*np.pi * w * tt) + np.random.normal(size = len(tt), scale = 1)

ff_welch, spectrum = welch(data, nperseg = 250, fs = 1/dt, scaling = 'spectrum')
ff_welch, density = welch(data, nperseg = 250, fs = 1/dt, scaling = 'density')
df_welch = (ff_welch[1]-ff_welch[0])
test =  spectrum / (density * df_welch)

变量 test 永远是 1.5...

Answer 1

密度是瓦特每赫兹，而频谱是瓦特每箱。如果你转换为 dB，你只会得到相同的图向上或向下移动。

在您的情况下（在 Welch 函数的输出中以 4 Hz 的采样率和 126 个 bin），每个 bin 得到 ~0.0159 Hz。直觉上它小于 1 赫兹，因此频谱值更小是有道理的。通常，您可以通过执行 1/0.0159 = 63 并使用 63 作为比例因子从频谱到密度。开窗使它变得更棘手，但你可以通过以下方式做到这一点：

• 频谱的比例因子是：

  sum(ff_welch)^2

• 密度的比例因子是：

  sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2)

这里有一些代码（抱歉：我在一个旧的 Python 2 盒子上）显示 Welch 输出的频谱和密度之间的比率。如果您使用自己的窗口和移位进行 FFT，然后使用比例因子，您可以获得与频谱的韦尔奇函数相同的值。我从来没有得到过精确的密度，所以它在小数点后第三位略有变化。

#!/usr/bin/env python
from __future__ import print_function

import numpy as np
from scipy.signal import welch
import matplotlib
import matplotlib.pyplot as plt

#Added '4.0'. This or 'float(1)/4' works. Python 2 doesn't know to make this a float.
dt = 1/4.0
tt = np.arange(0,1000,dt)
#Added '5.0'. This or 'float(1)/5' works. Python 2 doesn't know to make this a float.
w = 1/5.0
data = 2 * np.sin(2*np.pi * w * tt) + np.random.normal(size = len(tt), scale = 1)

ff_welch, spectrum = welch(data, nperseg = 250, fs = 1/dt, scaling = 'spectrum')
ff_welch, density = welch(data, nperseg = 250, fs = 1/dt, scaling = 'density')
#df_welch = (f1_welch[1]-f1_welch[0])
#test =  spectrum / (density * df_welch)
print ('density divided by spectrum is:', density/spectrum)

# Spectrum is for power per bin, density is for power per Hz.
# Scaling factor for spectrum is: 
#     sum(ff_welch)**2
# Scaling factor for density is:  
#     sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2) 
# If you divide density by spectrum, the sum(ff_welch)**2 cancels and you get: 
#    fs/NFFT * sum(ff_welch) / sum(ff_welch**2) 
fs = 1.0/dt
NFFT = len(density)
scaling_factor_density = sum(ff_welch)**2 * fs/NFFT * sum(ff_welch) / sum(ff_welch**2)
scaling_factor_spectrum = sum(ff_welch)**2
scaling_factor_comparison = 1 / (fs/NFFT * sum(ff_welch) / sum(ff_welch**2))
print ('scaling_factor_density:', scaling_factor_density)
print ('scaling_factor_spectrum:', scaling_factor_spectrum)
print ('spectrum scaling divided by density scaling:', scaling_factor_comparison)

freq_step_size = (fs/2.0) / NFFT
x_axis_Hz = np.arange(0,fs/2.0,freq_step_size)
plt.plot(x_axis_Hz, 10*np.log10(density), label='Density')
plt.plot(x_axis_Hz, 10*np.log10(spectrum), label='Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Output from scipy.signal.welch in dB')
plt.legend()
plt.show()