信号处理 (scipy.signal)

信号处理工具箱目前包含一些滤波函数、一套有限的滤波器设计工具，以及一些用于一维和二维数据的 B 样条插值算法。虽然 B 样条算法在技术上可以放在插值类别下，但它们之所以包含在这里，是因为它们仅适用于等间隔数据，并大量使用滤波理论和传递函数形式来提供快速的 B 样条变换。要理解本节，您需要理解 SciPy 中的信号是实数或复数数组。

B 样条

B 样条是用 B 样条系数和节点点来近似表示有限域上的连续函数。如果节点点以间隔 \(\Delta x\) 等间距分布，则一维函数的 B 样条近似是有限基展开。

\[y\left(x\right)\approx\sum_{j}c_{j}\beta^{o}\left(\frac{x}{\Delta x}-j\right).\]

在节点间距为 \(\Delta x\) 和 \(\Delta y\) 的二维情况下，函数表示为

\[z\left(x,y\right)\approx\sum_{j}\sum_{k}c_{jk}\beta^{o}\left(\frac{x}{\Delta x}-j\right)\beta^{o}\left(\frac{y}{\Delta y}-k\right).\]

在这些表达式中，\(\beta^{o}\left(\cdot\right)\) 是阶数为 \(o\) 的空间有限 B 样条基函数。等间距节点点和等间距数据点的要求允许开发快速（逆滤波）算法，用于从样本值 \(y_{n}\) 确定系数 \(c_{j}\)。与一般的样条插值算法不同，这些算法可以快速找到大型图像的样条系数。

通过 B 样条基函数表示一组样本的优点是，可以从样条系数相对容易地计算连续域运算符（导数、重采样、积分等），这些运算符假设数据样本是从底层连续函数中提取的。例如，样条的二阶导数为

\[y{}^{\prime\prime}\left(x\right)=\frac{1}{\Delta x^{2}}\sum_{j}c_{j}\beta^{o\prime\prime}\left(\frac{x}{\Delta x}-j\right).\]

利用 B 样条的性质：

\[\frac{d^{2}\beta^{o}\left(w\right)}{dw^{2}}=\beta^{o-2}\left(w+1\right)-2\beta^{o-2}\left(w\right)+\beta^{o-2}\left(w-1\right),\]

可以看出：

\[y^{\prime\prime}\left(x\right)=\frac{1}{\Delta x^{2}}\sum_{j}c_{j}\left[\beta^{o-2}\left(\frac{x}{\Delta x}-j+1\right)-2\beta^{o-2}\left(\frac{x}{\Delta x}-j\right)+\beta^{o-2}\left(\frac{x}{\Delta x}-j-1\right)\right].\]

如果 \(o=3\)，则在采样点处：

\begin{eqnarray*} \Delta x^{2}\left.y^{\prime}\left(x\right)\right|_{x=n\Delta x} & = & \sum_{j}c_{j}\delta_{n-j+1}-2c_{j}\delta_{n-j}+c_{j}\delta_{n-j-1},\\ & = & c_{n+1}-2c_{n}+c_{n-1}.\end{eqnarray*}

因此，可以很容易地从样条拟合中计算出二阶导数信号。如果需要，可以找到平滑样条，以使二阶导数对随机误差不太敏感。

精明的读者可能已经注意到，数据样本通过卷积运算符与节点系数相关，因此，只需将采样的 B 样条函数进行简单卷积，即可从样条系数中恢复原始数据。卷积的输出会因边界的处理方式而变化（随着数据集中维数的增加，这一点变得越来越重要）。信号处理子包中与 B 样条相关的算法假定镜像对称边界条件。因此，样条系数是基于该假设计算的，并且假设样条系数也是镜像对称的，则可以从样条系数中精确地恢复数据样本。

目前，该软件包提供了用于从一维和二维等间隔样本确定二阶和三阶三次样条系数的函数 (qspline1d, qspline2d, cspline1d, cspline2d)。对于较大的 \(o\)，B 样条基函数可以用均值为零且标准差等于 \(\sigma_{o}=\left(o+1\right)/12\) 的高斯函数很好地近似表示

