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')[source]#
使用 Welch 方法估计功率谱密度。
Welch 方法 [1] 通过将数据分成重叠的段，计算每个段的修正周期图并平均周期图来计算功率谱密度的估计值。
- 参数：
- xarray_like
测量值的时序
- fsfloat, 可选
x 时序的采样频率。默认为 1.0。
- windowstr 或 tuple 或 array_like, 可选
要使用的窗口。如果 window 是字符串或元组，则将其传递给
get_window以生成窗口值，这些值默认情况下为 DFT-even。有关窗口和所需参数的列表，请参见get_window。如果 window 是 array_like，它将直接用作窗口，并且其长度必须为 nperseg。默认为 Hann 窗口。- npersegint, 可选
每个段的长度。默认为 None，但如果 window 是 str 或 tuple，则设置为 256，如果 window 是 array_like，则设置为窗口的长度。
- 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 周期图
注释
适当的重叠量将取决于窗口的选择和你的要求。对于默认的 Hann 窗口，50% 的重叠是准确估计信号功率和不重复计算任何数据之间的合理权衡。更窄的窗口可能需要更大的重叠。
如果 noverlap 为 0，则此方法等效于 Bartlett 方法 [2]。
有关功率谱密度的缩放和（平方）幅度谱的讨论，请参阅 频谱分析 部分的 SciPy 用户指南。
在版本 0.12.0 中添加。
参考文献
[1]P. Welch，“使用快速傅里叶变换估计功率谱：一种基于对短的、修正的周期图进行时间平均的方法”，IEEE Trans. Audio Electroacoust. 第 15 卷，第 70-73 页，1967 年。
[2]M.S. Bartlett，“周期图分析和连续谱”，Biometrika，第 37 卷，第 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 正弦波，被 10 kHz 采样频率的 0.001 V**2/Hz 白噪声所破坏。
>>> 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()
