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')[源代码]#
使用 Welch 方法估计功率谱密度。
Welch 方法 [1] 通过将数据分成重叠的片段，计算每个片段的修正周期图并对周期图进行平均来估计功率谱密度。
- 参数:
- x数组状
测量值的时间序列
- fs浮点数，可选
x 时间序列的采样频率。默认为 1.0。
- window字符串、元组或数组状，可选
要使用的所需窗口。如果 window 是字符串或元组，它将传递给
get_window以生成窗口值，这些值默认是 DFT 偶数的。有关窗口和所需参数的列表，请参阅get_window。如果 window 是数组状，它将直接用作窗口，其长度必须为 nperseg。默认为 Hann 窗口。- nperseg整数，可选
每个片段的长度。默认为 None，但如果 window 是字符串或元组，则设置为 256；如果 window 是数组状，则设置为窗口的长度。
- noverlap整数，可选
片段之间的重叠点数。如果为 None，则
noverlap = nperseg // 2。默认为 None。- nfft整数，可选
使用的 FFT 长度，如果需要零填充 FFT。如果为 None，则 FFT 长度为 nperseg。默认为 None。
- detrend字符串、函数或 False，可选
指定如何对每个片段进行去趋势。如果
detrend是字符串，它将作为 type 参数传递给detrend函数。如果它是一个函数，它接受一个片段并返回一个去趋势的片段。如果detrend为 False，则不进行去趋势。默认为 'constant'。- return_onesided布尔值，可选
如果为 True，则返回实数数据的一个单边谱。如果为 False，则返回一个双边谱。默认为 True，但对于复数数据，始终返回一个双边谱。
- scaling{ ‘density’（密度）, ‘spectrum’（频谱） }，可选
选择计算功率谱密度（'density'），此时 Pxx 的单位是 V**2/Hz；或者计算平方幅值谱（'spectrum'），此时 Pxx 的单位是 V**2，如果 x 以 V 为单位测量，fs 以 Hz 为单位测量。默认为 'density'。
- axis整数，可选
计算周期图的轴；默认为最后一个轴（即
axis=-1）。- average{ ‘mean’（平均值）, ‘median’（中位数） }，可选
对周期图进行平均时使用的方法。默认为 'mean'。
在 1.2.0 版中新增。
- 返回:
- fndarray
采样频率数组。
- Pxxndarray
x 的功率谱密度或功率谱。
另请参阅
csd使用 Welch 方法的互功率谱密度
periodogram简单，可选修正的周期图
lombscargle用于非均匀采样数据的 Lomb-Scargle 周期图
注意
适当的重叠量将取决于窗口的选择和您的要求。对于默认的 Hann 窗口，50% 的重叠是准确估计信号功率与不过度统计任何数据之间的合理权衡。更窄的窗口可能需要更大的重叠。如果 noverlap 为 0，此方法等同于 Bartlett 方法 [2]。
平方幅值（
scaling='spectrum'）除以频谱功率密度（scaling='density'）的比率是常数因子sum(abs(window)**2)*fs / abs(sum(window))**2。如果 return_onesided 为True，则负频率的值将添加到相应正频率的值中。有关功率谱密度和（平方）幅值谱缩放的讨论，请参阅 频谱分析 部分，该部分位于 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()
生成一个测试信号，一个 2 Vrms 的 1234 Hz 正弦波，被 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()
