信号处理 (scipy.signal)

信号处理工具箱目前包含一些滤波函数、一组有限的滤波器设计工具以及一些用于一维和二维数据的 B 样条插值算法。虽然从技术上讲，B 样条算法可以归入插值类别，但它们包含在此处是因为它们仅适用于等距数据，并且大量使用滤波器理论和传递函数形式主义来提供快速 B 样条变换。要理解本部分，你需要了解 SciPy 中的信号是实数或复数数组。

B 样条

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

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

在结点间距为 \(\Delta x\) 和 \(\Delta y\) 的二维中，函数表示为

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

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

通过 B 样条基函数表示一组样本的优点在于，连续域算子（导数、重新采样、积分等）假定数据样本是从一个底层连续函数中提取的，可以从样条系数中相对轻松地计算出来。例如，样条的二阶导数是

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

使用 B 样条的性质

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

可以看出

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

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

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

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

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

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

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

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

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

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

或者，我们可以这样做

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

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

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

滤波

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

卷积/相关

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

设\(x\left[n\right]\)定义由整数\(n.\)索引的 1-D 信号。两个 1-D 信号的完全卷积可以表示为

\[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.\)的有限序列。

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

第一个可选标志“mode”允许指定要返回的输出信号的哪一部分。默认值“full”返回整个信号。如果标志的值为“same”，则仅返回从\(y\left[\left\lfloor \frac{M-1}{2}\right\rfloor \right]\)开始的中间\(K\)个值，以便输出的长度与第一个输入相同。如果标志的值为“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”强制使用其他两种算法进行计算。

以下代码显示了 2 个序列卷积的简单示例

>>> 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]\) 和 \(\left[0,M\right],\) 之外的 \(y\left[n\right]=0\) 和 \(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]\) 是输出序列。如果我们假设初始静止状态，使得 \(y\left[n\right]=0\)，对于 \(n<0\)，则可以使用卷积实现这种类型的滤波器。但是，如果对于 \(k\geq1\)，\(a_{k}\neq0\)，则卷积滤波器序列 \(h\left[n\right]\) 可能为无穷大。此外，这一类通用的线性滤波器允许对 \(y\left[n\right]\) 设置初始条件，对于 \(n<0\)，从而导致无法使用卷积表示的滤波器。

差分方程滤波器可以被认为是递归地根据其先前值找到\(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 根据窗口方法设计滤波器。根据提供的参数，该函数返回不同的滤波器类型（例如，低通、带通……）。

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

>>> 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+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)，表示形式为

\[\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) 的单 2-D 数组，它表示一系列二阶传递函数，当串联时，会实现具有最小数值误差的高阶滤波器。每一行对应一个二阶 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-D 数组的通用中值滤波器是 medfilt。仅适用于 2-D 数组的专用版本可用作 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 宏，使 LaTeX 公式更具可读性： \newcommand{\IC}{{\mathbb{C}}} % 复数集 \newcommand{\IN}{{\mathbb{N}}} % 自然数集 \newcommand{\IR}{{\mathbb{R}}} % 实数集 \newcommand{\IZ}{{\mathbb{Z}}} % 整数集 \newcommand{\jj}{{\mathbb{j}}} % 虚数单位 \newcommand{\e}{\operatorname{e}} % 欧拉数 \newcommand{\dd}{\operatorname{d}} % 无穷小算子 \newcommand{\abs}[1]{\left|#1\right|} % 绝对值 \newcommand{\conj}[1]{\overline{#1}} % 复共轭 \newcommand{\conjT}[1]{\overline{#1^T}} % 转置复共轭 \newcommand{\inv}[1]{\left(#1\right)^{\!-1}} % 逆 % 由于未加载物理学包，因此我们自己定义宏： \newcommand{\vb}[1]{\mathbf{#1}} % 向量和矩阵为粗体 % 新宏： \newcommand{\rect}{\operatorname{rect}} % 矩形或箱形函数 \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}\ ,\]

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

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

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

同样地。例如，可以通过直接计算证明方程式(4)中的\(X(f)\)的能量为\(|a|^2\tau/2\)。因此，信号在频带\([f_a, f_b]\)中的功率可以用以下公式确定

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

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

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

正弦信号 \(x(t)\) 的频谱表示比较，公式 (1)，单位 \(\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(-f)= \conj{X}(f)\)，如果 \(x(t)\in\IR\)）。有时将负频率的值添加到其正频率对应值中。然后，幅度频谱允许读取 \(x(t)\) 在 \(f_z\) 处的完整（不是一半）幅度正弦，并且 PSD 中间隔的面积表示其完整（不是一半）功率。请注意，对于幅度谱密度，正值不会加倍，而是乘以 \(\sqrt{2}\)，因为它是 PSD 的平方根。此外，没有规范的方法来命名加倍的频谱。

以下图表显示了公式 (1) 中四个正弦信号 \(x(t)\) 的三种不同的频谱表示，具有不同的幅度持续时间 \(a\) 和持续时间 \(\tau\)。为了减少混乱，频谱以 \(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\) 毫秒，持续时间为 \(\tau=1\) 秒，采样频率为 \(f_S=100\) 赫兹。

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 箱宽的倍数，而 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 频谱的最大值也十分接近一半的预期高度。

频谱能量和频谱功率可以类似于方程式 (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}\) 可以类似地定义。

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

公式 (6) 中窗 DFT \(X^w_l\) 的谱表示比较，对于单位为 \(\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\) 处的幅度值与连续时间情况不同，这是由于从积分变为求和来确定谱能量/功率。