\[\beta^{o}\left(x\right)\approx\frac{1}{\sqrt{2\pi\sigma_{o}^{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{o}}\right).\]

还提供了计算任意 \(x\) 和 \(o\) 的高斯函数的函数（ gauss_spline ）。以下代码和图使用样条滤波来计算浣熊脸部的边缘图像（平滑样条的二阶导数），这是一个由命令 scipy.datasets.face 返回的数组。命令 sepfir2d 用于将具有镜像对称边界条件的可分离二维 FIR 滤波器应用于样条系数。此函数非常适合从样条系数重建样本，并且比 convolve2d 快，后者会卷积任意二维滤波器并允许选择镜像对称边界条件。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = datasets.face(gray=True).astype(np.float32)
>>> derfilt = np.array([1.0, -2, 1.0], dtype=np.float32)
>>> ck = signal.cspline2d(image, 8.0)
>>> deriv = (signal.sepfir2d(ck, derfilt, [1]) +
...          signal.sepfir2d(ck, [1], derfilt))

或者，我们可以这样做：

laplacian = np.array([[0,1,0], [1,-4,1], [0,1,0]], dtype=np.float32)
deriv2 = signal.convolve2d(ck,laplacian,mode='same',boundary='symm')

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(deriv)
>>> plt.gray()
>>> plt.title('Output of spline edge filter')
>>> plt.show()

滤波

滤波是对输入信号进行某种修改的任何系统的通用名称。在 SciPy 中，信号可以被视为 NumPy 数组。有不同类型的滤波器用于不同类型的操作。滤波操作主要分为两种：线性和非线性。线性滤波器总是可以简化为将展平的 NumPy 数组乘以适当的矩阵，从而产生另一个展平的 NumPy 数组。当然，这通常不是计算滤波器的最佳方式，因为所涉及的矩阵和向量可能非常庞大。例如，使用此方法对 \(512 \times 512\) 图像进行滤波需要将 \(512^2 \times 512^2\) 矩阵乘以 \(512^2\) 向量。仅仅尝试使用标准 NumPy 数组存储 \(512^2 \times 512^2\) 矩阵就需要 \(68,719,476,736\) 个元素。每个元素 4 个字节，这将需要 \(256\textrm{GB}\) 的内存。在大多数应用程序中，此矩阵的大部分元素为零，并且采用不同的方法来计算滤波器的输出。

卷积/相关

许多线性滤波器也具有移位不变性。这意味着滤波操作在信号中的不同位置是相同的，并且这意味着可以仅从矩阵的一行（或一列）的知识来构造滤波矩阵。在这种情况下，可以使用傅里叶变换来完成矩阵乘法。

让 \(x\left[n\right]\) 定义一个由整数 \(n\) 索引的一维信号。两个一维信号的完全卷积可以表示为

\[y\left[n\right]=\sum_{k=-\infty}^{\infty}x\left[k\right]h\left[n-k\right].\]

只有当我们限制序列为可以在计算机中存储的有限支持序列时，才能直接实现此方程，选择 \(n=0\) 作为两个序列的起始点，让 \(K+1\) 作为 \(x\left[n\right]=0\) 对于所有 \(n\geq K+1\) 的值，让 \(M+1\) 作为 \(h\left[n\right]=0\) 对于所有 \(n\geq M+1\) 的值，则离散卷积表达式为

\[y\left[n\right]=\sum_{k=\max\left(n-M,0\right)}^{\min\left(n,K\right)}x\left[k\right]h\left[n-k\right].\]

为方便起见，假设 \(K\geq M.\) 那么，更明确地说，此操作的输出为：

\begin{eqnarray*} y\left[0\right] & = & x\left[0\right]h\left[0\right]\\ y\left[1\right] & = & x\left[0\right]h\left[1\right]+x\left[1\right]h\left[0\right]\\ y\left[2\right] & = & x\left[0\right]h\left[2\right]+x\left[1\right]h\left[1\right]+x\left[2\right]h\left[0\right]\\ \vdots & \vdots & \vdots\\ y\left[M\right] & = & x\left[0\right]h\left[M\right]+x\left[1\right]h\left[M-1\right]+\cdots+x\left[M\right]h\left[0\right]\\ y\left[M+1\right] & = & x\left[1\right]h\left[M\right]+x\left[2\right]h\left[M-1\right]+\cdots+x\left[M+1\right]h\left[0\right]\\ \vdots & \vdots & \vdots\\ y\left[K\right] & = & x\left[K-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[0\right]\\ y\left[K+1\right] & = & x\left[K+1-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[1\right]\\ \vdots & \vdots & \vdots\\ y\left[K+M-1\right] & = & x\left[K-1\right]h\left[M\right]+x\left[K\right]h\left[M-1\right]\\ y\left[K+M\right] & = & x\left[K\right]h\left[M\right].\end{eqnarray*}

因此，两个长度分别为 \(K+1\) 和 \(M+1\) 的有限序列的完整离散卷积，会产生一个长度为 \(K+M+1=\left(K+1\right)+\left(M+1\right)-1.\) 的有限序列。

一维卷积在 SciPy 中使用函数 convolve 实现。此函数将信号 \(x,\) \(h\) 以及两个可选标志 'mode' 和 'method' 作为输入，并返回信号 \(y.\)

第一个可选标志 'mode' 允许指定要返回的输出信号的哪个部分。默认值 'full' 返回整个信号。如果标志的值为 'same'，则仅返回中间的 \(K\) 个值，从 \(y\left[\left\lfloor \frac{M-1}{2}\right\rfloor \right]\) 开始，使得输出与第一个输入具有相同的长度。如果标志的值为 'valid'，则仅返回中间的 \(K-M+1=\left(K+1\right)-\left(M+1\right)+1\) 个输出值，其中 \(z\) 取决于从 \(h\left[0\right]\) 到 \(h\left[M\right]\) 的最小输入的所有值。换句话说，仅返回包括 \(y\left[M\right]\) 到 \(y\left[K\right]\) 的值。

第二个可选标志 'method' 确定如何计算卷积，可以通过使用 fftconvolve 的傅里叶变换方法，或者通过直接方法。默认情况下，它选择预期更快的计算方法。傅里叶变换方法的阶数为 \(O(N\log N)\)，而直接方法的阶数为 \(O(N^2)\)。根据大 O 常数和 \(N\) 的值，这两种方法中的一种可能会更快。默认值 'auto' 执行粗略的计算并选择预期更快的方法，而值 'direct' 和 'fft' 强制使用其他两种方法进行计算。

下面的代码显示了两个序列卷积的简单示例。

>>> x = np.array([1.0, 2.0, 3.0])
>>> h = np.array([0.0, 1.0, 0.0, 0.0, 0.0])
>>> signal.convolve(x, h)
array([ 0.,  1.,  2.,  3.,  0.,  0.,  0.])
>>> signal.convolve(x, h, 'same')
array([ 2.,  3.,  0.])

同一个函数 convolve 实际上可以接受 N 维数组作为输入，并将返回两个数组的 N 维卷积，如下面的代码示例所示。同样，在这种情况下也可以使用相同的输入标志。

>>> x = np.array([[1., 1., 0., 0.], [1., 1., 0., 0.], [0., 0., 0., 0.], [0., 0., 0., 0.]])
>>> h = np.array([[1., 0., 0., 0.], [0., 0., 0., 0.], [0., 0., 1., 0.], [0., 0., 0., 0.]])
>>> signal.convolve(x, h)
array([[ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])

相关性与卷积非常相似，只是减号变成了加号。因此，

\[w\left[n\right]=\sum_{k=-\infty}^{\infty}y\left[k\right]x\left[n+k\right],\]

是信号 \(y\) 和 \(x.\) 的（交叉）相关性。对于有限长度的信号，其中 \(y\left[n\right]=0\) 在 \(\left[0,K\right]\) 范围外，\(x\left[n\right]=0\) 在 \(\left[0,M\right],\) 范围外，求和可以简化为

\[w\left[n\right]=\sum_{k=\max\left(0,-n\right)}^{\min\left(K,M-n\right)}y\left[k\right]x\left[n+k\right].\]

再次假设 \(K\geq M\)，则得到

\begin{eqnarray*} w\left[-K\right] & = & y\left[K\right]x\left[0\right]\\ w\left[-K+1\right] & = & y\left[K-1\right]x\left[0\right]+y\left[K\right]x\left[1\right]\\ \vdots & \vdots & \vdots\\ w\left[M-K\right] & = & y\left[K-M\right]x\left[0\right]+y\left[K-M+1\right]x\left[1\right]+\cdots+y\left[K\right]x\left[M\right]\\ w\left[M-K+1\right] & = & y\left[K-M-1\right]x\left[0\right]+\cdots+y\left[K-1\right]x\left[M\right]\\ \vdots & \vdots & \vdots\\ w\left[-1\right] & = & y\left[1\right]x\left[0\right]+y\left[2\right]x\left[1\right]+\cdots+y\left[M+1\right]x\left[M\right]\\ w\left[0\right] & = & y\left[0\right]x\left[0\right]+y\left[1\right]x\left[1\right]+\cdots+y\left[M\right]x\left[M\right]\\ w\left[1\right] & = & y\left[0\right]x\left[1\right]+y\left[1\right]x\left[2\right]+\cdots+y\left[M-1\right]x\left[M\right]\\ w\left[2\right] & = & y\left[0\right]x\left[2\right]+y\left[1\right]x\left[3\right]+\cdots+y\left[M-2\right]x\left[M\right]\\ \vdots & \vdots & \vdots\\ w\left[M-1\right] & = & y\left[0\right]x\left[M-1\right]+y\left[1\right]x\left[M\right]\\ w\left[M\right] & = & y\left[0\right]x\left[M\right].\end{eqnarray*}

SciPy 函数 correlate 实现此操作。此操作可以使用等效的标志来返回完整的 \(K+M+1\) 长度序列 ('full')，或返回与最大序列大小相同的序列，从 \(w\left[-K+\left\lfloor \frac{M-1}{2}\right\rfloor \right]\) ('same') 开始，或者返回一个序列，其中值取决于最小序列 ('valid') 的所有值。最后一个选项返回 \(K-M+1\) 个值，包括 \(w\left[M-K\right]\) 到 \(w\left[0\right]\)。

函数 correlate 还可以将任意 N 维数组作为输入，并在输出时返回两个数组的 N 维卷积。

当 \(N=2,\) 时，correlate 和/或 convolve 可用于构建任意图像滤波器，以执行诸如模糊、增强和图像边缘检测等操作。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = datasets.face(gray=True)
>>> w = np.zeros((50, 50))
>>> w[0][0] = 1.0
>>> w[49][25] = 1.0
>>> image_new = signal.fftconvolve(image, w)

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

当其中一个信号远小于另一个信号时（\(K\gg M\)），主要使用如上所述的时域卷积来执行滤波。否则，使用函数 fftconvolve 提供的频域线性滤波效率更高。默认情况下，convolve 使用 choose_conv_method 估计最快的方法。

如果滤波器函数 \(w[n,m]\) 可以根据以下方式分解：

\[h[n, m] = h_1[n] h_2[m],\]

则可以通过函数 sepfir2d 计算卷积。例如，我们考虑一个高斯滤波器 gaussian

\[h[n, m] \propto e^{-x^2-y^2} = e^{-x^2} e^{-y^2},\]

该滤波器通常用于模糊处理。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt

>>> image = np.asarray(datasets.ascent(), np.float64)
>>> w = signal.windows.gaussian(51, 10.0)
>>> image_new = signal.sepfir2d(image, w, w)

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

差分方程滤波

线性一维滤波器（包括卷积滤波器）的一般类别是由差分方程描述的滤波器：

\[\sum_{k=0}^{N}a_{k}y\left[n-k\right]=\sum_{k=0}^{M}b_{k}x\left[n-k\right],\]

其中 \(x\left[n\right]\) 是输入序列，\(y\left[n\right]\) 是输出序列。如果我们假设初始静止状态，使得当 \(n<0\) 时 \(y\left[n\right]=0\)，则可以使用卷积来实现这种滤波器。但是，如果 \(k\geq1.\) 时 \(a_{k}\neq0\)，则卷积滤波器序列 \(h\left[n\right]\) 可能是无限的。此外，这种通用类别的线性滤波器允许对 \(n<0\) 的 \(y\left[n\right]\) 施加初始条件，从而导致无法使用卷积表示的滤波器。

差分方程滤波器可以被认为是根据其先前的值递归地寻找 \(y\left[n\right]\)：

\[a_{0}y\left[n\right]=-a_{1}y\left[n-1\right]-\cdots-a_{N}y\left[n-N\right]+\cdots+b_{0}x\left[n\right]+\cdots+b_{M}x\left[n-M\right].\]

通常，选择 \(a_{0}=1\) 进行归一化。在 SciPy 中实现此通用差分方程滤波器比前面的方程式所暗示的要复杂一些。它的实现方式是只需要延迟一个信号。实际的实现方程为（假设 \(a_{0}=1\) ）

\begin{eqnarray*} y\left[n\right] & = & b_{0}x\left[n\right]+z_{0}\left[n-1\right]\\ z_{0}\left[n\right] & = & b_{1}x\left[n\right]+z_{1}\left[n-1\right]-a_{1}y\left[n\right]\\ z_{1}\left[n\right] & = & b_{2}x\left[n\right]+z_{2}\left[n-1\right]-a_{2}y\left[n\right]\\ \vdots & \vdots & \vdots\\ z_{K-2}\left[n\right] & = & b_{K-1}x\left[n\right]+z_{K-1}\left[n-1\right]-a_{K-1}y\left[n\right]\\ z_{K-1}\left[n\right] & = & b_{K}x\left[n\right]-a_{K}y\left[n\right],\end{eqnarray*}

其中 \(K=\max\left(N,M\right).\) 请注意，如果 \(K>M\)，则 \(b_{K}=0\)，如果 \(K>N.\)，则 \(a_{K}=0\)。这样，时间 \(n\) 的输出仅取决于时间 \(n\) 的输入和前一时间的值 \(z_{0}\) 。只要在每个时间步计算并存储 \(K\) 个值 \(z_{0}\left[n-1\right]\ldots z_{K-1}\left[n-1\right]\)，就可以始终计算出此值。

差分方程滤波器通过 SciPy 中的 lfilter 命令调用。该命令的输入包括向量 \(b,\) 向量 \(a,\) 一个信号 \(x\)，并返回向量 \(y\)（与 \(x\) 长度相同），该向量使用上述方程计算得出。如果 \(x\) 是 N 维的，则沿提供的轴计算滤波器。如果需要，可以提供初始条件，即 \(z_{0}\left[-1\right]\) 到 \(z_{K-1}\left[-1\right]\) 的值，否则将假定它们都为零。如果提供了初始条件，则还会返回中间变量的最终条件。例如，这些条件可用于以相同的状态重新开始计算。

有时，用信号 \(x\left[n\right]\) 和 \(y\left[n\right]\) 来表示初始条件更方便。换句话说，也许您拥有 \(x\left[-M\right]\) 到 \(x\left[-1\right]\) 的值以及 \(y\left[-N\right]\) 到 \(y\left[-1\right]\) 的值，并且想确定应该将哪些 \(z_{m}\left[-1\right]\) 的值作为初始条件提供给差分方程滤波器。不难证明，对于 \(0\leq m<K,\)

\[z_{m}\left[n\right]=\sum_{p=0}^{K-m-1}\left(b_{m+p+1}x\left[n-p\right]-a_{m+p+1}y\left[n-p\right]\right).\]

使用此公式，我们可以在给定 \(y\)（和 \(x\)）的初始条件的情况下，找到初始条件向量 \(z_{0}\left[-1\right]\) 到 \(z_{K-1}\left[-1\right]\)。 lfiltic 命令执行此功能。

例如，考虑以下系统

\[y[n] = \frac{1}{2} x[n] + \frac{1}{4} x[n-1] + \frac{1}{3} y[n-1]\]

以下代码计算给定信号 \(x[n]\) 的信号 \(y[n]\)；首先，对于初始条件 \(y[-1] = 0\)（默认情况），然后对于 \(y[-1] = 2\)，通过 lfiltic 实现。

>>> import numpy as np
>>> from scipy import signal

>>> x = np.array([1., 0., 0., 0.])
>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.lfilter(b, a, x)
array([0.5, 0.41666667, 0.13888889, 0.0462963])
>>> zi = signal.lfiltic(b, a, y=[2.])
>>> signal.lfilter(b, a, x, zi=zi)
(array([ 1.16666667,  0.63888889,  0.21296296,  0.07098765]), array([0.02366]))

请注意，输出信号 \(y[n]\) 的长度与输入信号 \(x[n]\) 的长度相同。

线性系统分析

由线性差分方程描述的线性系统可以完全由系数向量 \(a\) 和 \(b\) 描述，如上所述；另一种表示方法是提供一个因子 \(k\)、\(N_z\) 个零点 \(z_k\) 和 \(N_p\) 个极点 \(p_k\)，分别通过传递函数 \(H(z)\) 来描述系统，根据

\[H(z) = k \frac{ (z-z_1)(z-z_2)...(z-z_{N_z})}{ (z-p_1)(z-p_2)...(z-p_{N_p})}.\]

这种替代表示可以使用 scipy 函数 tf2zpk 获得；其逆运算由 zpk2tf 提供。

对于上面的示例，我们有

>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.tf2zpk(b, a)
(array([-0.5]), array([ 0.33333333]), 0.5)

即，该系统在 \(z=-1/2\) 处有一个零点，在 \(z=1/3\) 处有一个极点。

scipy 函数 freqz 允许计算由系数 \(a_k\) 和 \(b_k\) 描述的系统的频率响应。 有关综合示例，请参阅 freqz 函数的帮助。

滤波器设计

时域离散滤波器可以分为有限响应 (FIR) 滤波器和无限响应 (IIR) 滤波器。 FIR 滤波器可以提供线性相位响应，而 IIR 滤波器则不能。 SciPy 提供了用于设计这两种类型滤波器的函数。

FIR 滤波器

firwin 函数根据窗口方法设计滤波器。根据提供的参数，该函数返回不同的滤波器类型（例如，低通、带通...）。

下面的示例分别设计了一个低通滤波器和一个带阻滤波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b1 = signal.firwin(40, 0.5)
>>> b2 = signal.firwin(41, [0.3, 0.8])
>>> w1, h1 = signal.freqz(b1)
>>> w2, h2 = signal.freqz(b2)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w1, 20*np.log10(np.abs(h1)), 'b')
>>> plt.plot(w2, 20*np.log10(np.abs(h2)), 'r')
>>> plt.ylabel('Amplitude Response (dB)')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

请注意，默认情况下，firwin 使用归一化频率，使得值 \(1\) 对应于奈奎斯特频率，而函数 freqz 的定义是使值 \(\pi\) 对应于奈奎斯特频率。

firwin2 允许通过分别指定拐角频率数组和相应的增益来设计几乎任意的频率响应。

下面的示例设计一个具有这种任意幅度响应的滤波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b = signal.firwin2(150, [0.0, 0.3, 0.6, 1.0], [1.0, 2.0, 0.5, 0.0])
>>> w, h = signal.freqz(b)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, np.abs(h))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

请注意，y 轴的线性缩放以及 firwin2 和 freqz 中奈奎斯特频率的不同定义（如上所述）。

IIR 滤波器

SciPy 提供了两个函数来直接设计 IIR： iirdesign 和 iirfilter，其中滤波器类型（例如，椭圆）作为参数传递，以及用于特定滤波器类型的更多滤波器设计函数，例如 ellip。

下面的示例设计了一个具有已定义的通带和阻带纹波的椭圆低通滤波器。请注意，与上面示例中的 FIR 滤波器相比，为了达到相同的 \(\approx 60\) dB 的阻带衰减，滤波器阶数要低得多（阶数为 4）。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b, a = signal.iirfilter(4, Wn=0.2, rp=5, rs=60, btype='lowpass', ftype='ellip')
>>> w, h = signal.freqz(b, a)

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

注意

重要的是要注意，firwin 和 iirfilter 的截止频率的定义不同。对于 firwin，截止频率位于半幅度处（即 -6dB）。对于 iirfilter，截止频率位于半功率处（即 -3dB）。

>>> import numpy as np
>>> from matplotlib import pyplot as plt
>>> from scipy import signal as sig

>>> fs = 16000
>>> b = sig.firwin(101, 2500, fs=fs)
>>> f, h_fft = sig.freqz(b, fs=fs)
>>> h_amp = 20 * np.log10(np.abs(h_fft))
>>> _, ax = plt.subplots(layout="constrained")
>>> ax.plot(f, h_amp, label="FIR")
>>> ax.grid(True)

>>> b, a = sig.iirfilter(15, 2500, btype="low", fs=fs)
>>> f, h_fft = sig.freqz(b, a, fs=fs)
>>> h_amp = 20 * np.log10(np.abs(h_fft))
>>> ax.plot(f, h_amp, label="IIR")
>>> ax.set(xlim=[2100, 2900], ylim=[-10, 2])
>>> ax.set(xlabel="Frequency (Hz)", ylabel="Amplitude Response [dB]")
>>> ax.legend()

滤波器系数

滤波器系数可以以几种不同的格式存储

• ‘ba’ 或 ‘tf’ = 传递函数系数

• ‘zpk’ = 零点、极点和整体增益

• ‘ss’ = 状态空间系统表示

• ‘sos’ = 二阶部分的传递函数系数

诸如 tf2zpk 和 zpk2ss 之类的函数可以在它们之间进行转换。

传递函数表示

ba 或 tf 格式是一个 2 元组 (b, a)，表示一个传递函数，其中 b 是长度为 M+1 的数组，包含 M 阶分子多项式的系数，a 是长度为 N+1 的数组，包含 N 阶分母的系数，为传递函数变量的正降幂。因此，\(b = [b_0, b_1, ..., b_M]\) 和 \(a =[a_0, a_1, ..., a_N]\) 的元组可以表示以下形式的模拟滤波器

\[H(s) = \frac {b_0 s^M + b_1 s^{(M-1)} + \cdots + b_M} {a_0 s^N + a_1 s^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i s^{(M-i)}} {\sum_{i=0}^N a_i s^{(N-i)}}\]

或以下形式的离散时间滤波器

\[H(z) = \frac {b_0 z^M + b_1 z^{(M-1)} + \cdots + b_M} {a_0 z^N + a_1 z^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i z^{(M-i)}} {\sum_{i=0}^N a_i z^{(N-i)}}.\]

这种“正幂”形式在控制工程中更常见。如果 M 和 N 相等（对于通过双线性变换生成的所有滤波器来说都是如此），那么这恰好等效于 DSP 中首选的“负幂”离散时间形式

\[H(z) = \frac {b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}} {a_0 + a_1 z^{-1} + \cdots + a_N z^{-N}} = \frac {\sum_{i=0}^M b_i z^{-i}} {\sum_{i=0}^N a_i z^{-i}}.\]

虽然这对于常见的滤波器是正确的，但请记住，这在一般情况下并不成立。如果 M 和 N 不相等，则在找到极点和零点之前，必须先将离散时间传递函数系数转换为“正幂”形式。

这种表示在高阶时会出现数值误差，因此在可能的情况下，最好使用其他格式。

零点和极点表示

zpk 格式是一个三元组 (z, p, k)，其中 z 是传递函数的复零点的 M 长度数组 \(z = [z_0, z_1, ..., z_{M-1}]\)，p 是传递函数的复极点的 N 长度数组 \(p = [p_0, p_1, ..., p_{N-1}]\)，而 k 是一个标量增益。它们表示数字传递函数

\[H(z) = k \cdot \frac {(z - z_0) (z - z_1) \cdots (z - z_{(M-1)})} {(z - p_0) (z - p_1) \cdots (z - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (z - z_i)} {\prod_{i=0}^{N-1} (z - p_i)}\]

或模拟传递函数

\[H(s) = k \cdot \frac {(s - z_0) (s - z_1) \cdots (s - z_{(M-1)})} {(s - p_0) (s - p_1) \cdots (s - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (s - z_i)} {\prod_{i=0}^{N-1} (s - p_i)}.\]

尽管根的集合存储为有序的 NumPy 数组，但它们的顺序并不重要：([-1, -2], [-3, -4], 1) 与 ([-2, -1], [-4, -3], 1) 是相同的滤波器。

状态空间系统表示

ss 格式是数组 (A, B, C, D) 的 4 元组，表示形式为 N 阶数字/离散时间系统的状态空间

\[\begin{split}\mathbf{x}[k+1] = A \mathbf{x}[k] + B \mathbf{u}[k]\\ \mathbf{y}[k] = C \mathbf{x}[k] + D \mathbf{u}[k]\end{split}\]

或形式为连续/模拟系统的状态空间

\[\begin{split}\dot{\mathbf{x}}(t) = A \mathbf{x}(t) + B \mathbf{u}(t)\\ \mathbf{y}(t) = C \mathbf{x}(t) + D \mathbf{u}(t),\end{split}\]

其中有 P 个输入，Q 个输出和 N 个状态变量，其中

• x 是状态向量

• y 是长度为 Q 的输出向量

• u 是长度为 P 的输入向量

• A 是状态矩阵，形状为 (N, N)

• B 是输入矩阵，形状为 (N, P)

• C 是输出矩阵，形状为 (Q, N)

• D 是前馈或直通矩阵，形状为 (Q, P)。（在系统没有直接直通的情况下，D 中的所有值均为零。）

状态空间是最通用的表示形式，也是唯一允许使用多输入多输出 (MIMO) 系统的表示形式。对于给定的传递函数，存在多种状态空间表示形式。具体而言，“可控规范型”和“可观测规范型”具有与 tf 表示形式相同的系数，因此，会遭受相同的数值误差。

二阶节表示

sos 格式是一个形状为 (n_sections, 6) 的二维数组，表示一系列二阶传递函数，当串联在一起时，可以实现具有最小数值误差的高阶滤波器。每行对应一个二阶 tf 表示形式，其中前三列提供分子系数，后三列提供分母系数

\[[b_0, b_1, b_2, a_0, a_1, a_2]\]

这些系数通常被归一化，使得 \(a_0\) 始终为 1。使用浮点数计算时，节的顺序通常并不重要；无论顺序如何，滤波器的输出都将相同。

滤波器变换

IIR 滤波器设计函数首先生成一个原型模拟低通滤波器，其归一化截止频率为 1 rad/sec。然后使用以下替换将其转换为其他频率和频带类型

类型	变换
lp2lp	\(s \rightarrow \frac{s}{\omega_0}\)
lp2hp	\(s \rightarrow \frac{\omega_0}{s}\)
lp2bp	\(s \rightarrow \frac{s^2 + {\omega_0}^2}{s \cdot \mathrm{BW}}\)
lp2bs	\(s \rightarrow \frac{s \cdot \mathrm{BW}}{s^2 + {\omega_0}^2}\)

此处，\(\omega_0\) 是新的截止或中心频率，而 \(\mathrm{BW}\) 是带宽。这些在对数频率轴上保持对称性。

要将变换后的模拟滤波器转换为数字滤波器，将使用 bilinear 变换，该变换进行以下替换

\[s \rightarrow \frac{2}{T} \frac{z - 1}{z + 1},\]

其中 T 是采样时间（采样频率的倒数）。

其他滤波器

信号处理包还提供了许多其他滤波器。

中值滤波器

当噪声明显是非高斯噪声或希望保留边缘时，通常会应用中值滤波器。中值滤波器的工作原理是对围绕感兴趣点的矩形区域中的所有数组像素值进行排序。此邻域像素值列表的样本中值用作输出数组的值。样本中值是邻域值的排序列表中间的数组值。如果邻域中存在偶数个元素，则中间两个值的平均值用作中值。可以在 N 维数组上工作的通用中值滤波器是 medfilt。仅适用于二维数组的专用版本可作为 medfilt2d 使用。

阶数滤波器

中值滤波器是称为阶数滤波器的更通用的一类滤波器的具体示例。为了计算特定像素处的输出，所有阶数滤波器都使用围绕该像素的区域中的数组值。这些数组值被排序，然后选择其中一个作为输出值。对于中值滤波器，数组值列表的样本中值用作输出。通用阶数滤波器允许用户选择排序值中的哪个值将用作输出。因此，例如，可以选择列表中的最大值或最小值。除了输入数组和区域掩码之外，阶数滤波器还采用一个额外的参数，用于指定排序的邻域数组值列表中应使用哪个元素作为输出。执行阶数滤波器的命令是 order_filter。

维纳滤波器

维纳滤波器是用于图像去噪的简单去模糊滤波器。这不是图像重建问题中通常描述的维纳滤波器，而是一个简单的局部均值滤波器。令 \(x\) 为输入信号，则输出为

\[\begin{split}y=\left\{ \begin{array}{cc} \frac{\sigma^{2}}{\sigma_{x}^{2}}m_{x}+\left(1-\frac{\sigma^{2}}{\sigma_{x}^{2}}\right)x & \sigma_{x}^{2}\geq\sigma^{2},\\ m_{x} & \sigma_{x}^{2}<\sigma^{2},\end{array}\right.\end{split}\]

其中 \(m_{x}\) 是均值的局部估计，\(\sigma_{x}^{2}\) 是方差的局部估计。这些估计的窗口是一个可选的输入参数（默认为 \(3\times3\)）。参数 \(\sigma^{2}\) 是阈值噪声参数。如果未给出 \(\sigma\)，则将其估计为局部方差的平均值。

希尔伯特滤波器

希尔伯特变换从实信号构建复值解析信号。例如，如果 \(x=\cos\omega n\)，那么 \(y=\textrm{hilbert}\left(x\right)\) 将返回（除了边缘附近） \(y=\exp\left(j\omega n\right).\) 在频域中，希尔伯特变换执行

\[Y=X\cdot H,\]

其中 \(H\) 对于正频率为 \(2\)，对于负频率为 \(0\)，对于零频率为 \(1\)。

模拟滤波器设计

函数 iirdesign，iirfilter，以及特定滤波器类型（例如，ellip）的滤波器设计函数都具有一个标志 analog，允许设计模拟滤波器。

下面的示例设计一个模拟（IIR）滤波器，通过 tf2zpk 获取极点和零点，并在复数 s 平面中绘制它们。在幅度响应中可以清楚地看到 \(\omega \approx 150\) 和 \(\omega \approx 300\) 处的零点。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt

>>> b, a = signal.iirdesign(wp=100, ws=200, gpass=2.0, gstop=40., analog=True)
>>> w, h = signal.freqs(b, a)

>>> plt.title('Analog filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency')
>>> plt.grid()
>>> plt.show()

>>> z, p, k = signal.tf2zpk(b, a)

>>> plt.plot(np.real(z), np.imag(z), 'ob', markerfacecolor='none')
>>> plt.plot(np.real(p), np.imag(p), 'xr')
>>> plt.legend(['Zeros', 'Poles'], loc=2)

>>> plt.title('Pole / Zero Plot')
>>> plt.xlabel('Real')
>>> plt.ylabel('Imaginary')
>>> plt.grid()
>>> plt.show()

\[% LaTeX Macros to make the LaTeX formulas more readable: \newcommand{\IC}{{\mathbb{C}}} % set of complex numbers \newcommand{\IN}{{\mathbb{N}}} % set of natural numbers \newcommand{\IR}{{\mathbb{R}}} % set of real numbers \newcommand{\IZ}{{\mathbb{Z}}} % set of integers \newcommand{\jj}{{\mathbb{j}}} % imaginary unit \newcommand{\e}{\operatorname{e}} % Euler's number \newcommand{\dd}{\operatorname{d}} % infinitesimal operator \newcommand{\abs}[1]{\left|#1\right|} % absolute value \newcommand{\conj}[1]{\overline{#1}} % complex conjugate \newcommand{\conjT}[1]{\overline{#1^T}} % transposed complex conjugate \newcommand{\inv}[1]{\left(#1\right)^{\!-1}} % inverse % Since the physics package is not loaded, we define the macros ourselves: \newcommand{\vb}[1]{\mathbf{#1}} % vectors and matrices are bold % new macros: \newcommand{\rect}{\operatorname{rect}} % rect or boxcar function \newcommand{\sinc}{\operatorname{sinc}} % sinc(t) := sin(pi*t) / (pi*t)\]

频谱分析

频谱分析是指研究信号的傅里叶变换 [1]。根据上下文，傅里叶变换的各种频谱表示存在诸如频谱、谱密度或周期图等各种名称。[2] 本节通过固定持续时间的连续时间正弦波信号的示例来说明最常见的表示。然后讨论了对该正弦波的采样版本使用离散傅里叶变换 [3]。

单独的小节专门介绍频谱的相位，估计没有 (periodogram) 和有平均 (welch) 的功率谱密度，以及非等间隔信号 (lombscargle)。

请注意，傅里叶级数的概念密切相关，但在一个关键点上有所不同：傅里叶级数具有由离散频率谐波组成的频谱，而本节中的频谱在频率上是连续的。

连续时间正弦信号

考虑一个具有幅度 \(a\)、频率 \(f_x\) 和持续时间 \(\tau\) 的正弦信号，即

(1)\[x(t) = a \sin(2 \pi f_x t)\, \rect(\frac{t}{\tau}-\frac{1}{2}) = \left(\frac{a}{2\jj} \e^{\jj 2 \pi f_x t} - \frac{a}{2\jj} \e^{-\jj 2 \pi f_x t} \right) \rect(\frac{t}{\tau}-\frac{1}{2})\ .\]

由于 \(\rect(t)\) 函数在 \(|t|<1/2\) 时为 1，在 \(|t|>1/2\) 时为 0，因此它将 \(x(t)\) 限制在区间 \([0, \tau]\)。用复指数表示正弦波显示了其频率为 \(\pm f_x\) 的两个周期性分量。我们假设 \(x(t)\) 是一个电压信号，因此其单位为 \(\text{V}\)。

在信号处理中，绝对值的平方 \(|x(t)|^2\) 的积分用于定义信号的能量和功率，即

(2)\[E_x := \int_0^\tau |x(t)|^2 \dd t\ = \frac{1}{2}|a|^2\tau\ , \qquad P_x := \frac{1}{\tau}E_x = \frac{1}{2}|a|^2\ .\]

功率 \(P_x\) 可以解释为单位时间间隔内的能量 \(E_x\)。按单位计算，对 \(t\) 积分会导致乘以秒。因此，\(E_x\) 的单位为 \(\text{V}^2\text{s}\)，\(P_x\) 的单位为 \(\text{V}^2\)。

对 \(x(t)\) 应用傅里叶变换，即

(3)\[X(f) = \int_{\IR} x(t)\, \e^{-2\jj\pi f t}\, \dd t = \frac{a \tau}{2\jj} \Big( \sinc\!\big(\tau (f-f_x)\big) - \sinc\!\big(\tau (f+f_x)\big) \Big)\e^{-\jj\pi\tau f}\ ,\]

得到两个以 \(\pm f_x\) 为中心的 \(\sinc(f) := \sin(\pi f) /(\pi f)\) 函数。幅度（绝对值）\(|X(f)|\) 在 \(\pm f_x\) 处有两个最大值，其值为 \(|a|\tau/2\)。可以在下面的图中看到，\(X(f)\) 并不集中在 \(\pm f_x\) 处的主瓣周围，而是包含高度与 \(1/(\tau f)\) 成比例减小的旁瓣。这种所谓的“频谱泄漏”[4] 是由将正弦波限制在有限区间内引起的。请注意，信号持续时间 \(\tau\) 越短，泄漏越大。为了独立于信号持续时间，可以使用所谓的“幅度谱”\(X(f)/\tau\) 来代替频谱 \(X(f)\)。它在 \(f\) 处的值对应于复指数 \(\exp(\jj2\pi f t)\) 的幅度。

根据帕塞瓦尔定理，也可以通过其傅里叶变换 \(X(f)\) 计算能量，即

(4)\[ E_X := \int_\IR \abs{X(f)}^2 \dd f = E_x\]

。例如，可以直接计算得出等式 (4) 中 \(X(f)\) 的能量为 \(|a|^2\tau/2\)。因此，可以使用以下公式确定频带 \([f_a, f_b]\) 中的信号功率

(5)\[P_X^{a,b} = \frac{1}{\tau} \int_{f_a}^{f_b} \abs{X(f)}^2 \dd f\ .\]

因此，函数 \(|X(f)|^2\) 可以定义为所谓的“能量谱密度”，\(S_{xx}(f) := |X(f)|^2 / \tau\) 定义为 \(x(t)\) 的“功率谱密度”（PSD）。除了 PSD 之外，还使用所谓的“幅度谱密度” \(X(f) / \sqrt{\tau}\)，它仍然包含相位信息。它的绝对平方是 PSD，因此它与信号的均方根 (RMS) 值 \(\sqrt{P_x}\) 的概念密切相关。

总而言之，本小节介绍了五种表示频谱的方法

正弦信号 \(x(t)\) 的频谱表示比较，其中 \(x(t)\) 的单位为 \(\text{V}\)，如公式 (1) 所示：

	频谱	幅度谱	能量谱密度	功率谱密度 (PSD)	幅度谱密度
定义	\(X(f)\)	\(X(f) / \tau\)	\(|X(f)|^2\)	\(|X(f)|^2 / \tau\)	\(X(f) / \sqrt{\tau}\)
在 \(\pm f_x\) 处的幅度	\(\frac{1}{2}|a|\tau\)	\(\frac{1}{2}|a|\)	\(\frac{1}{4}|a|^2\tau^2\)	\(\frac{1}{4}|a|^2\tau\)	\(\frac{1}{2}|a|\sqrt{\tau}\)
单位	\(\text{V} / \text{Hz}\)	\(\text{V}\)	\(\text{V}^2\text{s} / \text{Hz}\)	\(\text{V}^2 / \text{Hz}\)	\(\text{V} / \sqrt{\text{Hz}}\)

请注意，上表中给出的单位不是明确的，例如，\(\text{V}^2\text{s} / \text{Hz} = \text{V}^2\text{s}^2 = \text{V}^2/ \text{Hz}^2\)。当使用幅度谱的绝对值 \(|X(f) / \tau|\) 时，它被称为幅度谱。此外，请注意，表示的命名方案不一致，并且在文献中有所不同。

对于实值信号，通常使用所谓的“单边”频谱表示。它只使用非负频率（因为如果 \(x(t)\in\IR\)，则 \(X(-f)= \conj{X}(f)\)）。有时，负频率的值会加到它们对应的正频率值上。这样，幅度谱就可以读取 \(x(t)\) 在 \(f_x\) 处的完整（而不是一半）幅度正弦波，并且 PSD 中一个区间的面积表示其完整（而不是一半）的功率。请注意，对于幅度谱密度，正值不会加倍，而是乘以 \(\sqrt{2}\)，因为它是 PSD 的平方根。此外，没有规范的方法来命名加倍的频谱。

下图显示了四个正弦信号 \(x(t)\)（如公式 (1) 所示）的三种不同频谱表示，这些信号具有不同的振幅 \(a\) 和持续时间 \(\tau\)。为了减少混乱，频谱以 \(f_x\) 为中心，并彼此相邻绘制。

import matplotlib.pyplot as plt
import numpy as np

aa = [1, 1, 2, 2]  # amplitudes
taus = [1, 2, 1, 2]  # durations

fg0, axx = plt.subplots(3, 1, sharex='all', tight_layout=True, figsize=(6., 4.))
axx[0].set(title=r"Spectrum $|X(f)|$", ylabel="V/Hz")
axx[1].set(title=r"Magnitude Spectrum $|X(f)/\tau|$ ", ylabel=r"V")
axx[2].set(title=r"Amplitude Spectral Density $|X(f)/\sqrt{\tau}|$",
           ylabel=r"$\operatorname{V} / \sqrt{\operatorname{Hz}}$",
           xlabel="Frequency $f$ in Hertz",)

x_labels, x_ticks = [], []
f = np.linspace(-2.5, 2.5, 400)
for c_, (a_, tau_) in enumerate(zip(aa, taus), start=1):
    aZ_, f_ = abs(a_ * tau_ * np.sinc(tau_ * f) / 2), f + c_ * 5
    axx[0].plot(f_, aZ_)
    axx[1].plot(f_, aZ_ / tau_)
    axx[2].plot(f_, aZ_ / np.sqrt(tau_))
    x_labels.append(rf"$a={a_:g}$, $\tau={tau_:g}$")
    x_ticks.append(c_ * 5)

axx[2].set_xticks(x_ticks)
axx[2].set_xticklabels(x_labels)
plt.show()

请注意，根据表示形式的不同，峰值的高度会有所不同。只有幅度谱的解释是直接的：第二个图中 \(f_x\) 处的峰值表示正弦信号幅度 \(|a|\) 的一半。对于所有其他表示形式，需要考虑持续时间 \(\tau\) 才能提取有关信号振幅的信息。

采样正弦信号

在实践中，采样信号被广泛使用。即，信号由 \(n\) 个样本 \(x_k := x(kT)\) 表示，其中 \(k=0, \ldots, n-1\)，\(T\) 是采样间隔，\(\tau:=nT\) 是信号的持续时间，\(f_S := 1/T\) 是采样频率。请注意，连续信号需要被限制在 \([-f_S/2, f_S/2]\) 以避免混叠，其中 \(f_S/2\) 被称为奈奎斯特频率。 [5] 将积分替换为求和来计算信号的能量和功率，即：

\[E_x = T\sum_{k=0}^{n-1} \abs{x_k}^2 = \frac{1}{2}|a|^2\tau\ , \qquad P_x = \frac{1}{\tau}E_x = \frac{1}{2}|a|^2\ ,\]

与公式 (2) 中的连续时间情况产生相同的结果。离散傅里叶变换 (DFT) 及其逆变换（使用 scipy.fft 模块中的高效 FFT 计算实现）由下式给出：

\[X_l := \sum_{k=0}^{n-1} x_k \e^{-2\jj\pi k l / n}\ ,\qquad x_k = \frac{1}{n} \sum_{l=0}^{n-1} X_l \e^{2\jj\pi k l / n}\ .\]

DFT 可以解释为公式 (3) 中连续傅里叶变换的未缩放采样版本，即：

\[X(l\Delta f) = T X_l\ , \quad \Delta f := 1/\tau = 1/(nT)\ .\]

下图显示了两个具有单位振幅和频率分别为 20 Hz 和 20.5 Hz 的正弦信号的幅度谱。该信号由 \(n=100\) 个样本组成，采样间隔为 \(T=10\) ms，导致持续时间为 \(\tau=1\) s，采样频率为 \(f_S=100\) Hz。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq

n, T = 100, 0.01  # number of samples and sampling interval
fcc = (20, 20.5)  # frequencies of sines
t = np.arange(n) * T
xx = (np.sin(2 * np.pi * fx_ * t) for fx_ in fcc)  # sine signals

f = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_) / n for x_ in xx)  # one-sided magnitude spectrum

fg1, ax1 = plt.subplots(1, 1, tight_layout=True, figsize=(6., 3.))
ax1.set(title=r"Magnitude Spectrum (no window) of $x(t) = \sin(2\pi f_x t)$ ",
        xlabel=rf"Frequency $f$ in Hertz (bin width $\Delta f = {f[1]}\,$Hz)",
        ylabel=r"Magnitude $|X(f)|/\tau$", xlim=(f[0], f[-1]))
for X_, fc_, m_ in zip(XX, fcc, ('x-', '.-')):
    ax1.plot(f, abs(X_), m_, label=rf"$f_x={fc_}\,$Hz")

ax1.grid(True)
ax1.legend()
plt.show()

对 20 Hz 信号的解释似乎很简单：除了 20 Hz 处为 0.5 外，所有值都为零，这与公式 (1) 一致，对应于输入信号幅度的一半。另一方面，20.5 Hz 信号的峰值沿频率轴分散。公式 (3) 表明，这种差异的原因是 20 Hz 是 1 Hz 的频带宽度的倍数，而 20.5 Hz 不是。下图通过在采样频谱上叠加连续频谱来阐明这一点。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq

n, T = 100, 0.01  # number of samples and sampling interval
tau = n*T
t = np.arange(n) * T
fcc = (20, 20.5)  # frequencies of sines
xx = (np.sin(2 * np.pi * fc_ * t) for fc_ in fcc)  # sine signals

f = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_) / n for x_ in xx)  # one-sided FFT normalized to magnitude

i0, i1 = 15, 25
f_cont = np.linspace(f[i0], f[i1], 501)

fg1, axx = plt.subplots(1, 2, sharey='all', tight_layout=True,
                        figsize=(6., 3.))
for c_, (ax_, X_, fx_) in enumerate(zip(axx, XX, fcc)):
    Xc_ = (np.sinc(tau * (f_cont - fx_)) +
           np.sinc(tau * (f_cont + fx_))) / 2
    ax_.plot(f_cont, abs(Xc_), f'-C{c_}', alpha=.5, label=rf"$f_x={fx_}\,$Hz")
    m_line, _, _, = ax_.stem(f[i0:i1+1], abs(X_[i0:i1+1]), markerfmt=f'dC{c_}',
                             linefmt=f'-C{c_}', basefmt=' ')
    plt.setp(m_line, markersize=5)

    ax_.legend(loc='upper left', frameon=False)
    ax_.set(xlabel="Frequency $f$ in Hertz", xlim=(f[i0], f[i1]),
            ylim=(0, 0.59))

axx[0].set(ylabel=r'Magnitude $|X(f)/\tau|$')
fg1.suptitle("Continuous and Sampled Magnitude Spectrum ", x=0.55, y=0.93)
fg1.tight_layout()
plt.show()

频率的微小变化会产生外观明显不同的幅度谱，这显然是许多实际应用中不希望出现的行为。可以使用以下两种常用技术来改善频谱：

所谓的“零填充”通过在信号末尾附加零来减小 \(\Delta f\)。要对频率进行 q 倍的过采样，请将参数 n=q*n_x 传递给 fft / rfft 函数，其中 n_x 是输入信号的长度。

第二种技术称为加窗，即用合适的函数乘以输入信号，使得通常以加宽主瓣为代价来抑制次瓣。加窗 DFT 可以表示为：

(6)\[X^w_l := \sum_{k=0}^{n-1} x_k w_k\e^{-2\jj\pi k l / n}\ ,\]

其中 \(w_k\)，\(k=0,\ldots,n-1\) 是采样的窗口函数。为了计算上一小节中给出的频谱表示的采样版本，需要使用以下归一化常数：

\[c^\text{amp}:= \abs{\sum_{k=0}^{n-1} w_k}\ ,\qquad c^\text{den} := \sqrt{\sum_{k=0}^{n-1} \abs{w_k}^2}\]

第一个常数确保频谱中的峰值与该频率下的信号振幅一致。例如，幅度谱可以用 \(|X^w_l / c^\text{amp}|\) 表示。第二个常数保证公式 (5) 中定义的频率区间的功率是一致的。由于不禁止复值窗口，因此需要绝对值。

下图显示了对 \(x(t)\) 应用 hann 窗口和三倍过采样的结果。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq
from scipy.signal.windows import hann

n, T = 100, 0.01  # number of samples and sampling interval
tau = n*T
q = 3  # over-sampling factor
t = np.arange(n) * T
fcc = (20, 20.5)  # frequencies of sines
xx = [np.sin(2 * np.pi * fc_ * t) for fc_ in fcc]  # sine signals
w = hann(n)
c_w = abs(sum(w))  # normalize constant for window

f_X = rfftfreq(n, T)  # frequency bins range from 0 Hz to Nyquist freq.
XX = (rfft(x_ * w) / c_w for x_ in xx)  # one-sided amplitude spectrum
# Oversampled spectrum:
f_Y = rfftfreq(n*q, T)  # frequency bins range from 0 Hz to Nyquist freq.
YY = (rfft(x_ * w, n=q*n) / c_w for x_ in xx)  # one-sided magnitude spectrum

i0, i1 = 15, 25
j0, j1 = i0*q, i1*q

fg1, axx = plt.subplots(1, 2, sharey='all', tight_layout=True,
                        figsize=(6., 3.))
for c_, (ax_, X_, Y_, fx_) in enumerate(zip(axx, XX, YY, fcc)):
    ax_.plot(f_Y[j0:j1 + 1], abs(Y_[j0:j1 + 1]), f'.-C{c_}',
             label=rf"$f_x={fx_}\,$Hz")
    m_ln, s_ln, _, = ax_.stem(f_X[i0:i1 + 1], abs(X_[i0:i1 + 1]), basefmt=' ',
                              markerfmt=f'dC{c_}', linefmt=f'-C{c_}')
    plt.setp(m_ln, markersize=5)
    plt.setp(s_ln, alpha=0.5)

    ax_.legend(loc='upper left', frameon=False)
    ax_.set(xlabel="Frequency $f$ in Hertz", xlim=(f_X[15], f_X[25]),
            ylim=(0, 0.59))

axx[0].set(ylabel=r'Magnitude $|X(f)/\tau|$')
fg1.suptitle(r"Magnitude Spectrum (Hann window, $%d\times$oversampled)" % q,
             x=0.55, y=0.93)
plt.show()

现在，两个波瓣看起来几乎相同，并且旁瓣得到了很好的抑制。20.5 Hz 频谱的最大值也非常接近预期的一半高度。

频谱能量和频谱功率可以类似于公式 (4) 进行计算，得到相同的结果，即：

\[E^w_X = T\sum_{k=0}^{n-1} \abs{\frac{X^w_k}{c^\text{den}}}^2 = E_x\ ,\qquad P^w_X = \frac{1}{\tau} E^w_X = \frac{1}{n} \sum_{k=0}^{n-1} \abs{\frac{X^w_k}{c^\text{den}}}^2 = P_x\ .\]

此公式不要与矩形窗口（或无窗口）的特殊情况混淆，即 \(w_k = 1\)，\(X^w_l=X_l\)，\(c^\text{den}=\sqrt{n}\)，这会导致：

\[E_X = \frac{T}{n}\sum_{k=0}^{n-1} \abs{X_k}^2\ ,\qquad P_X = \frac{1}{n^2} \sum_{k=0}^{n-1} \abs{X_k}^2\ .\]

加窗的频域离散功率谱密度为：

\[S^w_{xx} := \frac{1}{f_S}\abs{\frac{X^w_l}{c^\text{den}}}^2 = T \abs{\frac{X^w_l}{c^\text{den}}}^2\]

定义在频率范围 \([0, f_S)\) 上，可以解释为每个频率间隔 \(\Delta f\) 的功率。像公式 (5) 中一样，对频带 \([l_a\Delta f, l_b\Delta f)\) 进行积分，将变为求和：

\[P_X^{a,b} = \Delta f\sum_{k=l_a}^{l_b-1} S^w_{xx} = \frac{1}{nT}\sum_{k=l_a}^{l_b-1} S^w_{xx}\ .\]

加窗的频域离散能量谱密度 \(\tau S^w_{xx}\) 可以类似地定义。

以上讨论表明，可以定义连续时间情况下频谱表示的采样版本。下表总结了这些：

采样信号的加窗 DFT \(X^w_l\)（如公式 (6) 所示）的频谱表示比较，单位为 \(\text{V}\)：

	频谱	幅度谱	能量谱密度	功率谱密度 (PSD)	幅度谱密度
定义	\(\tau X^w_l / c^\text{amp}\)	\(X^w_l / c^\text{amp}\)	\(\tau T |X^w_l / c^\text{den}|^2\)	\(T |X^w_l / c^\text{den}|^2\)	\(\sqrt{T} X^w_l / c^\text{den}\)
单位	\(\text{V} / \text{Hz}\)	\(\text{V}\)	\(\text{V}^2\text{s} / \text{Hz}\)	\(\text{V}^2 / \text{Hz}\)	\(\text{V} / \sqrt{\text{Hz}}\)

请注意，对于密度，由于确定谱能量/功率时从积分更改为求和，因此 \(\pm f_x\) 处的幅度值与连续时间情况下不同。

虽然汉宁窗是频谱分析中最常用的窗口函数，但还存在其他窗口。下图显示了 windows 子模块中各种窗口函数的幅度谱。它可以解释为单频输入信号的波瓣形状。请注意，仅显示右半部分，并且 \(y\) 轴以分贝为单位，即以对数方式缩放。

import matplotlib.pyplot as plt
import numpy as np

from scipy.fft import rfft, rfftfreq
from scipy.signal import get_window

n, n_zp = 128, 16384  # number of samples without and with zero-padding
t = np.arange(n)
f = rfftfreq(n_zp, 1 / n)

ww = ['boxcar', 'hann', 'hamming', 'tukey', 'blackman', 'flattop']
fg0, axx = plt.subplots(len(ww), 1, sharex='all', sharey='all', figsize=(6., 4.))
for c_, (w_name_, ax_) in enumerate(zip(ww, axx)):
    w_ = get_window(w_name_, n, fftbins=False)
    W_ = rfft(w_ / abs(sum(w_)), n=n_zp)
    W_dB = 20*np.log10(np.maximum(abs(W_), 1e-250))
    ax_.plot(f, W_dB, f'C{c_}-', label=w_name_)
    ax_.text(0.1, -50, w_name_, color=f'C{c_}', verticalalignment='bottom',
             horizontalalignment='left', bbox={'color': 'white', 'pad': 0})
    ax_.set_yticks([-20, -60])
    ax_.grid(axis='x')

axx[0].set_title("Spectral Leakage of various Windows")
fg0.supylabel(r"Normalized Magnitude $20\,\log_{10}|W(f)/c^\operatorname{amp}|$ in dB",
              x=0.04, y=0.5, fontsize='medium')
axx[-1].set(xlabel=r"Normalized frequency $f/\Delta f$ in bins",
            xlim=(0, 9), ylim=(-75, 3))

fg0.tight_layout(h_pad=0.4)
plt.show()

此图表明，窗口函数的选择通常是在主瓣的宽度和旁瓣的高度之间进行权衡。请注意，boxcar 窗口对应于 \(\rect\) 函数，即不加窗。此外，许多描绘的窗口在滤波器设计中比在频谱分析中更常用。

频谱的相位

傅里叶变换的相位（即，angle）通常用于研究信号通过滤波器等系统时频谱分量的时间延迟。在以下示例中，标准测试信号（具有单位功率的脉冲）通过一个简单的滤波器，该滤波器将输入延迟三个样本。输入由 \(n=50\) 个样本组成，采样间隔为 \(T = 1\) s。该图显示了输入和输出信号的幅度和相位随频率的变化。

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np

from scipy import signal
from scipy.fft import rfft, rfftfreq

# Create input signal:
n = 50
x = np.zeros(n)
x[0] = n

# Apply FIR filter which delays signal by 3 samples:
y = signal.lfilter([0, 0, 0, 1], 1, x)

X, Y = (rfft(z_) / n for z_ in (x, y))
f = rfftfreq(n, 1)  # sampling interval T = 1 s

fig = plt.figure(tight_layout=True, figsize=(6., 4.))
gs = gridspec.GridSpec(3, 1)
ax0 = fig.add_subplot(gs[0, :])
ax1 = fig.add_subplot(gs[1:, :], sharex=ax0)

for Z_, n_, m_ in zip((X, Y), ("Input $X(f)$", "Output $Y(f)$"), ('+-', 'x-')):
    ax0.plot(f, abs(Z_), m_, alpha=.5, label=n_)
    ax1.plot(f, np.unwrap(np.angle(Z_)), m_, alpha=.5, label=n_)

ax0.set(title="Frequency Response of 3 Sample Delay Filter (no window)",
        ylabel="Magnitude", xlim=(0, f[-1]), ylim=(0, 1.1))
ax1.set(xlabel=rf"Frequency $f$ in Hertz ($\Delta f = {1/n}\,$Hz)",
        ylabel=r"Phase in rad")
ax1.set_yticks(-np.arange(0, 7)*np.pi/2,
               ['0', '-π/2', '-π', '-3/2π', '-2π', '-4/3π', '-3π'])

ax2 = ax1.twinx()
ax2.set(ylabel=r"Phase in Degrees", ylim=np.rad2deg(ax1.get_ylim()),
        yticks=np.arange(-540, 90, 90))
for ax_ in (ax0, ax1):
    ax_.legend()
    ax_.grid()
plt.show()

输入信号具有单位幅度和零相位傅里叶变换，这是将其用作测试信号的原因。输出信号也具有单位幅度，但相位呈线性下降，斜率为 \(-6\pi\)。这是预期结果，因为将信号 \(x(t)\) 延迟 \(\Delta t\) 会在其傅里叶变换中产生一个额外的线性相位项，即：

\[\int_\IR x(t-\Delta t) \e^{-\jj2\pi f t}\dd t = X(f)\, e^{-\jj2\pi \Delta t f}\ .\]

请注意，在图中，相位并未限制在区间 \((+\pi, \pi]\) （angle 的输出）内，因此没有任何不连续性。这是通过利用 numpy.unwrap 函数实现的。如果已知滤波器的传递函数，则可以使用 freqz 直接确定滤波器的频谱响应。

带平均的频谱

periodogram 函数计算功率谱密度（scaling='density'）或平方幅度谱（scaling='spectrum'）。要获得平滑的周期图，可以使用 welch 函数。它通过将输入信号分成重叠的段来执行平滑，然后计算每个段的加窗 DFT。结果是这些 DFT 的平均值。

下面的示例显示了由一个幅度为 \(\sqrt{2}\,\text{V}\) 的 \(1.27\,\text{kHz}\) 正弦信号和平均频谱功率密度为 \(10^{-3}\,\text{V}^2/\text{Hz}\) 的加性高斯噪声组成的信号的平方幅度谱和功率谱密度。

import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal

rng = np.random.default_rng(73625)  # seeding for reproducibility

fs, n = 10e3, 10_000
f_x, noise_power = 1270, 1e-3 * fs / 2
t = np.arange(n) / fs
x = (np.sqrt(2) * np.sin(2 * np.pi * f_x * t) +
     rng.normal(scale=np.sqrt(noise_power), size=t.shape))

fg, axx = plt.subplots(1, 2, sharex='all', tight_layout=True, figsize=(7, 3.5))
axx[0].set(title="Squared Magnitude Spectrum", ylabel="Square of Magnitude in V²")
axx[1].set(title="Power Spectral Density", ylabel="Power Spectral Density in V²/Hz")
for ax_, s_ in zip(axx, ('spectrum', 'density')):
    f_p, P_p = signal.periodogram(x, fs, 'hann', scaling=s_)
    f_w, P_w = signal.welch(x, fs, scaling=s_)
    ax_.semilogy(f_p/1e3, P_p, label=f"Periodogram ({len(f_p)} bins)")
    ax_.semilogy(f_w/1e3, P_w, label=f"Welch's Method ({len(f_w)} bins)")
    ax_.set(xlabel="Frequency in kHz", xlim=(0, 2), ylim=(1e-7, 1.3))
    ax_.grid(True)
    ax_.legend(loc='lower center')
plt.show()

这些图显示，welch 函数以牺牲频率分辨率为代价，产生了更平滑的噪声基底。由于平滑，正弦波瓣的高度更宽，并且不如周期图中的高。左图可用于读取波瓣的高度，即正弦波平方幅度的一半，为 \(1\,\text{V}^2\)。右图可用于确定 \(10^{-3}\,\text{V}^2/\text{Hz}\) 的噪声基底。请注意，由于有限的频率分辨率，平均平方幅度谱的波瓣高度并非正好为 1。可以使用零填充（例如，将 nfft=4*len(x) 传递给 welch）或通过增加段长度（设置参数 nperseg）来减少段的数量，从而增加频率箱的数量。

Lomb-Scargle 周期图（lombscargle）

最小二乘谱分析 (LSSA) [6] [7] 是一种基于正弦曲线对数据样本的最小二乘拟合来估计频谱的方法，类似于傅里叶分析。傅里叶分析是科学中最常用的频谱方法，它通常会放大长间隔记录中的长周期噪声；而 LSSA 则可以减轻此类问题。

Lomb-Scargle 方法对不均匀采样的数据执行频谱分析，并且已知是查找和测试弱周期信号显著性的强大方法。

对于包含 \(N_{t}\) 个测量值 \(X_{j}\equiv X(t_{j})\) 的时间序列，这些测量值在时间 \(t_{j}\) 处采样，其中 \((j = 1, \ldots, N_{t})\)，假设已缩放和平移，使其均值为零且方差为 1，则频率 \(f\) 处的归一化 Lomb-Scargle 周期图为

\[P_{n}(f) = \frac{1}{2}\left\{\frac{\left[\sum_{j}^{N_{t}}X_{j}\cos\omega(t_{j}-\tau)\right]^{2}}{\sum_{j}^{N_{t}}\cos^{2}\omega(t_{j}-\tau)}+\frac{\left[\sum_{j}^{N_{t}}X_{j}\sin\omega(t_{j}-\tau)\right]^{2}}{\sum_{j}^{N_{t}}\sin^{2}\omega(t_{j}-\tau)}\right\}.\]

其中，\(\omega \equiv 2\pi f\) 是角频率。与频率相关的时间偏移 \(\tau\) 由下式给出：

\[\tan 2\omega\tau = \frac{\sum_{j}^{N_{t}}\sin 2\omega t_{j}}{\sum_{j}^{N_{t}}\cos 2\omega t_{j}}.\]

lombscargle 函数使用 Zechmeister 和 Kürster 创建的稍微修改过的算法 [8] 计算周期图，该算法允许对各个样本进行加权，并独立地计算每个频率的未知偏移（也称为“浮动平均值”）。

短时傅里叶变换

本节提供有关使用 ShortTimeFFT 类的一些背景信息：短时傅里叶变换 (STFT) 可用于分析信号随时间的频谱特性。它通过利用滑动窗口将信号分成重叠的块，并计算每个块的傅里叶变换。对于连续时间复值信号 \(x(t)\)，STFT 定义 [9] 为

\[S(f, t) := \int_\IR x(\xi)\, \conj{w(\xi-t)}\,\e^{-\jj2\pi f \xi}\dd\xi\ ,\]

其中，\(w(t)\) 是一个复值窗函数，其复共轭为 \(\conj{w(t)}\)。它可以解释为确定 \(x\) 与窗口 \(w\) 的标量积，该窗口通过时间 \(t\) 平移，然后通过频率 \(f\) 进行调制（即，频率偏移）。对于使用采样信号 \(x[k] := x(kT)\)，\(k\in\IZ\)，采样间隔为 \(T\)（采样频率 fs 的倒数），离散版本，即仅在离散网格点 \(S[q, p] := S(q \Delta f, p\Delta t)\)，\(q,p\in\IZ\) 处评估 STFT，需要使用。它可以表示为

(7)\[S[q,p] = \sum_{k=0}^{N-1} x[k]\,\conj{w[k-p h]}\, \e^{-\jj2\pi q k / N}\ , \quad q,p\in\IZ\ ,\]

其中，p 表示 \(S\) 的时间索引，时间间隔为 \(\Delta t := h T\)，\(h\in\IN\)（请参阅 delta_t），这可以表示为 \(h\) 个样本的 hop 大小。\(q\) 表示 \(S\) 的频率索引，步长为 \(\Delta f := 1 / (N T)\)（请参阅 delta_f），这使其与 FFT 兼容。\(w[m] := w(mT)\)，\(m\in\IZ\) 是采样的窗函数。

为了更符合 ShortTimeFFT 的实现，将公式 (7) 重述为两步过程是有意义的。

1. 通过使用由 \(M\) 个样本组成的窗口 \(w[m]\)（请参阅 m_num），提取以 \(t[p] := p \Delta t = h T\) 为中心的第 \(p\) 个切片（请参阅 delta_t），即：

   (8)\[x_p[m] = x\!\big[m - \lfloor M/2\rfloor + h p\big]\, \conj{w[m]}\ , \quad m = 0, \ldots M-1\ ,\]

   其中整数 \(\lfloor M/2\rfloor\) 表示 M//2，即，它是窗口的中点（m_num_mid）。为方便表示，假设 \(k\not\in\{0, 1, \ldots, N-1\}\) 时 \(x[k]:=0\)。在子章节 滑动窗口 中，将更详细地讨论切片的索引。

2. 然后对\(x_p[m]\)执行离散傅里叶变换（即FFT）。

   (9)\[S[q, p] = \sum_{m=0}^{M-1} x_p[m] \exp\!\big\{% -2\jj\pi (q + \phi_m)\, m / M\big\}\ .\]

   请注意，可以指定线性相位\(\phi_m\)（参见phase_shift），这相当于将输入移动\(\phi_m\)个采样点。默认值为\(\phi_m = \lfloor M/2\rfloor\)（根据定义，对应于phase_shift = 0），它会抑制未移位信号的线性相位分量。此外，可以通过用零填充\(w[m]\)来对FFT进行过采样。这可以通过指定mfft大于窗口长度m_num来实现——这将\(M\)设置为mfft（这意味着对于\(m\not\in\{0, 1, \ldots, M-1\}\)，\(w[m]:=0\)也成立）。

逆短时傅里叶变换（istft）是通过反转这两个步骤实现的

1. 执行逆离散傅里叶变换，即：

   \[x_p[m] = \frac{1}{M}\sum_{q=0}^M S[q, p]\, \exp\!\big\{ 2\jj\pi (q + \phi_m)\, m / M\big\}\ .\]

2. 将加权的移位切片与\(w_d[m]\)相加，以重建原始信号，即：

   \[x[k] = \sum_p x_p\!\big[\mu_p(k)\big]\, w_d\!\big[\mu_p(k)\big]\ ,\quad \mu_p(k) = k + \lfloor M/2\rfloor - h p\]

   对于\(k \in [0, \ldots, n-1]\)。\(w_d[m]\)是\(w[m]\)的所谓规范对偶窗，也由\(M\)个采样组成。

请注意，对于所有窗口和跳跃大小，并不一定存在逆STFT。对于给定的窗口\(w[m]\)，跳跃大小\(h\)必须足够小，以确保\(x[k]\)的每个采样至少被一个窗口切片的非零值触及。这有时被称为“非零重叠条件”（参见check_NOLA）。有关更多详细信息，请参阅子章节逆 STFT 和对偶窗口。

滑动窗口

本小节通过一个示例讨论如何在ShortTimeFFT中索引滑动窗口：考虑一个长度为6的窗口，其hop间隔为2，采样间隔T为1，例如ShortTimeFFT (np.ones(6), 2, fs=1)。下图示意性地描述了前四个窗口位置，也称为时间切片

x轴表示时间\(t\)，它对应于底部蓝色框行指示的采样索引k。y轴表示时间切片索引\(p\)。信号\(x[k]\)从索引\(k=0\)开始，并以浅蓝色背景标记。根据定义，第零个切片（\(p=0\)）以\(t=0\)为中心。每个切片的中心（m_num_mid），这里是采样点6//2=3，用文本“mid”标记。默认情况下，stft计算所有与信号有部分重叠的切片。因此，第一个切片位于p_min = -1，最低的采样索引为k_min = -5。第一个不受切片影响的采样索引，该切片不会突出到信号的左侧，是\(p_{lb} = 2\)，第一个不受边界效应影响的采样索引是\(k_{lb} = 5\)。属性lower_border_end返回元组\((k_{lb}, p_{lb})\)。

下图描绘了一个具有\(n=50\)个采样的信号在信号末端的行为，如蓝色背景所示

此处，最后一个切片的索引为\(p=26\)。因此，按照Python约定，结束索引超出范围，p_max = 27 表示第一个不接触信号的切片。相应的采样索引为k_max = 55。第一个突出到右侧的切片是\(p_{ub} = 24\)，其第一个采样位于\(k_{ub}=45\)。函数upper_border_begin返回元组\((k_{ub}, p_{ub})\)。

逆STFT和对偶窗口

对偶窗口一词源于框架理论[10]，其中框架是一种级数展开，可以表示给定希尔伯特空间中的任何函数。在那里，如果对于给定希尔伯特空间\(\mathcal{H}\)中的所有函数\(f\)，展开\(\{g_k\}\)和\(\{h_k\}\)是对偶框架，则满足以下条件：

\[f = \sum_{k\in\IN} \langle f, g_k\rangle h_k = \sum_{k\in\IN} \langle f, h_k\rangle g_k\ , \quad f \in \mathcal{H}\ ,\]

成立，其中\(\langle ., .\rangle\)表示\(\mathcal{H}\)的标量积。所有框架都有对偶框架[10]。

在文献中，仅在离散网格点\(S(q \Delta f, p\Delta t)\)上计算的STFT称为“Gabor框架”[9][10]。由于窗口\(w[m]\)的支持仅限于有限区间，因此ShortTimeFFT属于所谓的“无痛非正交展开”类[9]。在这种情况下，对偶窗口始终具有相同的支持，并且可以通过反转对角矩阵来计算。下面将简要概述一个粗略的推导，它只需要对操作矩阵有一些了解。

由于公式(7)中给出的STFT是\(x[k]\)中的线性映射，因此可以用向量矩阵符号表示。这允许我们通过线性最小二乘法的正式解（如lstsq中）表示逆，这导致了一个美观而简单的结果。

我们首先重新制定公式(8)的加窗操作

(10)\[\begin{split} \vb{x}_p = \vb{W}_{\!p}\,\vb{x} = \begin{bmatrix} \cdots & 0 & w[0] & 0 & \cdots&&&\\ & \cdots & 0 & w[1] & 0 & \cdots&&\\ & & & & \ddots&&&\\ &&\cdots & 0 & 0 & w[M-1] & 0 & \cdots \end{bmatrix}\begin{bmatrix} x[0]\\ x[1]\\ \vdots\\ x[N-1] \end{bmatrix}\ ,\end{split}\]

其中\(M\times N\)矩阵\(\vb{W}_{\!p}\)仅在第\((ph)\)个次对角线上具有非零项，即：

(11)\[\begin{split} W_p[m,k] = w[m]\, \delta_{m+ph,k}\ ,\quad \delta_{k,l} &= \begin{cases} 1 & \text{ 如果 } k=l\ ,\\ 0 & \text{ 其他情况 ,} \end{cases}\end{split}\]

其中\(\delta_{k,l}\)是克罗内克δ。公式(9)可以表示为

\[\vb{s}_p = \vb{F}\,\vb{x}_p \quad\text{，其中}\quad F[q,m] =\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,\]

这允许将第\(p\)个切片的STFT写为

(12)\[\vb{s}_p = \vb{F}\vb{W}_{\!p}\,\vb{x} =: \vb{G}_p\,\vb{x} \quad\text{，其中}\quad s_p[q] = S[p,q]\ .\]

请注意，\(\vb{F}\)是酉矩阵，即逆矩阵等于其共轭转置，这意味着\(\conjT{\vb{F}}\vb{F} = \vb{I}\)。

为了获得STFT的单个向量矩阵方程，将切片堆叠成一个向量，即：

\[\begin{split}\vb{s} := \begin{bmatrix} \vb{s}_0\\ \vb{s}_1\\ \vdots\\ \vb{s}_{P-1} \end{bmatrix} = \begin{bmatrix} \vb{G}_0\\ \vb{G}_1\\ \vdots\\ \vb{G}_{P-1} \end{bmatrix}\, \vb{x} =: \vb{G}\, \vb{x}\ ,\end{split}\]

其中 \(P\) 是结果 STFT 的列数。为了反转这个方程，可以使用 Moore-Penrose 逆 \(\vb{G}^\dagger\)

(13)\[\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}\, \conjT{\vb{G}} \vb{s} =: \vb{G}^\dagger \vb{s}\ ,\]

