信号处理 (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\) 的向量相乘。仅存储该 \(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\) 为使得对于所有 \(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\) 的（互）相关。对于在 \(\left[0,K\right]\) 范围之外 \(y\left[n\right]=0\) 且在 \(\left[0,M\right]\) 范围之外 \(x\left[n\right]=0\) 的有限长度信号，求和可以简化为

\[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\)（默认情况），然后通过 lfiltic 针对 \(y[-1] = 2\) 进行计算。

>>> 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 根据窗函数法设计滤波器。根据提供的参数，该函数返回不同的滤波器类型（例如，低通、带通等）。以下示例设计了一个低通滤波器和一个带阻滤波器，并绘制了它们的频率响应。0.5 Hz 的低通截止频率和 0.3 Hz 至 0.8 Hz 的阻带区间用垂直虚线表示。由于采样频率设置为 2 Hz，因此仅绘制 0 Hz 到 1 Hz 奈奎斯特频率之间的频率范围。

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

fs = 2  # sampling frequency
n0, f0_c = 40, 0.5  # number of taps and cut-off frequency for H_0(f)
n1, ff1_c = 41, [0.3, 0.8]  # number of taps and cut-off frequencies for H_1(f)

b0 = signal.firwin(n0, f0_c, fs=fs)  # design FIR filters
b1 = signal.firwin(n1, ff1_c, fs=fs)

f0, H0 = signal.freqz(b0, fs=fs)  # calculate frequency responses
f1, H1 = signal.freqz(b1, fs=fs)
H0_dB, H1_dB = (20 * np.log10(np.abs(H_)) for H_ in (H0, H1))  # convert to dB

# do the plotting:
fg0, (ax0, ax1) = plt.subplots(2, 1, sharex='all', layout="constrained",
                               figsize=(6, 4))
ax0.set(title=f"Frequency Response of {n0}-tap Low-pass FIR-Filter", xlim=(0, fs/2))
ax0.plot(f0, H0_dB, 'C0', label="Response $|H_0(f)|$")
ax0.axvline(f0_c, color='C2', linestyle='--', label=rf"$f_0={f0_c}\,$Hz")

ax1.set_title(f"Frequency Response of {n1}-tap Band-stop FIR-Filter")
ax1.set_xlabel(rf"Frequency $f\,$ in hertz (sampling frequency $f_S={fs}\,$Hz)")
ax1.plot(f1, H1_dB, 'C0', label="Response $|H_1(f)|$")
ax1.axvline(ff1_c[0], color='C2', linestyle='--', label=rf"$f_0={ff1_c[0]}\,$Hz")
ax1.axvline(ff1_c[1], color='C3', linestyle='--', label=rf"$f_1={ff1_c[1]}\,$Hz")

for ax_ in (ax0, ax1):
    ax_.set_ylabel("Gain in dB")
    ax_.legend()
    ax_.grid()
plt.show()

函数 firwin2 允许通过指定一组截止频率和相应的增益来设计几乎任意的频率响应。下面的示例设计了一个具有这种任意幅度响应的滤波器，并在线性幅度标度上绘制了其频率响应。四个指定的增益截止点由用虚线连接的灰色点表示，而滤波器的响应由连续的蓝线描绘。

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

fs, numtaps = 2, 41  # sampling frequency and number of taps
freqs = [0.0, 0.3, 0.6, 1.0]  # corner frequencies
gains = [1.0, 2.0, 0.5, 0.8]  # desired gains at corner frequencies

b = signal.firwin2(numtaps, freqs, gains, fs=fs)  # design filter
f, H = signal.freqz(b, fs=fs)  # calculate frequency response

fg, ax = plt.subplots(1, 1, layout="constrained")  # do the plotting
ax.set(title=f"Frequency Response of {numtaps}-tap Band-pass FIR-Filter",
       ylabel="Gain", xlim=(0, fs/2),
       xlabel=rf"Frequency $f\,$ in hertz (sampling frequency $f_S={fs}\,$Hz)")
ax.plot(freqs, gains, 'ko--', alpha=.5, label="Desired Gains")
ax.plot(f, np.abs(H), label="Response $|H(f)|$")
ax.grid()
ax.legend()
plt.show()

可以通过增加抽头（taps）的数量来减少所需增益与计算响应之间的偏差。

注意，函数 firwin、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 格式是一个 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 是采样时间（采样频率的倒数）。

其他滤波器

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

