python scipy.signal 包络_python scipy signal.welch用法及代码示例

使用Welch方法估算功率谱密度。

韦尔奇的方法[1]通过将数据划分为重叠的段，计算每个段的修改后的周期图，并对周期图取平均值，可以计算出功率谱密度的估计值。

参数：

x：array_like测量值的时间序列

fs：float, 可选参数x时间序列的采样频率。默认为1.0。

window：str 或 tuple 或 array_like, 可选参数希望使用的窗口。如果window是字符串或元组，则将其传递给get_window生成窗口值，默认情况下为DFT-even。参考get_window有关窗口和必需参数的列表。如果窗口是数组，它将直接用作窗口，并且其长度必须为nperseg。默认为Hann窗口。

nperseg：int, 可选参数每个段的长度。默认为无，但是如果window是str或tuple，则设置为256，如果window是数组，则设置为窗口的长度。

noverlap：int, 可选参数段之间重叠的点数。如果没有，noverlap = nperseg // 2。默认为没有。

nfft：int, 可选参数如果需要零填充的FFT，则使用的FFT的长度。如果为None，则FFT长度为nperseg。默认为无。

detrend：str 或 function 或 False, 可选参数指定如何使每个段趋势消失。如果detrend是一个字符串，它将作为类型参数传递给detrend函数。如果它是一个函数，它将采用一个段并返回一个去趋势的段。如果detrend是假，不进行趋势消除。默认为‘constant’。

return_onesided：bool, 可选参数如果为True，则返回真实数据的one-sided频谱。如果为False，则返回two-sided频谱。默认为True，但是对于复杂数据，始终返回two-sided频谱。

scaling：{ ‘density’, ‘spectrum’ }, 可选参数如果x的单位为V且fs的单位为fs，则在计算Pxx的单位为V ** 2 /Hz的功率谱密度(‘density’)和计算Pxx的单位为V ** 2的功率谱(‘spectrum’)之间进行选择。赫兹。默认为‘density’

axis：int, 可选参数计算周期图的轴；默认值位于最后一个轴上(即axis=-1)。

average：{ ‘mean’, ‘median’ }, 可选参数平均周期图时使用的方法。默认为‘mean’。

1.2.0版的新函数。

返回值：

f：ndarray采样频率数组。

Pxx：ndarrayx的功率谱密度或功率谱。

注意：

适当的重叠量将取决于窗口的选择和您的要求。对于默认的Hann窗口，重叠50％是在准确估计信号功率与不过度计算任何数据之间的合理折衷。较窄的窗口可能需要较大的重叠。

如果noverlap为0，则此方法等效于Bartlett的方法[2]。

版本0.12.0中的新函数。

参考文献：

P. Welch，“使用快速傅立叶变换估计功率谱：一种基于短修正周期图上时间平均的方法”，IEEE Trans。音频电声。卷15，第70-73页，1967年。

女士。 Bartlett，“周期图分析和连续光谱”，Biometrika，第一卷。 37卷，第1-16页，1950年。

例子：

>>> from scipy import signal

>>> import matplotlib.pyplot as plt

>>> np.random.seed(1234)

产生一个测试信号，即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 += np.random.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()