信号处理 (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\右)。\]

在具有结间距的二维中 \(\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}{\Deltay}-k\Right)。\]

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

通过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\右)。\]

使用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 )。该包还提供了一个函数( bspline )用于计算B样条基函数， \(\beta^{{o}}\left(x\right)\) 对于任意顺序和 \(x.\) 对于大型 \(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.misc.face 。该命令 sepfir2d 将具有镜像对称边界条件的可分离二维FIR过滤应用于样条系数。此函数非常适合从样条系数重建样本，并且速度比 convolve2d 它可以卷积任意的2-D滤波器，并允许选择镜像对称的边界条件。

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

>>> image = misc.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\) 使用标准NumPy数组的矩阵需要 \(68,719,476,736\) 元素。在每个元素4个字节的情况下，这将需要 \(256\textrm{{GB}}\) 记忆中的一部分。在大多数应用中，该矩阵的大多数元素为零，并且使用不同的方法来计算过滤的输出。

卷积/相关

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

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

\[Y\Left [n\right] =\sum_{k=-\infty}^{\infty}x\Left [k\right] H\Left [n-k\right] 。\]

只有当我们将序列限制为可以存储在计算机中的有限支持序列时，这个等式才能直接实现，选择 \(n=0\) 要成为这两个序列的起点，让 \(K+1\) 是其值 \(x\left[n\right]=0\) 为了所有人 \(n\geq K+1\) 和 \(M+1\) 是其值 \(h\left[n\right]=0\) 为了所有人 \(n\geq M+1\) ，则离散卷积表达式为

\[Y\Left [n\right] =\sum_{k=\max\left(n-M，0\right)}^{\min\left(n，K\右)}x\左 [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’使用其他两种方法进行力计算。

下面的代码显示了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-D数组作为输入，并将返回两个数组的N-D卷积，如下面的代码示例所示。同样的输入标志也可用于该情况。

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

相关与卷积非常相似，不同之处在于减号变成了加号。因此，

\[W\Left [n\right] =\sum_{k=-\infty}^{\infty}y\Left [k\right] X\Left [n+k\right] ，\]

是信号的(交叉)相关 \(y\) 和 \(x.\) 对于具有以下特性的有限长度信号 \(y\left[n\right]=0\) 超出范围 \(\left[0,K\right]\) 和 \(x\left[n\right]=0\) 超出范围 \(\left[0,M\right],\) 求和可以简化为

\[W\Left [n\right] =\sum_{k=\max\Left(0，-n\Right)}^{\min\Left(K，M-n\Right)}y\Left [k\right] X\Left [n+k\right] 。\]

再次假设 \(K\geq M\) ，这是

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

SciPy函数 correlate 实现此操作。此操作有等效的标志可用来返回完整的 \(K+M+1\) 长度序列(‘Full’)或与最大序列大小相同的序列，从 \(w\left[-K+\left\lfloor \frac{{M-1}}{{2}}\right\rfloor \right]\) (“相同”)或其中值取决于最小序列的所有值的序列(“有效”)。最后一个选项返回 \(K-M+1\) 值 \(w\left[M-K\right]\) 至 \(w\left[0\right]\) 包括在内。

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

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

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

>>> image = misc.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, misc
>>> import matplotlib.pyplot as plt

>>> image = misc.ascent()
>>> 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\) ，那么这种过滤就可以用卷积来实现。然而，卷积过滤序列 \(h\left[n\right]\) 可以是无穷大的，如果 \(a_{{k}}\neq0\) 为 \(k\geq1.\) 此外，这类一般的线性过滤允许将初始条件施加在 \(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\) 选择进行标准化。这个通用差分方程过滤在SCHPY中的实现比前面的方程所暗示的要稍微复杂一些。它的实施使得只需要延迟一个信号。实际实施方程式为(假设 \(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).\) 请注意， \(b_{{K}}=0\) 如果 \(K>M\) 和 \(a_{{K}}=0\) 如果 \(K>N.\) 这样，此时的输出 \(n\) 仅取决于此时的输入 \(n\) 以及它的价值 \(z_{{0}}\) 在之前的时间。只要 \(K\) 值 \(z_{{0}}\left[n-1\right]\ldots z_{{K-1}}\left[n-1\right]\) 在每个时间步长计算并存储。

使用以下命令调用差分方程过滤 lfilter 在本网站上。此命令将向量作为输入 \(b,\) 向量， \(a,\) 一个信号 \(x\) 并返回向量 \(y\) (长度与 \(x\) )使用上面给出的方程式计算。如果 \(x\) 为N-D，则沿提供的轴计算过滤。如果需要，提供以下值的初始条件 \(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] \右)。\]