中值滤波器

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

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

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

	频谱	幅度谱	能量谱密度	功率谱密度 (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 的平方根。此外，没有加倍频谱的标准命名方式。

下图显示了具有不同振幅 \(a\) 和持续时间 \(\tau\) 的方程 (1) 的四个正弦信号 \(x(t)\) 的三种不同频谱表示。为了减少混乱，频谱以 \(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 计算实现）由下式给出

(6)\[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 可以表示为

(7)\[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(rf"Magnitude Spectrum (Hann window, ${q}\times$oversampled)",
             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\) 的功率。在频率带 \([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}\) 可以类似地定义。

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

采样信号单位为 \(\text{V}\) 时方程 (7) 的加窗 DFT \(X^w_l\) 的频谱表示比较：

	频谱	幅度谱	能量谱密度	功率谱密度 (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\) 处的量值与连续时间情况不同，这是由于在确定频谱能量/功率时从积分变为求和。

虽然 hann 窗是频谱分析中最常用的窗函数，但还存在其他窗口。下图显示了 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。可以通过零填充（例如向 welch 传递 nfft=4*len(x)）或通过增加段长度（设置参数 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})\)，假设已对其进行了缩放和平移，使其均值为零且方差为单位值，则频率 \(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] 创建的稍作修改的算法计算周期图，该算法允许对单个样本进行加权，并为每个频率独立计算一个未知偏移量（也称为“浮动均值”）。

短时傅里叶变换 (STFT)

本节提供了关于使用 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。它可以表示为：

(8)\[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 的实现保持一致，将公式 (8) 重新表述为两个步骤是有意义的：

1. 提取第 \(p\) 个切片，方法是使用由 \(M\) 个采样点组成（参见 m_num）并以 \(t[p] := p \Delta t = h T\)（参见 delta_t）为中心的窗口 \(w[m]\) 进行加窗，即：

   (9)\[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）。

   (10)\[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]，其中框架是一种级数展开，可以表示给定希尔伯特空间中的任何函数。在其中，如果对于给定希尔伯特空间 \(\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]。在这种情况下，对偶窗口总是具有相同的支撑集，并且可以通过反转对角矩阵来计算。下面将勾勒出一个仅需要了解矩阵操作的粗略推导。

由于公式 (8) 中给出的 STFT 是对 \(x[k]\) 的线性映射，它可以以向量-矩阵表示法表示。这允许我们通过线性最小二乘法的形式解（如 lstsq 中）来表达逆运算，从而得出美观且简单的结果。

我们首先重新表述公式 (9) 的加窗：

(11)\[\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)\) 条次对角线上有非零条目，即：

(12)\[\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}\) 是克罗内克 Delta。公式 (10) 可以表示为：

(13)\[\vb{s}_p = \vb{F}\,\vb{x}_p \quad\text{with}\quad F[q,m] = \frac{1}{\sqrt{M}}\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,\]

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

(14)\[\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]\ .\]

由于 \(M^{-1/2}\) 的缩放因子，\(\vb{F}\) 是幺正的，即逆矩阵等于其共轭转置，这意味着 \(\conjT{\vb{F}}\vb{F} = \vb{I}\)。其他缩放比例（例如公式 (6) 中的缩放比例）也是允许的，但在本节中会使符号稍微复杂一些。

为了获得 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\)：

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

其存在前提是：

(16)\[\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}\) 进一步简化为：

(17)\[\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}\) 是幺正的）时，这一点变得很清楚。此外：

(18)\[\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\) 求和保留了该属性。这允许进一步简化公式 (15)，即：

(19)\[\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}\]

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

(20)\[\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\) 具有与公式 (12) 中 \(\vb{W}_p\) 相同的结构，即仅在第 \((ph)\) 条次对角线上有非零条目。逆矩阵中的和项可以解释为在 \(w[m]\) 上滑动 \(|w[\mu]|^2\)（并包含一个求逆过程），因此只有与 \(w[m]\) 重叠的分量才起作用。因此，所有距离边界足够远的 \(U_p[m, k]\) 都是相同的窗口。为了规避边界效应，用零填充 \(x[k]\)，从而扩大 \(\vb{U}\)，使得所有触及 \(x[k]\) 的切片都包含相同的对偶窗口：

(21)\[\begin{split} w_d[m] &= w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta\, h]\big|^2}\\ &=\ w[m] \inv{\sum_{l=0}^{M-1} \big|w[l]\big|^2 \delta_{l+m,\eta h}}\ , \quad \eta\in \IZ\ .\end{split}\]