当且仅当

\[\begin{split}\vb{D} := \conjT{\vb{G}}\vb{G} = \begin{bmatrix} \conjT{\vb{G}_0}& \conjT{\vb{G}_1}& \cdots & \conjT{\vb{G}_{P-1}} \end{bmatrix}^T \begin{bmatrix} \vb{G}_0\\ \vb{G}_1\\ \vdots\\ \vb{G}_{P-1} \end{bmatrix} = \sum_{p=0}^{P-1} \conjT{\vb{G}_p}\vb{G}_p \ .\end{split}\]

可逆时，该式成立。那么 \(\vb{x} = \vb{G}^\dagger\vb{G}\,\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}\,\conjT{\vb{G}}\vb{G}\,\vb{x}\) 显然成立。\(\vb{D}\) 始终是一个对角矩阵，具有非负对角线元素。当进一步简化 \(\vb{D}\) 到以下形式时，这一点会很清楚：

(14)\[\vb{D} = \sum_{p=0}^{P-1} \conjT{\vb{G}_p}\vb{G}_p = \sum_{p=0}^{P-1} \conjT{(\vb{F}\,\vb{W}_{\!p})}\, \vb{F}\,\vb{W}_{\!p} = \sum_{p=0}^{P-1} \conjT{\vb{W}_{\!p}}\vb{W}_{\!p} =: \sum_{p=0}^{P-1} \vb{D}_p\]

