信号处理 (scipy.signal)

信号处理工具箱目前包含一些滤波函数，有限的滤波器设计工具集以及一些用于 1 维和 2 维数据的 B 样条插值算法。虽然 B 样条算法在技术上可以归类为插值类别，但它们包含在这里是因为它们只对等间距数据有效，并且大量使用滤波器理论和传递函数形式来提供快速 B 样条变换。要理解本节内容，您需要了解 SciPy 中的信号是一个实数或复数数组。

B 样条

B 样条是用 B 样条系数和节点点来近似有限域上的连续函数。如果节点点等间距，间距为 \(\Delta x\)，则 1 维函数的 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.\) 索引的 1 维信号。两个 1 维信号的完全卷积可以表示为

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

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

\[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]\) 是输出序列。如果我们假设初始静止，使得 \(y\left[n\right]=0\) 对于 \(n<0\)，那么这种滤波器可以使用卷积来实现。但是，如果 \(a_{k}\neq0\) 对于 \(k\geq1.\)，则卷积滤波器序列 \(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).\]

使用此公式，我们可以找到初始条件向量 \(z_{0}\left[-1\right]\) 到 \(z_{K-1}\left[-1\right]\)，前提是在 \(y\)（以及 \(x\)）上给出初始条件。命令 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 滤波器相比，滤波器阶数（阶数 4）低得多，以达到相同的 \(\approx 60\) dB 阻带衰减。

>>> 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()

滤波器系数

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

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

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

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

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

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

传递函数表示

格式 ba 或 tf 是一个 2 元组 (b, a)，表示传递函数，其中 b 是 M 阶分子多项式的系数的长度为 M+1 的数组，而 a 是 N 阶分母的系数的长度为 N+1 的数组，作为传递函数变量的正的降幂。因此，\(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 是一个 3 元组 (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 是一个由数组组成的 4 元组 (A, B, C, D)，表示以下形式的 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 是采样时间（采样频率的倒数）。

其他滤波器

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

中值滤波器

当噪声明显是非高斯噪声或需要保留边缘时，通常会应用中值滤波器。中值滤波器的工作原理是，对围绕感兴趣点周围的矩形区域内的所有数组像素值进行排序。此邻域像素值列表的样本中值用作输出数组的值。样本中值是排序的邻域值列表中的中间值。如果邻域中元素数量为偶数，则使用中间两个值的平均值作为中值。 medfilt 是一个适用于 N 维数组的通用中值滤波器。一个专门针对二维数组的版本可作为 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}\) 的概念密切相关。

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

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

	频谱	幅度谱	能量谱密度	功率谱密度 (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_z\) 处的完整（不是一半）幅度正弦，而 PSD 中一个间隔的面积表示其完整（不是一半）功率。需要注意的是，对于幅度谱密度，正值不是加倍，而是乘以 \(\sqrt{2}\)，因为它是对 PSD 的平方根。此外，并没有对加倍的频谱进行规范命名的方法。

下图显示了四个正弦信号 \(x(t)\) 的三种不同的频谱表示，这些信号具有不同的幅度 \(a\) 和持续时间 \(\tau\)，它们都符合公式 (1)。为了减少混乱，频谱以 \(f_z\) 为中心，并并排绘制

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_z\) 处的峰值表示正弦信号幅度 \(|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)\ .\]

下图显示了两个正弦信号的幅度谱，这两个正弦信号的幅度为 1，频率分别为 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 处，所有值都为零。它在 20 Hz 处为 0.5，这与公式 (1) 中输入信号幅度的一半一致。另一方面，20.5 Hz 信号的峰值沿频率轴分散。公式 (3) 表明，这种差异是由于 20 Hz 是 1 Hz 频谱 bin 宽度 的倍数，而 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) 中定义的频率间隔的功率一致。由于复值窗是允许的，因此需要使用绝对值。

下图显示了将 hann 窗和三次过采样应用于 \(x(t)\) 的结果

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 频谱的最大值也非常接近预期的 0.5 高度。

频谱能量和频谱功率可以类似于公式 (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\) 的功率。对频率带 \([l_a\Delta f, l_b\Delta f)\) 进行积分，如公式 (5) 中所示，将变为求和

\[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 的平均值。

下面的示例显示了由 \(1.27\,\text{kHz}\) 正弦信号（幅度为 \(\sqrt{2}\,\text{V}\)）和具有 \(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 函数使用 Townsend [8] 的略微修改过的算法来计算周期图，该算法允许使用每个频率的输入数组的单次遍历来计算周期图。

该方程被重构为

\[P_{n}(f) = \frac{1}{2}\left[\frac{(c_{\tau}XC + s_{\tau}XS)^{2}}{c_{\tau}^{2}CC + 2c_{\tau}s_{\tau}CS + s_{\tau}^{2}SS} + \frac{(c_{\tau}XS - s_{\tau}XC)^{2}}{c_{\tau}^{2}SS - 2c_{\tau}s_{\tau}CS + s_{\tau}^{2}CC}\right]\]

和

\[\tan 2\omega\tau = \frac{2CS}{CC-SS}.\]

这里，

\[c_{\tau} = \cos\omega\tau,\qquad s_{\tau} = \sin\omega\tau,\]

而求和为

\[\begin{split}XC &= \sum_{j}^{N_{t}} X_{j}\cos\omega t_{j}\\ XS &= \sum_{j}^{N_{t}} X_{j}\sin\omega t_{j}\\ CC &= \sum_{j}^{N_{t}} \cos^{2}\omega t_{j}\\ SS &= \sum_{j}^{N_{t}} \sin^{2}\omega t_{j}\\ CS &= \sum_{j}^{N_{t}} \cos\omega t_{j}\sin\omega t_{j}.\end{split}\]

这需要 \(N_{f}(2N_{t}+3)\) 个三角函数计算，比直接实现快约 \(\sim 2\) 倍。

短时傅立叶变换

本节介绍使用 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）对信号进行加窗，提取第 \(p\) 个切片，该窗以 \(t[p] := p \Delta t = h T\)（见 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）。 为了便于表示，假设 \(x[k]:=0\) 对于 \(k\not\in\{0, 1, \ldots, N-1\}\)。 在 滑动窗口 小节中将更详细地讨论切片的索引。

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. 将加权后的移位切片求和，以重建原始信号，即

   \[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]，其中框架是希尔伯特空间中的一系列展开，可以表示该空间中的任何函数。在该理论中，展开 \(\{g_k\}\) 和 \(\{h_k\}\) 是对偶框架，如果对于给定希尔伯特空间 \(\mathcal{H}\) 中的所有函数 \(f\)

\[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{ for } k=l\ ,\\ 0 & \text{ elsewhere ,} \end{cases}\end{split}\]

其中 \(\delta_{k,l}\) 为克罗内克德尔塔。式 (9) 可以表示为

\[\vb{s}_p = \vb{F}\,\vb{x}_p \quad\text{with}\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{with}\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 的列数。为了对该方程求逆，可以使用摩尔-彭罗斯逆 \(\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[\mu]|^2\) 滑过 \(w[m]\)（包含反转操作），因此仅与 \(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]\)。注意，\(w_d[m]\) 不是唯一的对偶窗口，因为 \(\vb{s}\) 通常比 \(\vb{x}\) 具有更多元素。可以证明，\(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 相同的结果

>>> # 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_mode	scaling
psd	density	True		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 会切换到 two-sided 模式。 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]

R.H.D. Townsend，“使用图形处理单元快速计算 Lomb-Scargle 周期图。”，天体物理学杂志增刊，第 191 卷，第 247-253 页，2010 年

[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