由于对于 \(m \not\in\{0, \ldots, M-1\}\)，\(w[m] = 0\) 成立，因此只需要对满足 \(|\eta| < M/h\) 的索引 \(\eta\) 求和。第二个表达式是另一种形式，通过以 \(h\) 的增量从索引 \(l=0\) 到 \(M\) 进行求和得到。\(\delta_{l+m,\eta h}\) 是 (l+m) % h == 0 的克罗内克 Delta 符号。“对偶窗口”这一名称可以通过将公式 (14) 插入公式 (19) 来证明，即：

(22)\[ \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}\) 是可以互换的。由于 \(\vb{U}_p\) 和 \(\vb{W}_{\!p}\) 具有公式 (11) 中给出的相同结构，公式 (22) 包含了所有可能的对偶窗口的通用条件，即：

(23)\[\sum_{p=0}^{P-1} \conjT{\vb{W}_{\!p}}\,\vb{U}_p = \vb{I} \quad\Leftrightarrow\quad \sum_{p=0}^{P-1} \sum_{m=0}^{M-1} \conj{w[m]} u[m]\, \delta_{m+ph,i}\, \delta_{m+ph,j} = \delta_{i,j}\]

这可以重构为：

(24)\[\conjT{\vb{V}}\,\vb{u} = \vb{1}\ , \qquad V[i,j] := w[i]\, \delta_{i, j+ph}\ , \qquad \vb{V}\in \IC^{M\times h},\ p\in\IZ\ ,\]

其中 \(\vb{1}\in\IR^h\) 是全 1 向量。\(\vb{V}\) 只有 \(h\) 列的原因是公式 (22) 中的第 \(i\) 行和第 \((i+m)\) 行（\(i\in\IN\)）是相同的。因此只有 \(h\) 个不同的方程。

实际感兴趣的是找到最接近给定向量 \(\vb{d}\in\IC^M\) 的有效对偶窗口 \(\vb{u}_d\)。通过利用拉格朗日乘子的 \(h\) 维向量 \(\vb{\lambda}\)，我们得到凸二次规划问题：

\[\begin{split}\vb{u}_d = \min_{\vb{u}}{ \frac{1}{2}\lVert\vb{d} - \vb{u}\rVert^2 } \quad\text{w.r.t.}\quad \conjT{\vb{V}}\vb{u} = \vb{1} \qquad\Leftrightarrow\qquad \begin{bmatrix} \vb{I} & \vb{V}\\ \conjT{\vb{V}} & \vb{0} \end{bmatrix} \begin{bmatrix} \vb{u}_d \\ \vb{\lambda} \end{bmatrix} = \begin{bmatrix} \vb{d}\\ \vb{1} \end{bmatrix}\ .\end{split}\]

通过符号化反转 \(2\times2\) 分块矩阵，可以计算出闭式解，这给出：

(25)\[\begin{split}\vb{u}_d &= \underbrace{\vb{V}\inv{\conjT{\vb{V}}\vb{V}}}_{% =:\vb{Q}\in\IC^{M\times h}} \vb{1} + \big(\vb{I} - \vb{Q}\conjT{\vb{V}}\big)\vb{d} = \vb{w}_d + \vb{d} - \vb{q}_d\ ,\\ \vb{w}_d &= \vb{Q}\vb{1}\ ,\qquad \vb{q}_d = \vb{Q}\conjT{\vb{V}}\vb{d}\ ,\\ Q[i,j] &= \frac{ w[i] \delta_{i-j, \eta h} }{% \sum_{k=0}^M |w[k]|^2\,\delta_{i-k, \xi h} } ,\qquad w_d[i] = \frac{w[i] }{ \sum_{l=0}^M \conj{w[l]} w[l] \delta_{i-l, \xi h} }\ ,\\ q_d[i] &= w[i] \frac{\sum_{j=1}^h \delta_{i-j, \eta h} \sum_{l=1}^M \conj{w[l]} d[l] \delta_{j-l, \xi h} }{% \sum_{l=0}^M \conj{w[l]} w[l] \delta_{i-l, \zeta h} } = w_d[i] \sum_{l=1}^M \conj{w[l]} d[l] \delta_{i-l, \eta h}\ ,\end{split}\]