这是因为 \(\vb{F}\) 是酉矩阵。此外，

(15)\[\begin{split}D_p[r,s] &= \sum_{m=0}^{M-1} \conj{W_p^T[r,m]}\,W_p[m,s] = \sum_{m=0}^{M-1} \left(\conj{w[m]}\, \delta_{m+ph,r}\right) \left(w[m]\, \delta_{m+ph,s}\right)\\ &= \sum_{m=0}^{M-1} \big|w[m]\big|^2\, \delta_{r,s}\, \delta_{r,m+ph}\ .\end{split}\]

表明 \(\vb{D}_p\) 是一个具有非负元素的对角矩阵。因此，对 \(\vb{D}_p\) 求和会保留该属性。这允许进一步简化公式 (13)，即：

(16)\[\begin{split}\vb{x} &= \vb{D}^{-1} \conjT{\vb{G}}\vb{s} = \sum_{p=0}^{P-1} \vb{D}^{-1}\conjT{\vb{W}_{\!p}}\, \conjT{\vb{F}}\vb{s}_p = \sum_{p=0}^{P-1} (\conj{\vb{W}_{\!p}\vb{D}^{-1}})^T\, \conjT{\vb{F}}\vb{s}_p\\ &=: \sum_{p=0}^{P-1}\conjT{\vb{U}_p}\,\conjT{\vb{F}}\vb{s}_p\ .\end{split}\]

