SciPy-1-12-中文文档-十四-
SciPy 1.12 中文文档（十四）
scipy.signal.zpk2ss
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.zpk2ss.html#scipy.signal.zpk2ss
scipy.signal.zpk2ss(z, p, k)
零极点增益表示转换为状态空间表示
参数：
z, p 序列
零点和极点。
k 浮点数
系统增益。
返回：
A, B, C, D 数组
系统的状态空间表示，处于控制规范形式。
scipy.signal.ss2tf
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.ss2tf.html#scipy.signal.ss2tf
scipy.signal.ss2tf(A, B, C, D, input=0)
从状态空间到传递函数。
A、B、C、D 定义了一个具有p个输入、q个输出和n个状态变量的线性状态空间系统。
参数：
Aarray_like
形状为(n, n)的状态（或系统）矩阵。
Barray_like
形状为(n, p)的输入矩阵。
Carray_like
形状为(q, n)的输出矩阵。
Darray_like
形状为(q, p)的馈送矩阵（或前馈矩阵）。
inputint，可选
对于多输入系统，使用的输入索引。
返回：
num2-D ndarray
结果传递函数的分子。num每行对应系统输出。每行是分子多项式的序列表示。
den1-D ndarray
结果传递函数的分母。den是分母多项式的序列表示。
示例
转换状态空间表示：
[ \begin{align}\begin{aligned}\begin{split}\dot{\textbf{x}}(t) = \begin{bmatrix} -2 & -1 \ 1 & 0 \end{bmatrix} \textbf{x}(t) + \begin{bmatrix} 1 \ 0 \end{bmatrix} \textbf{u}(t) \\end{split}\\textbf{y}(t) = \begin{bmatrix} 1 & 2 \end{bmatrix} \textbf{x}(t) + \begin{bmatrix} 1 \end{bmatrix} \textbf{u}(t)\end{aligned}\end{align} ]
>>> A = [[-2, -1], [1, 0]]
>>> B = [[1], [0]] # 2-D column vector
>>> C = [[1, 2]] # 2-D row vector
>>> D = 1
到传递函数：
[H(s) = \frac{s² + 3s + 3}{s² + 2s + 1}]
>>> from scipy.signal import ss2tf
>>> ss2tf(A, B, C, D)
(array([[1., 3., 3.]]), array([ 1., 2., 1.]))
scipy.signal.ss2zpk
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.ss2zpk.html#scipy.signal.ss2zpk
scipy.signal.ss2zpk(A, B, C, D, input=0)
将状态空间表示转换为零极点增益表示。
A, B, C, D 定义了具有p个输入，q个输出和n个状态变量的线性状态空间系统。
参数：
Aarray_like
形状为(n, n)的状态（或系统）矩阵。
Barray_like
形状为(n, p)的输入矩阵。
Carray_like
形状为(q, n)的输出矩阵。
Darray_like
形状为(q, p)的递送（或前馈）矩阵。
inputint, optional
对于多输入系统，使用的输入索引。
返回：
z, psequence
零点和极点。
kfloat
系统增益。
scipy.signal.sos2zpk
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.sos2zpk.html#scipy.signal.sos2zpk
scipy.signal.sos2zpk(sos)
返回一系列第二阶段的零点、极点和增益
参数：
sosarray_like
第二阶滤波器系数的数组，必须具有形状(n_sections, 6)。参见sosfilt以获取 SOS 滤波器格式规范。
返回：
zndarray
传递函数的零点。
pndarray
传递函数的极点。
kfloat
系统增益。
注意事项
即使某些零点和极点（实际上）为零，返回的零点和极点数量将为n_sections * 2。
自版本 0.16.0 新增。
scipy.signal.sos2tf
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.sos2tf.html#scipy.signal.sos2tf
scipy.signal.sos2tf(sos)
从一系列二阶段得到单一传递函数
参数：
sosarray_like
必须具有形状(n_sections, 6)的二阶滤波器系数数组。请参阅sosfilt以获取 SOS 滤波器格式规范。
返回：
bndarray
分子多项式系数。
andarray
分母多项式系数。
注意事项
新版本 0.16.0 中的新增内容。
示例
找到椭圆滤波器的多项式表示，使用其‘sos’（二阶段形式）格式。
>>> from scipy.signal import sos2tf
>>> from scipy import signal
>>> sos = signal.ellip(1, 0.001, 50, 0.1, output='sos')
>>> sos2tf(sos)
( array([0.91256522, 0.91256522, 0\. ]),
array([1\. , 0.82513043, 0\. ]))
scipy.signal.cont2discrete
scipy.signal.cont2discrete(system, dt, method='zoh', alpha=None)
将连续状态空间系统转换为离散系统。
参数：
system：描述系统的元组或 lti 的实例
以下内容给出元组中的元素数量和解释：
- 1: （
lti的实例）- 2: （分子，分母）
- 3: （零点，极点，增益）
- 4: （A, B, C, D）
dt：浮点数
离散化时间步长。
method：字符串，可选
使用哪种方法：
- gbt：广义双线性变换
- bilinear：Tustin 逼近法（“gbt” with alpha=0.5）
- euler: 欧拉（或向前差分）方法（“gbt” with alpha=0）
- backward_diff：后向差分（“gbt” with alpha=1.0）
- zoh：零阶保持（默认）
- foh：一阶保持（versionadded: 1.3.0）
- impulse：等效冲激响应（versionadded: 1.3.0）
alpha：在 [0, 1] 范围内的浮点数，可选
应仅在 method=”gbt” 时指定的广义双线性变换加权参数，否则将被忽略
返回：
sysd：包含离散系统的元组
根据输入类型，输出将采用以下形式
-
(num, den, dt)：用于传递函数输入
-
(zeros, poles, gain, dt)：用于零点-极点-增益输入
-
(A, B, C, D, dt)：用于状态空间系统输入
注意事项
默认情况下，该例程使用零阶保持（zoh）方法执行转换。也可以使用广义双线性变换，其中包括常见的 Tustin 双线性逼近法、欧拉方法技术或后向差分技术。
零阶保持（zoh）方法基于 [1]，广义双线性逼近基于 [2] 和 [3]，一阶保持（foh）方法基于 [4]。
参考文献
[1]
en.wikipedia.org/wiki/Discretization#Discretization_of_linear_state_space_models
[2]
techteach.no/publications/discretetime_signals_systems/discrete.pdf
[3]
G. Zhang, X. Chen 和 T. Chen，《Digital redesign via the generalized bilinear transformation》，Int. J. Control，第 82 卷，第 4 期，2009 年，页码 741-754。(www.mypolyuweb.hk/~magzhang/Research/ZCC09_IJC.pdf)
[4]
G. F. Franklin, J. D. Powell 和 M. L. Workman，《Digital control of dynamic systems》，第 3 版，Menlo Park, Calif: Addison-Wesley，1998 年，页码 204-206。
示例
我们可以将连续状态空间系统转换为离散系统：
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import cont2discrete, lti, dlti, dstep
定义一个连续状态空间系统。
>>> A = np.array([[0, 1],[-10., -3]])
>>> B = np.array([[0],[10.]])
>>> C = np.array([[1., 0]])
>>> D = np.array([[0.]])
>>> l_system = lti(A, B, C, D)
>>> t, x = l_system.step(T=np.linspace(0, 5, 100))
>>> fig, ax = plt.subplots()
>>> ax.plot(t, x, label='Continuous', linewidth=3)
使用多种方法将其转换为离散状态空间系统。
>>> dt = 0.1
>>> for method in ['zoh', 'bilinear', 'euler', 'backward_diff', 'foh', 'impulse']:
... d_system = cont2discrete((A, B, C, D), dt, method=method)
... s, x_d = dstep(d_system)
... ax.step(s, np.squeeze(x_d), label=method, where='post')
>>> ax.axis([t[0], t[-1], x[0], 1.4])
>>> ax.legend(loc='best')
>>> fig.tight_layout()
>>> plt.show()
scipy.signal.place_poles
scipy.signal.place_poles(A, B, poles, method='YT', rtol=0.001, maxiter=30)
计算 K，使得特征值(A - dot(B, K))=poles。
K 是增益矩阵，使得由线性系统AX+BU描述的过程的闭环极点，即特征值A - B*K，尽可能接近所要求的极点。
支持 SISO、MISO 和 MIMO 系统。
参数：
A, Bndarray
线性系统AX + BU的状态空间表示。
polesarray_like
所需的实部极点和/或共轭复极点。仅支持method="YT"（默认）的复极点。
method: {‘YT’, ‘KNV0’}, optional
选择用于找到增益矩阵 K 的方法之一：
- ‘YT’：Yang Tits
- ‘KNV0’：Kautsky、Nichols、Van Dooren 更新方法 0
有关算法的详细信息，请参见参考文献和注释。
rtol: float, optional
每次迭代后，比较A - B*K的特征向量的行列式与其先前值，当这两个值之间的相对误差低于rtol时，算法停止。默认值为 1e-3。
maxiter: int, optional
计算增益矩阵的最大迭代次数。默认值为 30。
返回：
full_state_feedbackBunch 对象
full_state_feedback 由以下组成：
gain_matrix1-D ndarray
闭环矩阵 K，使得A-BK的特征值尽可能接近要求的极点。
computed_poles1-D ndarray
与A-BK对应的极点，首先按升序排列实部极点，然后按字典顺序排列共轭复极点。
requested_poles1-D ndarray
算法要求的极点如上所述排序，可能与实际实现的不同。
X2-D ndarray
传递矩阵如X * diag(poles) = (A - B*K)*X（参见注释）
rtolfloat
实现在det(X)上实现的相对容差（参见注释）。如果能够解决系统diag(poles) = (A - B*K)，则rtol将为 NaN；当优化算法无法做任何事情，即B.shape[1] == 1时，为 0。
nb_iterint
在收敛之前执行的迭代次数。如果能够解决系统diag(poles) = (A - B*K)，则nb_iter将为 NaN；当优化算法无法做任何事情，即B.shape[1] == 1时，为 0。
注释
Tits 和 Yang（YT）论文[2]是 Kautsky 等人（KNV）原始论文[1]的更新。KNV 依赖于秩-1 更新，以找到传递矩阵 X，使得X * diag(poles) = (A - B*K)*X，而 YT 使用秩-2 更新。这通常会产生更为健壮的解决方案（参见[2]第 21-22 页），此外 YT 算法支持复极点，而 KNV 原始版本不支持。此处仅实现了 KNV 提出的更新方法 0，因此命名为'KNV0'。
KNV 扩展到复极点被用于 Matlab 的place函数，YT 以非自由许可证由 Slicot 发布，名称为robpole。目前尚不清楚和未记录如何将 KNV0 扩展到复极点（Tits 和 Yang 在其论文第 14 页声称他们的方法不能用于扩展 KNV 到复极点），因此只有 YT 在这个实现中支持它们。
由于对于 MIMO 系统来说，极点配置问题的解并不唯一，因此两种方法都从一个试探性的传递矩阵开始，该矩阵以不同的方式改变以增加其行列式。已经证明这两种方法都会收敛到一个稳定的解，然而，根据初始传递矩阵的选择方式，它们将收敛到不同的解，因此使用'KNV0'不能保证产生与 Matlab 或任何其他实现这些算法的结果相似。
在大多数情况下，使用默认方法'YT'应该是可以的；'KNV0'仅提供是因为在某些特定情况下'YT'需要。此外，'YT'比使用绝对值abs(det(X))作为鲁棒性指标时，'KNV0'平均提供更稳健的结果。
[2] 可作为技术报告在以下网址获取：hdl.handle.net/1903/5598
参考文献
[1] (1,2)
J. Kautsky, N.K. Nichols 和 P. van Dooren，“线性状态反馈中的鲁棒极点分配”，《国际控制杂志》，第 41 卷，第 1129-1155 页，1985 年。
[2] (1,2,3)
A.L. Tits 和 Y. Yang，“用于鲁棒状态反馈极点分配的全局收敛算法”，《IEEE 自动控制杂志》，第 41 卷，第 1432-1452 页，1996 年。
示例
一个简单的示例，展示了使用 KNV 和 YT 算法进行实际极点配置的方法。这是参考 KNV 出版物第四部分中的示例 1（1）：
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> A = np.array([[ 1.380, -0.2077, 6.715, -5.676 ],
... [-0.5814, -4.290, 0, 0.6750 ],
... [ 1.067, 4.273, -6.654, 5.893 ],
... [ 0.0480, 4.273, 1.343, -2.104 ]])
>>> B = np.array([[ 0, 5.679 ],
... [ 1.136, 1.136 ],
... [ 0, 0, ],
... [-3.146, 0 ]])
>>> P = np.array([-0.2, -0.5, -5.0566, -8.6659])
现在使用 KNV 方法 0 计算 K，使用默认的 YT 方法以及使用强制算法 100 次迭代的 YT 方法，并在每次调用后打印一些结果。
>>> fsf1 = signal.place_poles(A, B, P, method='KNV0')
>>> fsf1.gain_matrix
array([[ 0.20071427, -0.96665799, 0.24066128, -0.10279785],
[ 0.50587268, 0.57779091, 0.51795763, -0.41991442]])
>>> fsf2 = signal.place_poles(A, B, P) # uses YT method
>>> fsf2.computed_poles
array([-8.6659, -5.0566, -0.5 , -0.2 ])
>>> fsf3 = signal.place_poles(A, B, P, rtol=-1, maxiter=100)
>>> fsf3.X
array([[ 0.52072442+0.j, -0.08409372+0.j, -0.56847937+0.j, 0.74823657+0.j],
[-0.04977751+0.j, -0.80872954+0.j, 0.13566234+0.j, -0.29322906+0.j],
[-0.82266932+0.j, -0.19168026+0.j, -0.56348322+0.j, -0.43815060+0.j],
[ 0.22267347+0.j, 0.54967577+0.j, -0.58387806+0.j, -0.40271926+0.j]])
X 的行列式的绝对值是检查结果鲁棒性的良好指标，'KNV0'和'YT'都旨在最大化它。以下是上述结果鲁棒性的比较：
>>> abs(np.linalg.det(fsf1.X)) < abs(np.linalg.det(fsf2.X))
True
>>> abs(np.linalg.det(fsf2.X)) < abs(np.linalg.det(fsf3.X))
True
现在是一个复极点的简单示例：
>>> A = np.array([[ 0, 7/3., 0, 0 ],
... [ 0, 0, 0, 7/9. ],
... [ 0, 0, 0, 0 ],
... [ 0, 0, 0, 0 ]])
>>> B = np.array([[ 0, 0 ],
... [ 0, 0 ],
... [ 1, 0 ],
... [ 0, 1 ]])
>>> P = np.array([-3, -1, -2-1j, -2+1j]) / 3.
>>> fsf = signal.place_poles(A, B, P, method='YT')
我们可以在复平面上绘制期望和计算的极点：
>>> t = np.linspace(0, 2*np.pi, 401)
>>> plt.plot(np.cos(t), np.sin(t), 'k--') # unit circle
>>> plt.plot(fsf.requested_poles.real, fsf.requested_poles.imag,
... 'wo', label='Desired')
>>> plt.plot(fsf.computed_poles.real, fsf.computed_poles.imag, 'bx',
... label='Placed')
>>> plt.grid()
>>> plt.axis('image')
>>> plt.axis([-1.1, 1.1, -1.1, 1.1])
>>> plt.legend(bbox_to_anchor=(1.05, 1), loc=2, numpoints=1)
scipy.signal.chirp
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.chirp.html#scipy.signal.chirp
scipy.signal.chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True)
频率扫描余弦发生器。
在下文中，“Hz”应理解为“每单位的循环次数”；这里没有要求单位必须为一秒。重要的区别在于旋转的单位是循环，而不是弧度。同样，t 可能是空间的度量，而不是时间的度量。
参数：
t：类似数组
评估波形的时间点。
f0：浮点数
时间 t=0 时的频率（例如 Hz）。
t1：浮点数
指定 f1 的时间。
f1：浮点数
时间 t1 处波形的频率（例如 Hz）。
method：{‘linear’, ‘quadratic’, ‘logarithmic’, ‘hyperbolic’}，可选
频率扫描类型。如果未给出，则假定为 linear。有关更多详细信息，请参见下面的注意事项。
phi：浮点数，可选
相位偏移，以度为单位。默认值为 0。
vertex_zero：布尔值，可选
该参数仅在 method 为 ‘quadratic’ 时使用。它决定频率图的抛物线顶点是否在 t=0 或 t=t1 处。
返回：
y：ndarray
包含在 t 上评估请求的时变频率信号的 numpy 数组。更精确地说，函数返回 cos(phase + (pi/180)*phi)，其中 phase 是 2*pi*f(t) 的积分（从 0 到 t）。f(t) 如下定义。
另请参见
sweep_poly
注意
method 有四个选项。以下公式给出了由 chirp() 生成的信号的瞬时频率（以 Hz 为单位）。为方便起见，下面显示的较短名称也可以使用。
linear、lin、li：
f(t) = f0 + (f1 - f0) * t / t1
quadratic、quad、q：
频率 f(t) 的图形是通过点 (0, f0) 和 (t1, f1) 的抛物线。默认情况下，抛物线顶点位于 (0, f0) 处。如果 vertex_zero 为 False，则顶点位于 (t1, f1) 处。公式如下：
如果
vertex_zero为 True：
f(t) = f0 + (f1 - f0) * t**2 / t1**2else：
f(t) = f1 - (f1 - f0) * (t1 - t)**2 / t1**2要使用更一般的二次函数或任意多项式，请使用函数
scipy.signal.sweep_poly。
logarithmic、log、lo：
f(t) = f0 * (f1/f0)**(t/t1)f0 和 f1 必须非零，并且符号相同。
该信号也称为几何或指数啁啾。
hyperbolic、hyp：
f(t) = f0*f1*t1 / ((f0 - f1)*t + f1*t1)f0 和 f1 必须非零。
示例
在示例中将使用以下内容：
>>> import numpy as np
>>> from scipy.signal import chirp, spectrogram
>>> import matplotlib.pyplot as plt
首个示例中，我们将绘制从 6 Hz 到 1 Hz 的线性啁啾波形，时长为 10 秒：
>>> t = np.linspace(0, 10, 1500)
>>> w = chirp(t, f0=6, f1=1, t1=10, method='linear')
>>> plt.plot(t, w)
>>> plt.title("Linear Chirp, f(0)=6, f(10)=1")
>>> plt.xlabel('t (sec)')
>>> plt.show()
对于其余示例，我们将使用更高的频率范围，并使用scipy.signal.spectrogram来展示结果。我们将使用 7200 Hz 采样的 4 秒间隔。
>>> fs = 7200
>>> T = 4
>>> t = np.arange(0, int(T*fs)) / fs
我们将使用此函数在每个示例中绘制频谱图。
>>> def plot_spectrogram(title, w, fs):
... ff, tt, Sxx = spectrogram(w, fs=fs, nperseg=256, nfft=576)
... fig, ax = plt.subplots()
... ax.pcolormesh(tt, ff[:145], Sxx[:145], cmap='gray_r',
... shading='gouraud')
... ax.set_title(title)
... ax.set_xlabel('t (sec)')
... ax.set_ylabel('Frequency (Hz)')
... ax.grid(True)
...
从 1500 Hz 到 250 Hz 的二次啁啾（频率抛物线曲线的顶点在 t=0）：
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic')
>>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250', w, fs)
>>> plt.show()
从 1500 Hz 到 250 Hz 的二次啁啾（频率抛物线曲线的顶点在 t=T）：
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic',
... vertex_zero=False)
>>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250\n' +
... '(vertex_zero=False)', w, fs)
>>> plt.show()
从 1500 Hz 到 250 Hz 的对数啁啾：
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='logarithmic')
>>> plot_spectrogram(f'Logarithmic Chirp, f(0)=1500, f({T})=250', w, fs)
>>> plt.show()
从 1500 Hz 到 250 Hz 的双曲线啁啾：
>>> w = chirp(t, f0=1500, f1=250, t1=T, method='hyperbolic')
>>> plot_spectrogram(f'Hyperbolic Chirp, f(0)=1500, f({T})=250', w, fs)
>>> plt.show()
scipy.signal.gausspulse
scipy.signal.gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False, retenv=False)
返回一个高斯调制的正弦波：
exp(-a t²) exp(1j*2*pi*fc*t).
如果 retquad 为 True，则返回实部和虚部（同相和象限）。如果 retenv 为 True，则返回包络（未调制信号）。否则，返回调制正弦波的实部。
参数：
tndarray 或字符串 'cutoff'
输入数组。
fcfloat，可选
中心频率（例如 Hz）。默认为 1000。
bwfloat，可选
脉冲在频率域的分数带宽（例如 Hz）。默认为 0.5。
bwrfloat，可选
用于计算分数带宽的参考级别（dB）。默认为 -6。
tprfloat，可选
如果 t 为 'cutoff'，则函数返回脉冲幅度下降至 tpr（以 dB 为单位）的截止时间。默认为 -60。
retquadbool，可选
如果为 True，则返回信号的象限（虚部）以及实部。默认为 False。
retenvbool，可选
如果为 True，则返回信号的包络。默认为 False。
返回：
yIndarray
信号的实部。始终返回。
yQndarray
信号的虚部。仅在 retquad 为 True 时返回。
yenvndarray
信号的包络。仅在 retenv 为 True 时返回。
另请参阅
scipy.signal.morlet
示例
绘制实部、虚部和 5 Hz 脉冲的包络，以 100 Hz 采样 2 秒：
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(-1, 1, 2 * 100, endpoint=False)
>>> i, q, e = signal.gausspulse(t, fc=5, retquad=True, retenv=True)
>>> plt.plot(t, i, t, q, t, e, '--')
scipy.signal.max_len_seq
scipy.signal.max_len_seq(nbits, state=None, length=None, taps=None)
最大长度序列（MLS）生成器。
参数：
nbitsint
要使用的位数。生成的序列长度将为 (2**nbits) - 1。请注意，生成长序列（例如大于 nbits == 16）可能需要很长时间。
statearray_like，可选
如果是数组，则必须是 nbits 长度，并将被转换为二进制（bool）表示。如果为 None，则使用全 1 的种子，生成可重复的表示。如果 state 全为零，则会引发错误，因为这是无效的。默认值：None。
lengthint，可选
要计算的样本数。如果为 None，则计算整个长度 (2**nbits) - 1。
tapsarray_like，可选
用于生成多项式 taps（例如 [7, 6, 1] 用于 8 位序列）。如果为 None，则会自动选择 taps（最多支持 nbits == 32）。
返回：
seqarray
结果的 MLS 序列，由 0 和 1 组成。
statearray
移位寄存器的最终状态。
笔记
MLS 生成算法在以下面描述：
taps 的默认值专门取自每个 nbits 值的第一个选项：
0.15.0 版中的新功能。
示例
MLS 使用二进制约定：
>>> from scipy.signal import max_len_seq
>>> max_len_seq(4)[0]
array([1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0], dtype=int8)
MLS 具有白色频谱（除了直流分量）：
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from numpy.fft import fft, ifft, fftshift, fftfreq
>>> seq = max_len_seq(6)[0]*2-1 # +1 and -1
>>> spec = fft(seq)
>>> N = len(seq)
>>> plt.plot(fftshift(fftfreq(N)), fftshift(np.abs(spec)), '.-')
>>> plt.margins(0.1, 0.1)
>>> plt.grid(True)
>>> plt.show()
MLS 的循环自相关是一个冲激：
>>> acorrcirc = ifft(spec * np.conj(spec)).real
>>> plt.figure()
>>> plt.plot(np.arange(-N/2+1, N/2+1), fftshift(acorrcirc), '.-')
>>> plt.margins(0.1, 0.1)
>>> plt.grid(True)
>>> plt.show()
MLS 的线性自相关大致上是一个冲激：
>>> acorr = np.correlate(seq, seq, 'full')
>>> plt.figure()
>>> plt.plot(np.arange(-N+1, N), acorr, '.-')
>>> plt.margins(0.1, 0.1)
>>> plt.grid(True)
>>> plt.show()
scipy.signal.sawtooth
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.sawtooth.html#scipy.signal.sawtooth
scipy.signal.sawtooth(t, width=1)
返回一个周期性的锯齿波或三角波形。
锯齿波形的周期是 2*pi，在区间 0 到 width*2*pi 上升从-1 到 1，然后在区间 width*2*pi 到 2*pi 下降从 1 到-1。width 必须在区间 [0, 1] 内。
请注意这不是带限制的。它产生无限多的谐波，这些谐波在频率谱上来回反射。
参数：
tarray_like
时间。
widtharray_like, 可选
上升斜坡的宽度，作为总周期的比例。默认为 1，生成上升斜坡，而 0 生成下降斜坡。width = 0.5 生成三角波。如果是一个数组，则导致波形随时间变化，并且必须与 t 具有相同的长度。
返回：
yndarray
包含锯齿波形的输出数组。
示例
以 500 Hz 对 1 秒钟进行采样的 5 Hz 波形：
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> t = np.linspace(0, 1, 500)
>>> plt.plot(t, signal.sawtooth(2 * np.pi * 5 * t))
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’），其中 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’。
新版本 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()
scipy.signal.coherence
scipy.signal.coherence(x, y, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, detrend='constant', axis=-1)
使用 Welch 方法估计离散时间信号 X 和 Y 的幅度平方相干性估计，Cxy。
Cxy = abs(Pxy)**2/(Pxx*Pyy)，其中 Pxx 和 Pyy 是 X 和 Y 的功率谱密度估计，Pxy 是 X 和 Y 的交叉谱密度估计。
参数：
xarray_like
测量值的时间序列
yarray_like
测量值的时间序列
fsfloat，可选
x 和 y 时间序列的采样频率。默认为 1.0。
windowstr 或者 tuple 或者 array_like，可选
所需使用的窗口。如果 window 是字符串或元组，则传递给 get_window 以生成窗口值，默认情况下为 DFT-even。参见 get_window 获取窗口列表和必需的参数。如果 window 是 array_like，则直接用作窗口，其长度必须为 nperseg。默认为汉宁窗口。
npersegint，可选
每个段的长度。默认为 None，但如果窗口是字符串或元组，则设置为 256，如果窗口是 array_like，则设置为窗口的长度。
noverlap: int, 可选
在分段之间重叠的点数。如果 None，则 noverlap = nperseg // 2。默认为 None。
nfftint，可选
如果需要零填充 FFT，则使用的 FFT 长度。如果 None，则 FFT 长度为 nperseg。默认为 None。
detrendstr 或者 函数 或者 False，可选
指定如何去趋势化每个段。如果 detrend 是一个字符串，则作为 type 参数传递给 detrend 函数。如果它是一个函数，则它接受一个段并返回去趋势化的段。如果 detrend 是 False，则不执行去趋势化。默认为 'constant'。
axisint，可选
计算两个输入信号的相干性的轴；默认为最后一个轴（即 axis=-1）。
返回：
fndarray
样本频率的数组。
Cxyndarray
x 和 y 的幅度平方相干性。
另请参阅
简单的，可选修改的周期图
不均匀采样数据的 Lomb-Scargle 周期图
Welch 方法计算的功率谱密度。
Welch 方法的交叉谱密度。
注意事项
适当的重叠量将取决于窗口的选择和您的要求。对于默认的 Hann 窗口，50%的重叠是在准确估计信号功率和不过多计算任何数据之间的合理折衷。更窄的窗口可能需要更大的重叠。
从版本 0.16.0 开始新增。
参考文献
[1]
P. Welch，“用于估计功率谱的快速傅立叶变换的使用：一种基于短期修改周期图平均的方法”，IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967 年
[2]
Stoica, Petre 和 Randolph Moses，“信号的频谱分析”，Prentice Hall，2005 年
示例
>>> 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, Cxy = signal.coherence(x, y, fs, nperseg=1024)
>>> plt.semilogy(f, Cxy)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('Coherence')
>>> plt.show()
scipy.signal.spectrogram
scipy.signal.spectrogram(x, fs=1.0, window=('tukey', 0.25), nperseg=None, noverlap=None, nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1, mode='psd')
使用连续的傅里叶变换计算频谱图。
频谱图可用作可视化非平稳信号频率内容随时间变化的一种方法。
遗留
此函数被视为遗留版本，将不再接收更新。这可能意味着在未来的 SciPy 版本中将被移除。ShortTimeFFT是一个更新的 STFT / ISTFT 实现，具有更多功能，还包括一个spectrogram方法。在SciPy 用户指南的Short-Time Fourier Transform部分中可以找到这些实现之间的比较。
参数：
xarray_like
测量值的时间序列
fsfloat，可选
x时间序列的采样频率。默认为 1.0。
windowstr 或元组或 array_like，可选
期望使用的窗口。如果window是字符串或元组，则会传递给get_window以生成窗口数值，默认情况下为 DFT 偶数。请参阅get_window获取窗口列表和所需参数。如果window是 array_like，则将直接使用作为窗口，并且其长度必须为nperseg。默认为 Tukey 窗口，形状参数为 0.25。
npersegint，可选
每个段的长度。默认为 None，但如果window是字符串或元组，则设置为 256，如果window是 array_like，则设置为窗口的长度。
noverlapint，可选
每个段之间重叠的点数。如果为None，则noverlap = nperseg // 8。默认为None。
nfftint，可选
所使用的 FFT 长度，如果需要零填充 FFT。如果为None，则 FFT 长度为nperseg。默认为None。
detrendstr 或函数或False，可选
指定如何去趋势化每个段。如果detrend是一个字符串，则传递为detrend函数的type参数。如果是一个函数，则接受一个段并返回去趋势化的段。如果detrend为False，则不进行去趋势化。默认为‘constant’。
return_onesidedbool，可选
如果True，返回实数据的单侧频谱。如果False，返回双侧频谱。默认为True，但对于复杂数据，始终返回双侧频谱。
scaling，可选
选择计算功率谱密度（‘density’）或功率谱（‘spectrum’），其中Sxx的单位为 V**2/Hz，如果x以 V 为单位，fs以 Hz 为单位。默认为‘density’。
axisint，可选
计算谱图的轴；默认为最后一个轴（即axis=-1）。
modestr，可选
定义预期的返回值类型。选项有[‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]。‘complex’等同于没有填充或边界扩展的stft的输出。‘magnitude’返回 STFT 的绝对幅度。‘angle’和‘phase’分别返回 STFT 的复角，带有和不带有展开。
返回：
fndarray
样本频率的数组。
tndarray
分段时间的数组。
Sxxndarray
x 的谱图。默认情况下，Sxx 的最后一个轴对应于段时间。
另请参阅
periodogram
简单的、可选修改后的周期图
lombscargle
Lomb-Scargle 不规则采样数据的周期图
welch
Welch 方法的功率谱密度。
csd
Welch 方法的交叉谱密度
ShortTimeFFT
提供更多功能的新 STFT/ISTFT 实现，其中还包括一个spectrogram方法。
注释
适当的重叠量取决于窗口的选择和您的需求。与 Welch 方法相反，在计算谱图时，人们可能希望使用较小的重叠（或者根本不重叠），以保持各个段的统计独立性。因此，默认窗口是 Tukey 窗口，每端重叠窗口长度的 1/8。
新版本 0.16.0 中引入。
参考文献
[1]
Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”，Prentice Hall，1999。
示例
>>> import numpy as np
>>> from scipy import signal
>>> from scipy.fft import fftshift
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
生成一个测试信号，幅值为 2 Vrms 的正弦波，其频率围绕 3kHz 缓慢调制，被指数衰减的白噪声污染，采样频率为 10 kHz。
>>> fs = 10e3
>>> N = 1e5
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.01 * fs / 2
>>> time = np.arange(N) / float(fs)
>>> mod = 500*np.cos(2*np.pi*0.25*time)
>>> carrier = amp * np.sin(2*np.pi*3e3*time + mod)
>>> noise = rng.normal(scale=np.sqrt(noise_power), size=time.shape)
>>> noise *= np.exp(-time/5)
>>> x = carrier + noise
计算并绘制谱图。
>>> f, t, Sxx = signal.spectrogram(x, fs)
>>> plt.pcolormesh(t, f, Sxx, shading='gouraud')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()
注意，如果使用的输出不是单边的话，请使用以下内容：
>>> f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False)
>>> plt.pcolormesh(t, fftshift(f), fftshift(Sxx, axes=0), shading='gouraud')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()
scipy.signal.lombscargle
scipy.signal.lombscargle(x, y, freqs)
计算 Lomb-Scargle 周期图。
Lomb-Scargle 周期图由 Lomb [1]开发，并由 Scargle [2]进一步扩展，用于发现和测试不均匀时间采样中弱周期信号的显著性。
当 normalize 设置为 False（默认值）时，计算得到的周期图未归一化，对于具有振幅 A 的谐波信号，对于足够大的 N，它取值为(A**2) * N/4。
当 normalize 设置为 True 时，计算得到的周期图将通过数据围绕常数参考模型（在零点）的残差进行归一化。
输入数组应为 1-D，并将转换为 float64 类型。
参数：
xarray_like
样本时间。
yarray_like
测量值。
freqsarray_like
输出周期图的角频率。
precenterbool, optional
通过减去均值预置测量值。
normalizebool, optional
计算归一化周期图。
返回：
pgramarray_like
Lomb-Scargle 周期图。
Raises：
ValueError
如果输入数组 x 和 y 的形状不同。
另见
istft
逆短时傅立叶变换
check_COLA
检查是否满足常数重叠加（COLA）约束
welch
Welch 方法的功率谱密度
spectrogram
Welch 方法的谱图
csd
Welch 方法的交叉谱密度
Notes
此子程序使用了由 Townsend 稍作修改的算法来计算周期图[3]，该算法允许在每个频率上仅通过输入数组的一次传递计算周期图。
算法运行时间大致按 O(x * freqs)或 O(N²)缩放，适用于大量样本和频率。
参考文献
[1]
N.R. Lomb，“不等间隔数据的最小二乘频率分析”，《天体物理学和空间科学》，第 39 卷，第 447-462 页，1976 年
[2]
J.D. Scargle，“天文时间序列分析研究 II - 不均匀间隔数据谱分析的统计方面”，《天体物理学期刊》，第 263 卷，第 835-853 页，1982 年
[3]
R.H.D. Townsend，“使用图形处理单元快速计算 Lomb-Scargle 周期图”，《天体物理学期刊增刊》，第 191 卷，第 247-253 页，2010 年
示例
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
首先为信号定义一些输入参数：
>>> A = 2.
>>> w0 = 1. # rad/sec
>>> nin = 150
>>> nout = 100000
随机生成样本时间：
>>> x = rng.uniform(0, 10*np.pi, nin)
绘制所选时间的正弦波：
>>> y = A * np.cos(w0*x)
定义用于计算周期图的频率数组：
>>> w = np.linspace(0.01, 10, nout)
计算 Lomb-Scargle 周期图：
>>> import scipy.signal as signal
>>> pgram = signal.lombscargle(x, y, w, normalize=True)
现在制作输入数据的图表：
>>> fig, (ax_t, ax_w) = plt.subplots(2, 1, constrained_layout=True)
>>> ax_t.plot(x, y, 'b+')
>>> ax_t.set_xlabel('Time [s]')
然后绘制归一化周期图：
>>> ax_w.plot(w, pgram)
>>> ax_w.set_xlabel('Angular frequency [rad/s]')
>>> ax_w.set_ylabel('Normalized amplitude')
>>> plt.show()
scipy.signal.vectorstrength
scipy.signal.vectorstrength(events, period)
确定与给定周期对应的事件的矢量强度。
矢量强度是相位同步的一个度量，表明事件的定时如何与周期信号的单个周期同步。
如果使用多个周期，计算每个的矢量强度。这称为“共振矢量强度”。
参数：
events1D 数组类似
包含事件时间点的时间点数组。
periodfloat 或 array_like
事件应该与之同步的信号周期。周期与 events 单位相同。它也可以是周期数组，此时输出也是相同长度的数组。
返回：
strengthfloat 或 1D 数组
同步的强度。1.0 是完美同步，0.0 是没有同步。如果 period 是一个数组，则这也是一个数组，其中每个元素包含相应周期的矢量强度。
phasefloat 或 array
事件与弧度最强同步的相位。如果 period 是一个数组，则这也是一个数组，其中每个元素包含相应周期的相位。
参考文献：
van Hemmen, JL, Longtin, A 和 Vollmayr, AN。测试共振矢量
strength：听觉系统、电鱼和噪声。混沌 21, 047508 (2011); DOI:10.1063/1.3670512.
van Hemmen, JL。Goldberg、Brown 和 von Mises 后的矢量强度：
生物和数学视角。生物控制。2013 年 8 月;107(4):385-96. DOI:10.1007/s00422-013-0561-7.
van Hemmen, JL 和 Vollmayr, AN。共振矢量强度：发生了什么
当我们改变“探测”频率但保持尖峰时间不变时。生物控制。2013 年 8 月;107(4):491-94。DOI:10.1007/s00422-013-0560-8.
scipy.signal.ShortTimeFFT
class scipy.signal.ShortTimeFFT(win, hop, fs, *, fft_mode='onesided', mfft=None, dual_win=None, scale_to=None, phase_shift=0)
提供参数化的离散短时傅里叶变换 (stft) 及其逆变换 (istft)。
stft 通过滑动窗口 (win) 在输入信号上以 hop 增量计算连续的 FFT，可用于量化频谱随时间的变化。
stft 由复数矩阵 S[q,p] 表示，其中第 p 列代表以时间 t[p] = p * delta_t = p * hop * T 居中的窗口的 FFT，其中 T 是输入信号的采样间隔。第 q 行表示在频率 f[q] = q * delta_f 处的值，其中 delta_f = 1 / (mfft * T) 是 FFT 的频率分辨率。
逆 STFT istft 通过逆转 STFT 步骤计算：取 S[q,p] 的第 p 切片的 IFFT，并与所谓的双窗口（参见 dual_win）结果相乘。将结果按 p * delta_t 移动，并将结果添加到先前移动的结果以重建信号。如果仅知道双窗口并且 STFT 可逆，则可以使用 from_dual 实例化此类。
由于时间 t = 0 约定为输入信号的第一个样本，STFT 值通常具有负时间槽。因此，像p_min或k_min这样的负索引不像标准 Python 索引中的倒数计数从数组末尾开始，而是位于 t = 0 的左侧。
更详细的信息可以在 SciPy 用户指南的短时傅里叶变换部分找到。
请注意，除了使用scaling的scale_to之外，初始化器的所有参数都具有相同的命名属性。
参数：
winnp.ndarray
窗口必须是一个实数或复数值的一维数组。
hopint
在每个步骤中窗口移动的样本增量。
fsfloat
输入信号和窗口的采样频率。其与采样间隔T的关系为T = 1 / fs。
fft_mode‘twosided’, ‘centered’, ‘onesided’, ‘onesided2X’
要使用的 FFT 模式（默认为'onesided'）。有关详细信息，请参见属性fft_mode。
mfft: int | None
如果需要零填充 FFT，则使用的 FFT 的长度。如果为None（默认），则使用窗口win的长度。
dual_winnp.ndarray | None
win的双重窗口。如果设置为None，则在需要时进行计算。
scale_to‘magnitude’, ‘psd’ | None
如果不为None（默认），则缩放窗口函数，使每个 STFT 列表示“幅度”或功率谱密度（'psd'）谱。此参数将属性scaling设置为相同值。有关详细信息，请参见方法scale_to。
phase_shiftint | None
如果设置，对每个频率 f 添加一个线性相位 phase_shift / mfft * f。默认值 0 确保在零切片上没有相位移（其中 t=0 居中）。有关详细信息，请参阅属性 phase_shift。
示例
以下示例显示了带有变频 (f_i(t)) 的正弦波的 STFT 幅度（在图中由红色虚线标记）：
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import ShortTimeFFT
>>> from scipy.signal.windows import gaussian
...
>>> T_x, N = 1 / 20, 1000 # 20 Hz sampling rate for 50 s signal
>>> t_x = np.arange(N) * T_x # time indexes for signal
>>> f_i = 1 * np.arctan((t_x - t_x[N // 2]) / 2) + 5 # varying frequency
>>> x = np.sin(2*np.pi*np.cumsum(f_i)*T_x) # the signal
使用的高斯窗口为 50 个样本或 2.5 秒长。参数 mfft=200 在 ShortTimeFFT 中导致频谱过采样 4 倍：
>>> g_std = 8 # standard deviation for Gaussian window in samples
>>> w = gaussian(50, std=g_std, sym=True) # symmetric Gaussian window
>>> SFT = ShortTimeFFT(w, hop=10, fs=1/T_x, mfft=200, scale_to='magnitude')
>>> Sx = SFT.stft(x) # perform the STFT
在图中，信号 x 的时间范围由垂直虚线标记。注意，SFT 产生的值超出 x 的时间范围。左侧和右侧的阴影区域表示由于窗口片段未完全位于 x 的时间范围内而导致的边界效应：
>>> fig1, ax1 = plt.subplots(figsize=(6., 4.)) # enlarge plot a bit
>>> t_lo, t_hi = SFT.extent(N)[:2] # time range of plot
>>> ax1.set_title(rf"STFT ({SFT.m_num*SFT.T:g}$\,s$ Gaussian window, " +
... rf"$\sigma_t={g_std*SFT.T}\,$s)")
>>> ax1.set(xlabel=f"Time $t$ in seconds ({SFT.p_num(N)} slices, " +
... rf"$\Delta t = {SFT.delta_t:g}\,$s)",
... ylabel=f"Freq. $f$ in Hz ({SFT.f_pts} bins, " +
... rf"$\Delta f = {SFT.delta_f:g}\,$Hz)",
... xlim=(t_lo, t_hi))
...
>>> im1 = ax1.imshow(abs(Sx), origin='lower', aspect='auto',
... extent=SFT.extent(N), cmap='viridis')
>>> ax1.plot(t_x, f_i, 'r--', alpha=.5, label='$f_i(t)$')
>>> fig1.colorbar(im1, label="Magnitude $|S_x(t, f)|$")
...
>>> # Shade areas where window slices stick out to the side:
>>> for t0_, t1_ in [(t_lo, SFT.lower_border_end[0] * SFT.T),
... (SFT.upper_border_begin(N)[0] * SFT.T, t_hi)]:
... ax1.axvspan(t0_, t1_, color='w', linewidth=0, alpha=.2)
>>> for t_ in [0, N * SFT.T]: # mark signal borders with vertical line:
... ax1.axvline(t_, color='y', linestyle='--', alpha=0.5)
>>> ax1.legend()
>>> fig1.tight_layout()
>>> plt.show()
使用 istft 重构信号很简单，但请注意应指定 x1 的长度，因为在 hop 步骤中 SFT 的长度会增加：
>>> SFT.invertible # check if invertible
True
>>> x1 = SFT.istft(Sx, k1=N)
>>> np.allclose(x, x1)
True
可以计算信号部分的 SFT：
>>> p_q = SFT.nearest_k_p(N // 2)
>>> Sx0 = SFT.stft(x[:p_q])
>>> Sx1 = SFT.stft(x[p_q:])
在组装连续的 STFT 部分时，需要考虑重叠：
>>> p0_ub = SFT.upper_border_begin(p_q)[1] - SFT.p_min
>>> p1_le = SFT.lower_border_end[1] - SFT.p_min
>>> Sx01 = np.hstack((Sx0[:, :p0_ub],
... Sx0[:, p0_ub:] + Sx1[:, :p1_le],
... Sx1[:, p1_le:]))
>>> np.allclose(Sx01, Sx) # Compare with SFT of complete signal
True
也可以计算信号部分的 itsft：
>>> y_p = SFT.istft(Sx, N//3, N//2)
>>> np.allclose(y_p, x[N//3:N//2])
True
属性：
T
输入信号和窗口的采样间隔。
delta_f
STFT 的频率箱宽度。
delta_t
STFT 的时间增量。
dual_win
规范双窗口。
f
STFT 的频率值。
f_pts
频率轴上的点数。
fac_magnitude
Factor to multiply the STFT values by to scale each frequency slice to a magnitude spectrum.
fac_psd
Factor to multiply the STFT values by to scale each frequency slice to a power spectral density (PSD).
fft_mode
Mode of utilized FFT (‘twosided’, ‘centered’, ‘onesided’ or ‘onesided2X’).
fs
Sampling frequency of input signal and of the window.
hop
Time increment in signal samples for sliding window.
invertible
Check if STFT is invertible.
k_min
The smallest possible signal index of the STFT.
lower_border_end
First signal index and first slice index unaffected by pre-padding.
m_num
Number of samples in window win.
m_num_mid
Center index of window win.
mfft
Length of input for the FFT used - may be larger than window length m_num.
onesided_fft
Return True if a one-sided FFT is used.
p_min
The smallest possible slice index.
phase_shift
如果设置，为每个 FFT 频率片段添加线性相位phase_shift / mfft * f。
正规化应用于窗口函数（‘magnitude’、‘psd’或None）。
窗口函数作为实值或复值 1 维数组。
方法
extent(n[, axes_seq, center_bins]) |
返回最小和最大值的时频值。 |
|---|---|
from_dual(dual_win, hop, fs, *[, fft_mode, ...]) |
仅通过提供双窗口实例化ShortTimeFFT。 |
from_window(win_param, fs, nperseg, noverlap, *) |
使用get_window实例化ShortTimeFFT。 |
istft(S[, k0, k1, f_axis, t_axis]) |
逆短时傅里叶变换。 |
k_max(n) |
信号结束后首个未触及时段的样本索引。 |
nearest_k_p(k[, left]) |
返回最接近的样本索引 k_p，其中 t[k_p] == t[p]成立。 |
p_max(n) |
第一个非重叠的上时段索引，用于n个样本输入。 |
p_num(n) |
n个样本输入信号的时段数。 |
p_range(n[, p0, p1]) |
确定和验证切片索引范围。 |
scale_to(scaling) |
缩放窗口以获得 STFT 的‘magnitude’或‘psd’缩放。 |
spectrogram(x[, y, detr, p0, p1, k_offset, ...]) |
计算频谱图或交叉谱图。 |
stft(x[, p0, p1, k_offset, padding, axis]) |
执行短时傅里叶变换。 |
stft_detrend(x, detr[, p0, p1, k_offset, ...]) |
在每个段之前从中减去趋势的短时傅里叶变换。 |
t(n[, p0, p1, k_offset]) |
用于具有n个样本的输入信号的 STFT 的时间。 |
upper_border_begin(n) |
受后填充影响的第一个信号索引和第一个切片索引。 |
scipy.signal.stft
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.stft.html#scipy.signal.stft
scipy.signal.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=-1, scaling='spectrum')
计算短时傅里叶变换（STFT）。
STFT 可用作量化非平稳信号随时间的频率和相位内容变化的一种方法。
Legacy
此函数被视为传统功能，将不再接收更新。这可能意味着它将在未来的 SciPy 版本中被移除。ShortTimeFFT是一种新的 STFT / ISTFT 实现，具有更多功能。在SciPy 用户指南的教程-STFT部分中可以找到这两种实现的比较。
参数：
x类数组
测量值的时间序列
fs浮点数，可选
x时间序列的采样频率。默认为 1.0。
window字符串或元组或类数组，可选
欲使用的窗口。如果window为字符串或元组，则将其传递给get_window以生成窗口值，默认情况下为 DFT-even。有关窗口和必需参数的列表，请参阅get_window。如果window为类数组，则直接使用它作为窗口，其长度必须为 nperseg。默认为 Hann 窗口。
nperseg整数，可选
每个片段的长度。默认为 256。
noverlap整数，可选
分段之间重叠的点数。如果为None，则noverlap = nperseg // 2。默认为None。当指定时，必须满足 COLA 约束（见下面的说明）。
nfft整数，可选
使用的 FFT 长度，如果需要零填充的 FFT。如果为None，则 FFT 长度为nperseg。默认为None。
detrend字符串或函数或False，可选
指定如何对每个片段进行去趋势化处理。如果detrend为字符串，则将其作为detrend函数的type参数传递。如果它是一个函数，则接受一个片段并返回一个去趋势化的片段。如果detrend为False，则不进行去趋势化处理。默认为False。
return_onesided布尔值，可选
如果为True，则为实数据返回单边谱。如果为False，则返回双边谱。默认为True，但对于复杂数据，始终返回双边谱。
boundary字符串或 None，可选
指定输入信号是否在两端进行扩展，以及如何生成新值，以便将第一个窗段居中在第一个输入点上。这样做有利于在使用窗函数从零开始时重构第一个输入点。有效选项为 ['even', 'odd', 'constant', 'zeros', None]。默认为‘zeros’，用于零填充扩展。例如，对于 nperseg=3， [1, 2, 3, 4] 扩展为 [0, 1, 2, 3, 4, 0]。
paddedbool, optional
指定输入信号是否在末尾进行零填充，以使信号恰好适合整数个窗段，以便所有信号都包含在输出中。默认为True。如果boundary不是None，并且padded为True（默认情况下是这样），填充将在边界扩展之后进行。
axisint, optional
计算 STFT 的轴；默认情况下是在最后一个轴上（即 axis=-1）。
scaling: {‘spectrum’, ‘psd’}
默认的“spectrum”缩放使得Zxx的每个频率线都可以解释为幅度谱。选项“psd”将每行缩放为功率谱密度 - 它允许通过数值积分计算信号的能量 abs(Zxx)**2。
自版本 1.9.0 起新增。
返回：
fndarray
采样频率的数组。
tndarray
段时间的数组。
Zxxndarray
x 的短时傅里叶变换（STFT）。默认情况下，Zxx 的最后一个轴对应于各段时间。
另请参阅
逆短时傅里叶变换
提供更多功能的新 STFT/ISTFT 实现。
检查是否满足恒定重叠添加（COLA）约束
检查是否满足非零重叠添加（NOLA）约束
Welch 方法的功率谱密度。
Welch 方法的谱图。
Welch 方法的交叉谱密度。
不均匀采样数据的 Lomb-Scargle 周期图
注释
为了通过istft中的逆短时傅里叶变换启用 STFT 的反演，信号窗必须遵守“非零重叠加”（NOLA）约束，并且输入信号必须具有完整的窗覆盖（即 (x.shape[axis] - nperseg) % (nperseg-noverlap) == 0）。padded 参数可用于实现此目的。
给定一个时域信号 (x[n])、一个窗口 (w[n]) 和一个跳跃大小 (H) = nperseg - noverlap，时间索引 (t) 处的窗口帧由以下公式给出
[x_{t}[n]=x[n]w[n-tH]]
重叠-添加 (OLA) 重构方程如下所示：
[x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]}]
NOLA 约束确保 OLA 重构方程分母中的每个归一化项都不为零。 可以使用 check_NOLA 来测试 window、nperseg 和 noverlap 是否满足此约束。
新版本 0.19.0 中的新功能。
参考文献
[1]
Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.
[2]
Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified Short-Time Fourier Transform”, IEEE 1984, 10.1109/TASSP.1984.1164317
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
生成一个测试信号，一个振幅为 2 Vrms 的正弦波，其频率围绕 3kHz 缓慢调制，同时受到以 10 kHz 采样的指数衰减幅度的白噪声的影响。
>>> fs = 10e3
>>> N = 1e5
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.01 * fs / 2
>>> time = np.arange(N) / float(fs)
>>> mod = 500*np.cos(2*np.pi*0.25*time)
>>> carrier = amp * np.sin(2*np.pi*3e3*time + mod)
>>> noise = rng.normal(scale=np.sqrt(noise_power),
... size=time.shape)
>>> noise *= np.exp(-time/5)
>>> x = carrier + noise
计算并绘制 STFT 的幅度。
>>> f, t, Zxx = signal.stft(x, fs, nperseg=1000)
>>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')
>>> plt.title('STFT Magnitude')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()
比较信号 x 的能量与其 STFT 的能量：
>>> E_x = sum(x**2) / fs # Energy of x
>>> # Calculate a two-sided STFT with PSD scaling:
>>> f, t, Zxx = signal.stft(x, fs, nperseg=1000, return_onesided=False,
... scaling='psd')
>>> # Integrate numerically over abs(Zxx)**2:
>>> df, dt = f[1] - f[0], t[1] - t[0]
>>> E_Zxx = sum(np.sum(Zxx.real**2 + Zxx.imag**2, axis=0) * df) * dt
>>> # The energy is the same, but the numerical errors are quite large:
>>> np.isclose(E_x, E_Zxx, rtol=1e-2)
True
scipy.signal.istft
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.istft.html#scipy.signal.istft
scipy.signal.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2, scaling='spectrum')
执行反短时傅立叶变换（iSTFT）。
传统
此函数被视为遗留版本，将不再接收更新。这也可能意味着它将在未来的 SciPy 版本中删除。ShortTimeFFT是一个新的 STFT / ISTFT 实现，具有更多功能。可以在SciPy 用户指南的短时傅立叶变换部分找到这些实现的比较。
参数：
Zxxarray_like
要重构的信号的 STFT。如果传递的是纯实数组，则将其转换为复杂数据类型。
fsfloat，可选
时间序列的采样频率。默认为 1.0。
windowstr 或 tuple 或 array_like，可选
所需使用的窗口。如果window是字符串或元组，则将其传递给get_window以生成窗口值，默认为 DFT-even。详见get_window获取窗口列表和所需参数。如果window是 array_like，则直接用作窗口，其长度必须为 nperseg。默认为 Hann 窗口。必须与用于生成 STFT 的窗口匹配，以确保忠实反演。
npersegint，可选
数据点数对应于每个 STFT 段。如果每段数据点数为奇数，或者 STFT 通过nfft > nperseg进行填充，则必须指定此参数。如果为None，则其值取决于Zxx和input_onesided的形状。如果input_onesided为 True，则nperseg=2*(Zxx.shape[freq_axis] - 1)。否则，nperseg=Zxx.shape[freq_axis]。默认为None。
noverlapint，可选
点之间重叠的点数。如果为None，则为段长度的一半。默认为None。在指定时，必须满足 COLA 约束（参见下面的注释），并且应与用于生成 STFT 的参数匹配。默认为None。
nfftint，可选
FFT 点数对应于每个 STFT 段。如果 STFT 通过nfft > nperseg进行填充，则必须指定此参数。如果为None，则默认值与nperseg相同，详见上文，但有一例外：如果input_onesided为 True 且nperseg==2*Zxx.shape[freq_axis] - 1，则nfft也取该值。这种情况允许使用nfft=None正确反演奇数长度未填充的 STFT。默认为None。
input_onesidedbool，可选
如果为True，将输入数组解释为单边 FFT，例如由stft返回的return_onesided=True和numpy.fft.rfft。如果为False，将输入解释为双边 FFT。默认为True。
boundarybool, 可选
指定输入信号是否通过向 stft 提供非None boundary 参数来在其边界上扩展。默认为True。
time_axisint, 可选
STFT 的时间段所在位置；默认为最后一轴（即axis=-1）。
freq_axisint, 可选
STFT 的频率轴所在位置；默认为倒数第二轴（即axis=-2）。
scaling: {‘spectrum’, ‘psd’}
默认的'spectrum'缩放允许解释Zxx的每个频率线为幅度谱。'psd'选项将每行缩放到功率谱密度 - 允许通过数值积分计算信号的能量 abs(Zxx)**2。
返回:
tndarray
输出数据数组的时间。
xndarray
Zxx的逆短时傅立叶变换。
另请参阅
stft
短时傅立叶变换
ShortTimeFFT
更多功能的新 STFT/ISTFT 实现。
check_COLA
检查是否满足 Constant OverLap Add (COLA)约束
check_NOLA
检查是否满足 Nonzero Overlap Add (NOLA)约束
注意事项
为了通过istft反转 STFT 以进行反 STFT，信号窗必须遵守“非零重叠添加”（NOLA）约束：
[\sum_{t}w^{2}[n-tH] \ne 0]
这确保了出现在重叠添加重建方程分母中的归一化因子
[x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]}]
不为零。使用check_NOLA函数可以检查 NOLA 约束。
已修改的 STFT（通过掩蔽或其他方式）不能保证与确切可实现信号对应。该函数通过最小二乘估计算法实现了 iSTFT，该算法详细说明见[2]，其生成的信号最小化了返回信号的 STFT 和修改后 STFT 之间的均方误差。
版本 0.19.0 中的新功能。
参考文献
[1]
Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “离散时间信号处理”，Prentice Hall，1999 年。
[2]
Daniel W. Griffin, Jae S. Lim “从修改后的短时傅里叶变换估计信号”, IEEE 1984, 10.1109/TASSP.1984.1164317
示例
>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()
生成一个测试信号，一个 2 Vrms 的 50Hz 正弦波，受 1024 Hz 采样的 0.001 V**2/Hz 白噪声的影响。
>>> fs = 1024
>>> N = 10*fs
>>> nperseg = 512
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / float(fs)
>>> carrier = amp * np.sin(2*np.pi*50*time)
>>> noise = rng.normal(scale=np.sqrt(noise_power),
... size=time.shape)
>>> x = carrier + noise
计算 STFT，并绘制其幅度
>>> f, t, Zxx = signal.stft(x, fs=fs, nperseg=nperseg)
>>> plt.figure()
>>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')
>>> plt.ylim([f[1], f[-1]])
>>> plt.title('STFT Magnitude')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.yscale('log')
>>> plt.show()
将幅度为载波幅度的 10%或更少的分量置零，然后通过逆 STFT 转换回时间序列
>>> Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0)
>>> _, xrec = signal.istft(Zxx, fs)
将清理后的信号与原始和真实的载波信号进行比较。
>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([2, 2.1])
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()
注意，清理后的信号并不像原始信号那样突然开始，因为某些瞬态的系数也被移除了：
>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([0, 0.1])
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()
scipy.signal.check_COLA
scipy.signal.check_COLA(window, nperseg, noverlap, tol=1e-10)
检查常量重叠添加（COLA）约束是否满足。
参数：
window字符串或元组或 array_like
所需使用的窗口。如果 window 是字符串或元组，则将其传递给 get_window 以生成窗口值，默认情况下为 DFT-even。有关窗口和所需参数的列表，请参见 get_window。如果 window 是 array_like，则将其直接用作窗口，其长度必须为 nperseg。
nperseg整数
每个片段的长度。
noverlap整数
段之间重叠的点数。
tol浮点数，可选
每个频段加权和与中位数频段和的允许方差。
返回：
verdict布尔值
True 如果选择的组合在 tol 范围内满足 COLA，否则 False
请参见
check_NOLA
检查是否满足非零重叠添加（NOLA）约束
stft
短时傅里叶变换
istft
逆短时傅里叶变换
注释
为了通过逆短时傅里叶变换中的逆 STFT 实现 STFT 的反演，在 istft 中，只需确保信号窗口符合“常数重叠添加”（COLA）的约束即可。这确保了输入数据中的每个点都被等权重，从而避免混叠，并允许完全重建。
满足 COLA 的一些窗口示例：
-
重叠为 0、1/2、2/3、3/4 等的矩形窗口
-
Bartlett 窗口在 1/2、3/4、5/6 等重叠时
-
Hann 窗口在 1/2、2/3、3/4 等重叠时
-
任何 Blackman 家族窗口的 2/3 重叠
-
任何具有
noverlap = nperseg-1的窗口
在[2]中可以找到其他窗口的非常全面的列表，在“幅度平坦度”为单位时满足 COLA 条件。
从版本 0.19.0 开始新增。
参考文献
[1]
Julius O. Smith III，《Spectral Audio Signal Processing》，W3K Publishing，2011 年，ISBN 978-0-9745607-3-1。
[2]
G. Heinzel, A. Ruediger and R. Schilling，《Spectrum and spectral density estimation by the Discrete Fourier transform (DFT)，including a comprehensive list of window functions and some new at-top windows》，2002 年，hdl.handle.net/11858/00-001M-0000-0013-557A-5
示例
>>> from scipy import signal
确认 75%（3/4）重叠的矩形窗口的 COLA 条件：
>>> signal.check_COLA(signal.windows.boxcar(100), 100, 75)
True
对于 25%（1/4）重叠，COLA 不成立：
>>> signal.check_COLA(signal.windows.boxcar(100), 100, 25)
False
“对称”Hann 窗口（用于滤波器设计）不满足 COLA：
>>> signal.check_COLA(signal.windows.hann(120, sym=True), 120, 60)
False
“周期性”或“DFT-even”的 Hann 窗口（用于 FFT 分析）在 1/2、2/3、3/4 等重叠情况下是 COLA 的：
>>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 60)
True
>>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 80)
True
>>> signal.check_COLA(signal.windows.hann(120, sym=False), 120, 90)
True
scipy.signal.check_NOLA
scipy.signal.check_NOLA(window, nperseg, noverlap, tol=1e-10)
检查是否满足非零重叠添加（NOLA）约束。
参数：
windowstr 或 tuple 或 array_like
要使用的期望窗口。如果window是字符串或元组，则将其传递给get_window以生成窗口值，默认情况下为 DFT-even。有关窗口和所需参数的列表，请参见get_window。如果window是 array_like，则将其直接用作窗口，其长度必须为 nperseg。
npersegint
每个段的长度。
noverlapint
分段之间重叠的点数。
tolfloat, 可选
一个频段的加权和与中位数频段的加权和的允许方差。
返回：
verdictbool
如果所选组合符合tol内的 NOLA 约束条件，则返回True，否则返回False
参见
check_COLA
检查是否满足恒定重叠添加（COLA）约束
stft
短时傅立叶变换
istft
逆短时傅立叶变换
注释
为了通过istft中的逆 STFT 启用 STFT 的反演，信号窗必须遵守“非零重叠添加”（NOLA）约束：
[\sum_{t}w^{2}[n-tH] \ne 0]
对于所有的(n)，其中(w)是窗口函数，(t)是帧索引，(H)是跨步大小（(H) = nperseg - noverlap）。
这确保了重叠添加反演方程中的归一化因子不为零。只有非常异常的窗口才会不满足 NOLA 约束。
1.2.0 版中的新功能。
参考文献
[1]
Julius O. Smith III，“音频信号谱分析”，W3K Publishing，2011 年，ISBN 978-0-9745607-3-1。
[2]
G. Heinzel, A. Ruediger and R. Schilling, “离散傅立叶变换（DFT）估计的频谱和谱密度，包括详细的窗函数列表和一些新的顶部窗口”，2002 年，hdl.handle.net/11858/00-001M-0000-0013-557A-5
示例
>>> import numpy as np
>>> from scipy import signal
确认 75%（3/4）重叠的矩形窗口的 NOLA 条件：
>>> signal.check_NOLA(signal.windows.boxcar(100), 100, 75)
True
对于 25%（1/4）重叠，NOLA 也成立：
>>> signal.check_NOLA(signal.windows.boxcar(100), 100, 25)
True
“对称”Hann 窗（用于滤波器设计）也满足 NOLA：
>>> signal.check_NOLA(signal.windows.hann(120, sym=True), 120, 60)
True
只要有重叠，就需要非常异常的窗口才能不满足 NOLA：
>>> w = np.ones(64, dtype="float")
>>> w[::2] = 0
>>> signal.check_NOLA(w, 64, 32)
False
如果重叠不足，末端带有零的窗口将无法工作：
>>> signal.check_NOLA(signal.windows.hann(64), 64, 0)
False
>>> signal.check_NOLA(signal.windows.hann(64), 64, 1)
False
>>> signal.check_NOLA(signal.windows.hann(64), 64, 2)
True
scipy.signal.czt
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/czt-function.html#scipy.signal.czt
scipy.signal.czt(x, m=None, w=None, a=1 + 0j, *, axis=-1)
计算 Z 平面中螺旋周围的频率响应。
参数：
x：array
要变换的信号。
m：int，可选
所需输出点的数量。默认为输入数据的长度。
w：complex，可选
在每个步骤中点之间的比率。这必须精确，否则累积误差将使输出序列的尾部退化。默认为整个单位圆周围均匀分布的点。
a：complex，可选
复平面中的起始点。默认为 1+0j。
axis：int，可选
计算 FFT 的轴。如果未给出，则使用最后一个轴。
返回：
out：ndarray
一个与 x 相同尺寸的数组，但是变换轴的长度设置为 m。
另见
创建可调用的啁啾 z 变换函数的类。
部分 FFT 计算的便捷函数。
注释
默认值选取为 signal.czt(x) 等同于 fft.fft(x)，如果 m > len(x)，则 signal.czt(x, m) 等同于 fft.fft(x, m)。
如果需要重复变换，请使用 CZT 来构建一个专门的变换函数，可以在不重新计算常量的情况下重复使用。
一个示例应用是在系统识别中，重复评估系统的 Z 变换的小片段，以精炼估计极点的真实位置。
参考文献
[1]
Steve Alan Shilling，《啁啾 z 变换及其应用研究》，第 20 页（1970）krex.k-state.edu/dspace/bitstream/handle/2097/7844/LD2668R41972S43.pdf
示例
生成一个正弦波：
>>> import numpy as np
>>> f1, f2, fs = 8, 10, 200 # Hz
>>> t = np.linspace(0, 1, fs, endpoint=False)
>>> x = np.sin(2*np.pi*t*f2)
>>> import matplotlib.pyplot as plt
>>> plt.plot(t, x)
>>> plt.axis([0, 1, -1.1, 1.1])
>>> plt.show()
其离散傅里叶变换的能量全集中在单一频率箱中：
>>> from scipy.fft import rfft, rfftfreq
>>> from scipy.signal import czt, czt_points
>>> plt.plot(rfftfreq(fs, 1/fs), abs(rfft(x)))
>>> plt.margins(0, 0.1)
>>> plt.show()
然而，如果正弦波是对数衰减的：
>>> x = np.exp(-t*f1) * np.sin(2*np.pi*t*f2)
>>> plt.plot(t, x)
>>> plt.axis([0, 1, -1.1, 1.1])
>>> plt.show()
DFT 将具有频谱泄漏：
>>> plt.plot(rfftfreq(fs, 1/fs), abs(rfft(x)))
>>> plt.margins(0, 0.1)
>>> plt.show()
尽管 DFT 总是在单位圆周围采样 Z 变换，啁啾 z 变换允许我们沿任何对数螺旋（例如半径小于单位的圆）采样 Z 变换：
>>> M = fs // 2 # Just positive frequencies, like rfft
>>> a = np.exp(-f1/fs) # Starting point of the circle, radius < 1
>>> w = np.exp(-1j*np.pi/M) # "Step size" of circle
>>> points = czt_points(M + 1, w, a) # M + 1 to include Nyquist
>>> plt.plot(points.real, points.imag, '.')
>>> plt.gca().add_patch(plt.Circle((0,0), radius=1, fill=False, alpha=.3))
>>> plt.axis('equal'); plt.axis([-1.05, 1.05, -0.05, 1.05])
>>> plt.show()
使用正确的半径，这将转换衰减正弦波（以及具有相同衰减率的其他波形），而不会出现频谱泄漏：
>>> z_vals = czt(x, M + 1, w, a) # Include Nyquist for comparison to rfft
>>> freqs = np.angle(points)*fs/(2*np.pi) # angle = omega, radius = sigma
>>> plt.plot(freqs, abs(z_vals))
>>> plt.margins(0, 0.1)
>>> plt.show()
scipy.signal.zoom_fft
原文链接：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.zoom_fft.html#scipy.signal.zoom_fft
scipy.signal.zoom_fft(x, fn, m=None, *, fs=2, endpoint=False, axis=-1)
仅计算范围在 fn 中的频率的 x 的 DFT。
参数：
x：数组
要变换的信号。
fn：类似数组
长度为 2 的序列 [f1, f2] 给出频率范围，或者一个标量，其中假设范围为 [0, fn]。
m：整数，可选
要评估的点数。默认为 x 的长度。
fs：浮点数，可选
采样频率。例如，如果 fs=10 表示 10 kHz，那么 f1 和 f2 也应该以 kHz 表示。默认采样频率为 2，因此 f1 和 f2 应在 [0, 1] 范围内以保持变换低于奈奎斯特频率。
endpoint：布尔值，可选
如果为 True，f2 是最后一个样本。否则，不包括它。默认为 False。
axis：整数，可选
计算 FFT 的轴。如果未给出，则使用最后一个轴。
返回：
out：ndarray
转换后的信号。傅里叶变换将在点 f1, f1+df, f1+2df, …, f2 处计算，其中 df=(f2-f1)/m。
另请参阅
ZoomFFT
创建一个可调用的部分 FFT 函数的类。
注意事项
默认选择这样，使得 signal.zoom_fft(x, 2) 等价于 fft.fft(x)，如果 m > len(x)，那么 signal.zoom_fft(x, 2, m) 等价于 fft.fft(x, m)。
要绘制结果变换的幅度图，请使用：
plot(linspace(f1, f2, m, endpoint=False), abs(zoom_fft(x, [f1, f2], m)))
如果需要重复变换，请使用 ZoomFFT 构建一个专门的变换函数，可以在不重新计算常数的情况下重复使用。
示例
要绘制变换结果，请使用类似以下的方法：
>>> import numpy as np
>>> from scipy.signal import zoom_fft
>>> t = np.linspace(0, 1, 1021)
>>> x = np.cos(2*np.pi*15*t) + np.sin(2*np.pi*17*t)
>>> f1, f2 = 5, 27
>>> X = zoom_fft(x, [f1, f2], len(x), fs=1021)
>>> f = np.linspace(f1, f2, len(x))
>>> import matplotlib.pyplot as plt
>>> plt.plot(f, 20*np.log10(np.abs(X)))
>>> plt.show()
scipy.signal.CZT
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.CZT.html#scipy.signal.CZT
class scipy.signal.CZT(n, m=None, w=None, a=1 + 0j)
创建一个可调用的啁啾变换函数。
转换以计算螺旋周围的频率响应。此类对象是可调用的，可以在其输入上计算啁啾变换。此对象预先计算给定变换中使用的恒定啁啾。
参数：
nint
信号的大小。
mint，可选
所需的输出点数。默认为n。
wcomplex，可选
每个步骤中点之间的比例。这必须是精确的，否则累积误差将降低输出序列的尾部。默认为整个单位圆周围均匀分布的点。
acomplex，可选
复平面中的起始点。默认为 1+0j。
返回：
fCZT
用于在x上计算啁啾变换的可调用对象f(x, axis=-1)。
另请参见
czt
用于快速计算 CZT 的便利函数。
ZoomFFT
创建可调用的部分 FFT 函数的类。
注意事项
默认值选为使f(x)等同于fft.fft(x)，如果m > len(x)，则使f(x, m)等同于fft.fft(x, m)。
如果w不位于单位圆上，则变换将围绕指数增长半径的螺旋进行。无论如何，角度将线性增加。
对于位于单位圆上的变换，当使用ZoomFFT时，精度更高，因为w中的任何数值误差在长数据长度上累积，偏离单位圆。
与等效零填充 FFT 相比，啁啾变换可能更快。尝试使用您自己的数组大小进行测试。
然而，啁啾变换的精度明显低于等效的零填充 FFT。
由于此 CZT 使用 Bluestein 算法实现，因此可以在 O(N log N)时间内计算大素数长度的傅里叶变换，而不是直接 DFT 计算所需的 O(N**2)时间。(scipy.fft也使用 Bluestein 的算法。)
（“啁啾变换”名称来自 Bluestein 算法中使用的啁啾。它不像其他带有“啁啾”名称的变换那样将信号分解为啁啾。）
参考文献
[1]
Leo I. Bluestein，“离散傅里叶变换计算的线性滤波方法”，东北电子研究与工程会议记录 10，218-219（1968）。
[2]
Rabiner、Schafer 和 Rader，“啁啾变换算法及其应用”，贝尔系统技术杂志 48，1249-1292（1969）。
示例
计算多个素数长度 FFT：
>>> from scipy.signal import CZT
>>> import numpy as np
>>> a = np.random.rand(7)
>>> b = np.random.rand(7)
>>> c = np.random.rand(7)
>>> czt_7 = CZT(n=7)
>>> A = czt_7(a)
>>> B = czt_7(b)
>>> C = czt_7(c)
显示计算 FFT 的点：
>>> czt_7.points()
array([ 1.00000000+0.j , 0.62348980+0.78183148j,
-0.22252093+0.97492791j, -0.90096887+0.43388374j,
-0.90096887-0.43388374j, -0.22252093-0.97492791j,
0.62348980-0.78183148j])
>>> import matplotlib.pyplot as plt
>>> plt.plot(czt_7.points().real, czt_7.points().imag, 'o')
>>> plt.gca().add_patch(plt.Circle((0,0), radius=1, fill=False, alpha=.3))
>>> plt.axis('equal')
>>> plt.show()
方法
__call__(x, *[, axis]) |
计算信号的奇异变换。 |
|---|---|
points() |
返回进行奇异变换的点的位置。 |
scipy.signal.ZoomFFT
原文：
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.ZoomFFT.html#scipy.signal.ZoomFFT
class scipy.signal.ZoomFFT(n, fn, m=None, *, fs=2, endpoint=False)
创建一个可调用的变焦 FFT 变换函数。
这是圆周单位周围等间距频率的啁啾变换（CZT）的特化，用于比计算整个 FFT 并截断更有效地计算 FFT 的一部分。
参数：
n整数
信号的大小。
fn类似数组
长度为 2 的序列[f1, f2]表示频率范围，或者标量，假定范围[0, fn]。
m整数，可选
评估点数。默认为n。
fs浮点数，可选
采样频率。例如，如果fs=10表示 10 kHz，则f1和f2也应以 kHz 为单位。默认的采样频率为 2，因此f1和f2的范围应在[0, 1]之间，以使变换保持在奈奎斯特频率以下。
endpoint布尔值，可选
如果为 True，则f2为最后一个样本。否则，不包括在内。默认为 False。
返回：
fZoomFFT
可调用对象f(x, axis=-1)用于计算x上的变焦 FFT。
另请参阅
用于计算变焦 FFT 的便捷函数。
注：
默认设置使得f(x, 2)等同于fft.fft(x)，如果m > len(x)，那么f(x, 2, m)等同于fft.fft(x, m)。
采样频率是信号x中样本之间的时间步长的倒数。单位圆对应从 0 到采样频率的频率。默认的采样频率为 2，因此f1和f2的值应在范围[0, 1)内，以保持变换在奈奎斯特频率以下。
请记住，变焦 FFT 只能插值现有 FFT 的点。它无法帮助解决两个分开的附近频率。只能通过增加采集时间来增加频率分辨率。
这些函数使用 Bluestein 算法实现（就像scipy.fft一样）。[2]
参考文献
[1]
Steve Alan Shilling，“啁啾变换及其应用研究”，第 29 页（1970 年）krex.k-state.edu/dspace/bitstream/handle/2097/7844/LD2668R41972S43.pdf
Leo I. Bluestein，“离散傅立叶变换的线性滤波方法”，东北电子研究与工程会议记录第 10 卷，218-219 页（1968 年）。
示例
要绘制变换结果，请使用类似以下的内容：
>>> import numpy as np
>>> from scipy.signal import ZoomFFT
>>> t = np.linspace(0, 1, 1021)
>>> x = np.cos(2*np.pi*15*t) + np.sin(2*np.pi*17*t)
>>> f1, f2 = 5, 27
>>> transform = ZoomFFT(len(x), [f1, f2], len(x), fs=1021)
>>> X = transform(x)
>>> f = np.linspace(f1, f2, len(x))
>>> import matplotlib.pyplot as plt
>>> plt.plot(f, 20*np.log10(np.abs(X)))
>>> plt.show()
方法
__call__ (x, *[, axis]) |
计算信号的奇异变换。 |
|---|---|
points() |
返回进行奇异变换计算的点。 |
scipy.signal.czt_points
scipy.signal.czt_points(m, w=None, a=1 + 0j)
返回进行啁啾变换的点。
参数：
mint
所需点的数量。
w复数，可选
每个步骤中点的比率。默认为均匀分布在整个单位圆周围的点。
a复数，可选
复平面中的起始点。默认为 1+0j。
返回：
outndarray
Z 平面中的点，CZT 在调用时以复数m、w和a作为参数进行 z 变换采样。
另请参阅：
CZT
创建一个可调用的啁啾变换函数的类。
czt
用于快速计算 CZT 的便捷函数。
示例
绘制 16 点 FFT 的点：
>>> import numpy as np
>>> from scipy.signal import czt_points
>>> points = czt_points(16)
>>> import matplotlib.pyplot as plt
>>> plt.plot(points.real, points.imag, 'o')
>>> plt.gca().add_patch(plt.Circle((0,0), radius=1, fill=False, alpha=.3))
>>> plt.axis('equal')
>>> plt.show()
和一个穿过单位圆的 91 点对数螺旋：
>>> m, w, a = 91, 0.995*np.exp(-1j*np.pi*.05), 0.8*np.exp(1j*np.pi/6)
>>> points = czt_points(m, w, a)
>>> plt.plot(points.real, points.imag, 'o')
>>> plt.gca().add_patch(plt.Circle((0,0), radius=1, fill=False, alpha=.3))
>>> plt.axis('equal')
>>> plt.show()
稀疏矩阵（scipy.sparse）
SciPy 二维稀疏数组包，用于数值数据。
注意
此软件包正在切换到与 NumPy 数组兼容的数组接口，而不再使用旧的矩阵接口。我们建议您对所有新工作使用数组对象（bsr_array, coo_array 等）。
在使用数组接口时，请注意：
-
x * y现在不再执行矩阵乘法，而是执行元素级乘法（与 NumPy 数组类似）。为了使代码同时适用于数组和矩阵，使用x @ y来进行矩阵乘法。 -
诸如 sum 的操作，原先生成密集矩阵，现在生成数组，其乘法行为类似但有所不同。
-
稀疏数组目前必须是二维的。这也意味着这些对象上的所有 切片 操作必须产生二维结果，否则将导致错误。这将在未来版本中解决。
构造实用程序（eye, kron, random, diags 等）尚未移植完成，但可以将其结果封装成数组：
A = csr_array(eye(3))
内容
稀疏数组类
bsr_array(arg1[, shape, dtype, copy, blocksize]) |
块稀疏行（Block Sparse Row）格式的稀疏数组。 |
|---|---|
coo_array(arg1[, shape, dtype, copy]) |
COOrdinate 格式的稀疏数组。 |
csc_array(arg1[, shape, dtype, copy]) |
压缩稀疏列（Compressed Sparse Column）数组。 |
csr_array(arg1[, shape, dtype, copy]) |
压缩稀疏行（Compressed Sparse Row）数组。 |
dia_array(arg1[, shape, dtype, copy]) |
带有对角线存储的稀疏数组。 |
dok_array(arg1[, shape, dtype, copy]) |
基于键的字典（Dictionary Of Keys）稀疏数组。 |
lil_array(arg1[, shape, dtype, copy]) |
基于行的列表（LIst of Lists）稀疏数组。 |
sparray() |
该类为所有稀疏数组提供基类。 |
稀疏矩阵类：
bsr_matrix(arg1[, shape, dtype, copy, blocksize]) |
块稀疏行格式的稀疏矩阵。 |
|---|---|
coo_matrix(arg1[, shape, dtype, copy]) |
COO 格式的稀疏矩阵。 |
csc_matrix(arg1[, shape, dtype, copy]) |
压缩稀疏列矩阵。 |
csr_matrix(arg1[, shape, dtype, copy]) |
压缩稀疏行矩阵。 |
dia_matrix(arg1[, shape, dtype, copy]) |
带有对角线存储的稀疏矩阵。 |
dok_matrix(arg1[, shape, dtype, copy]) |
基于键的字典稀疏矩阵。 |
lil_matrix(arg1[, shape, dtype, copy]) |
基于行的链表稀疏矩阵。 |
spmatrix() |
该类为所有稀疏矩阵类提供基类。 |
函数：
构建稀疏数组：
diags_array(diagonals, /, *[, offsets, ...]) |
从对角线构造稀疏数组。 |
|---|---|
eye_array(m[, n, k, dtype, format]) |
稀疏数组格式中的单位矩阵 |
random_array(shape, *[, density, format, ...]) |
返回一个 0, 1) 范围内均匀随机数的稀疏数组 |
[block_array(blocks, *[, format, dtype]) |
从稀疏子块构建稀疏数组 |
构建稀疏矩阵：
eye(m[, n, k, dtype, format]) |
对角线上有 1 的稀疏矩阵 |
|---|---|
identity(n[, dtype, format]) |
稀疏格式中的单位矩阵 |
diags(diagonals[, offsets, shape, format, dtype]) |
从对角线构造稀疏矩阵。 |
spdiags(data, diags[, m, n, format]) |
从对角线返回稀疏矩阵。 |
bmat(blocks[, format, dtype]) |
从稀疏子块构建稀疏数组或矩阵 |
random(m, n[, density, format, dtype, ...]) |
生成给定形状和密度的稀疏矩阵，值为随机分布。 |
rand(m, n[, density, format, dtype, ...]) |
生成给定形状和密度的稀疏矩阵，值均匀分布。 |
从更小的结构（数组或矩阵）构建更大的结构
kron(A, B[, format]) |
稀疏矩阵 A 和 B 的 Kronecker 乘积 |
|---|---|
kronsum(A, B[, format]) |
方阵稀疏矩阵 A 和 B 的 Kronecker 和 |
block_diag(mats[, format, dtype]) |
从提供的矩阵构建块对角稀疏矩阵或数组 |
tril(A[, k, format]) |
返回稀疏数组或矩阵的下三角部分 |
triu(A[, k, format]) |
返回稀疏数组或矩阵的上三角部分 |
hstack(blocks[, format, dtype]) |
水平堆叠稀疏矩阵（按列） |
vstack(blocks[, format, dtype]) |
垂直堆叠稀疏数组（按行） |
保存和加载稀疏矩阵：
save_npz(file, matrix[, compressed]) |
使用.npz格式将稀疏矩阵或数组保存到文件中。 |
|---|---|
load_npz(file) |
使用.npz格式从文件加载稀疏数组/矩阵。 |
稀疏工具：
find(A) |
返回矩阵非零元素的索引和值 |
|---|
辨识稀疏数组：
-
使用 isinstance(A, sp.sparse.sparray) 检查是否为数组或矩阵。
-
使用 A.format == ‘csr’ 来检查稀疏格式
辨识稀疏矩阵：
issparse(x) |
x 是否为稀疏数组或稀疏矩阵类型？ |
|---|---|
isspmatrix(x) |
x 是否为稀疏矩阵类型？ |
isspmatrix_csc(x) |
x 是否为 csc_matrix 类型？ |
isspmatrix_csr(x) |
x 是否为 csr_matrix 类型？ |
isspmatrix_bsr(x) |
x 是否为 bsr_matrix 类型？ |
isspmatrix_lil(x) |
x 是否为 lil_matrix 类型？ |
isspmatrix_dok(x) |
x 是否为 dok_array 类型？ |
isspmatrix_coo(x) |
x 是否为 coo_matrix 类型？ |
isspmatrix_dia(x) |
x 是否为 dia_matrix 类型？ |
子模块
csgraph |
压缩稀疏图例程 (scipy.sparse.csgraph) |
|---|---|
linalg |
稀疏线性代数 (scipy.sparse.linalg) |
异常情况
SparseEfficiencyWarning |
|
|---|---|
SparseWarning |
使用信息
有七种可用的稀疏数组类型：
csc_array: 压缩稀疏列格式csr_array: 压缩稀疏行格式bsr_array: 块稀疏行格式lil_array: 列表列表格式dok_array: 键字典格式coo_array: COO 格式（即 IJV，三元组格式）dia_array: 对角线格式
要高效构造数组，请使用dok_array或者lil_array。lil_array类支持基本切片和与 NumPy 数组类似语法的花式索引。正如下文所示，COO 格式也可用于高效构造数组。尽管它们与 NumPy 数组相似，强烈不建议直接在这些数组上使用 NumPy 函数，因为 NumPy 可能无法正确转换它们以进行计算，导致意外（和错误）的结果。如果确实要在这些数组上应用 NumPy 函数，请首先检查 SciPy 是否有适用于给定稀疏数组类的自己的实现，或者在应用方法之前将稀疏数组转换为 NumPy 数组（例如，使用类的toarray方法）。
要执行诸如乘法或求逆之类的操作，首先将数组转换为 CSC 或 CSR 格式。lil_array格式是基于行的，因此转换为 CSR 是有效的，而转换为 CSC 则不太有效。
在 CSR、CSC 和 COO 格式之间的所有转换都是高效的、线性时间的操作。
矩阵向量乘积
要在稀疏数组和向量之间进行向量乘积，简单地使用数组的dot方法，如其文档字符串中所述：
>>> import numpy as np
>>> from scipy.sparse import csr_array
>>> A = csr_array([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
>>> v = np.array([1, 0, -1])
>>> A.dot(v)
array([ 1, -3, -1], dtype=int64)
警告
从 NumPy 1.7 版本开始，np.dot不知道稀疏数组，因此使用它将导致意外的结果或错误。应该首先获得相应的密集数组：
>>> np.dot(A.toarray(), v)
array([ 1, -3, -1], dtype=int64)
但这样一来，所有的性能优势都会丧失。
CSR 格式特别适合快速矩阵向量乘积。
示例 1
构造一个 1000x1000 的lil_array并给它添加一些值：
>>> from scipy.sparse import lil_array
>>> from scipy.sparse.linalg import spsolve
>>> from numpy.linalg import solve, norm
>>> from numpy.random import rand
>>> A = lil_array((1000, 1000))
>>> A[0, :100] = rand(100)
>>> A[1, 100:200] = A[0, :100]
>>> A.setdiag(rand(1000))
现在将其转换为 CSR 格式并解决 A x = b 得到 x：
>>> A = A.tocsr()
>>> b = rand(1000)
>>> x = spsolve(A, b)
将其转换为密集数组并求解，并检查结果是否相同：
>>> x_ = solve(A.toarray(), b)
现在我们可以计算误差的范数：
>>> err = norm(x-x_)
>>> err < 1e-10
True
应该很小 😃
示例 2
在 COO 格式中构造一个数组：
>>> from scipy import sparse
>>> from numpy import array
>>> I = array([0,3,1,0])
>>> J = array([0,3,1,2])
>>> V = array([4,5,7,9])
>>> A = sparse.coo_array((V,(I,J)),shape=(4,4))
注意索引不需要排序。
在转换为 CSR 或 CSC 时，重复的(i,j)条目将被求和。
>>> I = array([0,0,1,3,1,0,0])
>>> J = array([0,2,1,3,1,0,0])
>>> V = array([1,1,1,1,1,1,1])
>>> B = sparse.coo_array((V,(I,J)),shape=(4,4)).tocsr()
这对于构造有限元刚度和质量矩阵非常有用。
进一步细节
CSR 列索引不一定排序。同样适用于 CSC 行索引。当需要排序索引时，请使用.sorted_indices()和.sort_indices()方法（例如，将数据传递给其他库时）。