利用这个公式，我们可以求出初始条件向量 \(z_{{0}}\left[-1\right]\) 至 \(z_{{K-1}}\left[-1\right]\) 在给定初始条件的情况下 \(y\) (及 \(x\) )。该命令 lfiltic 执行此功能。

例如，请考虑以下系统：

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

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

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

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

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

线性系统的分析

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

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

此替代表示可以使用scipy函数获得 tf2zpk ；反转由 zpk2tf 。

对于上面的示例，我们有

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

即，系统在以下位置具有零 \(z=-1/2\) 一根杆子放在 \(z=1/3\) 。

Scipy函数 freqz 允许计算由系数描述的系统的频率响应 \(a_k\) 和 \(b_k\) 。请参阅 freqz 函数提供了一个全面的示例。

过滤设计

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

冷杉过滤

该函数 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 (如上所述)。

国际IR过滤

SciPy提供了两个功能来直接设计IIR iirdesign 和 iirfilter 其中过滤类型(例如，椭圆)作为自变量传递，并且用于特定过滤类型的若干更多过滤设计函数，例如， ellip 。

下面的示例分别设计了具有定义的通带和阻带纹波的椭圆低通过滤。请注意，与上述示例中的FIR滤波器相比，过滤阶数(阶数为4)要低得多，以便达到相同的阻带衰减 \(\approx 60\) 分贝。

>>> 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_1z^{(M-1)}+\cdots+b_M} {a_0 z^N+a_1z^{(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+b1 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}点{\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 is the output vector of length Q

• u is the input vector of length P

• A is the state matrix, with shape (N, N)

• B is the input matrix with shape (N, P)

• C is the output matrix with shape (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。使用浮点计算时，节顺序通常并不重要；无论顺序如何，过滤输出都是相同的。

过滤转型

红外过滤设计函数首先产生一个归一化截止频率为1rad/s的模拟低通过滤原型。然后使用以下替换将其转换为其他频率和频段类型：

类型	转型
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 使用了Transform，它进行了以下替换：

\[s\right tarrow\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}\右。\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), 'xb')
>>> plt.plot(np.real(p), np.imag(p), 'or')
>>> plt.legend(['Zeros', 'Poles'], loc=2)

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

光谱分析

周期图测量

Scipy函数 periodogram 提出了一种利用周期图法估计谱密度的方法。

下面的例子计算了高斯白噪声中正弦信号的周期图。

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

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*np.sqrt(2)
>>> freq = 1270.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> x = amp*np.sin(2*np.pi*freq*time)
>>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

>>> f, Pper_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')

>>> plt.semilogy(f, Pper_spec)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD')
>>> plt.grid()
>>> plt.show()

用韦尔奇方法进行光谱分析

一种改进的方法，特别是在抗噪性方面，是由Scipy函数实现的Welch方法 welch 。

下面的示例使用Welch的方法估计频谱，并使用与上面示例相同的参数。请注意，频谱图中的噪音底板要平滑得多。

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

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*np.sqrt(2)
>>> freq = 1270.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> x = amp*np.sin(2*np.pi*freq*time)
>>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

>>> f, Pwelch_spec = signal.welch(x, fs, scaling='spectrum')

>>> plt.semilogy(f, Pwelch_spec)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD')
>>> plt.grid()
>>> plt.show()

Lomb-Scarger周期图 (lombscargle )

最小二乘谱分析(LSSA) 1 2 是一种基于正弦对数据样本的最小二乘拟合来估计频谱的方法，类似于傅立叶分析。傅立叶分析是科学上最常用的频谱分析方法，它通常会增强长间隔记录中的长周期噪声；LSSA可以缓解这种问题。

Lomb-Scarger方法对不均匀采样的数据进行频谱分析，被认为是发现和测试弱周期信号重要性的一种强有力的方法。

对于包括以下内容的时间序列 \(N_{{t}}\) 测量结果 \(X_{{j}}\equiv X(t_{{j}})\) 每隔一段时间采样 \(t_{{j}}\) ，在哪里 \((j = 1, \ldots, N_{{t}})\) ，假设其均值为零，方差为1，则得到归一化的频域Lomb-Scargl型周期图，并对其进行了定标和移位，得到了归一化的频域Lomb-Scargl型周期图 \(f\) 是

\[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而稍作修改的算法计算周期图 3, 这允许仅使用通过每个频率的输入阵列的单次通过来计算周期图。

方程被重构为：

\[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\) 在直接实施的基础上提高了速度。

解除趋势

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

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

2

J.D.斯卡格尔：“天文时间序列分析研究。II-非均匀分布数据的谱分析的统计方面”，“天体物理杂志”，第263卷，第835-853页，1982年。

3

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

上一页

傅立叶变换 (scipy.fft )

下一页

线性代数 (scipy.linalg )