利用公式 (11)、(14)、(15)，\(\vb{U}_p=\vb{W}_{\!p}\vb{D}^{-1}\) 可以表示为：

(17)\[\begin{split}U_p[m, k] &= W[m,k]\, D^{-1}[k,k] = \left(w[m] \delta_{m+ph,k}\right) \inv{\sum_{\eta=0}^{P-1} \vb{D}_\eta[k,k]} \delta_{m+ph,k}\\ &= w[m] \inv{\sum_{\eta=0}^{P-1}\sum_{\mu=0}^{M-1} \big|w[\mu]\big|^2\,\delta_{m+ph, \mu+\eta h}} \delta_{m+ph,k}\\ &= w[m] \inv{\sum_{\eta=0}^{P-1} \big|w[m+(p-\eta)h]\big|^2} \delta_{m+ph,k} \ .\end{split}\]

这表明 \(\vb{U}_p\) 具有与公式 (11) 中 \(\vb{W}_p\) 相同的结构，即仅在第 \((ph)\) 个次对角线上具有非零条目。逆中的求和项可以解释为在 \(w[m]\) 上滑动 \(|w[\mu]|^2\) (并进行反转)，因此只有与 \(w[m]\) 重叠的分量才会产生影响。因此，所有远离边界的 \(U_p[m, k]\) 都是相同的窗口。为了规避边界效应，\(x[k]\) 用零填充，扩大了 \(\vb{U}\)，从而使所有接触 \(x[k]\) 的切片都包含相同的对偶窗口。

