SciPy 1.12 中文文档(三十三)
scipy.signal.square
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.square.html#scipy.signal.square
scipy.signal.square(t, duty=0.5)
返回周期性方波波形。
方波的周期为2*pi,在0到2*pi*duty之间取值为+1,在2*pi*duty到2*pi之间取值为-1。duty必须在区间[0,1]内。
请注意,此波形不是带限制的。它产生无限多个谐波,这些谐波在频谱上来回混叠。
参数:
tarray_like
输入时间数组。
占空比array_like,可选
占空比。默认为 0.5(50%占空比)。如果是数组,则导致波形随时间变化,并且必须与 t 具有相同的长度。
返回:
yndarray
输出包含方波波形的数组。
示例
一个 5 Hz 波形,以 500 Hz 采样,持续 1 秒钟:
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(0, 1, 500, endpoint=False)
>>> plt.plot(t, signal.square(2 * np.pi * 5 * t))
>>> plt.ylim(-2, 2)
一个脉宽调制的正弦波:
>>> plt.figure()
>>> sig = np.sin(2 * np.pi * t)
>>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
>>> plt.subplot(2, 1, 1)
>>> plt.plot(t, sig)
>>> plt.subplot(2, 1, 2)
>>> plt.plot(t, pwm)
>>> plt.ylim(-1.5, 1.5)
scipy.signal.sweep_poly
scipy.signal.sweep_poly(t, poly, phi=0)
频率扫描余弦生成器,带有时间依赖的频率。
此函数生成一个正弦函数,其即时频率随时间变化。时间 t 处的频率由多项式 poly 给出。
参数:
tndarray
评估波形的时间点。
poly1-D 数组或者是 numpy.poly1d 的实例
所需频率表示为一个多项式。如果 poly 是长度为 n 的列表或 ndarray,则 poly 的元素为多项式的系数,即时频率为
f(t) = poly[0]*t**(n-1) + poly[1]*t**(n-2) + ... + poly[n-1]
如果 poly 是 numpy.poly1d 的实例,则即时频率为
f(t) = poly(t)
phi浮点数,可选
相位偏移,以度为单位,默认为 0。
返回:
sweep_polyndarray
包含信号在 t 处评估的 numpy 数组,具有请求的时间变化频率。更确切地说,函数返回 cos(phase + (pi/180)*phi),其中 phase 是积分(从 0 到 t)的 2 * pi * f(t);f(t) 如上所定义。
另请参阅
chirp
注意事项
自 0.8.0 版本开始引入。
如果 poly 是长度为 n 的列表或 ndarray,则 poly 的元素为多项式的系数,即时频率为:
f(t) = poly[0]*t**(n-1) + poly[1]*t**(n-2) + ... + poly[n-1]
如果 poly 是 numpy.poly1d 的实例,则即时频率为:
f(t) = poly(t)
最后,输出 s 为:
cos(phase + (pi/180)*phi)
其中 phase 是从 0 到 t 的积分,式子为 2 * pi * f(t),其中 f(t) 如上所定义。
示例
计算具有即时频率的波形:
f(t) = 0.025*t**3 - 0.36*t**2 + 1.25*t + 2
在 0 <= t <= 10 的区间内。
>>> import numpy as np
>>> from scipy.signal import sweep_poly
>>> p = np.poly1d([0.025, -0.36, 1.25, 2.0])
>>> t = np.linspace(0, 10, 5001)
>>> w = sweep_poly(t, p)
绘制它:
>>> import matplotlib.pyplot as plt
>>> plt.subplot(2, 1, 1)
>>> plt.plot(t, w)
>>> plt.title("Sweep Poly\nwith frequency " +
... "$f(t) = 0.025t³ - 0.36t² + 1.25t + 2$")
>>> plt.subplot(2, 1, 2)
>>> plt.plot(t, p(t), 'r', label='f(t)')
>>> plt.legend()
>>> plt.xlabel('t')
>>> plt.tight_layout()
>>> plt.show()
scipy.signal.unit_impulse
scipy.signal.unit_impulse(shape, idx=None, dtype=<class 'float'>)
单位脉冲信号(离散δ函数)或单位基向量。
参数:
shape整数或整数元组
输出中的样本数量(1 维),或者表示输出形状的元组(N 维)。
idxNone 或整数或整数元组或‘mid’,可选
值为 1 的索引位置。如果为 None,则默认为第 0 个元素。如果idx='mid',则脉冲信号将在所有维度上居中于shape // 2。如果为整数,则脉冲信号将在所有维度上位于idx。
dtype数据类型,可选
数组的期望数据类型,例如,numpy.int8。默认为numpy.float64。
返回:
yndarray
输出数组,包含脉冲信号。
注意
1 维情况也称为 Kronecker delta。
新版本 0.19.0 中新增。
示例
一个在第 0 个元素处的脉冲信号((\delta[n])):
>>> from scipy import signal
>>> signal.unit_impulse(8)
array([ 1., 0., 0., 0., 0., 0., 0., 0.])
脉冲信号偏移了 2 个样本((\delta[n-2])):
>>> signal.unit_impulse(7, 2)
array([ 0., 0., 1., 0., 0., 0., 0.])
二维脉冲信号,居中:
>>> signal.unit_impulse((3, 3), 'mid')
array([[ 0., 0., 0.],
[ 0., 1., 0.],
[ 0., 0., 0.]])
在(2, 2)处的脉冲信号,使用广播:
>>> signal.unit_impulse((4, 4), 2)
array([[ 0., 0., 0., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 1., 0.],
[ 0., 0., 0., 0.]])
绘制 4 阶 Butterworth 低通滤波器的脉冲响应:
>>> imp = signal.unit_impulse(100, 'mid')
>>> b, a = signal.butter(4, 0.2)
>>> response = signal.lfilter(b, a, imp)
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> plt.plot(np.arange(-50, 50), imp)
>>> plt.plot(np.arange(-50, 50), response)
>>> plt.margins(0.1, 0.1)
>>> plt.xlabel('Time [samples]')
>>> plt.ylabel('Amplitude')
>>> plt.grid(True)
>>> plt.show()
scipy.signal.get_window
scipy.signal.get_window(window, Nx, fftbins=True)
返回给定长度和类型的窗口。
参数:
window字符串、浮点数或元组
要创建的窗口类型。详见下文。
Nxint
窗口中的样本数。
fftbinsbool,可选
如果为 True(默认),创建一个“周期性”窗口,准备用于 ifftshift 并乘以 FFT 的结果(还请参阅 fftfreq)。如果为 False,创建一个“对称”窗口,用于滤波器设计。
返回:
get_windowndarray
返回长度为 Nx、类型为 window 的窗口
注意
窗口类型:
-
kaiser(需要 beta) -
kaiser_bessel_derived(需要 beta) -
gaussian(需要标准差) -
general_cosine(需要加权系数) -
general_gaussian(需要功率、宽度) -
general_hamming(需要窗口系数) -
dpss(需要归一化半带宽) -
chebwin(需要衰减)
如果窗口不需要参数,则window可以是一个字符串。
如果窗口需要参数,则window必须是一个元组,第一个参数是窗口的字符串名称,后续参数是所需的参数。
如果window是一个浮点数,则被解释为kaiser窗口的β参数。
上述每种窗口类型也是一个可以直接调用以创建该类型窗口的函数名称。
示例
>>> from scipy import signal
>>> signal.get_window('triang', 7)
array([ 0.125, 0.375, 0.625, 0.875, 0.875, 0.625, 0.375])
>>> signal.get_window(('kaiser', 4.0), 9)
array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093,
0.97885093, 0.82160913, 0.56437221, 0.29425961])
>>> signal.get_window(('exponential', None, 1.), 9)
array([ 0.011109 , 0.03019738, 0.082085 , 0.22313016, 0.60653066,
0.60653066, 0.22313016, 0.082085 , 0.03019738])
>>> signal.get_window(4.0, 9)
array([ 0.08848053, 0.29425961, 0.56437221, 0.82160913, 0.97885093,
0.97885093, 0.82160913, 0.56437221, 0.29425961])
scipy.signal.cascade
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.cascade.html#scipy.signal.cascade
scipy.signal.cascade(hk, J=7)
从滤波器系数计算出在二分点K/2**J处的(x, phi, psi)。
自 1.12.0 版本起已弃用:scipy.signal.cascade 在 SciPy 1.12 中已弃用,将在 SciPy 1.15 中移除。我们建议改用 PyWavelets。
参数:
hk数组型
低通滤波器的系数。
J 整型,可选
值将在网格点K/2**J处计算。默认值为 7。
返回:
x 数组型
对于K=0...N * (2**J)-1,K/2**J是二分点,其中len(hk) = len(gk) = N+1。
phi 数组型
缩放函数phi(x)在x处的定义为:phi(x) = sum(hk * phi(2x-k)),其中k的取值范围是从 0 到N。
psi 数组型,可选
小波函数psi(x)在x处的定义为:phi(x) = sum(gk * phi(2x-k)),其中k的取值范围是从 0 到N。当gk不为 None 时才返回psi。
注解
算法使用 Strang 和 Nguyen 在《小波与滤波器组》中描述的向量级联算法。它构建一个值和切片的字典以便快速重用。然后在最后将向量插入到最终向量中。
scipy.signal.daub
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.daub.html#scipy.signal.daub
scipy.signal.daub(p)
生成 Daubechies 小波的 FIR 低通滤波器的系数。
自版本 1.12.0 起已弃用:scipy.signal.daub 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中移除。我们建议改用 PyWavelets。
当 p>=1 时,这是在 f=1/2 处的零点的阶数。有 2p 个滤波器系数。
参数:
pint
在 f=1/2 处的零点的阶数,可以取 1 到 34 的值。
返回:
daubndarray
返回:
scipy.signal.morlet
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.morlet.html#scipy.signal.morlet
scipy.signal.morlet(M, w=5.0, s=1.0, complete=True)
复数 Morlet 小波。
自版本 1.12.0 起不建议使用:scipy.signal.morlet 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中移除。我们建议改用 PyWavelets。
参数:
M整数
小波的长度。
w浮点数,可选
Omega0. 默认值为 5。
s浮点数,可选
缩放因子,从-s*2*pi到+s*2*pi窗口化。默认值为 1。
complete布尔值,可选
是否使用完整版或标准版。
返回:
morlet(M,) ndarray
另请参阅
morlet2
Morlet 小波的实现,与cwt兼容。
scipy.signal.gausspulse
注意事项
标准版本:
pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))
这种常用的小波通常被简称为 Morlet 小波。请注意,这个简化版本在w的低值时可能会导致可接受性问题。
完整版本:
pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))
此版本具有改正项以改善可接受性。对于w大于 5,改正项可忽略不计。
注意,返回小波的能量未根据s进行标准化。
此小波的基本频率(以 Hz 为单位)由f = 2*s*w*r / M给出,其中r是采样率。
注意:此函数在cwt之前创建,与其不兼容。
示例
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> M = 100
>>> s = 4.0
>>> w = 2.0
>>> wavelet = signal.morlet(M, s, w)
>>> plt.plot(wavelet.real, label="real")
>>> plt.plot(wavelet.imag, label="imag")
>>> plt.legend()
>>> plt.show()
scipy.signal.qmf
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.qmf.html#scipy.signal.qmf
scipy.signal.qmf(hk)
返回高通 qmf 滤波器自低通
自 1.12.0 版本起不推荐使用:scipy.signal.qmf 在 SciPy 1.12 中已被弃用,并将在 SciPy 1.15 中移除。我们建议改用 PyWavelets 替代。
参数:
hkarray_like
高通滤波器的系数。
返回:
array_like
高通滤波器系数。
scipy.signal.ricker
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.ricker.html#scipy.signal.ricker
scipy.signal.ricker(points, a)
返回一个 Ricker 小波,也称为“墨西哥帽小波”。
自 SciPy 1.12 版本起已弃用:scipy.signal.ricker 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中移除。我们建议改用 PyWavelets。
它模拟函数:
A * (1 - (x/a)**2) * exp(-0.5*(x/a)**2),
其中 A = 2/(sqrt(3*a)*(pi**0.25))。
参数:
points整数
vector 中的点数。将以 0 为中心。
a标量
小波的宽度参数。
返回值:
向量(N,) ndarray
形状为 ricker 曲线的长度为 points 的数组。
示例
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> points = 100
>>> a = 4.0
>>> vec2 = signal.ricker(points, a)
>>> print(len(vec2))
100
>>> plt.plot(vec2)
>>> plt.show()
scipy.signal.morlet2
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.morlet2.html#scipy.signal.morlet2
scipy.signal.morlet2(M, s, w=5)
复杂的莫尔雷特小波,设计用于与cwt配合使用。
自 SciPy 1.12 版本起弃用:scipy.signal.morlet2 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中移除。我们建议改用 PyWavelets。
返回归一化后的完整莫尔雷特小波,根据s进行归一化:
exp(1j*w*x/s) * exp(-0.5*(x/s)**2) * pi**(-0.25) * sqrt(1/s)
参数:
Mint
小波的长度。
sfloat
小波的宽度参数。
wfloat, optional
Omega0. 默认值为 5
返回:
morlet(M,) ndarray
另请参阅
morlet
莫尔雷特小波的实现,与cwt不兼容
注意事项
新功能 1.4.0 版。
此函数设计用于与cwt配合使用。因为morlet2返回一个复数数组,所以最好将cwt的dtype参数设置为complex128以获得最佳结果。
注意与morlet实现上的差异。该小波的基频(单位:Hz)由以下公式给出:
f = w*fs / (2*s*np.pi)
其中fs为采样率,s为小波宽度参数。类似地,我们可以在f处得到小波宽度参数:
s = w*fs / (2*f*np.pi)
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> M = 100
>>> s = 4.0
>>> w = 2.0
>>> wavelet = signal.morlet2(M, s, w)
>>> plt.plot(abs(wavelet))
>>> plt.show()
此示例展示了在时间频率分析中使用morlet2与cwt的基本用法:
>>> t, dt = np.linspace(0, 1, 200, retstep=True)
>>> fs = 1/dt
>>> w = 6.
>>> sig = np.cos(2*np.pi*(50 + 10*t)*t) + np.sin(40*np.pi*t)
>>> freq = np.linspace(1, fs/2, 100)
>>> widths = w*fs / (2*freq*np.pi)
>>> cwtm = signal.cwt(sig, signal.morlet2, widths, w=w)
>>> plt.pcolormesh(t, freq, np.abs(cwtm), cmap='viridis', shading='gouraud')
>>> plt.show()
scipy.signal.cwt
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.cwt.html#scipy.signal.cwt
scipy.signal.cwt(data, wavelet, widths, dtype=None, **kwargs)
连续小波变换。
自版本 1.12.0 起弃用:scipy.signal.cwt 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中删除。我们建议改用 PyWavelets。
对data执行连续小波变换,使用wavelet函数。 CWT 使用wavelet函数对data进行卷积,该函数以宽度参数和长度参数为特征。 wavelet函数允许是复数。
参数:
data(N,) ndarray
执行变换的数据。
wavelet函数
小波函数,应该接受 2 个参数。第一个参数是返回的向量将具有的点数(len(wavelet(length,width)) == length)。第二个是宽度参数,定义小波的大小(例如,高斯标准差)。参见ricker,满足这些要求。
widths(M,) 序列
用于变换的宽度。
dtype数据类型,可选
输出的期望数据类型。如果wavelet的输出是实数,则默认为float64,如果是复数,则为complex128。
新版本 1.4.0 中新增。
kwargs
传递给小波函数的关键字参数。
新版本 1.4.0 中新增。
返回:
cwt:(M, N) ndarray
将具有形状(len(widths), len(data))。
注意事项
新版本 1.4.0 中新增。
对于非对称的复数值小波,输入信号与小波数据的时间反转共轭卷积[1]。
length = min(10 * width[ii], len(data))
cwt[ii,:] = signal.convolve(data, np.conj(wavelet(length, width[ii],
**kwargs))[::-1], mode='same')
参考资料
[1]
S. Mallat,“信号处理的小波之旅(第 3 版)”,Academic Press,2009。
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(-1, 1, 200, endpoint=False)
>>> sig = np.cos(2 * np.pi * 7 * t) + signal.gausspulse(t - 0.4, fc=2)
>>> widths = np.arange(1, 31)
>>> cwtmatr = signal.cwt(sig, signal.ricker, widths)
注意
对于 cwt 矩阵绘图,建议翻转 y 轴
>>> cwtmatr_yflip = np.flipud(cwtmatr)
>>> plt.imshow(cwtmatr_yflip, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
... vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
>>> plt.show()
scipy.signal.cwt
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.cwt.html#scipy.signal.cwt
scipy.signal.cwt(data, wavelet, widths, dtype=None, **kwargs)
连续小波变换。
自版本 1.12.0 起已弃用:scipy.signal.cwt 在 SciPy 1.12 中已弃用,并将在 SciPy 1.15 中删除。 我们建议改用 PyWavelets。
使用data和wavelet函数执行连续小波变换。 CWT 通过wavelet函数对data进行卷积,wavelet函数具有宽度参数和长度参数。 wavelet函数可以是复数。
参数:
data(N,) ndarray
执行变换的数据。
wavelet函数
小波函数,应该接受 2 个参数。 第一个参数是返回向量将具有的点数(len(wavelet(length,width)) == length)。 第二个是宽度参数,定义小波的大小(例如高斯分布的标准偏差)。 参见ricker,它满足这些要求。
widths(M,) 序列
用于变换的宽度。
dtype数据类型,可选
输出的期望数据类型。 如果wavelet的输出是实数,则默认为float64;如果是复数,则为complex128。
新版本 1.4.0 中新增。
kwargs
传递给小波函数的关键字参数。
新版本 1.4.0 中新增。
返回:
cwt:(M, N) ndarray
形状为(len(widths), len(data))。
注意
新版本 1.4.0 中新增。
对于非对称的复数值小波,输入信号与小波数据的时间反转复共轭进行卷积 [1]。
length = min(10 * width[ii], len(data))
cwt[ii,:] = signal.convolve(data, np.conj(wavelet(length, width[ii],
**kwargs))[::-1], mode='same')
参考文献
[1]
S. Mallat,《信号处理的小波之旅(第 3 版)》,Academic Press,2009 年。
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(-1, 1, 200, endpoint=False)
>>> sig = np.cos(2 * np.pi * 7 * t) + signal.gausspulse(t - 0.4, fc=2)
>>> widths = np.arange(1, 31)
>>> cwtmatr = signal.cwt(sig, signal.ricker, widths)
注意
对于 cwt 矩阵绘图,建议翻转 y 轴。
>>> cwtmatr_yflip = np.flipud(cwtmatr)
>>> plt.imshow(cwtmatr_yflip, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
... vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
>>> plt.show()
scipy.signal.argrelmin
scipy.signal.argrelmin(data, axis=0, order=1, mode='clip')
计算data的相对最小值。
参数:
datandarray
用于查找相对最小值的数组。
axisint,可选
选择从data中选取的轴。默认为 0。
orderint,可选
在每一侧用于比较的点数以便认为comparator(n, n+x)为 True。
modestr,可选
指定向量边缘的处理方式。可用选项为'wrap'(环绕)或'clip'(将溢出视为最后(或第一个)元素)。默认为'clip'。参见 numpy.take。
返回:
extremandarray 的元组
整数数组中的最小值的索引。extrema[k]是data的轴k的索引数组。请注意,即使data是 1-D,返回值也是元组。
另请参阅
argrelextrema,argrelmax,find_peaks
注意
此函数使用argrelextrema作为比较器的 np.less。因此,它要求在值的两侧都严格使用不等号才能将其视为最小值。这意味着平坦的最小值(多于一个样本宽度)不会被检测到。在 1-D data的情况下,可以通过使用反向的data调用find_peaks来检测所有本地最小值,包括平坦的最小值。
0.11.0 版本中新增。
示例
>>> import numpy as np
>>> from scipy.signal import argrelmin
>>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelmin(x)
(array([1, 5]),)
>>> y = np.array([[1, 2, 1, 2],
... [2, 2, 0, 0],
... [5, 3, 4, 4]])
...
>>> argrelmin(y, axis=1)
(array([0, 2]), array([2, 1]))
scipy.signal.argrelmax
scipy.signal.argrelmax(data, axis=0, order=1, mode='clip')
计算 data 的相对最大值。
参数:
datandarray
要在其中查找相对最大值的数组。
axisint,可选
用于从 data 中选择的轴。默认为 0。
orderint,可选
每侧使用多少点进行比较,以确定 comparator(n, n+x) 是否为真。
modestr,可选
如何处理向量的边缘。可用选项为‘wrap’(环绕)或‘clip’(将溢出视为与最后(或第一个)元素相同)。默认为‘clip’。参见numpy.take。
返回:
extremandarray 的元组
整数数组中极大值的索引。extrema[k] 是 data 的轴 k 的索引数组。注意,即使 data 是 1-D,返回值也是元组。
另见
argrelextrema,argrelmin,find_peaks
注意
此函数使用 argrelextrema 作为 np.greater 的比较器。因此,它要求在值的两侧都有严格的不等式才能将其视为最大值。这意味着平坦的最大值(多于一个样本宽度)不会被检测到。在 1-D data 的情况下,可以使用 find_peaks 来检测所有本地最大值,包括平坦的最大值。
从版本 0.11.0 开始。
示例
>>> import numpy as np
>>> from scipy.signal import argrelmax
>>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelmax(x)
(array([3, 6]),)
>>> y = np.array([[1, 2, 1, 2],
... [2, 2, 0, 0],
... [5, 3, 4, 4]])
...
>>> argrelmax(y, axis=1)
(array([0]), array([1]))
scipy.signal.argrelextrema
scipy.signal.argrelextrema(data, comparator, axis=0, order=1, mode='clip')
计算data的相对极值。
参数:
data:ndarray
要查找相对极值的数组。
comparator:callable
用于比较两个数据点的函数。应接受两个数组作为参数。
axis:int,可选
选择data的轴。默认为 0。
order:int,可选
用于比较comparator(n, n+x)是否为 True 时每侧要使用的点数。
mode:str,可选
向量边缘的处理方式。‘wrap’(环绕)或‘clip’(将溢出视为与最后(或第一个)元素相同)。默认为‘clip’。参见numpy.take.
返回值:
extrema:ndarrays 的元组
整数数组中的极大值的索引。extrema[k]是data的轴k的索引数组。请注意,即使data是 1-D,返回值也是元组。
参见
argrelmin, argrelmax
注意事项
自版本 0.11.0 新增。
示例
>>> import numpy as np
>>> from scipy.signal import argrelextrema
>>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
>>> argrelextrema(x, np.greater)
(array([3, 6]),)
>>> y = np.array([[1, 2, 1, 2],
... [2, 2, 0, 0],
... [5, 3, 4, 4]])
...
>>> argrelextrema(y, np.less, axis=1)
(array([0, 2]), array([2, 1]))
scipy.signal.find_peaks
scipy.signal.find_peaks(x, height=None, threshold=None, distance=None, prominence=None, width=None, wlen=None, rel_height=0.5, plateau_size=None)
根据峰值属性查找信号内的峰值。
此函数接受一个一维数组,并通过简单比较相邻值来找到所有局部最大值。可选地,可以通过指定峰值属性的条件来选择其中的一部分峰值。
参数:
x序列
带有峰值的信号。
height数字或数组或序列,可选
峰值的所需高度。可以是一个数字、None、与x匹配的数组或前述的两个元素的序列。第一个元素始终解释为最小值,如果提供第二个元素,则为最大所需高度。
threshold数字或数组或序列,可选
峰值的所需阈值,与其相邻样本的垂直距离。可以是一个数字、None、与x匹配的数组或前述的两个元素的序列。第一个元素始终解释为最小值,如果提供第二个元素,则为最大所需阈值。
distance数字,可选
相邻峰值之间的必需最小水平距离(>= 1)(以样本为单位)。直到所有剩余的峰值满足条件之前,较小的峰值会被首先移除。
prominence数字或数组或序列,可选
峰值的所需显著性。可以是一个数字、None、与x匹配的数组或前述的两个元素的序列。第一个元素始终解释为最小值,如果提供第二个元素,则为最大所需显著性。
width数字或数组或序列,可选
峰值的所需宽度(以样本为单位)。可以是一个数字、None、与x匹配的数组或前述的两个元素的序列。第一个元素始终解释为最小值,如果提供第二个元素,则为最大所需宽度。
wlen整数,可选
用于计算峰值显著性,因此只有在给定prominence或width之一的参数时才会使用。有关其效果的详细描述,请参见peak_prominences中的参数wlen。
rel_height浮点数,可选
用于计算峰值宽度,因此只有在给定width参数时才会使用。有关其效果的详细描述,请参见peak_widths中的参数rel_height。
plateau_size数字或数组或序列,可选
峰顶的所需平坦顶部大小(以样本为单位)。可以是一个数字、None、与x匹配的数组或前述的两个元素的序列。第一个元素始终解释为最小值,如果提供第二个元素,则为最大所需平顶大小。
1.2.0 版本中的新功能。
返回:
peaks数组
满足所有给定条件的x中的峰值的索引。
propertiesdict
包含在指定条件评估过程中计算的返回峰值的属性的字典:
-
‘peak_heights’
如果给定height,则为x中每个峰的高度。
-
‘left_thresholds’、‘right_thresholds’
如果给定threshold,则这些键包含峰值与其相邻样本的垂直距离。
-
‘prominences’、‘right_bases’、‘left_bases’
如果给定prominence,则可以访问这些键。详见
peak_prominences以获取其内容的描述。 -
‘width_heights’、‘left_ips’、‘right_ips’
如果给定width,则可以访问这些键。详见
peak_widths以获取其内容的描述。 -
‘plateau_sizes’、‘left_edges’、‘right_edges’
如果给定plateau_size,则可以访问这些键,并包含峰的边缘(边缘仍然是平台的一部分)的索引和计算的平台大小。
新版本 1.2.0 中提供。
若要计算并返回不排除峰值的属性,请将开放区间(None, None)作为适当参数的值(不包括distance)。
警告:
PeakPropertyWarning
如果峰值的属性具有意外的值(参见peak_prominences和peak_widths),则会引发此警告。
警告
对于包含 NaN 的数据,此函数可能返回意外结果。为避免此问题,应删除或替换 NaN。
另见
find_peaks_cwt
使用小波变换查找峰值。
peak_prominences
直接计算峰的显著性。
peak_widths
直接计算峰的宽度。
注释
在此函数的上下文中,峰值或局部最大值定义为任何两个直接相邻样本其振幅较小。对于平顶峰(多于一个相等振幅的样本宽度),返回中间样本的索引(如果样本数为偶数则向下取整)。对于噪声信号,峰位置可能会偏移,因为噪声可能会改变局部最大值的位置。在这些情况下,考虑在搜索峰值之前对信号进行平滑处理或使用其他峰值查找和拟合方法(如find_peaks_cwt)。
关于指定条件的一些额外评论:
-
几乎所有条件(除了 distance)都可以给出半开或闭区间,例如,
1或(1, None)定义了半开区间 ([1, \infty]),而(None, 1)定义了区间 ([-\infty, 1])。开区间(None, None)也可以被指定,返回匹配的属性而不排除峰值。 -
边界始终包含在用于选择有效峰值的区间中。
-
对于几个条件,区间边界可以用与 x 形状匹配的数组指定,从而基于样本位置实现动态约束。
-
条件按以下顺序进行评估:plateau_size、height、threshold、distance、prominence、width。在大多数情况下,这个顺序是最快的,因为会先应用更快的操作,以减少后续需要评估的峰值数量。
-
虽然 peaks 中的索引保证至少相隔 distance 个样本,但平坦峰的边缘可能比允许的 distance 更近。
-
如果 x 较大或有许多局部最大值,可以使用 wlen 减少评估 prominence 或 width 条件所需的时间(参见
peak_prominences)。
新功能在版本 1.1.0 中引入。
示例
为了演示这个函数的使用,我们使用了 SciPy 提供的信号 x(参见scipy.datasets.electrocardiogram)。让我们找出所有幅度大于 0 的 x 中的峰值(局部最大值)。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.datasets import electrocardiogram
>>> from scipy.signal import find_peaks
>>> x = electrocardiogram()[2000:4000]
>>> peaks, _ = find_peaks(x, height=0)
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.plot(np.zeros_like(x), "--", color="gray")
>>> plt.show()
我们可以使用 height=(None, 0) 或使用与 x 大小匹配的数组来反映不同信号部分的变化条件。
>>> border = np.sin(np.linspace(0, 3 * np.pi, x.size))
>>> peaks, _ = find_peaks(x, height=(-border, border))
>>> plt.plot(x)
>>> plt.plot(-border, "--", color="gray")
>>> plt.plot(border, ":", color="gray")
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()
对于周期信号,另一个有用的条件可以用 distance 参数给出。在这种情况下,我们可以通过要求至少 150 个样本的距离轻松地选择心电图(ECG)中的 QRS 复合体的位置。
>>> peaks, _ = find_peaks(x, distance=150)
>>> np.diff(peaks)
array([186, 180, 177, 171, 177, 169, 167, 164, 158, 162, 172])
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()
特别是对于嘈杂信号,可以通过它们的显著性轻松地分组峰值(参见peak_prominences)。例如,我们可以通过将允许的显著性限制为 0.6 来选择除了上述 QRS 复合体之外的所有峰值。
>>> peaks, properties = find_peaks(x, prominence=(None, 0.6))
>>> properties["prominences"].max()
0.5049999999999999
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()
最后,让我们检查包含不同形状节拍的 ECG 的不同部分。为了仅选择非典型心跳,我们结合两个条件:至少 1 的最小显著性和至少 20 个样本的宽度。
>>> x = electrocardiogram()[17000:18000]
>>> peaks, properties = find_peaks(x, prominence=1, width=20)
>>> properties["prominences"], properties["widths"]
(array([1.495, 2.3 ]), array([36.93773946, 39.32723577]))
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.vlines(x=peaks, ymin=x[peaks] - properties["prominences"],
... ymax = x[peaks], color = "C1")
>>> plt.hlines(y=properties["width_heights"], xmin=properties["left_ips"],
... xmax=properties["right_ips"], color = "C1")
>>> plt.show()
scipy.signal.find_peaks_cwt
scipy.signal.find_peaks_cwt(vector, widths, wavelet=None, max_distances=None, gap_thresh=None, min_length=None, min_snr=1, noise_perc=10, window_size=None)
使用小波变换在一维数组中找到峰值。
一般方法是通过将向量与每个宽度中的小波(width)卷积来平滑向量。足够长的多尺度上出现的相对最大值,并且具有足够高的信噪比,将被接受。
参数:
向量 ndarray
要在其中找到峰值的一维数组。
宽度浮点数或序列
单个宽度或用于计算 CWT 矩阵的一维类似宽度数组。一般来说,这个范围应该覆盖感兴趣峰值的预期宽度。
小波可调用函数,可选项
应接受两个参数并返回与向量卷积的一维数组。第一个参数确定返回的小波数组的点数,第二个参数是小波的尺度(宽度)。应该是归一化和对称的。默认为里克小波。
最大距离 ndarray,可选项
在每一行,只有当在row[n]处的相对最大值与row[n+1]处的相对最大值在max_distances[n]内时,才连接一条脊线。默认值为widths/4。
间隙阈值浮点数,可选项
如果在max_distances内找不到相对最大值,则会有一个间隙。如果有超过gap_thresh个点而不连接新的相对最大值,则脊线被中断。默认值是宽度数组的第一个值,即 widths[0]。
最小长度整数,可选项
脊线需要接受的最小长度。默认为cwt.shape[0] / 4,即宽度的四分之一。
最小信噪比浮点数,可选项
最小信噪比。默认值为 1。信号是最大的 CWT 系数在最大脊线上。噪声是noise_perc百分位数的数据点,这些数据点包含在同一脊线内。
噪声百分比浮点数,可选项
在计算噪声底线时,百分位数的数据点低于这个值被认为是噪声。使用stats.scoreatpercentile计算。默认值为 10。
窗口大小整数,可选项
用于计算噪声底线的窗口大小。默认值为cwt.shape[1] / 20。
返回:
峰值索引 ndarray
找到峰值的向量中的位置的索引。列表已排序。
另见
cwt
连续小波变换。
find_peaks
根据峰值属性在信号内部找到峰值。
笔记
此方法旨在从嘈杂数据中找出尖峰,但通过适当的参数选择,它应该能够很好地适应不同的峰形状。
算法如下:
-
对 vector 执行连续小波变换,使用提供的 widths。这是 vector 与每个 widths 中的 wavelet(width) 的卷积。参见
cwt。 -
在 cwt 矩阵中识别“脊线”。这些是每行的相对最大值,在相邻行之间连接。参见 identify_ridge_lines
-
使用
filter_ridge_lines过滤脊线。
新功能在版本 0.11.0 中引入。
参考
[1]
生物信息学(2006)22(17):2059-2065。DOI:10.1093/bioinformatics/btl355
示例
>>> import numpy as np
>>> from scipy import signal
>>> xs = np.arange(0, np.pi, 0.05)
>>> data = np.sin(xs)
>>> peakind = signal.find_peaks_cwt(data, np.arange(1,10))
>>> peakind, xs[peakind], data[peakind]
([32], array([ 1.6]), array([ 0.9995736]))
scipy.signal.peak_prominences
scipy.signal.peak_prominences(x, peaks, wlen=None)
计算信号中每个峰值的显著性。
峰值的显著性衡量了峰值在信号周围基线的突出程度,并定义为峰值与其最低轮廓线之间的垂直距离。
参数:
x序列
一个具有峰值的信号。
peaks序列
x中的峰值索引。
wlenint,可选
采样中的窗口长度,可选地限制每个峰值的评估区域为x的子集。峰值始终位于窗口的中间,因此给定的长度会向上舍入为下一个奇数整数。此参数可以加速计算(见注释)。
返回:
prominencesndarray
每个peaks中的峰值的计算显著性。
left_bases, right_basesndarray
峰值的基准作为x中每个峰值左右的索引。每对中较高的基准是峰值的最低轮廓线。
Raises:
ValueError
如果peaks中的值是x的无效索引。
警告:
PeakPropertyWarning
对于peaks中不指向x中有效局部最大值的索引,返回的显著性将为 0,并引发此警告。如果wlen小于峰值的平台大小,则也会发生这种情况。
警告
对于包含 NaN 的数据,此函数可能返回意外的结果。为避免此情况,应移除或替换 NaN。
另请参阅
find_peaks
根据峰值属性在信号内查找峰值。
peak_widths
计算峰值的宽度。
注释
计算峰值显著性的策略:
-
从当前峰值水平线向左右延伸,直到该线到达窗口边界(参见wlen)或再次在较高峰值的斜率处与信号相交。与相同高度的峰值的交叉将被忽略。
-
在每侧找到定义的区间内的最小信号值。这些点是峰值的基准。
-
两个基准中较高的一个标记了峰值的最低轮廓线。然后可以计算显著性,作为峰值本身高度与其最低轮廓线之间的垂直差异。
对于具有周期性行为的大x,寻找峰值基线可能会很慢,因为需要评估大块或甚至整个信号作为第一个算法步骤。可以使用参数wlen来限制评估区域,该参数将算法限制在当前峰值周围的窗口内,并且如果窗口长度相对于x较短,则可以缩短计算时间。然而,如果峰值的真实基线超出此窗口,则可能阻止算法找到真正的全局等高线。相反,会在限制的窗口内找到一个更高的等高线,导致计算出的突出度较小。实际上,这仅对x中最高一组峰值相关。此行为甚至可能会被有意用来计算“局部”突出度。
新版本 1.1.0 中的内容。
参考资料
[1]
维基百科关于地形突出度的文章: zh.wikipedia.org/wiki/地形突出度
示例
>>> import numpy as np
>>> from scipy.signal import find_peaks, peak_prominences
>>> import matplotlib.pyplot as plt
创建一个带有两个重叠谐波的测试信号
>>> x = np.linspace(0, 6 * np.pi, 1000)
>>> x = np.sin(x) + 0.6 * np.sin(2.6 * x)
查找所有峰值并计算突出度
>>> peaks, _ = find_peaks(x)
>>> prominences = peak_prominences(x, peaks)[0]
>>> prominences
array([1.24159486, 0.47840168, 0.28470524, 3.10716793, 0.284603 ,
0.47822491, 2.48340261, 0.47822491])
计算每个峰值的等高线高度并绘制结果
>>> contour_heights = x[peaks] - prominences
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.vlines(x=peaks, ymin=contour_heights, ymax=x[peaks])
>>> plt.show()
让我们评估第二个示例,该示例演示了索引为 5 的一个峰值的几种边缘情况。
>>> x = np.array([0, 1, 0, 3, 1, 3, 0, 4, 0])
>>> peaks = np.array([5])
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()
>>> peak_prominences(x, peaks) # -> (prominences, left_bases, right_bases)
(array([3.]), array([2]), array([6]))
注意在寻找左边界时,同高度的索引为 3 的峰值不被视为边界。相反,在寻找左边界时找到了两个最小值,索引为 0 和 2,在这种情况下,总是选择离评估峰值更近的那个。然而,在右侧,基线必须放在 6 处,因为更高的峰值代表了评估区域的右边界。
>>> peak_prominences(x, peaks, wlen=3.1)
(array([2.]), array([4]), array([6]))
在这里,我们将算法限制在从 3 到 7 的窗口范围内(长度为 5 个样本,因为wlen被四舍五入到下一个奇整数)。因此,在评估区域内只有两个候选样本,即两个相邻的样本和一个较小的突出度被计算。
scipy.signal.peak_widths
scipy.signal.peak_widths(x, peaks, rel_height=0.5, prominence_data=None, wlen=None)
计算信号中每个峰的宽度。
此函数计算峰的宽度,以样本为单位,相对于峰的高度和显著性的相对距离。
参数:
xsequence
一个带有峰的信号。
peakssequence
x 中峰值的索引。
rel_heightfloat,可选
选择峰宽度被测量的相对高度,作为其显著性的百分比。1.0 计算峰在其最低等高线处的宽度,而 0.5 在其显著性高度的一半处进行评估。必须至少为 0。详见注释以进一步解释。
prominence_datatuple,可选
一个元组,包含三个数组,与使用相同参数 x 和 peaks 调用 peak_prominences 时的输出相匹配。如果未提供此数据,则在内部计算。
wlenint,可选
作为内部计算 prominence_data 可选参数传递给 peak_prominences 的样本窗口长度。如果提供了 prominence_data,则忽略此参数。
返回:
widthsndarray
每个峰的宽度(以样本为单位)。
width_heightsndarray
widths 被评估的等高线高度。
left_ips, right_ipsndarray
在左右交点的插值位置,水平线分别在相应的评估高度。
引发:
ValueError
如果提供了 prominence_data 但不满足每个峰的条件 0 <= left_base <= peak <= right_base < x.shape[0],具有错误的 dtype,不是 C 连续的或形状不同。
警告:
PeakPropertyWarning
如果计算得到的任何宽度为 0,则引发此错误。这可能源于提供的 prominence_data 或如果 rel_height 设置为 0。
警告
此函数可能对包含 NaN 的数据返回意外结果。为了避免这种情况,应删除或替换 NaN。
另请参见
find_peaks
基于峰的特性在信号内找到峰。
peak_prominences
计算峰的显著性。
注释
计算峰宽度的基本算法如下:
-
使用公式计算评估高度 (h_{eval} = h_{Peak} - P \cdot R),其中 (h_{Peak}) 是峰本身的高度,(P) 是峰的显著性,(R) 是用参数 rel_height 指定的正比例。
-
在评估高度处绘制水平线,向两侧延伸,从峰值当前的垂直位置开始,直到这些线与斜坡、信号边界相交或越过峰值基底的垂直位置(详见
peak_prominences的定义)。对于第一种情况,与信号相交,使用线性插值估算真实的交点。 -
将宽度计算为两侧选择的端点之间的水平距离。因此,每个峰值的最大可能宽度是其基底之间的水平距离。
如上所示,要计算峰值的宽度,必须了解其突出和基底。您可以通过参数prominence_data自行提供这些数据。否则,它们将内部计算(详见peak_prominences)。
新版本 1.1.0 中引入。
示例
>>> import numpy as np
>>> from scipy.signal import chirp, find_peaks, peak_widths
>>> import matplotlib.pyplot as plt
创建一个包含两个重叠谐波的测试信号
>>> x = np.linspace(0, 6 * np.pi, 1000)
>>> x = np.sin(x) + 0.6 * np.sin(2.6 * x)
找到所有峰值,并计算它们在相对高度为 0.5(在高度的一半处的等高线)和 1(在完全突出高度处的最低等高线)时的宽度。
>>> peaks, _ = find_peaks(x)
>>> results_half = peak_widths(x, peaks, rel_height=0.5)
>>> results_half[0] # widths
array([ 64.25172825, 41.29465463, 35.46943289, 104.71586081,
35.46729324, 41.30429622, 181.93835853, 45.37078546])
>>> results_full = peak_widths(x, peaks, rel_height=1)
>>> results_full[0] # widths
array([181.9396084 , 72.99284945, 61.28657872, 373.84622694,
61.78404617, 72.48822812, 253.09161876, 79.36860878])
绘制信号、峰值和计算宽度的等高线
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.hlines(*results_half[1:], color="C2")
>>> plt.hlines(*results_full[1:], color="C3")
>>> plt.show()
scipy.signal.periodogram
scipy.signal.periodogram(x, fs=1.0, window='boxcar', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1)
使用周期图估计功率谱密度。
参数:
xarray_like
测量值的时间序列
fsfloat,可选项
x 时间序列的采样频率。默认为 1.0。
windowstr 或 元组 或 类似数组,可选项
需要使用的期望窗口。如果 window 是字符串或元组,则传递给 get_window 以生成窗口值,默认情况下是 DFT-偶数。参见 get_window 获取窗口列表及所需参数。如果 window 是类似数组,则直接用作窗口,并且其长度必须等于计算周期图的轴的长度。默认为 ‘boxcar’。
nfftint,可选项
使用的 FFT 长度。如果 None,则使用 x 的长度。
detrendstr 或 函数 或 False,可选项
指定如何去趋势化每个段。如果 detrend 是字符串,则传递给 detrend 函数的 type 参数。如果是函数,则接受一个段并返回去趋势化的段。如果 detrend 是 False,则不进行去趋势化。默认为 ‘constant’。
return_onesidedbool,可选项
如果 True,则为实数数据返回单边谱。如果 False,则返回双边谱。默认为 True,但对于复杂数据,始终返回双边谱。
scaling{ ‘density’, ‘spectrum’ },可选项
选择计算功率谱密度(‘density’),其中 Pxx 的单位为 V2/Hz,或计算功率谱(‘spectrum’),其中 Pxx 的单位为 V2,如果 x 单位为 V,fs 单位为 Hz。默认为 ‘density’。
axisint,可选项
计算周期图的轴;默认为最后一个轴(即 axis=-1)。
返回:
fndarray
样本频率数组。
Pxxndarray
x 的功率谱密度或功率谱。
参见
使用 Welch 方法估计功率谱密度
用于不均匀采样数据的 Lomb-Scargle 周期图
注意事项
版本 0.12.0 中的新功能。
示例
>>> 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.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
scipy.signal.welch
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.welch.html#scipy.signal.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')
使用韦尔奇方法估计功率谱密度。
韦尔奇方法[1]通过将数据分成重叠的段,计算每个段的修改周期图,并平均周期图来计算功率谱密度的估计。
参数:
xarray_like
测量值的时间序列
fs浮点数,可选项
x时间序列的采样频率。默认为 1.0。
window字符串或元组或 array_like,可选项
所用的期望窗口。如果window是字符串或元组,则传递给get_window以生成窗口值,默认情况下为 DFT-even。有关窗口和所需参数的列表,请参见get_window。如果window是 array_like,则直接用作窗口,其长度必须为 nperseg。默认为汉宁窗口。
nperseg整数,可选项
每个段的长度。默认为 None,但如果窗口是 str 或 tuple,则设置为 256,如果窗口是 array_like,则设置为窗口的长度。
noverlap整数,可选项
点数,用于段之间的重叠。如果为None,则noverlap = nperseg // 2。默认为None。
nfft整数,可选项
如果需要零填充的 FFT,则使用的 FFT 长度。如果为None,FFT 长度为nperseg。默认为None。
detrend字符串或函数或False,可选项
指定如何去趋势化每个段。如果detrend是一个字符串,则传递为detrend函数的type参数。如果它是一个函数,则取一个段并返回一个去趋势化的段。如果detrend是False,则不进行去趋势化。默认为'constant'。
return_onesided布尔值,可选项
如果为True,则针对实数数据返回单侧频谱。如果为False,则返回双侧频谱。默认为True,但对于复杂数据,始终返回双侧频谱。
scaling{ ‘密度’, ‘频谱’ },可选项
选择计算功率谱密度(‘密度’)还是计算功率谱(‘频谱’),其中Pxx的单位为 V**2/Hz,如果x以 V 测量,fs以 Hz 测量。默认为‘密度’
axis整数,可选项
计算周期图的轴;默认为最后一个轴(即axis=-1)。
average{ ‘mean’, ‘median’ },可选项
在平均周期图时使用的方法。默认为‘mean’。
新版本 1.2.0 中引入。
返回:
fndarray
采样频率阵列。
Pxxndarray
Power spectral density or power spectrum of x.
See also
Simple, optionally modified periodogram
Lomb-Scargle periodogram for unevenly sampled data
Notes
An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.
If noverlap is 0, this method is equivalent to Bartlett’s method [2].
New in version 0.12.0.
References
[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.
Examples
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 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)
Compute and plot the power spectral density.
>>> 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()
If we average the last half of the spectral density, to exclude the peak, we can recover the noise power on the signal.
>>> np.mean(Pxx_den[256:])
0.0009924865443739191
Now compute and plot the power spectrum.
>>> 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()
The peak height in the power spectrum is an estimate of the RMS amplitude.
>>> np.sqrt(Pxx_spec.max())
2.0077340678640727
If we now introduce a discontinuity in the signal, by increasing the amplitude of a small portion of the signal by 50, we can see the corruption of the mean average power spectral density, but using a median average better estimates the normal behaviour.
>>> 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()
scipy.signal.csd
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.csd.html#scipy.signal.csd
scipy.signal.csd(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, average='mean')
使用 Welch 方法估算交叉功率谱密度 Pxy。
参数:
x:类数组
测量值的时间序列
y:类数组
测量值的时间序列
fs:浮点数,可选
x 和 y 时间序列的采样频率。默认为 1.0。
window:字符串、元组或类数组,可选
所需使用的窗口。如果 window 是字符串或元组,则将其传递给 get_window 以生成窗口值,默认情况下为 DFT-even。有关窗口列表和所需参数,请参见 get_window。如果 window 是类数组,则直接使用作为窗口,其长度必须为 nperseg。默认为汉宁窗口。
nperseg:整数,可选
每个段的长度。默认为 None,但如果窗口为字符串或元组,则设置为 256;如果窗口为类数组,则设置为窗口的长度。
noverlap:整数,可选
分段之间重叠的点数。如果为 None,则 noverlap = nperseg // 2。默认为 None。
nfft:整数,可选
FFT 使用 FFT 使用的长度,如果需要进行零填充的 FFT。如果为 None,则 FFT 长度为 nperseg。默认为 None。
detrend:字符串、函数或 False,可选
指定如何去趋势每个段。如果 detrend 是字符串,则作为 detrend 函数的 type 参数传递。如果是函数,则接受一个段并返回去趋势后的段。如果 detrend 为 False,则不进行去趋势处理。默认为 ‘constant’。
return_onesided:布尔值,可选
如果为 True,则返回实数据的单边谱;如果为 False,则返回双边谱。默认为 True,但对于复杂数据,始终返回双边谱。
scaling:{‘density’, ‘spectrum’},可选
选择计算交叉功率谱密度(‘density’)还是交叉谱(‘spectrum’),其中 Pxy 的单位为 V2/Hz 或 V2,如果 x 和 y 分别以 V 和 Hz 计量,fs 以 Hz 计量。默认为 ‘density’。
axis:整数,可选
计算两个输入的 CSD 的轴;默认为最后一个轴(即 axis=-1)。
average:{‘mean’, ‘median’},可选
平均周期图的方法。如果频谱是复数,则分别计算实部和虚部的平均值。默认为 ‘mean’。
1.2.0 版新增功能。
返回:
f:ndarray
样本频率的数组。
Pxy:ndarray
x, y 的交叉谱密度或交叉功率谱。
另请参阅
简单的、可选修改后的周期图
不均匀采样数据的 Lomb-Scargle 周期图
使用威尔奇方法计算功率谱密度。[等同于 csd(x,x)]
威尔奇方法计算的幅度平方相干性。
注释
按照惯例,Pxy 是通过 X 的共轭 FFT 乘以 Y 的 FFT 来计算的。
如果输入序列长度不同,则较短的序列将被零填充以匹配。
适当的重叠量将取决于窗口的选择和您的需求。对于默认的 Hann 窗口,50%的重叠是在准确估计信号功率和不过度计数任何数据之间的合理折中。较窄的窗口可能需要更大的重叠。
0.16.0 版本的新增内容。
参考文献
[1]
P. Welch,“利用快速傅立叶变换估计功率谱的方法:基于短时平均和修改后的周期图”,IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967。
[2]
Rabiner, Lawrence R. 和 B. Gold。“数字信号处理的理论与应用” Prentice-Hall, pp. 414-419, 1975
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
生成两个具有一些共同特征的测试信号。
>>> fs = 10e3
>>> N = 1e5
>>> amp = 20
>>> freq = 1234.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> b, a = signal.butter(2, 0.25, 'low')
>>> x = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
>>> y = signal.lfilter(b, a, x)
>>> x += amp*np.sin(2*np.pi*freq*time)
>>> y += rng.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)
计算并绘制交叉谱密度的幅度。
>>> f, Pxy = signal.csd(x, y, fs, nperseg=1024)
>>> plt.semilogy(f, np.abs(Pxy))
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('CSD [V**2/Hz]')
>>> plt.show()