尽管 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\) 秒。该图显示了输入和输出信号的幅度和频率上的相位

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})\)，假设已经过缩放和移位，使得其均值为零，方差为一，则频率 \(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. 对加权的平移切片求和 \(w_d[m]\) 以重建原始信号，即

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

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

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

滑动窗口

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

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

信号末尾的行为对于具有 \(n=50\) 个样本的信号（如下面用蓝色背景表示）进行描述

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

逆 STFT 和对偶窗口

对偶窗口一词源自框架理论 [10]，其中框架是一个级数展开，可以表示给定希尔伯特空间中的任何函数。此处，展开 \(\{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]。

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

由于方程式中给出的 STFT (7) 是 \(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{其中}\quad F[q,m] =\exp\!\big\{-2\jj\pi (q + \phi_m)\, m / M\big\}\ ,\]

这允许将 \(p\) 片的 STFT 写成

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

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

为了获得 STFT 的单个向量矩阵方程式，将片叠加到一个向量中，即

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

其中 \(P\) 是结果 STFT 的列数。要反转此方程式，可以使用摩尔-彭若斯逆 \(\vb{G}^\dagger\)

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

如果

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

可逆，则存在。那么 \(\vb{x} = \vb{G}^\dagger\vb{G}\,\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}\,\conjT{\vb{G}}\vb{G}\,\vb{x}\) 显然成立。\(\vb{D}\) 始终是对角线元素为非负的对角矩阵。当进一步简化 \(\vb{D}\) 为

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

由于 \(\vb{F}\) 是酉矩阵。此外

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

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

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

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

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

这表明 \(\vb{U}_p\) 具有与方程式 (11) 中的 \(\vb{W}_p\) 相同的结构，即在 \((ph)\) 次对角线上只有非零项。逆项中的和项可以解释为在 \(w[m]\) 上滑动 \(|w[\mu]|^2\)（并进行求逆），因此只有与 \(w[m]\) 重叠的组件才有效。因此，所有远离边界的 \(U_p[m, k]\) 都是相同的窗口。为了规避边界效应，\(x[k]\) 用零填充，从而扩大 \(\vb{U}\)，使得所有与 \(x[k]\) 相交的切片都包含相同的对偶窗口

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

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

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

表明 \(\vb{U}_p\) 和 \(\vb{W}_{\!p}\) 是可以互换的。因此，\(w_d[m]\) 也是一个有效的窗口，其双窗口为 \(w[m]\)。请注意，\(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 还可以返回“角度”、“相位”、“psd”或“幅度”。旧版 spectrogram 的 scaling 行为并不直接，因为它取决于参数 mode、scaling 和 return_onesided。在 ShortTimeFFT 中，所有组合都没有直接对应关系，因为它仅提供“幅度”、“psd”或根本不提供窗口的 scaling。下表显示了这些对应关系

旧版 spectrogram				ShortTimeFFT
模式	缩放	返回单边		fft 模式	缩放
psd	密度	真		onesided2X	psd
psd	密度	假		双边	psd
幅度	频谱	真		单边	幅度
幅度	频谱	假		双边	幅度
复数	频谱	真		单边	幅度
复数	频谱	假		双边	幅度
psd	频谱	真		—	—
psd	频谱	假		—	—
复数	密度	真		—	—
复数	密度	假		—	—
幅度	密度	真		—	—
幅度	密度	假		—	—
—	—	—		*	无

在复值输入信号上使用 onesided 输出时，旧的 spectrogram 切换到 two-sided 模式。 ShortTimeFFT 引发 TypeError，因为所使用的 rfft 函数仅接受实值输入。请参阅上面的 频谱分析 部分，了解由各种参数化引起的各种频谱表示形式的讨论。

去趋势

SciPy 提供函数 detrend 以移除数据序列中的常量或线性趋势，以便查看高阶效应。

以下示例移除了二阶多项式时间序列的常量和线性趋势，并绘制了剩余信号分量。

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

>>> t = np.linspace(-10, 10, 20)
>>> y = 1 + t + 0.01*t**2
>>> yconst = signal.detrend(y, type='constant')
>>> ylin = signal.detrend(y, type='linear')

>>> plt.plot(t, y, '-rx')
>>> plt.plot(t, yconst, '-bo')
>>> plt.plot(t, ylin, '-k+')
>>> plt.grid()
>>> plt.legend(['signal', 'const. detrend', 'linear detrend'])
>>> plt.show()

参考文献

一些进一步的阅读和相关软件

[1]

“傅里叶变换”，维基百科，https://zh.wikipedia.org/wiki/%E5%82%85%E9%87%8C%E5%85%8B%E5%8F%98%E6%8D%A2

[2]

“谱密度”，维基百科，https://zh.wikipedia.org/wiki/%E7%B3%BB%E5%AF%86%E5%BC%8F

[3]

“离散傅里叶变换”，维基百科，https://zh.wikipedia.org/wiki/%E7%A6%BB%E6%95%A3%E5%82%85%E9%87%8C%E5%85%8B%E5%8F%98%E6%8D%A2

[4]

“频谱泄漏”，维基百科，https://zh.wikipedia.org/wiki/%E9%A2%91%E7%A2%B3%E6%B3%84%E6%BC%8F

[5]

“奈奎斯特-香农采样定理”，维基百科，https://zh.wikipedia.org/wiki/%E5%A5%88%E5%8D%B0%E6%96%AF%E5%9D%82-%E9%A6%99%E5%AE%89%E9%87%87%E6%A0%B7%E5%AE%9A%E7%90%86

[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