\[w_d[m] = w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta\, h]\big|^2}\ .\]

由于 \(w[m] = 0\) 对 \(m \not\in\{0, \ldots, M-1\}\) 成立，所以只需要对满足 \(|\eta| < M/h\) 的索引 \(\eta\) 求和。可以通过将公式 (12) 插入公式 (16) 来证明对偶窗口的名称是合理的，即：

\[\vb{x} = \sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\conjT{\vb{F}}\, \vb{F}\,\vb{W}_{\!p}\,\vb{x} = \left(\sum_{p=0}^{P-1} \conjT{\vb{U}_p}\,\vb{W}_{\!p}\right)\vb{x}\ ,\]

表明 \(\vb{U}_p\) 和 \(\vb{W}_{\!p}\) 是可以互换的。因此，\(w_d[m]\) 也是一个有效的窗口，其对偶窗口为 \(w[m]\)。请注意，由于 \(\vb{s}\) 通常比 \(\vb{x}\) 具有更多的条目，因此 \(w_d[m]\) 不是唯一的对偶窗口。可以证明，\(w_d[m]\) 具有最小的能量（或 \(L_2\) 范数）[9]，这也是它被称为“规范对偶窗口”的原因。

与传统实现的比较

函数 stft、istft 和 spectrogram 早于 ShortTimeFFT 实现。本节讨论较旧的“传统”实现和较新的 ShortTimeFFT 实现之间的主要区别。重写的主要动机是认识到，如果不破坏兼容性，就无法以合理的方式集成对偶窗口。这为重新思考代码结构和参数化提供了机会，从而使一些隐式行为更加明确。

