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