IPython-交互式计算和可视化秘籍第二版-四-IPython 交互式计算和可视化秘籍第二版（四）第九章：数值优化

IPython 交互式计算和可视化秘籍第二版（四）

原文：annas-archive.org/md5/62516af4e05f6a3b0523d8aa07fd5acb

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：数值优化

在本章中，我们将讨论以下主题：

寻找数学函数的根
最小化数学函数
使用非线性最小二乘法拟合数据
通过最小化潜在能量找到物理系统的平衡状态

引言

数学优化是应用数学的一个广泛领域，它涉及到寻找给定问题的最佳解。许多现实世界的问题可以用优化框架来表示。比如，从 A 点到 B 点的最短路径是什么？解决一个难题的最佳策略是什么？汽车的最节能形状是什么（汽车空气动力学）？数学优化在许多领域都有应用，包括工程学、经济学、金融学、运筹学、图像处理、数据分析等。

从数学角度看，优化问题通常包括找到一个函数的最大值或最小值。根据函数变量是实值的还是离散的，我们有时会使用连续优化或离散优化这两个术语。

在本章中，我们将重点讨论解决连续优化问题的数值方法。许多优化算法都在scipy.optimize模块中实现。我们将在本书的其他几章中遇到其他类型的优化问题。例如，我们将在第十四章中看到离散优化问题，图形、几何学与地理信息系统。在本引言中，我们将介绍一些与数学优化相关的重要定义和关键概念。

目标函数

我们将研究找到实值函数f的根或极值的方法，称为目标函数。极值可以是函数的最大值或最小值。这个数学函数通常在 Python 函数中实现。它可以接受一个或多个变量，可以是连续的或不连续的，等等。我们对函数的假设越多，优化起来就越容易。

注意

f的最大值即为*-f的最小值，因此任何最小化算法都可以通过考虑该函数的对立面来实现函数的最大化。因此，从现在开始，当我们谈论最小化时，实际上指的是最小化或最大化*。

凸函数通常比非凸函数更容易优化，因为它们满足某些有用的性质。例如，任何局部最小值必然是全局最小值。凸优化领域处理的是专门用于优化凸函数在凸领域上算法的研究。凸优化是一个高级主题，我们在这里不能深入讨论。

可微函数 具有梯度，这些梯度在优化算法中尤为有用。同样，连续函数 通常比非连续函数更容易优化。

此外，单变量函数比多变量函数更容易优化。

选择最合适的优化算法取决于目标函数所满足的属性。

局部和全局最小值

函数 f 的最小值是一个点 x[0]，满足 f(x) 局部和全局最小值 f(x[0]**)，对于 E 中的特定点集 x。当这个不等式在整个 E 集合上满足时，我们称 x[0] 为 全局最小值。当仅在局部（围绕点 x[0]）满足时，我们称 x[0] 为 局部最小值。最大值 的定义类似。

如果 f 可微，则极值 x[0] 满足：

局部和全局最小值

因此，寻找目标函数的极值与寻找导数的根密切相关。然而，满足此性质的点 x[0] 不一定是极值。

找到全局最小值比找到局部最小值更难。通常，当一个算法找到局部最小值时，并不能保证它也是全局最小值。经常有算法在寻找全局最小值时会被卡住在局部最小值。这一问题需要特别在全局最小化算法中考虑。然而，对于凸函数，情况较为简单，因为这些函数没有严格的局部最小值。此外，很多情况下找到局部最小值已经足够好（例如，在寻找一个问题的良好解决方案，而不是绝对最优解时）。最后，需要注意的是，全局最小值或最大值不一定存在（函数可能趋向无穷大）。在这种情况下，可能需要约束搜索空间；这就是 约束优化 的主题。

局部和全局最小值

局部和全局极值

约束与无约束优化

无约束优化：在函数 f 定义的整个集合 E 上寻找最小值
约束优化：在 E 的子集 E' 上寻找函数 f 的最小值；该集合通常通过等式和不等式来描述：

这里，g[i] 和 h[j] 是定义约束的任意函数。

例如，优化汽车的空气动力学形状需要对汽车的体积、质量、生产过程的成本等参数进行约束。

约束优化通常比无约束优化更难。

确定性算法和随机算法

一些全局优化算法是 确定性的，而其他则是 随机的。通常，确定性方法适用于表现良好的函数，而随机方法则可能适用于高度不规则和嘈杂的函数。

原因在于确定性算法可能会陷入局部最小值，尤其是在存在多个非全局局部最小值的情况下。通过花时间探索空间 E，随机算法有可能找到全局最小值。

参考文献

SciPy 讲义是一个关于使用 SciPy 进行数学优化的极好参考，访问 scipy-lectures.github.io/advanced/mathematical_optimization/index.html
scipy.optimize 的参考手册，访问 docs.scipy.org/doc/scipy/reference/optimize.html
维基百科上的数学优化概述，访问 en.wikipedia.org/wiki/Mathematical_optimization
维基百科上的极值、最小值和最大值，访问 en.wikipedia.org/wiki/Maxima_and_minima
维基百科上的凸优化，访问 en.wikipedia.org/wiki/Convex_optimization
Gabriel Peyré 提供的图像处理高级优化方法，访问 github.com/gpeyre/numerical-tours

寻找数学函数的根

在这个简短的教程中，我们将展示如何使用 SciPy 寻找单一实数变量的简单数学函数的根。

如何做…

让我们导入 NumPy、SciPy、scipy.optimize 和 matplotlib：

In [1]: import numpy as np
        import scipy as sp
        import scipy.optimize as opt
        import matplotlib.pyplot as plt
        %matplotlib inline

我们在 Python 中定义数学函数 f(x)=cos(x)-x，并尝试通过数值方法寻找该函数的根。这里，根对应于余弦函数的固定点：
```
In [2]: f = lambda x: np.cos(x) - x
```

让我们在区间 [-5, 5] 上绘制该函数（使用 1000 个样本）：

In [3]: x = np.linspace(-5, 5, 1000)
        y = f(x)
        plt.plot(x, y)
        plt.axhline(0, color='k')
        plt.xlim(-5,5)

如何做…

我们看到该函数在此区间内有唯一根（这是因为函数在该区间内的符号发生了变化）。scipy.optimize 模块包含几个根寻找函数，这里进行了相应的适配。例如，bisect() 函数实现了 二分法（也称为 二分法法）。它的输入是要寻找根的函数和区间：
```
In [4]: opt.bisect(f, -5, 5)
Out[4]: 0.7390851332155535
```
让我们在图表上可视化根的位置：
```
In [5]: plt.plot(x, y)
        plt.axhline(0, color='k')
        plt.scatter([_], [0], c='r', s=100)
        plt.xlim(-5,5)
```
更快且更强大的方法是brentq()（布伦特法）。该算法同样要求 f 为连续函数，并且 f(a) 与 f(b) 具有不同的符号：
```
In [6]: opt.brentq(f, -5, 5)
Out[6]: 0.7390851332151607
```
brentq() 方法比 bisect() 更快。如果条件满足，首先尝试布伦特法是一个好主意：
```
In [7]: %timeit opt.bisect(f, -5, 5)
        %timeit opt.brentq(f, -5, 5)
1000 loops, best of 3: 331 µs per loop
10000 loops, best of 3: 71 µs per loop
```

它是如何工作的…

二分法通过反复将区间一分为二，选择一个必定包含根的子区间来进行。该方法基于这样一个事实：如果 f 是一个单一实变量的连续函数，且 f(a)>0 且 f(b)<0，则 f 在 (a,b) 区间内必有根（中值定理）。

Brent 方法 是一种流行的混合算法，结合了根的括起来、区间二分法和反向二次插值。它是一个默认方法，在许多情况下都能工作。

让我们也提一下牛顿法。其基本思想是通过切线近似 f(x)（由 f'(x) 求得），然后找到与 y=0 线的交点。如果 f 足够规则，那么交点会更接近 f 的实际根。通过反复执行此操作，算法通常会收敛到所寻求的解。

还有更多……

下面是一些参考文献：

scipy.optimize 的文档，地址为 docs.scipy.org/doc/scipy/reference/optimize.html#root-finding
一个关于 SciPy 根查找的课程，地址为 quant-econ.net/scipy.html#roots-and-fixed-points
二分法的维基百科页面，地址为 en.wikipedia.org/wiki/Bisection_method
中值定理的维基百科页面，地址为 en.wikipedia.org/wiki/Intermediate_value_theorem
Brent 方法的维基百科页面，地址为 en.wikipedia.org/wiki/Brent%27s_method
牛顿法的维基百科页面，地址为 en.wikipedia.org/wiki/Newton%27s_method

另见

最小化数学函数 的教程

最小化数学函数

数学优化主要涉及寻找数学函数的最小值或最大值的问题。现实世界中的许多数值问题可以表达为函数最小化问题。这类问题可以在统计推断、机器学习、图论等领域中找到。

尽管有许多函数最小化算法，但并没有一个通用且普适的方法。因此，理解现有算法类别之间的差异、它们的特点以及各自的使用场景非常重要。我们还应该对问题和目标函数有清晰的了解；它是连续的、可微的、凸的、多维的、规则的，还是有噪声的？我们的优化问题是约束的还是无约束的？我们是在寻找局部最小值还是全局最小值？

在本教程中，我们将演示在 SciPy 中实现的几种函数最小化算法的使用示例。

如何实现……

我们导入库：

In [1]: import numpy as np
        import scipy as sp
        import scipy.optimize as opt
        import matplotlib.pyplot as plt
        %matplotlib inline

首先，我们定义一个简单的数学函数（基准正弦的反函数）。这个函数有许多局部最小值，但只有一个全局最小值(en.wikipedia.org/wiki/Sinc_function)：
```
In [2]: f = lambda x: 1-np.sin(x)/x
```

让我们在区间*[-20, 20]*上绘制这个函数（使用 1000 个样本）：

In [3]: x = np.linspace(-20., 20., 1000)
        y = f(x)
In [4]: plt.plot(x, y)

如何做…

scipy.optimize模块包含许多函数最小化的例程。minimize()函数提供了一个统一的接口，适用于多种算法。Broyden–Fletcher–Goldfarb–Shanno（BFGS）算法（minimize()中的默认算法）通常能给出良好的结果。minimize()函数需要一个初始点作为参数。对于标量一元函数，我们还可以使用minimize_scalar()：
```
In [5]: x0 = 3
        xmin = opt.minimize(f, x0).x
```
从x[0]**=3开始，算法能够找到实际的全局最小值，如下图所示：
```
In [6]: plt.plot(x, y)
        plt.scatter(x0, f(x0), marker='o', s=300)
        plt.scatter(xmin, f(xmin), marker='v', s=300)
        plt.xlim(-20, 20)
```

现在，如果我们从一个更远离实际全局最小值的初始点开始，算法只会收敛到一个局部最小值：

In [7]: x0 = 10
        xmin = opt.minimize(f, x0).x
In [8]: plt.plot(x, y)
        plt.scatter(x0, f(x0), marker='o', s=300)
        plt.scatter(xmin, f(xmin), marker='v', s=300)
        plt.xlim(-20, 20)

如何做…

像大多数函数最小化算法一样，BFGS 算法在找到局部最小值时效率很高，但不一定能找到全局最小值，尤其是在复杂或嘈杂的目标函数上。克服这个问题的一个通用策略是将此类算法与初始点的探索性网格搜索相结合。另一个选择是使用基于启发式和随机方法的不同类型算法。一个流行的例子是模拟退火方法：
```
In [9]: xmin = opt.minimize(f, x0, method='Anneal').x
In [10]: plt.plot(x, y)
         plt.scatter(x0, f(x0), marker='o', s=300)
         plt.scatter(xmin, f(xmin), marker='v', s=300)
         plt.xlim(-20, 20)
```
这次，算法成功找到了全局最小值。

现在，让我们定义一个新的函数，这次是二维的，称为Lévi 函数：如何做…

这个函数非常不规则，通常可能难以最小化。它是许多优化测试函数之一，研究人员为研究和基准测试优化算法而开发的(en.wikipedia.org/wiki/Test_functions_for_optimization)：

In [11]: def g(X):
             # X is a 2*N matrix, each column contains
             # x and y coordinates.
             x, y = X
             return (np.sin(3*np.pi*x)**2 +
                     (x-1)**2 * (1+np.sin(3*np.pi*y)**2) +
                     (y-1)**2 * (1+np.sin(2*np.pi*y)**2))

让我们使用imshow()在正方形区域*[-10,10]²*上显示这个函数：

In [12]: n = 200
         k = 10
         X, Y = np.mgrid[-k:k:n*1j,-k:k:n*1j]
In [13]: Z = g(np.vstack((X.ravel(),
                          Y.ravel()))).reshape(n,n)
In [14]: # We use a logarithmic scale for the color here.
         plt.imshow(np.log(Z), cmap=plt.cm.hot_r)
         plt.xticks([]); plt.yticks([])

如何做…

BFGS 算法也适用于多维：

