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, optional
x 時間序列的取樣頻率。預設值為 1.0。
- windowstr or tuple or array_like, optional
要使用的期望視窗。如果 window 是字串或元組，則會將其傳遞給
get_window以產生視窗值，預設情況下這些值是 DFT 偶數。 有關視窗列表和所需參數，請參閱get_window。如果 window 是 array_like，它將直接用作視窗，並且其長度必須為 nperseg。預設值為 Hann 視窗。- npersegint, optional
每個區段的長度。預設值為 None，但如果 window 是 str 或 tuple，則設定為 256；如果 window 是 array_like，則設定為視窗的長度。
- noverlapint, optional
區段之間重疊的點數。如果 None，
noverlap = nperseg // 2。預設值為 None。- nfftint, optional
使用的 FFT 長度，如果需要零填充 FFT。如果 None，則 FFT 長度為 nperseg。預設值為 None。
- detrendstr or function or False, optional
指定如何對每個區段進行去趨勢。如果
detrend是字串，則會將其作為 type 參數傳遞給detrend函數。如果它是一個函數，它會接收一個區段並返回一個去趨勢的區段。如果detrend是 False，則不進行去趨勢。預設值為 'constant'。- return_onesidedbool, optional
如果 True，則針對實數資料傳回單邊頻譜。如果 False，則傳回雙邊頻譜。預設值為 True，但對於複數資料，始終傳回雙邊頻譜。
- scaling{ ‘density’, ‘spectrum’ }, optional
選擇計算功率譜密度 ('density') (其中 Pxx 的單位為 V**2/Hz) 和計算平方幅度譜 ('spectrum') (其中 Pxx 的單位為 V**2，如果 x 以 V 為單位測量且 fs 以 Hz 為單位測量)。預設值為 'density'
- axisint, optional
計算週期圖所沿著的軸；預設值是最後一個軸 (即
axis=-1)。- average{ ‘mean’, ‘median’ }, optional
平均週期圖時使用的方法。預設值為 'mean'。
在版本 1.2.0 中新增。
- 返回:
- fndarray
取樣頻率的陣列。
- Pxxndarray
x 的功率譜密度或功率譜。
另請參閱
periodogram簡單、可選擇修改的週期圖
lombscargle用於不均勻取樣資料的 Lomb-Scargle 週期圖
註解
適當的重疊量將取決於視窗的選擇和您的需求。對於預設的 Hann 視窗，50% 的重疊是在準確估計訊號功率，同時不過度計算任何資料之間合理的權衡。較窄的視窗可能需要更大的重疊。
如果 noverlap 為 0，則此方法等效於 Bartlett 氏方法 [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()