其中 \(\eta,\xi,\zeta\in\IZ\)。注意，第一项 \(\vb{w}_d\) 等于公式 (21) 中给出的解，并且 \(\conjT{\vb{V}}\vb{V}\) 的逆必须存在，否则 STFT 不可逆。当 \(\vb{d}=\vb{0}\) 时，获得解 \(\vb{w}_d\)。因此，\(\vb{w}_d\) 最小化了 \(L^2\)-范数 \(\lVert\vb{u}\rVert\)，这就是它被称为“规范对偶窗口”的原因。有时，找到方向上最接近且忽略向量长度的向量更可取。这可以通过引入比例因子 \(\alpha\in\IC\) 来最小化 \(\lVert\alpha\vb{d} - \vb{u}\rVert^2\) 来实现。由于公式 (25) 已经提供了通用解，我们可以写成：

\[\begin{split} \alpha_{\min} &= \min_{\alpha}{ \frac{1}{2}\big\lVert\alpha\vb{d} - \vb{u}_d(\alpha)\big\rVert^2 }\ ,\qquad \vb{u}_d(\alpha) = \vb{w}_d + \alpha\vb{d} - \alpha\vb{q}_d \qquad\Leftrightarrow\\ \alpha_{\min} &= \conjT{\vb{q}_d}\vb{w}_d \big/\, \conjT{\vb{q}_d}\vb{q}_d \ .\end{split}\]

窗口 \(w[m]\) 和对偶窗口 \(u[m]\) 相等的情况可以很容易地从公式 (23) 推导出来，得到 \(h\) 个如下形式的条件：

(26)\[\sum_{p=0}^{\lfloor M / h \rfloor} \big|w[m+ph]\big|^2 = 1\ , \qquad m \in \{0, \ldots, h-1\}\ .\]

注意，每个窗口采样点 \(w[m]\) 在 \(h\) 个方程中仅出现一次。为给定的窗口 \(\vb{d}\) 找到最接近的窗口 \(\vb{w}\) 很简单：根据公式 (26) 对 \(\vb{d}\) 进行分区，并将每个分区的长度归一化为单位值。在这种情况下，\(w[m]\) 也是一个规范对偶窗口，可以通过识别出在公式 (24) 中设置 \(u[m]=w[m]\) 等同于公式 (21) 中的分母为单位值来看出。

此外，如果公式 (26) 成立，公式 (16) 的矩阵 \(\vb{D}\) 就是单位矩阵，使 STFT \(\vb{G}\) 成为幺正映射，即 \(\conjT{\vb{G}}\vb{G}=\vb{I}\)。注意，这仅在使用公式 (13) 的幺正 DFT 时才成立。ShortTimeFFT 实现使用公式 (6) 的标准 DFT。因此，需要将 STFT 空间中的标量积缩放 \(1/M\)，以确保幺正映射的关键属性（即标量积相等）成立。即：

(27)\[\langle x, y\rangle = \sum_k x[k]\, \conj{y[k]} \stackrel{\stackrel{\text{unitary}}{\downarrow}}{=} \frac{1}{M}\sum_{q,p} S_x[q,p]\, \conj{S_y[q,p]}\ ,\]

其中 \(S_{x,y}\) 分别是 \(x,y\) 的 STFT。或者，窗口可以缩放 \(1/\sqrt{M}\)，对偶窗口缩放 \(\sqrt{M}\) 以获得幺正映射，这在 from_win_equals_dual 中实现。

与旧版实现的比较

函数 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 s，这使其比默认的 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]

“Fourier transform”, Wikipedia, https://en.wikipedia.org/wiki/Fourier_transform

[2]

“Spectral density”, Wikipedia, https://en.wikipedia.org/wiki/Spectral_density

[3]

“Discrete Fourier transform”, Wikipedia, https://en.wikipedia.org/wiki/Discrete_Fourier_transform

[4]

“Spectral Leakage”, Wikipedia, https://en.wikipedia.org/wiki/Spectral_leakage

[5]

“Nyquist–Shannon sampling theorem”, Wikipedia, https://en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem

[6]

N.R. Lomb “Least-squares frequency analysis of unequally spaced data”, Astrophysics and Space Science, vol 39, pp. 447-462, 1976

[7]

J.D. Scargle “Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data”, The Astrophysical Journal, vol 263, pp. 835-853, 1982

[8]

M. Zechmeister and M. Kürster, “The generalised Lomb-Scargle periodogram. A new formalism for the floating-mean and Keplerian periodograms,” Astronomy and Astrophysics, vol. 496, pp. 577-584, 2009

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

Karlheinz Gröchenig: “Foundations of Time-Frequency Analysis”, Birkhäuser Boston 2001, DOI:10.1007/978-1-4612-0003-1

[10] (1,2,3)

Ole Christensen: “An Introduction to Frames and Riesz Bases”, Birkhäuser Boston 2016, DOI:10.1007/978-3-319-25613-9