In [15]: x0, y0 = opt.minimize(g, (8, 3)).x
In [16]: plt.imshow(np.log(Z), cmap=plt.cm.hot_r,
                    extent=(-k, k, -k, k), origin=0)
         plt.scatter([x0], [y0], c=['r'], s=100)
         plt.xticks([]); plt.yticks([])

如何做…

它是如何工作的…

许多函数最小化算法基于梯度下降的基本思想。如果一个函数f是可微的，那么在每个点，其梯度的相反方向指向函数下降速率最大的方向。沿着这个方向前进，我们可以预期找到一个局部最小值。

这个操作通常是通过迭代进行的，沿着梯度方向以小步长进行。这个步长的计算方法取决于优化方法。

牛顿法也可以在函数最小化的上下文中使用。其思路是利用牛顿法寻找 f' 的根，从而利用二阶导数 f''。换句话说，我们用一个二次函数来逼近 f，而不是用线性函数。在多维情况下，通过计算 f 的 Hessian（二阶导数），我们可以实现这一过程。通过迭代执行此操作，我们期望算法能够收敛到局部最小值。

当计算 Hessian 矩阵的代价过高时，我们可以计算 Hessian 的近似值。此类方法称为 准牛顿法。BFGS 算法属于这一类算法。

这些算法利用目标函数的梯度。如果我们能够计算梯度的解析表达式，应当将其提供给最小化程序。否则，算法将计算一个近似的梯度，可能不可靠。

模拟退火算法是一种通用的概率元启发式算法，用于全局优化问题。它基于热力学系统的类比：通过升高和降低温度，系统的配置可能会收敛到一个低能量状态。

有许多基于元启发式算法的随机全局优化方法。它们通常没有之前描述的确定性优化算法那样理论扎实，并且并不总是能保证收敛。然而，在目标函数非常不规则且含有许多局部最小值的情况下，这些方法可能会非常有用。协方差矩阵适应进化策略（CMA-ES）算法是一种在许多情况下表现良好的元启发式算法。它目前在 SciPy 中没有实现，但在稍后给出的参考文献中有 Python 实现。

SciPy 的 minimize() 函数接受一个 method 关键字参数，用于指定要使用的最小化算法。该函数返回一个包含优化结果的对象。x 属性是达到最小值的点。

还有更多内容……

下面是一些进一步的参考资料：

scipy.optimize 参考文档，链接地址：docs.scipy.org/doc/scipy/reference/optimize.html
一堂关于使用 SciPy 进行数学优化的精彩讲座，链接地址：scipy-lectures.github.io/advanced/mathematical_optimization/
维基百科上关于梯度的定义，链接地址：en.wikipedia.org/wiki/Gradient
维基百科上的牛顿法，链接地址：en.wikipedia.org/wiki/Newton%27s_method_in_optimization
维基百科上的准牛顿法，链接地址：en.wikipedia.org/wiki/Quasi-Newton_method
维基百科上的函数最小化元启发式方法，链接：en.wikipedia.org/wiki/Metaheuristic
维基百科上的模拟退火，链接：en.wikipedia.org/wiki/Simulated_annealing
维基百科上的 CMA-ES 算法描述，链接：en.wikipedia.org/wiki/CMA-ES
可在 www.lri.fr/~hansen/cmaes_inmatlab.html#python 获取 CMA-ES 的 Python 实现

参见其他资料

求解数学函数的根教程

用非线性最小二乘法拟合数据函数

在这个教程中，我们将展示数值优化应用于非线性最小二乘曲线拟合的一个例子。目标是根据多个参数拟合一个函数到数据点上。与线性最小二乘法不同，这个函数在这些参数上不必是线性的。

我们将用人工数据来演示这种方法。

如何实现…

让我们导入常用的库：

In [1]: import numpy as np
        import scipy.optimize as opt
        import matplotlib.pyplot as plt
        %matplotlib inline

我们定义了一个具有四个参数的逻辑斯蒂函数：如何实现…

In [2]: def f(x, a, b, c, d):
            return a/(1 + np.exp(-c * (x-d))) + b

让我们定义四个随机参数：

In [3]: a, c = np.random.exponential(size=2)
        b, d = np.random.randn(2)

现在，我们通过使用 sigmoid 函数并添加一点噪声来生成随机数据点：

In [4]: n = 100
        x = np.linspace(-10., 10., n)
        y_model = f(x, a, b, c, d)
        y = y_model + a * .2 * np.random.randn(n)

这里是数据点的图示，图中显示了用于生成数据的特定 sigmoid（用虚线黑色表示）：
```
In [5]: plt.plot(x, y_model, '--k')
        plt.plot(x, y, 'o')
```
我们现在假设我们只能访问数据点，而无法访问底层的生成函数。这些点可能是在实验中获得的。从数据来看，这些点似乎大致符合一个 sigmoid 曲线，因此我们可能希望尝试将这样的曲线拟合到这些点上。这就是曲线拟合的含义。SciPy 的 curve_fit() 函数允许我们将由任意 Python 函数定义的曲线拟合到数据上：
```
In [6]: (a_, b_, c_, d_), _ = opt.curve_fit(f, x, y,
                                            (a, b, c, d))
```
现在，让我们来看一下拟合后的 sigmoid 曲线：
```
In [7]: y_fit = f(x, a_, b_, c_, d_)
In [8]: plt.plot(x, y_model, '--k')
        plt.plot(x, y, 'o')
        plt.plot(x, y_fit, '-')
```
拟合后的 sigmoid 曲线似乎与用于数据生成的原始 sigmoid 曲线非常接近。

它是如何工作的…

在 SciPy 中，非线性最小二乘曲线拟合是通过最小化以下代价函数来实现的：

它是如何工作的…

这里，它是如何工作的… 是参数的向量（在我们的示例中， =(a,b,c,d)）。

非线性最小二乘法实际上与线性回归中的线性最小二乘法非常相似。在线性最小二乘法中，函数 f 在参数上是线性的，而在这里则不是线性的。因此，S( 它是如何工作的… ) 的最小化不能通过解析地解出 S 对的导数来完成。SciPy 实现了一种称为 Levenberg-Marquardt 算法 的迭代方法（高斯-牛顿算法的扩展）。

还有更多…

这里是更多的参考资料：

curvefit的参考文档，访问链接 docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
Wikipedia 上的非线性最小二乘法，访问链接 en.wikipedia.org/wiki/Non-linear_least_squares
Wikipedia 上的 Levenberg-Marquardt 算法，访问链接 en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm

另见

最小化数学函数的食谱

通过最小化潜在能量找到物理系统的平衡状态

在这个例子中，我们将给出前面描述的函数最小化算法的应用实例。我们将尝试通过最小化物理系统的潜在能量来数值地寻找其平衡状态。

更具体地说，我们将考虑一个由质量和弹簧构成的结构，弹簧固定在垂直墙面上，并受到重力作用。从初始位置开始，我们将寻找重力和弹性力相互平衡的平衡配置。

如何操作…

让我们导入 NumPy、SciPy 和 matplotlib：

In [1]: import numpy as np
        import scipy.optimize as opt
        import matplotlib.pyplot as plt
        %matplotlib inline

我们在国际单位制中定义一些常数：

In [2]: g = 9.81  # gravity of Earth
        m = .1  # mass, in kg
        n = 20  # number of masses
        e = .1  # initial distance between the masses
        l = e  # relaxed length of the springs
        k = 10000  # spring stiffness

我们定义质量的初始位置。它们排列在一个二维网格上，具有两行和n/2列：

