scipy.signal.

welch

scipy.signal.welch(x, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')

使用韦尔奇方法估计功率谱密度。

韦尔奇方法 [1] 通过将数据划分为重叠的段，计算每个段的修改后的周期图并对周期图进行平均，来计算功率谱密度的估计值。

参数:

xarray_like

测量值的时间序列

fsfloat, 可选

x 时间序列的采样频率。默认为 1.0。

windowstr 或 tuple 或 array_like, 可选

要使用的期望窗口。如果 window 是字符串或元组，则将其传递给 get_window 以生成窗口值，默认情况下这些窗口值是 DFT 偶数的。有关窗口和所需参数的列表，请参阅 get_window。如果 window 是类似数组的对象，则将其直接用作窗口，并且其长度必须为 nperseg。默认为汉宁窗。

npersegint, 可选

每个段的长度。默认为 None，但如果 window 是 str 或 tuple，则设置为 256，如果 window 是类似数组的对象，则设置为窗口的长度。

noverlapint, 可选

段之间重叠的点数。如果 None, noverlap = nperseg // 2. 默认为 None.

nfftint, 可选

如果需要零填充 FFT，则使用的 FFT 长度。如果 None，则 FFT 长度为 nperseg。默认为 None。

detrendstr 或 function 或 False, 可选

指定如何去除每个段的趋势。如果 detrend 是一个字符串，则将其作为 type 参数传递给 detrend 函数。如果它是一个函数，它会接受一个段并返回一个去趋势的段。如果 detrend 是 False，则不执行去趋势。默认为 “constant”。

return_onesidedbool, 可选

如果 True，则返回实数据的单边频谱。如果 False，则返回双边频谱。默认为 True，但对于复数数据，始终返回双边频谱。

scaling{ ‘density’, ‘spectrum’ }, 可选

选择计算功率谱密度（‘density’），其中 Pxx 的单位为 V**2/Hz；或计算平方幅度谱（‘spectrum’），其中 Pxx 的单位为 V**2，如果 x 的单位为 V，并且 fs 的单位为 Hz。默认为 ‘density’

axisint, 可选

计算周期图的轴；默认是在最后一个轴上计算（即 axis=-1）。

average{ ‘mean’, ‘median’ }, 可选

平均周期图时使用的方法。默认为 ‘mean’。

1.2.0 版本中新增。

返回:

fndarray

采样频率数组。

Pxxndarray

x 的功率谱密度或功率谱。

另请参阅

periodogram

简单、可选的修改周期图

lombscargle

针对不均匀采样数据的 Lomb-Scargle 周期图

注释

适当的重叠量取决于窗口的选择以及您的要求。对于默认的汉宁窗，50% 的重叠是在准确估计信号功率的同时，不会过度计算任何数据的合理折衷方案。较窄的窗口可能需要更大的重叠。

如果 noverlap 为 0，则此方法等效于巴特利特方法 [2]。

请参阅 频谱分析 部分的 SciPy 用户指南，以了解有关功率谱密度和（平方）幅度谱的缩放的讨论。

0.12.0 版本中新增。

参考

[1]

P. Welch, “The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms”, IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967.

[2]

M.S. Bartlett, “Periodogram Analysis and Continuous Spectra”, Biometrika, vol. 37, pp. 1-16, 1950.

示例

在您的浏览器中尝试！

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()

生成测试信号，一个 1234 Hz 的 2 Vrms 正弦波，被 0.001 V**2/Hz 的白噪声污染，采样频率为 10 kHz。

>>> fs = 10e3
>>> 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 += rng.normal(scale=np.sqrt(noise_power), size=time.shape)

计算并绘制功率谱密度。

>>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
>>> plt.semilogy(f, Pxx_den)
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.show()

如果我们对频谱密度的后半部分进行平均（不包括峰值），我们可以恢复信号上的噪声功率。

>>> np.mean(Pxx_den[256:])
0.0009924865443739191

现在计算并绘制功率谱。

>>> f, Pxx_spec = signal.welch(x, fs, 'flattop', 1024, scaling='spectrum')
>>> plt.figure()
>>> plt.semilogy(f, np.sqrt(Pxx_spec))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Linear spectrum [V RMS]')
>>> plt.show()

功率谱中的峰值高度是对 RMS 幅度的估计。

>>> np.sqrt(Pxx_spec.max())
2.0077340678640727

如果我们现在在信号中引入不连续性，将信号的一小部分幅度增加 50，我们可以看到平均功率谱密度的损坏，但使用中值平均值可以更好地估计正常行为。

>>> x[int(N//2):int(N//2)+10] *= 50.
>>> f, Pxx_den = signal.welch(x, fs, nperseg=1024)
>>> f_med, Pxx_den_med = signal.welch(x, fs, nperseg=1024, average='median')
>>> plt.semilogy(f, Pxx_den, label='mean')
>>> plt.semilogy(f_med, Pxx_den_med, label='median')
>>> plt.ylim([0.5e-3, 1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.legend()
>>> plt.show()

返回在新标签页中打开