Python SciPy signal.periodogram用法及代碼示例

本文簡要介紹 python 語言中 scipy.signal.periodogram 的用法。

用法:

scipy.signal.periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1)#

使用周期圖估計功率譜密度。

參數 ：：

x： array_like

測量值的時間序列

fs： 浮點數，可選

x 時間序列的采樣頻率。默認為 1.0。

window： str 或 tuple 或 數組，可選

想要使用的窗口。如果窗戶是一個字符串或元組，它被傳遞給scipy.signal.get_window生成窗口值，默認為DFT-even。看scipy.signal.get_window獲取窗口列表和所需參數。如果窗戶與 數組 類似，它將直接用作窗口，並且其長度必須等於計算周期圖的軸的長度。默認為‘boxcar’。

nfft： 整數，可選

使用的 FFT 的長度。如果 None 將使用 x 的長度。

detrend： str 或函數 或False， 可選的

指定如何去除每個段的趨勢。如果scipy.signal.detrend是一個字符串，它作為類型參數scipy.signal.detrend函數。如果它是一個函數，它接受一個段並返回一個去趨勢的段。如果scipy.signal.detrend是False，沒有去趨勢。默認為‘constant’。

return_onesided： 布爾型，可選

如果為 True，則返回真實數據的單側頻譜。如果為 False，則返回兩側頻譜。默認為 True，但對於複雜數據，始終返回兩側頻譜。

scaling： { ‘density’, ‘spectrum’ }，可選

在計算 Pxx 的單位為 V**2/Hz 的功率譜密度 (‘density’) 和計算 Pxx 的單位為 V**2 的功率譜 (‘spectrum’) 之間進行選擇，如果 x 以 V 為單位測量並且fs 以赫茲為單位。默認為‘density’

axis： 整數，可選

計算周期圖的軸；默認值在最後一個軸上(即 axis=-1 )。

返回 ：：

f： ndarray

采樣頻率數組。

Pxx： ndarray

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

注意：

例子：

>>> 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.periodogram(x, fs)
>>> plt.semilogy(f, Pxx_den)
>>> plt.ylim([1e-7, 1e2])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD [V**2/Hz]')
>>> plt.show()

如果我們平均頻譜密度的後半部分，以排除峰值，我們可以恢複信號上的噪聲功率。

>>> np.mean(Pxx_den[25000:])
0.000985320699252543

現在計算並繪製功率譜。

>>> f, Pxx_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')
>>> plt.figure()
>>> plt.semilogy(f, np.sqrt(Pxx_spec))
>>> plt.ylim([1e-4, 1e1])
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Linear spectrum [V RMS]')
>>> plt.show()

功率譜中的峰值高度是 RMS 幅度的估計值。

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