In [3]: P0 = np.zeros((n, 2))
        P0[:,0] = np.repeat(e*np.arange(n//2), 2)
        P0[:,1] = np.tile((0,-e), n//2)

现在，让我们定义质量之间的连接矩阵。系数*(i,j)为 1 表示质量i和j*之间通过弹簧连接，若没有连接则为 0：
```
In [4]: A = np.eye(n, n, 1) + np.eye(n, n, 2)
```

我们还指定了每个弹簧的刚度。它是l，除了对角线上的弹簧，其刚度为 How to do it… ：

In [5]: L = l * (np.eye(n, n, 1) + np.eye(n, n, 2))
        for i in range(n//2-1):
            L[2*i+1,2*i+2] *= np.sqrt(2)

我们获取弹簧连接的索引：
```
In [6]: I, J = np.nonzero(A)
```

这个dist函数计算距离矩阵（任何一对质量之间的距离）：

In [7]: dist = lambda P: np.sqrt(
                    (P[:,0]-P[:,0][:, np.newaxis])**2 + 
                    (P[:,1]-P[:,1][:, np.newaxis])**2)

我们定义一个函数来显示系统。弹簧根据其张力着色：

In [8]: def show_bar(P):
            # Wall.
            plt.axvline(0, color='k', lw=3)
            # Distance matrix.
            D = dist(P)
            # We plot the springs.
            for i, j in zip(I, J):
                # The color depends on the spring tension,
                # which is proportional to the spring 
                # elongation.
                c = D[i,j] - L[i,j]
                plt.plot(P[[i,j],0], P[[i,j],1], 
                         lw=2, color=plt.cm.copper(c*150))
            # We plot the masses.
            plt.plot(P[[I,J],0], P[[I,J],1], 'ok',)
            # We configure the axes.
            plt.axis('equal')
            plt.xlim(P[:,0].min()-e/2, P[:,0].max()+e/2)
            plt.ylim(P[:,1].min()-e/2, P[:,1].max()+e/2)
            plt.xticks([]); plt.yticks([])

这里是系统的初始配置：

In [9]: show_bar(P0)
        plt.title("Initial configuration")

How to do it…

要找到平衡状态，我们需要最小化系统的总潜在能量。以下函数根据质量的位置计算系统的能量。此函数在*如何工作…*部分中进行了说明：

In [10]: def energy(P):
             # The argument P is a vector (flattened 
             # matrix). We convert it to a matrix here.
             P = P.reshape((-1, 2))
             # We compute the distance matrix.
             D = dist(P)
             # The potential energy is the sum of the
             # gravitational and elastic potential 
             # energies.
             return (g * m * P[:,1].sum() + 
                     .5 * (k * A * (D - L)**2).sum())

让我们计算初始配置的潜在能量：

In [11]: energy(P0.ravel())
Out[11]: -0.98099

现在，让我们使用函数最小化方法来最小化潜在能量。我们需要一个约束优化算法，因为我们假设前两个质量被固定在墙壁上。因此，它们的位置不能改变。L-BFGS-B算法是 BFGS 算法的一个变种，支持边界约束。在这里，我们强制前两个点保持在初始位置，而其他点没有约束。minimize()函数接受一个包含每个维度的[min, max]值对的bounds列表：
```
In [12]: bounds = np.c_[P0[:2,:].ravel(),
                        P0[:2,:].ravel()].tolist() + \
                              [[None, None]] * (2*(n-2))
In [13]: P1 = opt.minimize(energy, P0.ravel(),
                           method='L-BFGS-B',
                           bounds=bounds).x \
                                     .reshape((-1, 2))
```
让我们显示稳定的配置：
```
In [14]: show_bar(P1)
         plt.title("Equilibrium configuration")
```
这个配置看起来很逼真。张力似乎在靠近墙壁的顶部弹簧上达到了最大值。

它是如何工作的……

这个例子在概念上很简单。系统的状态仅由质量的位置描述。如果我们能写出一个返回系统总能量的 Python 函数，那么找到平衡状态就只需要最小化这个函数。这就是最小总势能原理，源于热力学第二定律。

这里，我们给出了系统总能量的表达式。由于我们只关心平衡状态，我们省略了任何动能方面的内容，只考虑由重力（重力作用）和弹簧力（弹性势能）引起的势能。

令 U 为系统的总势能，U 可以表示为质量的重力势能与弹簧的弹性势能之和。因此：

它是如何工作的……

这里：

m 是质量
g 是地球的重力
k 是弹簧的刚度
p[i] = (x[i], y[i]) 是质量 i 的位置
a[ij] 如果质量 i 和 j 通过弹簧连接，则为 1，否则为 0
l[ij] 是弹簧 (i,j) 的放松长度，如果质量 i 和 j 没有连接，则为 0

energy() 函数使用 NumPy 数组上的向量化计算来实现这个公式。

还有更多……

以下参考资料包含关于这个公式背后物理学的详细信息：

维基百科上的势能，详情请见 en.wikipedia.org/wiki/Potential_energy
维基百科上的弹性势能，详情请见 en.wikipedia.org/wiki/Elastic_potential_energy
胡克定律是弹簧响应的线性近似，详情请见 en.wikipedia.org/wiki/Hooke%27s_law
维基百科上的最小能量原理，详情请见 en.wikipedia.org/wiki/Minimum_total_potential_energy_principle

这是关于优化算法的参考资料：

维基百科上的 L-BFGS-B 算法，详情请见 en.wikipedia.org/wiki/Limited-memory_BFGS#L-BFGS-B

另见

最小化数学函数的操作步骤

第十章：信号处理

在本章中，我们将涵盖以下主题：

使用快速傅里叶变换分析信号的频率成分
应用线性滤波器于数字信号
计算时间序列的自相关

引言

信号是描述一个量在时间或空间上变化的数学函数。依赖于时间的信号通常被称为时间序列。时间序列的例子包括股价，它们通常呈现为在均匀时间间隔内的连续时间点。在物理学或生物学中，实验设备记录如电磁波或生物过程等变量的演变。

在信号处理中，一般目标是从原始的、有噪声的测量数据中提取有意义和相关的信息。信号处理的主题包括信号获取、变换、压缩、滤波和特征提取等。当处理复杂数据集时，先清理数据可能对应用更先进的数学分析方法（例如机器学习）有帮助。

在本简短的章节中，我们将阐明并解释信号处理的主要基础。在下一章，第十一章中，我们将看到针对图像和声音的特定信号处理方法。

首先，在本引言中我们将给出一些重要的定义。

模拟信号与数字信号

信号可以是时间依赖的或空间依赖的。在本章中，我们将专注于时间依赖的信号。

设*x(t)*为时变信号。我们可以说：

如果t是连续变量，且*x(t)*是实数，则该信号为模拟信号
如果t是离散变量（离散时间信号），且*x(t)*只能取有限个值（量化信号），则该信号为数字信号

下图展示了模拟信号（连续曲线）与数字信号（点）的区别：

模拟和数字信号

模拟信号与数字（量化）信号之间的区别

模拟信号存在于数学中以及大多数物理系统中，如电路。然而，由于计算机是离散的机器，它们只能理解数字信号。这就是计算科学特别处理数字信号的原因。

由实验设备记录的数字信号通常由两个重要的量来描述：

采样率：每秒记录的值（或样本）数量（以赫兹为单位）
分辨率：量化的精度，通常以每个样本的位数（也称为位深度）表示

具有高采样率和位深度的数字信号更为精确，但它们需要更多的内存和处理能力。这两个参数受到记录信号的实验设备的限制。

奈奎斯特–香农采样定理

假设我们考虑一个连续（模拟）时变信号x(t)。我们使用实验设备记录这个物理信号，并以采样率*f[s]*获取数字信号。由于原始模拟信号具有无限精度，而记录的信号具有有限精度，我们预计在模拟到数字的过程中会丢失信息。

奈奎斯特-香农采样定理指出，在某些条件下，模拟信号和采样率可以确保采样过程中不丢失任何信息。换句话说，在这些条件下，我们可以从采样的数字信号中恢复原始的连续信号。更多细节，请参见en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem。

让我们定义这些条件。傅里叶变换 奈奎斯特-香农采样定理是通过以下公式定义的：

奈奎斯特-香农采样定理

在这里，傅里叶变换是时间依赖信号在频域中的表示。奈奎斯特准则指出：

奈奎斯特-香农采样定理

换句话说，信号必须是带限的，这意味着它不能包含任何高于某个截止频率B的频率。此外，采样率f[s]需要至少是该频率B的两倍。以下是几个定义：

奈奎斯特率是2B。对于给定的带限模拟信号，它是无损采样该信号所需的最小采样率。
奈奎斯特频率是f[s]/2。对于给定的采样率，它是信号可以包含的最大频率，以便无损采样。

在这些条件下，我们理论上可以从采样的数字信号中重建原始的模拟信号。

压缩感知

压缩感知是一种现代且重要的信号处理方法。它承认许多现实世界的信号本质上是低维的。例如，语音信号具有非常特定的结构，取决于人类发音道的普遍物理约束。

即使一个语音信号在傅里叶域中有许多频率，它也可以通过在合适的基（字典）上进行稀疏分解来很好地逼近。根据定义，稀疏分解是指大多数系数为零。如果字典选择得当，每个信号都是少数基信号的组合。

这本词典包含了特定于给定问题中考虑的信号的基本信号。这与将信号分解为通用正弦函数基础的傅里叶变换不同。换句话说，通过稀疏表示，可以规避奈奎斯特条件。我们可以从包含比奈奎斯特条件要求的样本更少的稀疏表示中精确重建连续信号。

稀疏分解可以通过复杂的算法找到。特别是，这些问题可以转化为凸优化问题，可以用特定的数值优化方法解决。

压缩感知在信号压缩、图像处理、计算机视觉、生物医学成像以及许多其他科学和工程领域中有许多应用。

关于压缩感知的更多参考资料：

参考文献

这里有几个参考资料：

理解数字信号处理，作者 Richard G. Lyons，出版社 Pearson Education，（2010）。
要详细了解压缩感知，请参阅《信号处理的小波之旅：稀疏方法》一书，作者 Mallat Stéphane，出版社 Academic Press，（2008）。
《Python 信号处理》一书由 Jose Unpingco 编写，比我们在本章中涵盖的内容要详细得多。代码可以在 GitHub 上的 IPython 笔记本上找到（python-for-signal-processing.blogspot.com）。
在 WikiBooks 上提供的数字信号处理可参考数字信号处理。

用快速傅里叶变换分析信号的频率成分

在这个案例中，我们将展示如何使用快速傅里叶变换（FFT）计算信号的频谱密度。频谱表示与频率相关联的能量（编码信号中的周期波动）。它是通过傅里叶变换获得的，这是信号的频率表示。信号可以在两种表示之间来回转换而不丢失信息。

在这个案例中，我们将阐述傅里叶变换的几个方面。我们将应用这一工具于从美国国家气候数据中心获取的法国 20 年的天气数据。

准备工作

从书籍的 GitHub 存储库中下载Weather数据集（github.com/ipython-books/cookbook-data），并将其提取到当前目录中。

数据来自www.ncdc.noaa.gov/cdo-web/dat…。

如何做……

让我们导入包，包括scipy.fftpack，其中包括许多与 FFT 相关的例程：

In [1]: import datetime
        import numpy as np
        import scipy as sp
        import scipy.fftpack
        import pandas as pd
        import matplotlib.pyplot as plt
        %matplotlib inline

我们从 CSV 文件中导入数据。数字-9999用于表示 N/A 值，pandas 可以轻松处理。此外，我们告诉 pandas 解析DATE列中的日期：

In [2]: df0 = pd.read_csv('data/weather.csv', 
                          na_values=(-9999), 
                          parse_dates=['DATE'])
In [3]: df = df0[df0['DATE']>='19940101']
In [4]: df.head()
Out[4]:  STATION                DATE  PRCP  TMAX  TMIN
365  FR013055001 1994-01-01 00:00:00     0   104    72
366  FR013055001 1994-01-02 00:00:00     4   128    49

每一行包含法国一个气象站每天记录的降水量和极端温度。对于日历中的每个日期，我们想得到整个国家的单一平均温度。pandas 提供的groupby()方法可以轻松实现这一点。我们还使用dropna()去除任何 N/A 值：
```
In [5]: df_avg = df.dropna().groupby('DATE').mean()
In [6]: df_avg.head()
Out[6]:
      DATE        PRCP        TMAX       TMIN
1994-01-01  178.666667  127.388889  70.333333
1994-01-02  122.000000  152.421053  81.736842
```
现在，我们得到日期列表和相应的温度列表。单位是十分之一度，我们计算的是最小和最大温度之间的平均值，这也解释了为什么我们要除以 20。
```
In [7]: date = df_avg.index.to_datetime()
        temp = (df_avg['TMAX'] + df_avg['TMIN']) / 20.
        N = len(temp)
```

让我们来看一下温度的变化：

In [8]: plt.plot_date(date, temp, '-', lw=.5)
        plt.ylim(-10, 40)
        plt.xlabel('Date')
        plt.ylabel('Mean temperature')

如何操作...

现在，我们计算信号的傅里叶变换和谱密度。第一步是使用fft()函数计算信号的 FFT：
```
In [9]: temp_fft = sp.fftpack.fft(temp)
```
一旦得到 FFT，我们需要取其绝对值的平方来获得功率谱密度（PSD）：
```
In [10]: temp_psd = np.abs(temp_fft) ** 2
```
下一步是获取与 PSD 值对应的频率。fftfreq()工具函数正是用来做这个的。它以 PSD 向量的长度和频率单位作为输入。在这里，我们选择年度单位：频率为 1 对应 1 年（365 天）。我们提供1/365，因为原始单位是天数。
```
In [11]: fftfreq = sp.fftpack.fftfreq(len(temp_psd), 
                                      1./365)
```
fftfreq()函数返回正频率和负频率。这里我们只关心正频率，因为我们有一个真实信号（这将在本配方的*它是如何工作的...*部分中解释）。
```
In [12]: i = fftfreq>0
```
现在我们绘制信号的功率谱密度，作为频率的函数（单位为1 年）。我们选择对y轴使用对数尺度（分贝）。
```
In [13]: plt.plot(fftfreq[i], 10*np.log10(temp_psd[i]))
         plt.xlim(0, 5)
         plt.xlabel('Frequency (1/year)')
         plt.ylabel('PSD (dB)')
```
因为信号的基频是温度的年度变化，我们观察到在f=1时有一个峰值。

现在，我们去除高于基频的频率：

In [14]: temp_fft_bis = temp_fft.copy()
         temp_fft_bis[np.abs(fftfreq) > 1.1] = 0

接下来，我们执行逆 FFT将修改后的傅里叶变换转换回时间域。这样，我们恢复了一个主要包含基频的信号，如下图所示：

In [15]: temp_slow = np.real(sp.fftpack.ifft(temp_fft_bis))
In [16]: plt.plot_date(date, temp, '-', lw=.5)
         plt.plot_date(date, temp_slow, '-')
         plt.xlim(datetime.date(1994, 1, 1),
                  datetime.date(2000, 1, 1))
         plt.ylim(-10, 40)
         plt.xlabel('Date')
         plt.ylabel('Mean temperature')

如何操作...

我们得到信号的平滑版本，因为当我们去除傅里叶变换中的高频时，快速变化的部分已经丢失。

它是如何工作的...

广义而言，傅立叶变换是信号的周期分量的超级位置的替代表示。这是一个重要的数学结果，任何良好行为的函数都可以用这种形式表示。而时间变化的信号最自然的考虑方式是作为时间的函数，傅立叶变换则将其表示为频率的函数。每个频率都与一个大小和相位相关联，这两者都被编码为单一复数。

离散傅立叶变换

让我们考虑一个由向量*(x[0], ..., x[(N-1)])表示的数字信号x*。我们假设这个信号是定期采样的。x的离散傅立叶变换（DFT）被定义为X = (X[0], ..., X[(N-1)])：

离散傅立叶变换

FFT 可以高效地计算离散傅立叶变换（DFT），这是一种利用该定义中的对称性和冗余性以显著加快计算速度的算法。FFT 的复杂度是O(N log N)，而朴素 DFT 的复杂度是O(N²)。FFT 是数字宇宙中最重要的算法之一。

这是关于 DFT 描述的直观解释。我们不是在实线上表示我们的信号，而是在圆上表示它。我们可以通过在圆上进行 1、2 或任意数量k的圈来播放整个信号。因此，当k固定时，我们用一个角度离散傅立叶变换和与原始距离相等的x[n]表示信号的每个值x[n]。

如果信号显示出某种周期性的k圈，意味着许多相关值将在该确切频率上叠加，以致系数X[k]将会很大。换句话说，第k个系数的模*|X[k]|表示与该频率相关的信号的能量*。

在下面的图中，信号是频率f=3 Hz的正弦波。该信号的点以蓝色显示，位于一个角度的位置离散傅立叶变换。它们在复平面上的代数和以红色显示。这些向量表示信号 DFT 的不同系数。

离散傅立叶变换

DFT 的插图

下一张图表示前一个信号的功率谱密度（PSD）：

离散傅立叶变换

前述示例中信号的 PSD

逆傅立叶变换

通过考虑所有可能的频率，我们在频率域中对我们的数字信号进行了精确表示。我们可以通过计算逆快速傅立叶变换来恢复初始信号，这计算了逆离散傅立叶变换。这个公式与 DFT 非常相似：

逆傅立叶变换

当需要寻找周期性模式时，DFT 非常有用。然而，通常来说，傅里叶变换无法检测到特定频率下的瞬态变化。此时需要更局部的谱方法，例如小波变换。

还有更多…

以下链接包含有关傅里叶变换的更多详细信息：

SciPy 中的 FFT 简介，链接地址为scipy-lectures.github.io/intro/scipy.html#fast-fourier-transforms-scipy-fftpack
SciPy 中 fftpack 的参考文档，链接地址为docs.scipy.org/doc/scipy/reference/fftpack.html
维基百科上的傅里叶变换，链接地址为en.wikipedia.org/wiki/Fourier_transform
维基百科上的离散傅里叶变换，链接地址为en.wikipedia.org/wiki/Discrete_Fourier_transform
维基百科上的快速傅里叶变换，链接地址为en.wikipedia.org/wiki/Fast_Fourier_transform
维基百科上的分贝，链接地址为en.wikipedia.org/wiki/Decibel

另见

将线性滤波器应用于数字信号方法
计算时间序列的自相关方法

将线性滤波器应用于数字信号

线性滤波器在信号处理中发挥着基础作用。通过线性滤波器，可以从数字信号中提取有意义的信息。

在本方法中，我们将展示两个使用股市数据（NASDAQ 股票交易所）的例子。首先，我们将使用低通滤波器平滑一个非常嘈杂的信号，以提取其慢速变化。我们还将对原始时间序列应用高通滤波器，以提取快速变化。这只是线性滤波器应用中常见的两个例子，实际应用有很多种。

做好准备

从本书的 GitHub 仓库下载纳斯达克数据集，链接地址为github.com/ipython-books/cookbook-data，并将其解压到当前目录。

数据来自finance.yahoo.com/q/hp?s=^IXIC&a=00&b=1&c=1990&d=00&e=1&f=2014&g=d。

如何操作…

让我们导入所需的包：

In [1]: import numpy as np
        import scipy as sp
        import scipy.signal as sg
        import pandas as pd
        import matplotlib.pyplot as plt
        %matplotlib inline

我们使用 pandas 加载 NASDAQ 数据：

In [2]: nasdaq_df = pd.read_csv('data/nasdaq.csv')
In [3]: nasdaq_df.head()
Out[3]:  Date     Open     High      Low    Close
0  2013-12-31  4161.51  4177.73  4160.77  4176.59 
1  2013-12-30  4153.58  4158.73  4142.18  4154.20

让我们提取两列数据：日期和每日收盘值：

In [4]: date = pd.to_datetime(nasdaq_df['Date'])
        nasdaq = nasdaq_df['Close']

让我们来看一下原始信号：

In [5]: plt.plot_date(date, nasdaq, '-')

如何操作...

现在，我们将采用第一种方法来获取信号的慢速变化。我们将信号与三角窗进行卷积，这相当于FIR 滤波器。我们将在本食谱的*它是如何工作的...*部分解释该方法背后的原理。现在，暂且说，我们将每个值替换为该值周围信号的加权平均值：

In [6]: # We get a triangular window with 60 samples.
        h = sg.get_window('triang', 60)
        # We convolve the signal with this window.
        fil = sg.convolve(nasdaq, h/h.sum())
In [7]: # We plot the original signal...
        plt.plot_date(date, nasdaq, '-', lw=1)
        # ... and the filtered signal.
        plt.plot_date(date, fil[:len(nasdaq)-1], '-')

如何操作...

现在，让我们使用另一种方法。我们创建一个 IIR 巴特沃斯低通滤波器来提取信号的慢速变化。filtfilt()方法允许我们前后应用滤波器，以避免相位延迟：

In [8]: plt.plot_date(date, nasdaq, '-', lw=1)
        # We create a 4-th order Butterworth low-pass
        # filter.
        b, a = sg.butter(4, 2./365)
        # We apply this filter to the signal.
        plt.plot_date(date, sg.filtfilt(b, a, nasdaq), 
                      '-')

如何操作...

最后，我们使用相同的方法来创建一个高通滤波器，并提取信号的快速变化：
```
In [9]: plt.plot_date(date, nasdaq, '-', lw=1)
        b, a = sg.butter(4, 2*5./365, btype='high')
        plt.plot_date(date, sg.filtfilt(b, a, nasdaq),
                      '-', lw=.5)
```
2000 年左右的快速变化对应于互联网泡沫的破裂，反映了当时股市的高波动性和股市指数的快速波动。更多细节，请参见en.wikipedia.org/wiki/Dot-com_bubble。

它是如何工作的...

在本节中，我们将解释在数字信号背景下线性滤波器的基本原理。

数字信号是一个离散的序列*(x[n])，由n索引！它是如何工作的... * 0。尽管我们常假设信号是无限序列，但在实际应用中，信号通常由有限大小N的向量表示。

在连续情况下，我们更倾向于操作依赖时间的函数f(t)。宽泛地说，我们可以通过离散化时间并将积分转化为求和，从连续信号转化为离散信号。

什么是线性滤波器？

线性滤波器 F 将输入信号 x = (x[n]) 转换为输出信号 y = (y[n])。这种转换是线性的——两个信号之和的转换是转换后的信号之和：F(x+y) = F(x)+F(y)。

除此之外，将输入信号乘以常数！什么是线性滤波器？会产生与将原始输出信号乘以相同常数相同的输出：。

线性时不变 (LTI) 滤波器有一个额外的性质：如果信号*(x[n])被转换为(y[n])，那么移位信号(x[(n-k)])将被转换为(y[(n-k)])，对于任何固定的k*。换句话说，系统是时不变的，因为输出不依赖于输入应用的具体时间。

注意

从现在开始，我们只考虑 LTI 滤波器。

线性滤波器与卷积

LTI 系统理论中的一个非常重要的结果是，LTI 滤波器可以通过一个单一的信号来描述：冲激响应h。这是滤波器对冲激信号的响应输出。对于数字滤波器，冲激信号为*(1, 0, 0, 0, ...)*。

可以证明，x = (x[n])通过卷积与冲激响应h及信号x变换为y = (y[n])：

线性滤波器和卷积

卷积是信号处理中一个基本的数学运算。从直觉上讲，考虑到一个在零点附近有峰值的卷积函数，卷积相当于对信号（这里是x）进行加权的局部平均，权重由给定的窗口（这里是h）决定。

根据我们的符号，隐含着我们将自己限制在因果滤波器（h[n] = 0 当n < 0）中。这个特性意味着信号的输出仅依赖于输入的当前和过去，而不是未来。这在许多情况下是一个自然的特性。

FIR 和 IIR 滤波器

一个信号*(h[n])的支持是满足！FIR 和 IIR 滤波器的n*集合。LTI 滤波器可以分为两类：

一个有限冲激响应（FIR）滤波器具有有限支持的冲激响应
一个无限冲激响应（IIR）滤波器具有无限支持的冲激响应

一个 FIR 滤波器可以通过大小为N（一个向量）的有限冲激响应来描述。它通过将信号与其冲激响应进行卷积来工作。我们定义b[n] = h[n]，当n满足！FIR 和 IIR 滤波器 * N*时。然后，y[n]是输入信号最后N+1个值的线性组合：

FIR 和 IIR 滤波器

另一方面，IIR 滤波器通过具有无限冲激响应来描述，这种形式下无法精确表示。出于这个原因，我们通常使用替代表示：

FIR 和 IIR 滤波器

这个差分方程将y[n]表示为输入信号的最后N+1个值的线性组合（前馈项，类似于 FIR 滤波器），以及输出信号的最后M个值的线性组合（反馈项）。反馈项使得 IIR 滤波器比 FIR 滤波器更复杂，因为输出不仅依赖于输入，还依赖于输出的先前值（动态性）。

频域中的滤波器

我们只描述了时域中的滤波器。其他域中的替代表示方法如拉普拉斯变换、Z 变换和傅里叶变换等也存在。

特别地，傅里叶变换具有一个非常方便的特性：它将卷积转化为频域中的乘法。换句话说，在频域中，一个 LTI 滤波器将输入信号的傅里叶变换与冲激响应的傅里叶变换相乘。

低通、高通和带通滤波器

滤波器可以通过它们对输入信号频率幅度的影响来表征。具体如下：

低通滤波器衰减高于截止频率的信号分量
高通滤波器衰减低于截止频率的信号成分，低频部分
带通滤波器通过特定频率范围内的信号成分，并衰减范围外的信号成分

在这个配方中，我们首先将输入信号与一个三角形窗口进行卷积（窗口有限支持）。可以证明，这个操作相当于一个低通 FIR 滤波器。这是移动平均方法的一个特殊案例，该方法通过计算每个值的局部加权平均来平滑信号。

接着，我们应用了两次巴特沃斯滤波器，这是一种特定类型的 IIR 滤波器，能够作为低通、高通或带通滤波器。在这个配方中，我们首先将其用作低通滤波器来平滑信号，接着再用作高通滤波器来提取信号中的快速变化部分。

还有更多...

这里有一些关于数字信号处理和线性滤波器的常规参考资料：

维基百科上的数字信号处理，详见en.wikipedia.org/wiki/Digital_signal_processing
维基百科上的线性滤波器，详见en.wikipedia.org/wiki/Linear_filter
维基百科上的 LTI 滤波器，详见en.wikipedia.org/wiki/LTI_system_theory

这里有一些关于冲激响应、卷积以及 FIR/IIR 滤波器的参考资料：

冲激响应在en.wikipedia.org/wiki/Impulse_response中有描述
卷积在en.wikipedia.org/wiki/Convolution中有描述
FIR 滤波器在en.wikipedia.org/wiki/Finite_impulse_response中有描述
IIR 滤波器在en.wikipedia.org/wiki/Infinite_impulse_response中有描述
低通滤波器在en.wikipedia.org/wiki/Low-pass_filter中有描述
高通滤波器在en.wikipedia.org/wiki/High-pass_filter中有描述
带通滤波器在en.wikipedia.org/wiki/Band-pass_filter中有描述

另请参见

使用快速傅里叶变换分析信号的频率成分 配方

计算时间序列的自相关

时间序列的自相关可以帮助我们了解重复的模式或序列相关性。后者是指信号在某一时刻与稍后时刻之间的相关性。自相关分析可以告诉我们波动的时间尺度。在这里，我们利用这一工具分析美国婴儿名字的变化，数据来源于美国社会保障管理局提供的数据。

准备工作

从本书的 GitHub 仓库下载Babies数据集，链接：github.com/ipython-books/cookbook-data，并将其解压到当前目录。

数据来自于www.data.gov（catalog.data.gov/dataset/baby-names-from-social-security-card-applications-national-level-data-6315b）。

如何操作...

我们导入以下包：

In [1]: import os
        import numpy as np
        import pandas as pd
        import matplotlib.pyplot as plt
        %matplotlib inline

我们使用 pandas 读取数据。数据集每年包含一个 CSV 文件。每个文件包含该年所有婴儿名字及其相应的频率。我们将数据加载到一个字典中，字典中每年对应一个DataFrame：

In [2]: files = [file for file in os.listdir('data/') 
                 if file.startswith('yob')]
In [3]: years = np.array(sorted([int(file[3:7]) 
                                 for file in files]))
In [4]: data = {year: 
                pd.read_csv(
                 'data/yob{y:d}.txt'.format(y=year),
                  index_col=0, header=None, 
                  names=['First name', 'Gender', 'Number']) 
                                         for year in years}
In [5]: data[2012].head()
Out[5]:    Gender  Number
First name               
Sophia          F   22158
Emma            F   20791
Isabella        F   18931
Olivia          F   17147
Ava             F   15418

我们编写函数以根据名字、性别和出生年份检索婴儿名字的频率：

In [6]: def get_value(name, gender, year):
            """Return the number of babies born a given 
            year, with a given gender and a given name."""
            try:
                return data[year] \
                       [data[year]['Gender'] == gender] \
                       ['Number'][name]
            except KeyError:
                return 0
In [7]: def get_evolution(name, gender):
            """Return the evolution of a baby name over the
            years."""
            return np.array([get_value(name, gender, year) 
                             for year in years])

让我们定义一个计算信号自相关的函数。这个函数本质上是基于 NumPy 的correlate()函数。

In [8]: def autocorr(x):
            result = np.correlate(x, x, mode='full')
            return result[result.size/2:]

现在，我们创建一个显示婴儿名字及其（归一化）自相关演变的函数：

In [9]: def autocorr_name(name, gender):
            x = get_evolution(name, gender)
            z = autocorr(x)
            # Evolution of the name.
            plt.subplot(121)
            plt.plot(years, x, '-o', label=name)
            plt.title("Baby names")
            # Autocorrelation.
            plt.subplot(122)
            plt.plot(z / float(z.max()), '-', label=name)
            plt.legend()
            plt.title("Autocorrelation")

让我们看一下两个女性名字：
```
In [10]: autocorr_name('Olivia', 'F')
         autocorr_name('Maria', 'F')
```
Olivia 的自相关衰减速度比 Maria 的要快得多。这主要是因为 Olivia 在 20 世纪末的急剧增加。相比之下，Maria 的变化较为缓慢，且其自相关衰减也较为缓慢。

它是如何工作的...

时间序列是按时间索引的序列。其重要应用包括股市、产品销售、天气预报、生物信号等。时间序列分析是统计数据分析、信号处理和机器学习的重要部分。

自相关有多种定义。在这里，我们将时间序列*(x[n])*的自相关定义为：

它是如何工作的...

在之前的图表中，我们通过最大值归一化自相关，以便比较两个信号的自相关。自相关量化了信号与其自身平移版本之间的平均相似性，这一相似性是延迟时间的函数。换句话说，自相关可以为我们提供有关重复模式以及信号波动时间尺度的信息。自相关衰减至零的速度越快，信号变化的速度越快。

还有更多...

以下是一些参考资料：

NumPy 的相关函数文档，见docs.scipy.org/doc/numpy/reference/generated/numpy.correlate.html
statsmodels 中的自相关函数，文档可见statsmodels.sourceforge.net/stable/tsa.html
维基百科上的时间序列，见en.wikipedia.org/wiki/Time_series
Wikipedia 上的序列依赖，详情请见 en.wikipedia.org/wiki/Serial_dependence
Wikipedia 上的自相关，详情请见 en.wikipedia.org/wiki/Autocorrelation

参见

使用快速傅里叶变换分析信号频率成分 的方法

第十一章：图像与音频处理

在本章中，我们将讨论以下主题：

操作图像的曝光度
对图像应用滤波器
对图像进行分割
在图像中找到兴趣点
使用 OpenCV 检测图像中的人脸
应用数字滤波器于语音声音
在笔记本中创建一个声音合成器

介绍

在前一章中，我们讨论了针对一维时间依赖信号的信号处理技术。在本章中，我们将看到针对图像和声音的信号处理技术。

通用信号处理技术可以应用于图像和声音，但许多图像或音频处理任务需要专门的算法。例如，我们将看到用于图像分割、检测图像中的兴趣点或检测人脸的算法。我们还将听到线性滤波器对语音声音的影响。

scikit-image是 Python 中的主要图像处理包之一。在本章的大多数图像处理实例中，我们将使用它。有关 scikit-image 的更多信息，请参阅scikit-image.org。

我们还将使用OpenCV（opencv.org），这是一个 C++计算机视觉库，具有 Python 包装器。它实现了专门的图像和视频处理任务的算法，但使用起来可能有些困难。一个有趣的（且更简单的）替代方案是SimpleCV（simplecv.org）。

在本介绍中，我们将从信号处理的角度讨论图像和声音的特点。

图像

灰度图像是一个二维信号，由一个函数f表示，该函数将每个像素映射到一个强度。强度可以是一个在 0（暗）和 1（亮）之间的实数值。在彩色图像中，该函数将每个像素映射到强度的三元组，通常是红色、绿色和蓝色（RGB）分量。

在计算机上，图像是数字化采样的。其强度不再是实数值，而是整数或浮动点数。一方面，连续函数的数学公式使我们能够应用诸如导数和积分之类的分析工具。另一方面，我们需要考虑我们处理的图像的数字特性。

声音

从信号处理的角度来看，声音是一个时间依赖的信号，在听觉频率范围内（约 20 Hz 到 20 kHz）具有足够的功率。然后，根据奈奎斯特-香农定理（在第十章，信号处理中介绍），数字声音信号的采样率需要至少为 40 kHz。44100 Hz 的采样率是常选的采样率。

参考资料

以下是一些参考资料：

维基百科上的图像处理，网址：en.wikipedia.org/wiki/Image_processing
由 Gabriel Peyré 撰写的高级图像处理算法，网址为github.com/gpeyre/numerical-tours
维基百科上的音频信号处理，网址为en.wikipedia.org/wiki/Audio_signal_processing
44100 Hz 采样率的特殊性解释在en.wikipedia.org/wiki/44,100_Hz

操纵图像的曝光

图像的曝光告诉我们图像是太暗、太亮还是平衡的。可以通过所有像素的强度值直方图来衡量。改善图像的曝光是一项基本的图像编辑操作。正如我们将在本篇中看到的，使用 scikit-image 可以轻松实现。

准备工作

您需要 scikit-image 来完成这个步骤。您可以在scikit-image.org/download.html找到安装说明。使用 Anaconda，您只需在终端中输入conda install scikit-image。

您还需要从该书的 GitHub 仓库下载Beach数据集，网址为github.com/ipython-books/cookbook-data。

操作步骤...

让我们导入包：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        import skimage.exposure as skie
        %matplotlib inline

我们使用 matplotlib 打开一幅图像。我们只取一个 RGB 分量，以获得灰度图像（有更好的方法将彩色图像转换为灰度图像）：
```
In [2]: img = plt.imread('data/pic1.jpg')[...,0]
```

我们创建一个函数，显示图像及其强度值直方图（即曝光）：

In [3]: def show(img):
            # Display the image.
            plt.subplot(121)
            plt.imshow(img, cmap=plt.cm.gray)
            plt.axis('off')
            # Display the histogram.
            plt.subplot(122)
            plt.hist(img.ravel(), lw=0, bins=256)
            plt.xlim(0, img.max())
            plt.yticks([])
            plt.show()

让我们显示图像及其直方图：
```
In [4]: show(img)
```
一幅图像及其直方图

直方图不平衡，图像看起来过曝（许多像素太亮）。
现在，我们使用 scikit-image 的rescale_intensity函数重新调整图像的强度。in_range和out_range参数定义了从原始图像到修改后图像的线性映射。超出in_range范围的像素被剪切到out_range的极值。在这里，最暗的像素（强度小于 100）变为完全黑色（0），而最亮的像素（>240）变为完全白色（255）：
```
In [5]: show(skie.rescale_intensity(img,
             in_range=(100, 240), out_range=(0, 255)))
```
一个粗糙的曝光操作技术

直方图中似乎缺少许多强度值，这反映了这种基本曝光校正技术的质量较差。
现在我们使用一种更高级的曝光校正技术，称为对比有限自适应直方图均衡化（CLAHE）：
```
In [6]: show(skie.equalize_adapthist(img))
```
对曝光校正的对比有限自适应直方图均衡化方法的结果

直方图看起来更平衡，图像现在显得更加对比。

工作原理...

图像的直方图代表像素强度值的分布。它是图像编辑、图像处理和计算机视觉中的一个核心工具。

rescale_intensity() 函数可以伸缩图像的强度级别。一个使用案例是确保图像使用数据类型允许的整个值范围。

equalize_adapthist() 函数的工作原理是将图像分割成矩形区域，并计算每个区域的直方图。然后，像素的强度值被重新分配，以改善对比度并增强细节。

还有更多内容...

这里是一些参考资料：

图像直方图相关内容可以在 Wikipedia 上找到。
直方图均衡化相关内容可以在 Wikipedia 上找到。
自适应直方图均衡化相关内容可以在 Wikipedia 上找到。
对比度相关内容可以在 Wikipedia 上找到。

参见

在图像上应用滤波器 示例

在图像上应用滤波器

在这个示例中，我们对图像应用了多种滤波器，以实现不同的目的：模糊、去噪和边缘检测。

如何运作...

让我们导入相关包：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        import skimage
        import skimage.filter as skif
        import skimage.data as skid
        %matplotlib inline

我们创建一个函数来显示灰度图像：

In [2]: def show(img):
            plt.imshow(img, cmap=plt.cm.gray)
            plt.axis('off')
            plt.show()

现在，我们加载 Lena 图像（包含在 scikit-image 中）。我们选择一个单一的 RGB 组件以获取灰度图像：
```
In [3]: img = skimage.img_as_float(skid.lena())[...,0]
In [4]: show(img)
```
让我们对图像应用模糊的高斯滤波器：
```
In [5]: show(skif.gaussian_filter(img, 5.))
```
现在，我们应用一个Sobel 滤波器，它增强了图像中的边缘：
```
In [6]: sobimg = skif.sobel(img)
        show(sobimg)
```
我们可以对过滤后的图像进行阈值处理，得到素描效果。我们得到一个只包含边缘的二值图像。我们使用笔记本小部件来找到适当的阈值值；通过添加 @interact 装饰器，我们在图像上方显示一个滑块。这个小部件让我们可以动态控制阈值。
```
In [7]: from IPython.html import widgets
        @widgets.interact(x=(0.01, .4, .005))
        def edge(x):
            show(sobimg<x)
```

最后，我们向图像中添加一些噪声，以展示去噪滤波器的效果：

In [8]: img = skimage.img_as_float(skid.lena())
        # We take a portion of the image to show the
        # details.
        img = img[200:-100, 200:-150]
        # We add Gaussian noise.
        img = np.clip(img + 0.3*np.random.rand(*img.shape),
                      0, 1)
In [9]: show(img)

如何运作...

denoise_tv_bregman() 函数实现了使用 Split Bregman 方法的全变差去噪：

In [10]: show(skimage.restoration.denoise_tv_bregman(img,
                                                     5.))

如何运作...

图像处理中使用的许多滤波器都是线性滤波器。这些滤波器与第十章中的滤波器非常相似，信号处理；唯一的区别是它们在二维中工作。对图像应用线性滤波器等同于对图像与特定函数进行离散卷积。高斯滤波器通过与高斯函数卷积来模糊图像。

Sobel 滤波器计算图像梯度的近似值。因此，它能够检测图像中快速变化的空间变化，通常这些变化对应于边缘。

图像去噪是指从图像中去除噪声的过程。总变差去噪通过找到一个与原始（有噪声）图像接近的规则图像来工作。规则性由图像的总变差来量化：

它是如何工作的...

Split Bregman 方法是基于 L¹范数的变种。它是压缩感知的一个实例，旨在找到真实世界中有噪声测量的规则和稀疏近似值。

还有更多内容...

以下是一些参考资料：

skimage.filter 模块的 API 参考，链接：scikit-image.org/docs/dev/api/skimage.filter.html
噪声去除，Wikipedia 上有介绍，链接：en.wikipedia.org/wiki/Noise_reduction
Wikipedia 上的高斯滤波器，链接：en.wikipedia.org/wiki/Gaussian_filter
Sobel 滤波器，Wikipedia 上有介绍，链接：en.wikipedia.org/wiki/Sobel_operator
图像去噪，Wikipedia 上有介绍，链接：en.wikipedia.org/wiki/Noise_reduction
总变差去噪，Wikipedia 上有介绍，链接：en.wikipedia.org/wiki/Total_variation_denoising
Split Bregman 算法的解释，链接：www.ece.rice.edu/~tag7/Tom_G…

另请参见

图像曝光调整的配方

图像分割

图像分割包括将图像分割成具有某些特征的不同区域。这是计算机视觉、面部识别和医学成像中的一项基本任务。例如，图像分割算法可以自动检测医学图像中器官的轮廓。

scikit-image 提供了几种分割方法。在这个配方中，我们将演示如何分割包含不同物体的图像。

如何操作...

让我们导入相关包：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        from skimage.data import coins
        from skimage.filter import threshold_otsu
        from skimage.segmentation import clear_border
        from skimage.morphology import closing, square
        from skimage.measure import regionprops, label
        from skimage.color import lab2rgb
        %matplotlib inline

我们创建一个显示灰度图像的函数：

In [2]: def show(img, cmap=None):
            cmap = cmap or plt.cm.gray
            plt.imshow(img, cmap=cmap)
            plt.axis('off')
            plt.show()

我们使用 scikit-image 中捆绑的测试图像，展示了几枚硬币放置在简单背景上的样子：
```
In [3]: img = coins()
In [4]: show(img)
```
分割图像的第一步是找到一个强度阈值，将（明亮的）硬币与（暗色的）背景分开。Otsu 方法定义了一个简单的算法来自动找到这个阈值。
```
In [5]: threshold_otsu(img)
Out[5]: 107
In [6]: show(img>107)
```
使用 Otsu 方法得到的阈值图像
图像的左上角似乎存在问题，背景的一部分过于明亮。让我们使用一个笔记本小部件来找到更好的阈值：
```
In [7]: from IPython.html import widgets
        @widgets.interact(t=(10, 240))
        def threshold(t):
            show(img>t)
```
使用手动选择的阈值的阈值化图像
阈值 120 看起来更好。下一步是通过平滑硬币并去除边界来清理二值图像。scikit-image 提供了一些功能来实现这些目的。
```
In [8]: img_bin = clear_border(closing(img>120, square(5)))
        show(img_bin)
```
清理过边界的阈值化图像
接下来，我们使用label()函数执行分割任务。此函数检测图像中的连接组件，并为每个组件分配唯一标签。在这里，我们在二值图像中为标签上色：
```
In [9]: labels = label(img_bin)
        show(labels, cmap=plt.cm.rainbow)
```
分割后的图像

图像中的小伪影会导致虚假的标签，这些标签并不对应硬币。因此，我们只保留大于 100 像素的组件。regionprops()函数允许我们检索组件的特定属性（在这里是面积和边界框）：

In [10]: regions = regionprops(labels, 
                               ['Area', 'BoundingBox'])
         boxes = np.array([label['BoundingBox']
                           for label in regions 
                               if label['Area'] > 100])
         print("There are {0:d} coins.".format(len(boxes)))
There are 24 coins.

最后，我们在原始图像中每个组件上方显示标签号：

In [11]: plt.imshow(img, cmap=plt.cm.gray)
         plt.axis('off')
         xs = boxes[:,[1,3]].mean(axis=1)
         ys = boxes[:,[0,2]].mean(axis=1)
         for i, box in enumerate(boxes):
             plt.text(xs[i]-5, ys[i]+5, str(i))

如何操作...

它是如何工作的...

为了清理阈值化图像中的硬币，我们使用了数学形态学技术。这些方法基于集合理论、几何学和拓扑学，使我们能够操控形状。

例如，首先让我们解释膨胀和腐蚀。假设 A 是图像中的一组像素，b 是一个二维向量，我们表示 A[b] 为通过 b 平移的 A 集合，如下所示：

它是如何工作的...

设 B 为一个整数分量的向量集合。我们称 B 为结构元素（在这里我们使用了方形结构元素）。此集合表示一个像素的邻域。A 通过 B 的膨胀操作是：

它是如何工作的...

A通过 B 的腐蚀操作是：

它是如何工作的...

膨胀操作通过在边界附近添加像素来扩展集合。腐蚀操作移除集合中与边界过于接近的像素。闭合操作是先膨胀后腐蚀。这一操作可以去除小的黑点并连接小的亮裂缝。在本配方中，我们使用了一个方形结构元素。

还有更多内容...

以下是一些参考资料：

有关图像处理的 SciPy 讲义，链接：scipy-lectures.github.io/packages/scikit-image/
维基百科上的图像分割，链接：en.wikipedia.org/wiki/Image_segmentation
Otsu 方法用于寻找阈值，详细解释见：en.wikipedia.org/wiki/Otsu's_method
使用 scikit-image 进行分割教程（本配方的灵感来源）可见：scikit-image.org/docs/dev/user_guide/tutorial_segmentation.html
维基百科上的数学形态学，查看 en.wikipedia.org/wiki/Mathematical_morphology
skimage.morphology模块的 API 参考，查看 scikit-image.org/docs/dev/api/skimage.morphology.html

另见

第十四章中的计算图像的连通分量食谱，图形、几何学与地理信息系统

在图像中找到兴趣点

在图像中，兴趣点是可能包含边缘、角点或有趣物体的位置。例如，在一幅风景画中，兴趣点可能位于房屋或人物附近。检测兴趣点在图像识别、计算机视觉或医学影像中非常有用。

在本食谱中，我们将使用 scikit-image 在图像中找到兴趣点。这将使我们能够围绕图像中的主题裁剪图像，即使该主题不在图像的中心。

准备就绪

从本书的 GitHub 仓库下载Child数据集，链接：github.com/ipython-books/cookbook-data，并将其解压到当前目录。

如何操作...

我们导入所需的包：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        import skimage
        import skimage.feature as sf
        %matplotlib inline

我们创建一个函数来显示彩色或灰度图像：

In [2]: def show(img, cmap=None):
            cmap = cmap or plt.cm.gray
            plt.imshow(img, cmap=cmap)
            plt.axis('off')

我们加载一张图像：

In [3]: img = plt.imread('data/pic2.jpg')
In [4]: show(img)

如何操作...

让我们使用哈里斯角点法在图像中找到显著点。第一步是使用corner_harris()函数计算哈里斯角点响应图像（我们将在*如何工作...*中解释这个度量）。该函数需要一个灰度图像，因此我们选择第一个 RGB 分量：
```
In [5]: corners = sf.corner_harris(img[:,:,0])
In [6]: show(corners)
```
我们看到这个算法能很好地检测到孩子外套上的图案。

下一步是使用corner_peaks()函数从这个度量图像中检测角点：

In [7]: peaks = sf.corner_peaks(corners)
In [8]: show(img)
        plt.plot(peaks[:,1], peaks[:,0], 'or', ms=4)

如何操作...

最后，我们在角点周围创建一个框，定义我们的兴趣区域：

In [9]: ymin, xmin = peaks.min(axis=0)
        ymax, xmax = peaks.max(axis=0)
        w, h = xmax-xmin, ymax-ymin
In [10]: k = .25
         xmin -= k*w
         xmax += k*w
         ymin -= k*h
         ymax += k*h
In [11]: show(img[ymin:ymax,xmin:xmax])

如何操作...

如何工作...

让我们解释本食谱中使用的方法。第一步是计算图像的结构张量（或哈里斯矩阵）：

如何工作...

这里，*I(x,y)*是图像，*I[x]和I[y]*是偏导数，括号表示围绕邻近值的局部空间平均。

这个张量在每个点上关联一个*(2,2)*的正对称矩阵。该矩阵用于计算图像在每个点上的自相关。

让它是如何工作的... 和成为这个矩阵的两个特征值（该矩阵是可对角化的，因为它是实数且对称的）。大致上，角点是通过各个方向的自相关变化很大来表征的，或者通过较大的正特征值和。角点度量图像定义为：

它是如何工作的...

在这里，k 是一个可调节的参数。当存在角点时，M 会很大。最后，corner_peaks() 通过查找角点度量图像中的局部最大值来找到角点。

使用 OpenCV 检测图像中的人脸

OpenCV（开放计算机视觉）是一个开源的 C++ 库，用于计算机视觉。它包含图像分割、物体识别、增强现实、人脸检测以及其他计算机视觉任务的算法。

在这个教程中，我们将使用 Python 中的 OpenCV 来检测图片中的人脸。

准备工作

你需要安装 OpenCV 和 Python 的包装器。你可以在 OpenCV 的官方网站找到安装说明：docs.opencv.org/trunk/doc/py_tutorials/py_tutorials.html。

在 Windows 上，你可以安装 Chris Gohlke 的包，链接地址：www.lfd.uci.edu/~gohlke/pyt…。

你还需要从书本的 GitHub 仓库下载 Family 数据集，链接地址：github.com/ipython-books/cookbook-data。

注意

在写这篇文章时，OpenCV 尚不兼容 Python 3。因此，本教程要求使用 Python 2。

如何操作...

首先，我们导入所需的包：

In [1]: import numpy as np
        import cv2
        import matplotlib.pyplot as plt
        %matplotlib inline

我们用 OpenCV 打开 JPG 图像：

In [2]: img = cv2.imread('data/pic3.jpg')

然后，我们使用 OpenCV 的 cvtColor() 函数将其转换为灰度图像。对于人脸检测，使用灰度图像已经足够且更快速。
```
In [3]: gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
```
为了检测人脸，我们将使用Viola–Jones 物体检测框架。一个 Haar-like 分类器级联已经在大量图像上训练，以检测人脸（更多细节将在下一节提供）。训练的结果存储在一个 XML 文件中（该文件属于Family数据集，数据集可在本书的 GitHub 仓库中找到）。我们使用 OpenCV 的CascadeClassifier类从这个 XML 文件中加载该级联：
```
In [4]: face_cascade = cv2.CascadeClassifier(
              'data/haarcascade_frontalface_default.xml')
```
最后，分类器的detectMultiScale()方法在灰度图像上检测物体，并返回围绕这些物体的矩形框列表：
```
In [5]: for x,y,w,h in \
                 face_cascade.detectMultiScale(gray, 1.3):
            cv2.rectangle(gray, (x,y), (x+w,y+h),
                          (255,0,0), 2)
        plt.imshow(gray, cmap=plt.cm.gray)
        plt.axis('off')
```
我们可以看到，尽管所有检测到的物体确实是人脸，但每四张脸中就有一张没有被检测到。这可能是因为这张脸并没有完全正对摄像头，而训练集中的人脸则都是正对着摄像头的。这表明该方法的有效性受限于训练集的质量和通用性。

它是如何工作的...

Viola–Jones 物体检测框架通过训练一系列使用 Haar-like 特征的提升分类器来工作。首先，我们考虑一组特征：

How it works...

Haar-like 特征

一个特征被定位在图像中的特定位置和大小。它覆盖了图像中的一个小窗口（例如 24 x 24 像素）。黑色区域的所有像素和白色区域的所有像素之和相减。这一操作可以通过积分图像高效地完成。

然后，所有分类器集通过提升技术进行训练；在训练过程中，只有最好的特征会被保留下来进入下一阶段。训练集包含正负图像（有脸和没有脸的图像）。尽管每个分类器单独的表现较差，但是这些提升分类器的级联方式既高效又快速。因此，这种方法非常适合实时处理。

XML 文件已从 OpenCV 包中获得。该文件对应多个训练集。你也可以使用自己的训练集训练自己的级联分类器。

还有更多...

以下是一些参考资料：

一个关于级联教程的 OpenCV (C++)教程可在docs.opencv.org/doc/tutorials/objdetect/cascade_classifier/cascade_classifier.html找到
训练级联的文档可在docs.opencv.org/doc/user_guide/ug_traincascade.html找到
Haar 级联库可在github.com/Itseez/opencv/tree/master/data/haarcascades找到
OpenCV 的级联分类 API 参考可在docs.opencv.org/modules/objdetect/doc/cascade_classification.html找到
维基百科上的 Viola–Jones 目标检测框架，网址为en.wikipedia.org/wiki/Viola%E2%80%93Jones_object_detection_framework
提升或如何从许多弱分类器创建一个强分类器，解释在en.wikipedia.org/wiki/Boosting_(machine_learning)

将数字滤波器应用于语音

在这个示例中，我们将展示如何在笔记本中播放声音。我们还将说明简单数字滤波器对语音的影响。

准备工作

您需要pydub包。您可以使用pip install pydub安装它，或从github.com/jiaaro/pydub/下载。

该软件包需要开源多媒体库 FFmpeg 来解压缩 MP3 文件，网址为www.ffmpeg.org。

这里给出的代码适用于 Python 3。您可以在书的 GitHub 存储库中找到 Python 2 版本。

如何做…

让我们导入这些包：

In [1]: import urllib
        from io import BytesIO
        import numpy as np
        import scipy.signal as sg
        import pydub
        import matplotlib.pyplot as plt
        from IPython.display import Audio, display
        %matplotlib inline

我们创建一个 Python 函数，从英文句子生成声音。这个函数使用 Google 的文本转语音（TTS）API。我们以 MP3 格式检索声音，并用 pydub 将其转换为 Wave 格式。最后，我们通过删除 NumPy 中的波形头检索原始声音数据：

In [2]: def speak(sentence):
            url = ("http://translate.google.com/"
                    "translate_tts?tl=en&q=") + 
                        urllib.parse.quote_plus(sentence)
            req = urllib.request.Request(url,
                             headers={'User-Agent': ''}) 
            mp3 = urllib.request.urlopen(req).read()
            # We convert the mp3 bytes to wav.
            audio = pydub.AudioSegment.from_mp3(
                                            BytesIO(mp3))
            wave = audio.export('_', format='wav')
            wave.seek(0)
            wave = wave.read()
            # We get the raw data by removing the 24 
            # first bytes of the header.
            x = np.frombuffer(wave, np.int16)[24:] / 2.**15

我们创建一个函数，在笔记本中播放声音（由 NumPy 向量表示），使用 IPython 的Audio类：
```
In [3]: def play(x, fr, autoplay=False):
            display(Audio(x, rate=fr, autoplay=autoplay))
```

让我们播放声音“Hello world.”，并使用 matplotlib 显示波形：

In [4]: x, fr = speak("Hello world")
        play(x, fr)
        t = np.linspace(0., len(x)/fr, len(x))
        plt.plot(t, x, lw=1)

如何做…

现在，我们将听到应用于此声音的 Butterworth 低通滤波器的效果（500 Hz 截止频率）：

In [5]: b, a = sg.butter(4, 500./(fr/2.), 'low')
        x_fil = sg.filtfilt(b, a, x)
In [6]: play(x_fil, fr)
        plt.plot(t, x, lw=1)
        plt.plot(t, x_fil, lw=1)

如何做…

我们听到了一个沉闷的声音。

现在，使用高通滤波器（1000 Hz 截止频率）：

In [7]: b, a = sg.butter(4, 1000./(fr/2.), 'high')
        x_fil = sg.filtfilt(b, a, x)
In [8]: play(x_fil, fr)
        plt.plot(t, x, lw=1)
        plt.plot(t, x_fil, lw=1)

如何做…

听起来像一个电话。

最后，我们可以创建一个简单的小部件，快速测试高通滤波器的效果，带有任意截止频率：

In [9]: from IPython.html import widgets
        @widgets.interact(t=(100., 5000., 100.))
        def highpass(t):
            b, a = sg.butter(4, t/(fr/2.), 'high')
            x_fil = sg.filtfilt(b, a, x)
            play(x_fil, fr, autoplay=True)

我们得到一个滑块，让我们改变截止频率并实时听到效果。

工作原理…

人耳可以听到高达 20 kHz 的频率。人声频带范围大约从 300 Hz 到 3000 Hz。

数字滤波器在第十章中有描述，信号处理。这里给出的示例允许我们听到低通和高通滤波器对声音的影响。

还有更多…

这里有一些参考资料：

维基百科上的音频信号处理，网址为en.wikipedia.org/wiki/Audio_signal_processing
维基百科上的音频滤波器，网址为en.wikipedia.org/wiki/Audio_filter
维基百科上的声音频率，网址为en.wikipedia.org/wiki/Voice_frequency
PyAudio，一个使用 PortAudio 库的音频 Python 包，链接在people.csail.mit.edu/hubert/pyaudio/

另见

在笔记本中创建声音合成器食谱

在笔记本中创建声音合成器

在这个食谱中，我们将在笔记本中创建一个小型电子钢琴。我们将使用 NumPy 合成正弦声音，而不是使用录制的音频。

如何操作...

我们导入模块：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        from IPython.display import (Audio, display,
                                     clear_output)
        from IPython.html import widgets
        from functools import partial
        %matplotlib inline

我们定义音符的采样率和持续时间：

In [2]: rate = 16000.
        duration = 0.5
        t = np.linspace(0., duration, rate * duration)

我们创建一个函数，使用 NumPy 和 IPython 的 Audio 类生成并播放指定频率的音符（正弦函数）声音：

In [3]: def synth(f):
            x = np.sin(f * 2\. * np.pi * t)
            display(Audio(x, rate=rate, autoplay=True))

这是基础的 440 Hz 音符：
```
In [4]: synth(440)
```

现在，我们生成钢琴的音符频率。十二平均律是通过一个公比为 2^(1/12) 的几何级数得到的：

In [5]: notes = zip(('C,C#,D,D#,E,F,F#,G,G#,'
                     'A,A#,B,C').split(','),
                     440\. * 2 ** (np.arange(3, 17) / 12.))

最后，我们使用笔记本小部件创建钢琴。每个音符是一个按钮，所有按钮都包含在一个水平框容器中。点击某个音符会播放对应频率的声音。钢琴布局与《第三章》中的使用交互式小部件——笔记本中的钢琴食谱相同，掌握笔记本一书中的布局相同。

In [6]: container = widgets.ContainerWidget()
        buttons = []
        for note, f in notes:
            button = widgets.ButtonWidget(description=note)
            def on_button_clicked(f, b):
                clear_output()
                synth(f)
            button.on_click(partial(on_button_clicked, f))
            button.set_css({...})
            buttons.append(button)
        container.children = buttons
        display(container)
        container.remove_class('vbox')
        container.add_class('hbox')

注意

此处使用的 IPython API 设计布局基于 IPython 2.x；在 IPython 3.0 中会有所不同。

它是如何工作的...

纯音是具有正弦波形的音调。它是表示音乐音符的最简单方式。由乐器发出的音符通常更加复杂，尽管声音包含多个频率，但我们通常感知到的是音乐音调（基频）。

通过生成另一个周期性函数代替正弦波形，我们会听到相同的音调，但音色（timbre）不同。电子音乐合成器基于这个原理。

还有更多...

以下是一些参考资料：

维基百科上的合成器，链接在en.wikipedia.org/wiki/Synthesizer
维基百科上的均等律，链接在en.wikipedia.org/wiki/Equal_temperament
维基百科上的音阶，链接在en.wikipedia.org/wiki/Chromatic_scale
维基百科上的纯音，链接在en.wikipedia.org/wiki/Pure_tone
维基百科上的音色，链接在en.wikipedia.org/wiki/Timbre

另见

应用数字滤波器于语音声音食谱
《第三章》中的使用交互式小部件——笔记本中的钢琴食谱，掌握笔记本。

第十二章：确定性动力学系统

在本章中，我们将讨论以下主题：

绘制混沌动力学系统的分叉图
模拟一个初级细胞自动机
使用 SciPy 模拟常微分方程
模拟偏微分方程——反应扩散系统与图灵模式

介绍

前几章讨论了数据科学中的经典方法：统计学、机器学习和信号处理。在本章和下一章中，我们将介绍另一种方法。我们不会直接分析数据，而是模拟代表数据生成方式的数学模型。一个具有代表性的模型能为我们提供数据背后现实世界过程的解释。

具体而言，我们将介绍几个动力学系统的例子。这些数学方程描述了量随时间和空间的演变。它们可以表示物理、化学、生物学、经济学、社会科学、计算机科学、工程学等学科中各种现实世界的现象。

在本章中，我们将讨论确定性动力学系统。该术语与随机系统相对，后者的规则中包含了随机性。我们将在下一章讨论随机系统。

动力学系统的类型

我们将在这里讨论的确定性动力学系统类型包括：

离散时间动力学系统（迭代函数）
细胞自动机
常微分方程（ODEs）
偏微分方程（PDEs）

在这些模型中，感兴趣的量依赖于一个或多个独立变量。这些变量通常包括时间和/或空间。独立变量可以是离散的或连续的，从而导致不同类型的模型和不同的分析与仿真技术。

离散时间动力学系统是通过在初始点上迭代应用一个函数来描述的：f(x), f(f(x)), f(f(f(x)))，以此类推。这种类型的系统可能会导致复杂且混沌的行为。

细胞自动机由一个离散的单元格网格表示，每个单元格可以处于有限个状态之一。规则描述了单元格状态如何根据相邻单元格的状态演化。这些简单的模型可以导致极其复杂的行为。

一个常微分方程描述了一个连续函数如何依赖于其相对于独立变量的导数。在微分方程中，未知变量是一个函数而不是数字。常微分方程特别出现在量的变化率依赖于该量的当前值的情况。例如，在经典力学中，运动定律（包括行星和卫星的运动）可以通过常微分方程来描述。

PDEs 与 ODEs 相似，但它们涉及多个独立变量（例如时间和空间）。这些方程包含关于不同独立变量的偏导数。例如，PDEs 描述波的传播（声波、电磁波或机械波）和流体（流体力学）。它们在量子力学中也很重要。

微分方程

ODE 和 PDE 可以是单维或多维的，取决于目标空间的维度。多个微分方程的系统可以看作是多维方程。

ODE 或 PDE 的阶数指的是方程中最大导数的阶数。例如，一阶方程仅涉及简单的导数，二阶方程则涉及二阶导数（导数的导数），以此类推。

常微分方程或偏微分方程有额外的规则：初始和 边界条件。这些公式描述了所求函数在空间和时间域边界上的行为。例如，在经典力学中，边界条件包括物体在力作用下的初始位置和初速度。

动态系统通常根据规则是否线性，分为线性和非线性系统（相对于未知函数）。非线性方程通常比线性方程在数学和数值上更加难以研究。它们可能会导致极其复杂的行为。

例如，Navier–Stokes 方程，一组描述流体物质运动的非线性 PDEs，可能会导致湍流，这是一种在许多流体流动中出现的高度混乱的行为。尽管在气象学、医学和工程学中具有重要意义，但 Navier-Stokes 方程的基本性质目前仍未为人所知。例如，三维中的存在性和平滑性问题是七个克雷数学研究所千年奖问题之一。对于任何能够提出解决方案的人，奖励为一百万美元。

参考文献

以下是一些参考文献：

维基百科上的动态系统概述，链接：en.wikipedia.org/wiki/Dynamical_system
动态系统的数学定义，链接：en.wikipedia.org/wiki/Dynamical_system_%28definition%29
动态系统主题列表，链接：en.wikipedia.org/wiki/List_of_dynamical_systems_and_differential_equations_topics
维基百科上的 Navier-Stokes 方程，链接：en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations
Prof. Lorena Barba 的《计算流体动力学》课程，使用 IPython 笔记本编写，课程内容可以在 github.com/barbagroup/CFDPython 获取

绘制混沌动力系统的分叉图

混沌动力系统对初始条件非常敏感；在任何给定时刻的微小扰动都会产生完全不同的轨迹。混沌系统的轨迹通常具有复杂且不可预测的行为。

许多现实世界的现象是混沌的，特别是那些涉及多个主体之间非线性相互作用的现象（复杂系统）。气象学、经济学、生物学和其他学科中都有著名的例子。

在这个示例中，我们将模拟一个著名的混沌系统：逻辑映射。这是一个经典的例子，展示了如何从一个非常简单的非线性方程中产生混沌。逻辑映射模型描述了一个种群的演化，考虑到繁殖和密度依赖性死亡（饥饿）。

我们将绘制系统的分叉图，它展示了作为系统参数函数的可能长期行为（平衡点、固定点、周期轨道和混沌轨迹）。我们还将计算系统的李雅普诺夫指数的近似值，以表征模型对初始条件的敏感性。

如何实现...

让我们导入 NumPy 和 matplotlib：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        %matplotlib inline

我们通过以下方式定义逻辑函数：

我们的离散动力系统由逻辑函数的递归应用定义：

这是该函数在 Python 中的实现：

In [2]: def logistic(r, x):
            return r*x*(1-x)

我们为 10000 个 r 值进行模拟，这些值在线性间隔的 2.5 和 4 之间，并通过 NumPy 向量化模拟，考虑到一个独立系统的向量（每个参数值对应一个动力系统）：
```
In [3]: n = 10000
        r = np.linspace(2.5, 4.0, n)
```
让我们模拟 1000 次逻辑映射的迭代，并保留最后 100 次迭代，以显示分叉图：
```
In [4]: iterations = 1000
        last = 100
```
我们使用相同的初始条件 x[0] = 0.00001 初始化我们的系统：
```
In [5]: x = 1e-5 * np.ones(n)
```
我们还计算了每个 r 值的李雅普诺夫指数近似值。李雅普诺夫指数的定义是：
我们首先初始化 lyapunov 向量：
```
In [6]: lyapunov = np.zeros(n)
```

现在，我们模拟系统并绘制分叉图。该模拟仅涉及在我们的向量 x 上迭代评估 logistic() 函数。然后，为了显示分叉图，我们在最后 100 次迭代期间每个点 x[n]^((r)) 绘制一个像素：

In [7]: plt.subplot(211)
        for i in range(iterations):
            x = logistic(r, x)
            # We compute the partial sum of the
            # Lyapunov exponent.
            lyapunov += np.log(abs(r-2*r*x))
            # We display the bifurcation diagram.
            if i >= (iterations - last):
                plt.plot(r, x, ',k', alpha=.02)
        plt.xlim(2.5, 4)
        plt.title("Bifurcation diagram")

        # We display the Lyapunov exponent.
        plt.subplot(212)
        plt.plot(r[lyapunov<0],
                 lyapunov[lyapunov<0] / iterations,
                 ',k', alpha=0.1)
        plt.plot(r[lyapunov>=0],
                 lyapunov[lyapunov>=0] / iterations,
                 ',r', alpha=0.25)
        plt.xlim(2.5, 4)
        plt.ylim(-2, 1)
        plt.title("Lyapunov exponent")

如何实现...

逻辑映射的分叉图和李雅普诺夫指数

分叉图揭示了在 r<3 时存在一个固定点，然后是两个和四个平衡点，当 r 属于参数空间的某些区域时，出现混沌行为。

我们观察到李雅普诺夫指数的一个重要特性：当系统处于混沌状态时，它是正值（此处为红色）。

还有更多…

以下是一些参考资料：

维基百科上的混沌理论，可以通过en.wikipedia.org/wiki/Chaos_theory查看。
维基百科上的复杂系统，可以通过en.wikipedia.org/wiki/Complex_system查看。
维基百科上的逻辑映射，可以通过en.wikipedia.org/wiki/Logistic_map查看。
维基百科上的迭代函数（离散动力学系统），可以通过en.wikipedia.org/wiki/Iterated_function查看。
维基百科上的分叉图，可以通过en.wikipedia.org/wiki/Bifurcation_diagram查看。
维基百科上的李雅普诺夫指数，可以通过en.wikipedia.org/wiki/Lyapunov_exponent查看。

另见

使用 SciPy 模拟常微分方程示例

模拟基础细胞自动机

细胞自动机是离散的动力学系统，在一个单元格网格上演化。这些单元格可以处于有限的状态（例如，开/关）。细胞自动机的演化遵循一组规则，描述了每个单元格的状态如何根据其邻居的状态发生变化。

尽管这些模型非常简单，但它们可以引发高度复杂和混乱的行为。细胞自动机可以模拟现实世界中的现象，例如汽车交通、化学反应、森林中的火灾传播、流行病传播等。细胞自动机也存在于自然界中。例如，一些海洋贝壳的图案就是由自然细胞自动机生成的。

基础细胞自动机是一种二进制的一维自动机，其中规则涉及每个单元格的直接左右邻居。

在这个实例中，我们将使用 NumPy 模拟基础细胞自动机，并使用它们的沃尔夫拉姆代码。

如何操作…

我们导入 NumPy 和 matplotlib：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        %matplotlib inline

我们将使用以下向量来获得二进制表示的数字：
```
In [2]: u = np.array([[4], [2], [1]])
```

让我们编写一个函数，执行网格上的迭代，根据给定的规则一次性更新所有单元格，规则以二进制表示（我们将在*它是如何工作的...*部分进行详细解释）。第一步是通过堆叠循环移位版本的网格，得到每个单元格的 LCR（左，中，右）三元组(y)。然后，我们将这些三元组转换为 3 位二进制数（z）。最后，我们使用指定的规则计算每个单元格的下一个状态：

In [3]: def step(x, rule_binary):
            """Compute a single stet of an elementary
            cellular automaton."""
            # The columns contain the L, C, R values
            # of all cells.
            y = np.vstack((np.roll(x, 1), x,
                           np.roll(x, -1))).astype(np.int8)
            # We get the LCR pattern numbers 
            # between 0 and 7.
            z = np.sum(y * u, axis=0).astype(np.int8)
            # We get the patterns given by the rule.
            return rule_binary[7-z]

现在，我们编写一个函数来模拟任何基础细胞自动机。首先，我们计算规则的二进制表示（沃尔夫拉姆代码）。然后，我们用随机值初始化网格的第一行。最后，我们在网格上迭代应用函数step()：

In [4]: def generate(rule, size=80, steps=80):
            """Simulate an elementary cellular automaton
            given its rule (number between 0 and 255)."""
            # Compute the binary representation of the
            # rule.
            rule_binary = np.array(
                [int(x) for x in np.binary_repr(rule, 8)],
                dtype=np.int8)
            x = np.zeros((steps, size), dtype=np.int8)
            # Random initial state.
            x[0,:] = np.random.rand(size) < .5
            # Apply the step function iteratively.
            for i in range(steps-1):
                x[i+1,:] = step(x[i,:], rule_binary)
            return x

现在，我们模拟并显示九种不同的自动机：

In [5]: rules = [  3,  18,  30, 
                  90, 106, 110, 
                 158, 154, 184]
        for i, rule in enumerate(rules):
            x = generate(rule)
            plt.subplot(331+i)
            plt.imshow(x, interpolation='none',
                       cmap=plt.cm.binary)
            plt.xticks([]); plt.yticks([])
            plt.title(str(rule))

如何操作…

它是如何工作的…

让我们考虑一个一维的初等细胞自动机。每个细胞C有两个邻居（L和R），并且它可以是关闭的（0）或开启的（1）。因此，一个细胞的未来状态依赖于其邻居 L、C 和 R 的当前状态。这个三元组可以编码为一个 0 到 7 之间的数字（二进制表示为三位数）。

一个特定的初等细胞自动机完全由这些八种配置的结果决定。因此，存在 256 种不同的初等细胞自动机（2⁸）。每一个这样的自动机由一个介于 0 和 255 之间的数字表示。

我们按顺序考虑所有八种 LCR 状态：111、110、101、...、001、000。自动机数字的二进制表示中的每一位对应一个 LCR 状态（使用相同的顺序）。例如，在规则 110 自动机（其二进制表示为01101110）中，状态 111 会产生中心细胞的 0，状态 110 产生 1，状态 101 产生 1，依此类推。已有研究表明，这种特定的自动机是图灵完备的（或称通用的）；理论上，它能够模拟任何计算机程序。

还有更多内容...

其他类型的细胞自动机包括康威的生命游戏，这是一个二维系统。这个著名的系统可以产生各种动态模式。它也是图灵完备的。

以下是一些参考资料：

维基百科上的细胞自动机，可通过en.wikipedia.org/wiki/Cellular_automaton访问
维基百科上的初等细胞自动机，可通过en.wikipedia.org/wiki/Elementary_cellular_automaton访问
规则 110，描述见en.wikipedia.org/wiki/Rule_110
解释见en.wikipedia.org/wiki/Wolfram_code的 Wolfram 代码，将一个一维初等细胞自动机分配给 0 到 255 之间的任何数字
维基百科上的《康威的生命游戏》，可通过en.wikipedia.org/wiki/Conway's_Game_of_Life访问

使用 SciPy 模拟常微分方程

常微分方程（ODEs）描述了一个系统在内外部动态影响下的演变。具体来说，ODE 将依赖于单一自变量（例如时间）的量与其导数联系起来。此外，系统还可能受到外部因素的影响。一个一阶 ODE 通常可以写成：

使用 SciPy 模拟常微分方程

更一般地说，一个n阶 ODE 涉及到y的连续导数直到阶数n。ODE 被称为线性或非线性，取决于f是否在y中是线性的。

当一个量的变化率依赖于其值时，常微分方程自然出现。因此，ODE 在许多科学领域都有应用，如力学（受动力学力作用的物体演化）、化学（反应产物的浓度）、生物学（流行病的传播）、生态学（种群的增长）、经济学和金融等。

虽然简单的常微分方程（ODE）可以通过解析方法求解，但许多 ODE 需要数值处理。在这个例子中，我们将模拟一个简单的线性二阶自治 ODE，描述一个在重力和粘性阻力作用下的空气中粒子的演化。虽然这个方程可以通过解析方法求解，但在这里我们将使用 SciPy 进行数值模拟。

如何实现...

让我们导入 NumPy、SciPy（integrate包）和 matplotlib：

In [1]: import numpy as np
        import scipy.integrate as spi
        import matplotlib.pyplot as plt
        %matplotlib inline

我们定义了一些在模型中出现的参数：

In [2]: m = 1\.  # particle's mass
        k = 1\.  # drag coefficient
        g = 9.81  # gravity acceleration

我们有两个变量：x和y（二维）。我们记u=(x,y)。我们将要模拟的 ODE 是：

这里，g是重力加速度向量。

提示

时间导数用变量上的点表示（一个点表示一阶导数，两个点表示二阶导数）。

为了使用 SciPy 模拟这个二阶 ODE，我们可以将其转换为一阶 ODE（另一种选择是先求解u'再进行积分）。为此，我们考虑两个二维变量：u和u'。我们记v = (u, u')。我们可以将v'表示为v的函数。现在，我们创建初始向量v[0]，其时间为t=0，它有四个分量。
```
In [3]: # The initial position is (0, 0).
        v0 = np.zeros(4)
        # The initial speed vector is oriented
        # to the top right.
        v0[2] = 4.
        v0[3] = 10.
```

让我们创建一个 Python 函数f，该函数接受当前向量v(t[0])和时间t[0]作为参数（可以有可选参数），并返回导数v'(t[0])：

In [4]: def f(v, t0, k):
            # v has four components: v=[u, u'].
            u, udot = v[:2], v[2:]
            # We compute the second derivative u'' of u.
            udotdot = -k/m * udot
            udotdot[1] -= g
            # We return v'=[u', u''].
            return np.r_[udot, udotdot]

现在，我们使用不同的k值来模拟系统。我们使用 SciPy 的odeint()函数，它定义在scipy.integrate包中。

In [5]: # We want to evaluate the system on 30 linearly
        # spaced times between t=0 and t=3.
        t = np.linspace(0., 3., 30)
        # We simulate the system for different values of k.
        for k in np.linspace(0., 1., 5):
            # We simulate the system and evaluate $v$ on
            # the given times.
            v = spi.odeint(f, v0, t, args=(k,))
            # We plot the particle's trajectory.
            plt.plot(v[:,0], v[:,1], 'o-',
                     mew=1, ms=8, mec='w',
                     label='k={0:.1f}'.format(k))
        plt.legend()
        plt.xlim(0, 12)
        plt.ylim(0, 6)

如何实现...

在上图中，最外层的轨迹（蓝色）对应于无阻力运动（没有空气阻力）。它是一条抛物线。在其他轨迹中，我们可以观察到空气阻力的逐渐增加，它通过k来参数化。

它是如何工作的...

让我们解释一下如何从我们的模型中获得微分方程。设* u = (x,y)* 表示粒子在二维空间中的位置，粒子质量为m。该粒子受到两个力的作用：重力 g = (0, -9.81)（单位：m/s）和空气阻力 F = -ku'。最后一项依赖于粒子的速度，并且只在低速下有效。对于更高的速度，我们需要使用更复杂的非线性表达式。

现在，我们使用牛顿第二定律来处理经典力学中的运动。该定律表明，在惯性参考系中，粒子的质量乘以其加速度等于作用于粒子的所有力的合力。在这里，我们得到了：

它是如何工作的...

我们立即得到我们的二阶 ODE：

它是如何工作的...

我们将其转化为一个一阶常微分方程系统，定义 v=(u, u')：

它是如何工作的...

最后一项可以仅以 v 的函数表示。

SciPy 的 odeint() 函数是一个黑箱求解器；我们只需指定描述系统的函数，SciPy 会自动求解。

该函数利用了 FORTRAN 库 ODEPACK，它包含了经过充分测试的代码，已经被许多科学家和工程师使用了几十年。

一个简单的数值求解器示例是 欧拉法。为了数值求解自治常微分方程 y'=f(y)，该方法通过用时间步长 dt 离散化时间，并将 y' 替换为一阶近似：

它是如何工作的...

然后，从初始条件 y[0] = y(t[0]) 开始，方法通过以下递推关系依次计算 y：

它是如何工作的...

还有更多...

下面是一些参考资料：

SciPy 中 integrate 包的文档，可以在 docs.scipy.org/doc/scipy/reference/integrate.html 查阅
维基百科上的常微分方程（ODEs）条目，可以在 en.wikipedia.org/wiki/Ordinary_differential_equation 查阅
维基百科上的牛顿运动定律条目，可以在 en.wikipedia.org/wiki/Newton's_laws_of_motion 查阅
维基百科上的空气阻力条目，可以在 en.wikipedia.org/wiki/Drag_%28physics%29 查阅
描述常微分方程的数值方法，可以在 en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations 查阅
维基百科上的欧拉法条目，可以在 en.wikipedia.org/wiki/Euler_method 查阅
FORTRAN 中 ODEPACK 包的文档，可以在 www.netlib.org/odepack/opk… 查阅

另见：

绘制混沌动力系统的分叉图 这个配方

模拟偏微分方程 —— 反应-扩散系统和图灵模式

偏微分方程（PDEs）描述了包含时间和空间的动力系统的演化。物理学中的例子包括声音、热、电子磁学、流体流动和弹性等。生物学中的例子包括肿瘤生长、种群动态和疫情传播。

偏微分方程（PDEs）很难通过解析方法求解。因此，偏微分方程通常通过数值模拟来研究。

在这个教程中，我们将演示如何模拟由称为FitzHugh–Nagumo 方程的偏微分方程（PDE）描述的反应扩散系统。反应扩散系统模拟了一个或多个变量的演变，这些变量受到两种过程的影响：反应（变量之间的相互转化）和扩散（在空间区域中的传播）。一些化学反应可以用这种模型来描述，但在物理学、生物学、生态学以及其他学科中也有其他应用。

在这里，我们模拟了艾伦·图灵提出的一个系统，用作动物皮毛图案形成的模型。两种影响皮肤着色的化学物质根据反应扩散模型相互作用。这个系统负责形成类似斑马、美洲豹和长颈鹿皮毛的图案。

我们将使用有限差分法模拟这个系统。该方法通过离散化时间和空间，并将导数替换为其离散等价物来实现。

如何实现...

让我们导入必要的包：

In [1]: import numpy as np
        import matplotlib.pyplot as plt
        %matplotlib inline

我们将在域E=[-1,1]²上模拟以下偏微分方程系统：

变量u表示有助于皮肤着色的物质浓度，而v表示另一种与第一种物质反应并抑制着色的物质。

在初始化时，我们假设每个网格点上的u和v包含独立的随机数。我们还采用Neumann 边界条件：要求变量相对于法向量的空间导数在域的边界上为零。

让我们定义模型的四个参数：

In [2]: a = 2.8e-4
        b = 5e-3
        tau = .1
        k = -.005

我们离散化时间和空间。以下条件确保我们在此使用的离散化方案是稳定的：如何实现...

In [3]: size = 80  # size of the 2D grid
        dx = 2./size  # space step
In [4]: T = 10.0  # total time
        dt = .9 * dx**2/2  # time step
        n = int(T/dt)

我们初始化变量u和v。矩阵U和V包含这些变量在二维网格顶点上的值。这些变量被初始化为介于 0 和 1 之间的均匀噪声：
```
In [5]: U = np.random.rand(size, size)
        V = np.random.rand(size, size)
```

现在，我们定义一个函数，使用五点模板有限差分法计算网格上二维变量的离散拉普拉斯算符。这个算符定义为：如何实现...

我们可以通过向量化的矩阵操作计算网格上该算符的值。由于矩阵边缘的副作用，我们需要在计算中去除网格的边界：

In [6]: def laplacian(Z):
            Ztop = Z[0:-2,1:-1]
            Zleft = Z[1:-1,0:-2]
            Zbottom = Z[2:,1:-1]
            Zright = Z[1:-1,2:]
            Zcenter = Z[1:-1,1:-1]
            return (Ztop + Zleft + Zbottom + Zright \
                    - 4 * Zcenter) / dx**2

现在，我们使用有限差分法模拟方程系统。在每个时间步长上，我们使用离散空间导数（拉普拉斯算符）计算网格上两个方程的右侧。然后，我们使用离散时间导数更新变量：

In [7]: for i in range(n):
            # We compute the Laplacian of u and v.
            deltaU = laplacian(U)
            deltaV = laplacian(V)
            # We take the values of u and v
            # inside the grid.
            Uc = U[1:-1,1:-1]
            Vc = V[1:-1,1:-1]
            # We update the variables.
            U[1:-1,1:-1], V[1:-1,1:-1] = (
              Uc + dt * (a*deltaU + Uc - Uc**3 - Vc + k), 
              Vc + dt * (b*deltaV + Uc - Vc) / tau)
            # Neumann conditions: derivatives at the edges
            # are null.
            for Z in (U, V):
                Z[0,:] = Z[1,:]
                Z[-1,:] = Z[-2,:]
                Z[:,0] = Z[:,1]
                Z[:,-1] = Z[:,-2]