下面的示例比较了具有负斜率的复值线性调频信号的两个 STFT

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.fft import fftshift
>>> from scipy.signal import stft, istft, spectrogram, ShortTimeFFT
...
>>> fs, N = 200, 1001  # 200 Hz sampling rate for 5 s signal
>>> t_z = np.arange(N) / fs  # time indexes for signal
>>> z = np.exp(2j*np.pi*70 * (t_z - 0.2*t_z**2))  # complex-valued chirp
...
>>> nperseg, noverlap = 50, 40
>>> win = ('gaussian', 1e-2 * fs)  # Gaussian with 0.01 s standard dev.
...
>>> # Legacy STFT:
>>> f0_u, t0, Sz0_u = stft(z, fs, win, nperseg, noverlap,
...                        return_onesided=False, scaling='spectrum')
>>> f0, Sz0 = fftshift(f0_u), fftshift(Sz0_u, axes=0)
...
>>> # New STFT:
>>> SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap, fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz1 = SFT.stft(z)
...
>>> # Plot results:
>>> fig1, axx = plt.subplots(2, 1, sharex='all', sharey='all',
...                          figsize=(6., 5.))  # enlarge figure a bit
>>> t_lo, t_hi, f_lo, f_hi = SFT.extent(N, center_bins=True)
>>> axx[0].set_title(r"Legacy stft() produces $%d\times%d$ points" % Sz0.T.shape)
>>> axx[0].set_xlim(t_lo, t_hi)
>>> axx[0].set_ylim(f_lo, f_hi)
>>> axx[1].set_title(r"ShortTimeFFT produces $%d\times%d$ points" % Sz1.T.shape)
>>> axx[1].set_xlabel(rf"Time $t$ in seconds ($\Delta t= {SFT.delta_t:g}\,$s)")
...
>>> # Calculate extent of plot with centered bins since
>>> # imshow does not interpolate by default:
>>> dt2 = (nperseg-noverlap) / fs / 2  # equals SFT.delta_t / 2
>>> df2 = fs / nperseg / 2  # equals SFT.delta_f / 2
>>> extent0 = (-dt2, t0[-1] + dt2, f0[0] - df2, f0[-1] - df2)
>>> extent1 = SFT.extent(N, center_bins=True)
...
>>> kw = dict(origin='lower', aspect='auto', cmap='viridis')
>>> im1a = axx[0].imshow(abs(Sz0), extent=extent0, **kw)
>>> im1b = axx[1].imshow(abs(Sz1), extent=extent1, **kw)
>>> fig1.colorbar(im1b, ax=axx, label="Magnitude $|S_z(t, f)|$")
>>> _ = fig1.supylabel(r"Frequency $f$ in Hertz ($\Delta f = %g\,$Hz)" %
...                    SFT.delta_f, x=0.08, y=0.5, fontsize='medium')
>>> plt.show()

ShortTimeFFT 生成的时间片比传统版本多 3 个，这是主要区别。正如滑动窗口部分所述，所有接触信号的切片都包含在新版本中。这样做的好处是，STFT 可以像 ShortTimeFFT 代码示例中所示那样进行切片和重组。此外，在某些地方为零的窗口的情况下，使用所有接触的切片可以使 ISTFT 更加健壮。

请注意，具有相同时间戳的切片会产生相同的结果（在数值精度范围内），即：

>>> np.allclose(Sz0, Sz1[:, 2:-1])
True

通常，这些额外的切片包含非零值。由于我们的示例中存在较大的重叠，因此它们很小。例如：

>>> abs(Sz1[:, 1]).min(), abs(Sz1[:, 1]).max()
(6.925060911593139e-07, 8.00271269218721e-07)

ISTFT 可用于重建原始信号

>>> t0_r, z0_r = istft(Sz0_u, fs, win, nperseg, noverlap,
...                    input_onesided=False, scaling='spectrum')
>>> z1_r = SFT.istft(Sz1, k1=N)
...
>>> len(z0_r), len(z)
(1010, 1001)
>>> np.allclose(z0_r[:N], z)
True
>>> np.allclose(z1_r, z)
True

请注意，传统实现返回的信号比原始信号长。另一方面，新的 istft 允许显式指定重建信号的起始索引 k0 和结束索引 k1。旧实现中的长度差异是由于信号长度不是切片的倍数造成的。

此示例中新版本和旧版本之间的其他差异是：

• 参数 fft_mode='centered' 确保在图中，对于双边 FFT，零频率在垂直方向上居中。对于传统实现，需要使用 fftshift。fft_mode='twosided' 产生与旧版本相同的行为。

• 参数 phase_shift=None 确保两个版本的相位相同。ShortTimeFFT 的默认值 0 会生成带有附加线性相位项的 STFT 切片。

频谱图定义为 STFT 的绝对平方值[9]。spectrogram 由 ShortTimeFFT 提供，它坚持该定义，即：

>>> np.allclose(SFT.spectrogram(z), abs(Sz1)**2)
True

另一方面，传统的 spectrogram 提供了另一个 STFT 实现，主要区别在于对信号边界的不同处理。以下示例演示如何使用 ShortTimeFFT 来获得与传统 spectrogram 生成的 SFT 相同的 SFT

>>> # Legacy spectrogram (detrending for complex signals not useful):
>>> f2_u, t2, Sz2_u = spectrogram(z, fs, win, nperseg, noverlap,
...                               detrend=None, return_onesided=False,
...                               scaling='spectrum', mode='complex')
>>> f2, Sz2 = fftshift(f2_u), fftshift(Sz2_u, axes=0)
...
>>> # New STFT:
... SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap,
...                                fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz3 = SFT.stft(z, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
>>> t3 = SFT.t(N, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
...
>>> np.allclose(t2, t3)
True
>>> np.allclose(f2, SFT.f)
True
>>> np.allclose(Sz2, Sz3)
True

与其他 STFT 的区别在于，时间片不是从 0 开始，而是从 nperseg//2 开始，即：

>>> t2
array([0.125, 0.175, 0.225, 0.275, 0.325, 0.375, 0.425, 0.475, 0.525,
       ...
       4.625, 4.675, 4.725, 4.775, 4.825, 4.875])

此外，只返回不向右突出的切片，将最后一个切片居中在 4.875 秒，这使其比默认的 stft 参数化更短。

使用 mode 参数，旧版的 spectrogram 也可以返回 ‘angle’，‘phase’，‘psd’ 或 ‘magnitude’。旧版 spectrogram 的 scaling 行为并不直接，因为它取决于参数 mode，scaling 和 return_onesided。ShortTimeFFT 中并非所有组合都有直接对应关系，因为它只提供 ‘magnitude’，‘psd’ 或根本不进行窗口 scaling。下表显示了这些对应关系。

旧版 spectrogram				ShortTimeFFT
模式（mode）	缩放（scaling）	返回单边（return_onesided）		fft_模式（fft_mode）	缩放（scaling）
功率谱密度（psd）	密度（density）	True		单边2X（onesided2X）	功率谱密度（psd）
功率谱密度（psd）	密度（density）	False		双边（twosided）	功率谱密度（psd）
幅度（magnitude）	频谱（spectrum）	True		单边（onesided）	幅度（magnitude）
幅度（magnitude）	频谱（spectrum）	False		双边（twosided）	幅度（magnitude）
复数（complex）	频谱（spectrum）	True		单边（onesided）	幅度（magnitude）
复数（complex）	频谱（spectrum）	False		双边（twosided）	幅度（magnitude）
功率谱密度（psd）	频谱（spectrum）	True		—	—
功率谱密度（psd）	频谱（spectrum）	False		—	—
复数（complex）	密度（density）	True		—	—
复数（complex）	密度（density）	False		—	—
幅度（magnitude）	密度（density）	True		—	—
幅度（magnitude）	密度（density）	False		—	—
—	—	—		*	None

当在复数值输入信号上使用 onesided 输出时，旧版 spectrogram 会切换到 双边 模式。ShortTimeFFT 会引发 TypeError，因为所使用的 rfft 函数只接受实数值输入。请参阅上面的 频谱分析 部分，了解有关各种参数化引起的各种频谱表示的讨论。

参考文献

一些深入阅读和相关软件

[1]

“傅里叶变换”，维基百科，https://en.wikipedia.org/wiki/Fourier_transform

[2]

“频谱密度”，维基百科，https://en.wikipedia.org/wiki/Spectral_density

[3]

“离散傅里叶变换”，维基百科，https://en.wikipedia.org/wiki/Discrete_Fourier_transform

[4]

“频谱泄漏”，维基百科，https://en.wikipedia.org/wiki/Spectral_leakage

[5]

“奈奎斯特-香农采样定理”，维基百科，https://en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem

[6]

N.R. Lomb，“非均匀间隔数据的最小二乘频率分析”，《天体物理学和空间科学》，第 39 卷，第 447-462 页，1976 年

[7]

J.D. Scargle，“天文时间序列分析研究。II - 不均匀间隔数据频谱分析的统计方面”，《天体物理学杂志》，第 263 卷，第 835-853 页，1982 年

[8]

M. Zechmeister 和 M. Kürster，“广义的隆布-斯卡格勒周期图。浮动平均和开普勒周期图的新形式主义”，《天文学和天体物理学》，第 496 卷，第 577-584 页，2009 年

[9] (1,2,3,4,5)

Karlheinz Gröchenig：“时频分析基础”，Birkhäuser Boston 2001，DOI:10.1007/978-1-4612-0003-1

[10] (1,2,3)

Ole Christensen：“框架和 Riesz 基简介”，Birkhäuser Boston 2016，DOI:10.1007/978-3-319-25613-9