Python 数学应用（一）（4）

Python 数学应用（一）（3）https://developer.aliyun.com/article/1506374

准备工作

我们将演示在 Python 中使用先前描述的冷却方程数值地解决微分方程的技术，因为在这种情况下我们可以计算真实解。我们将初始温度取为T[0]= 50和k = 0.2。让我们也找出t值在 0 到 5 之间的解。

一般（一阶）微分方程的形式为

其中f是t（自变量）和y（因变量）的某个函数。在这个公式中，T是因变量，f(t, T) = -kt。SciPy 包中用于解决微分方程的例程需要函数f和初始值y[0][以及我们需要计算解的t值范围。要开始，我们需要在 Python 中定义我们的函数f并创建变量y[0][和t范围，准备提供给 SciPy 例程：]]

def f(t, y):
    return -0.2*y
t_range = (0, 5)

接下来，我们需要定义应从中找到解的初始条件。出于技术原因，初始y值必须指定为一维 NumPy 数组：

T0 = np.array([50.])

由于在这种情况下，我们已经知道真实解，我们也可以在 Python 中定义这个解，以便与我们将计算的数值解进行比较：

def true_solution(t):
    return 50.*np.exp(-0.2*t)

如何做到…

按照以下步骤数值求解微分方程并绘制解以及误差：

1. 我们使用 SciPy 中的integrate模块中的solve_ivp例程来数值求解微分方程。我们添加了一个最大步长的参数，值为0.1，这样解就可以在合理数量的点上计算出来：

sol = integrate.solve_ivp(f, t_range, T0, max_step=0.1)

1. 接下来，我们从solve_ivp方法返回的sol对象中提取解的值：

t_vals = sol.t
T_vals = sol.y[0, :]

1. 接下来，我们按如下方式在一组坐标轴上绘制解。由于我们还将在同一图中绘制近似误差，我们使用subplots例程创建两个子图：

fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True)
ax1.plot(t_vals, T_vals)
ax1.set_xlabel("$t$")
ax1.set_ylabel("$T$")
ax1.set_title("Solution of the cooling equation")

这在图 3.1的左侧显示了解决方案的图。

1. 为此，我们需要计算从solve_ivp例程中获得的点处的真实解，然后计算真实解和近似解之间的差的绝对值：

err = np.abs(T_vals - true_solution(t_vals))

1. 最后，在图 3.1的右侧，我们使用y轴上的对数刻度绘制近似误差。然后，我们可以使用semilogy绘图命令在右侧使用对数刻度y轴，就像我们在第二章中看到的那样，使用 Matplotlib 进行数学绘图：

ax2.semilogy(t_vals, err)
ax2.set_xlabel("$t$")
ax2.set_ylabel("Error")
ax2.set_title("Error in approximation")

图 3.1中的左侧图显示随时间降低的温度，而右侧图显示随着我们远离初始条件给出的已知值，误差增加：

图 3.1：使用默认设置使用 solve_ivp 例程获得冷却方程的数值解的绘图

它是如何工作的…

解决微分方程的大多数方法都是“时间步进”方法。(t[i], y[i])对是通过采取小的t步骤并逼近函数y的值来生成的。这可能最好地通过欧拉方法来说明，这是最基本的时间步进方法。固定一个小的步长h > 0，我们使用以下公式在第i步形成近似值

从已知的初始值*y[0]*开始。我们可以轻松地编写一个执行欧拉方法的 Python 例程如下（当然，实现欧拉方法有许多不同的方法；这是一个非常简单的例子）：

1. 首先，通过创建将存储我们将返回的t值和y值的列表来设置方法：

def euler(func, t_range, y0, step_size):
    """Solve a differential equation using Euler's method"""
    t = [t_range[0]]
    y = [y0]
    i = 0

1. 欧拉方法一直持续到我们达到t范围的末尾。在这里，我们使用while循环来实现这一点。循环的主体非常简单；我们首先递增计数器i，然后将新的t和y值附加到它们各自的列表中。

while t[i] < t_range[1]:
        i += 1
        t.append(t[i-1] + step_size)  # step t
        y.append(y[i-1] + step_size*func(t[i-1], y[i-1]))  # step y
    return t, y

solve_ivp例程默认使用龙格-库塔-费尔伯格（RK45）方法，该方法能够自适应步长，以确保近似误差保持在给定的容差范围内。这个例程期望提供三个位置参数：函数f，应找到解的t范围，以及初始y值（在我们的例子中为T[0]）。可以提供可选参数来更改求解器、要计算的点数以及其他几个设置。

传递给solve_ivp例程的函数必须有两个参数，就像准备就绪部分中描述的一般微分方程一样。函数可以有额外的参数，可以使用args关键字为solve_ivp例程提供这些参数，但这些参数必须位于两个必要参数之后。将我们之前定义的euler例程与步长为 0.1 的solve_ivp例程进行比较，我们发现solve_ivp解的最大真实误差在 10^(-6)数量级，而euler解只能达到 31 的误差。euler例程是有效的，但步长太大，无法克服累积误差。

solve_ivp例程返回一个存储已计算解的信息的解对象。这里最重要的是t和y属性，它们包含计算解的t值和解y本身。我们使用这些值来绘制我们计算的解。y值存储在形状为(n, N)的 NumPy 数组中，其中n是方程的分量数（这里是 1），N是计算的点数。sol中的y值存储在一个二维数组中，在这种情况下有 1 行和许多列。我们使用切片y[0, :]来提取这个第一行作为一维数组，可以用来在步骤 4中绘制解。

我们使用对数缩放的y轴来绘制误差，因为有趣的是数量级。在非缩放的y轴上绘制它会得到一条非常靠近x轴的线，这不会显示出随着t值的变化误差的增加。对数缩放的y轴清楚地显示了这种增加。

还有更多…

solve_ivp例程是微分方程的多个求解器的便捷接口，默认为龙格-库塔-费尔伯格（RK45）方法。不同的求解器有不同的优势，但 RK45 方法是一个很好的通用求解器。

另请参阅

有关如何在 Matplotlib 中的图中添加子图的更详细说明，请参阅第二章中的添加子图示例，使用 Matplotlib 进行数学绘图。

解决微分方程系统

微分方程有时出现在由两个或更多相互关联的微分方程组成的系统中。一个经典的例子是竞争物种的简单模型。这是一个由以下方程给出的竞争物种的简单模型，标记为P（猎物）和W（捕食者）：*

*

第一个方程规定了猎物物种P的增长，如果没有任何捕食者，它将是指数增长。第二个方程规定了捕食者物种W的增长，如果没有任何猎物，它将是指数衰减。当然，这两个方程是耦合的；每种群体的变化都取决于两种群体。捕食者以与两种群体的乘积成比例的速率消耗猎物，并且捕食者以与猎物相对丰富度成比例的速率增长（再次是两种群体的乘积）。

在本示例中，我们将分析一个简单的微分方程系统，并使用 SciPy 的integrate模块来获得近似解。

准备就绪

使用 Python 解决常微分方程组的工具与解决单个方程的工具相同。我们再次使用 SciPy 中的integrate模块中的solve_ivp例程。然而，这只会给我们一个在给定起始种群下随时间预测的演变。因此，我们还将使用 Matplotlib 中的一些绘图工具来更好地理解这种演变。

如何做到…

以下步骤介绍了如何分析一个简单的常微分方程组：

1. 我们的第一个任务是定义一个包含方程组的函数。这个函数需要接受两个参数，就像单个方程一样，除了依赖变量y（在Solving simple differential equations numerically配方中的符号）现在将是一个包含与方程数量相同的元素的数组。在这里，将有两个元素。在这个配方中，我们需要的函数如下：

def predator_prey_system(t, y):
    return np.array([5*y[0] - 0.1*y[0]*y[1], 0.1*y[1]*y[0] -
       6*y[1]])

1. 现在我们已经在 Python 中定义了系统，我们可以使用 Matplotlib 中的quiver例程来生成一个图表，描述种群将如何演变——由方程给出——在许多起始种群中。我们首先设置一个网格点，我们将在这些点上绘制这种演变。选择相对较少的点数作为quiver例程的一个好主意，否则在图表中很难看到细节。对于这个例子，我们绘制种群值在 0 到 100 之间：

p = np.linspace(0, 100, 25)
w = np.linspace(0, 100, 25)
P, W = np.meshgrid(p, w)

1. 现在，我们计算每对这些点的系统值。请注意，系统中的任何一个方程都不是时间相关的（它们是自治的）；时间变量t在计算中并不重要。我们为t参数提供值0：

dp, dw = predator_prey_system(0, np.array([P, W]))

1. 现在变量dp和dw分别保存了种群P和W在我们的网格中每个点开始时将如何演变的“方向”。我们可以使用matplotlib.pyplot中的quiver例程一起绘制这些方向：

fig, ax = plt.subplots()
ax.quiver(P, W, dp, dw)
ax.set_title("Population dynamics for two competing species")
ax.set_xlabel("P")
ax.set_ylabel("W")

现在绘制这些命令的结果给出了图 3.2，它给出了解决方案演变的“全局”图像：

图 3.2：显示由常微分方程组模拟的两个竞争物种的种群动态的 quiver 图

为了更具体地理解解决方案，我们需要一些初始条件，这样我们就可以使用前面配方中描述的solve_ivp例程。

1. 由于我们有两个方程，我们的初始条件将有两个值。（回想一下，在Solving simple differential equations numerically配方中，我们看到提供给solve_ivp的初始条件需要是一个 NumPy 数组。）让我们考虑初始值P(0) = 85和W(0) = 40。我们在一个 NumPy 数组中定义这些值，小心地按正确的顺序放置它们：

initial_conditions = np.array([85, 40])

1. 现在我们可以使用scipy.integrate模块中的solve_ivp。我们需要提供max_step关键字参数，以确保我们在解决方案中有足够的点来得到平滑的解曲线：

from scipy import integrate
sol = integrate.solve_ivp(predator_prey_system, (0., 5.),
   initial_conditions, max_step=0.01)

1. 让我们在现有的图上绘制这个解，以展示这个特定解与我们已经生成的方向图之间的关系。我们同时也绘制初始条件：

ax.plot(initial_conditions[0], initial_conditions[1], "ko")
ax.plot(sol.y[0, :], sol.y[1, :], "k", linewidth=0.5)

这个结果显示在图 3.3中：

图 3.3：在显示一般行为的 quiver 图上绘制的解轨迹

它是如何工作的…

用于一组常微分方程的方法与单个常微分方程完全相同。我们首先将方程组写成一个单一的向量微分方程，

然后可以使用时间步进方法来解决，就好像y是一个简单的标量值一样。

使用quiver例程在平面上绘制方向箭头的技术是学习系统如何从给定状态演变的一种快速简单的方法。函数的导数代表曲线的梯度（x，u（x）），因此微分方程描述了解决方案函数在位置y和时间t的梯度。一组方程描述了在给定位置y和时间t的单独解决方案函数的梯度。当然，位置现在是一个二维点，所以当我们在一个点上绘制梯度时，我们将其表示为从该点开始的箭头，指向梯度的方向。箭头的长度表示梯度的大小；箭头越长，解决方案曲线在该方向上移动得越“快”。

当我们在这个方向场上绘制解的轨迹时，我们可以看到曲线（从该点开始）遵循箭头指示的方向。解的轨迹所显示的行为是一个极限环，其中每个变量的解都是周期性的，因为两种物种的人口增长或下降。如果我们将每个人口随时间的变化绘制成图，从图 3.4中可以看出，这种行为描述可能更清晰。从图 3.3中并不立即明显的是解的轨迹循环几次，但这在图 3.4中清楚地显示出来：

图 3.4：人口P和W随时间的变化图。两种人口都表现出周期性行为。

还有更多…

通过在变量之间绘制相空间（平面）分析系统的普通微分方程组的技术称为相空间（平面）分析。在这个示例中，我们使用quiver绘图例程快速生成了微分方程系统的相平面的近似值。通过分析微分方程系统的相平面，我们可以识别解的不同局部和全局特征，如极限环。

数值求解偏微分方程

偏微分方程是涉及函数在两个或多个变量中的偏导数的微分方程，而不是仅涉及单个变量的普通导数。偏微分方程是一个广泛的主题，可以轻松填满一系列书籍。偏微分方程的典型例子是（一维）热方程。

其中α是一个正常数，f（t，x）是一个函数。这个偏微分方程的解是一个函数u（t，x），它表示了在给定时间t>0 时，占据x范围 0≤x≤L的杆的温度。为了简单起见，我们将取f（t，x）=0，这相当于说系统没有加热/冷却，α=1，L=2。在实践中，我们可以重新调整问题来修复常数α，所以这不是一个限制性的问题。在这个例子中，我们将使用边界条件

这相当于说杆的两端保持在恒定温度 0。我们还将使用初始温度剖面

这个初始温度剖面描述了 0 和 2 之间的值之间的平滑曲线，峰值为 3，这可能是将杆在中心加热到 3 度的结果。

我们将使用一种称为有限差分的方法，将杆分成若干相等的段，并将时间范围分成若干离散步骤。然后我们计算每个段和每个时间步长的解的近似值。

在这个示例中，我们将使用有限差分来解一个简单的偏微分方程。

准备工作

对于这个示例，我们将需要 NumPy 包和 Matplotlib 包，通常导入为np和plt。我们还需要从mpl_toolkits导入mplot3d模块，因为我们将生成一个 3D 图：

from mpl_toolkits import mplot3d

我们还需要一些来自 SciPy 包的模块。

如何做…

在接下来的步骤中，我们将通过有限差分来解决热方程：

1. 首先创建代表系统物理约束的变量：杆的范围和α的值：

alpha = 1
x0 = 0 # Left hand x limit
xL = 2 # Right hand x limit

1. 我们首先将x范围分成N个相等的间隔—我们在这个例子中取N = 10*—使用N+1个点。我们可以使用 NumPy 中的linspace例程生成这些点。我们还需要每个间隔的公共长度h:

N = 10
x = np.linspace(x0, xL, N+1)
h = (xL - x0) / N

1. 接下来，我们需要设置时间方向上的步长。我们在这里采取了稍微不同的方法；我们设置时间步长k和步数（隐含地假设我们从时间 0 开始）：

k = 0.01
steps = 100
t = np.array([i*k for i in range(steps+1)])

1. 为了使方法正常运行，我们必须满足：

否则系统可能变得不稳定。我们将左侧存储在一个变量中，以便在步骤 4中使用，并使用断言来检查这个不等式是否成立：

r = alpha*k / h**2
assert r < 0.5, f"Must have r < 0.5, currently r={r}"

1. 现在我们可以构建一个矩阵，其中包含来自有限差分方案的系数。为此，我们使用scipy.sparse模块中的diags例程创建一个稀疏的三对角矩阵：

from scipy import sparse
diag = [1, *(1-2*r for _ in range(N-1)), 1]
abv_diag = [0, *(r for _ in range(N-1))]
blw_diag = [*(r for _ in range(N-1)), 0]
A = sparse.diags([blw_diag, diag, abv_diag], (-1, 0, 1), shape=(N+1,
      N+1), dtype=np.float64, format="csr")

1. 接下来，我们创建一个空白矩阵来保存解决方案：

u = np.zeros((steps+1, N+1), dtype=np.float64)

1. 我们需要将初始配置添加到第一行。这样做的最佳方法是创建一个保存初始配置的函数，并将在刚刚创建的矩阵u上的x数组上评估这个函数的结果：

def initial_profile(x):
    return 3*np.sin(np.pi*x/2)
u[0, :] = initial_profile(x)

1. 现在我们可以简单地循环每一步，通过将A和前一行相乘来计算矩阵u的下一行：

for i in range(steps):
    u[i+1, :] = A @ u[i, :]

1. 最后，为了可视化我们刚刚计算的解，我们可以使用 Matplotlib 将解作为曲面绘制出来：

X, T = np.meshgrid(x, t)
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
ax.plot_surface(T, X, u, cmap="hot")
ax.set_title("Solution of the heat equation")
ax.set_xlabel("t")
ax.set_ylabel("x")
ax.set_zlabel("u")

这导致了图 3.5中显示的曲面图：

图 3.5：使用有限差分方法计算的热方程解在范围 0 ≤ x ≤ 2 上的曲面图，使用了 10 个网格点

它是如何工作的…

有限差分方法通过用仅涉及函数值的简单分数替换每个导数来工作，我们可以估计这些分数。要实现这种方法，我们首先将空间范围和时间范围分解为若干离散间隔，由网格点分隔。这个过程称为离散化。然后我们使用微分方程和初始条件以及边界条件来形成连续的近似，这与solve_ivp例程在Solving differential equations numerically示例中使用的时间步进方法非常相似。

为了解决热方程这样的偏微分方程，我们至少需要三个信息。通常，对于热方程，这将以空间维度的边界条件的形式出现，告诉我们在杆的两端行为如何，以及时间维度的初始条件，即杆上的初始温度分布。

前面描述的有限差分方案通常被称为前向时间中心空间（FTCS）方案，因为我们使用前向有限差分来估计时间导数，使用中心有限差分来估计（二阶）空间导数。这些有限差分的公式如下：

和

将这些近似代入热方程，并使用近似值u[i]^(j)来表示在i空间点上经过j时间步后u(t[j], x[i])的值，我们得到

可以重新排列以获得公式

粗略地说，这个方程表示给定点的下一个温度取决于以前时间的周围温度。这也显示了为什么r值的条件是必要的；如果条件不成立，右侧的中间项将是负的。

我们可以将这个方程组写成矩阵形式，

其中u*^{j*是包含近似*u[i]}j和矩阵A的向量，矩阵A在步骤 4中定义。这个矩阵是三对角的，这意味着非零条目出现在或邻近主对角线上。我们使用 SciPy sparse模块中的diag例程，这是一种定义这种矩阵的实用程序。这与本章中解方程配方中描述的过程非常相似。这个矩阵的第一行和最后一行都是零，除了在左上角和右下角，分别代表（不变的）边界条件。其他行的系数由微分方程两侧的有限差分近似给出。我们首先创建对角线条目和对角线上下方的条目，然后我们使用diags例程创建稀疏矩阵。矩阵应该有N+1*行和列，以匹配网格点的数量，并且我们将数据类型设置为双精度浮点数和 CSR 格式。

初始配置给我们向量u⁰，从这一点开始，我们可以通过简单地执行矩阵乘法来计算每个后续时间步骤，就像我们在步骤 7中看到的那样。

还有更多…

我们在这里描述的方法相当粗糙，因为近似可能变得不稳定，正如我们提到的，如果时间步长和空间步长的相对大小没有得到仔细控制。这种方法是显式的，因为每个时间步骤都是显式地使用来自上一个时间步骤的信息来计算的。还有隐式方法，它给出了一个可以求解以获得下一个时间步骤的方程组。不同的方案在解的稳定性方面具有不同的特性。

当函数f(t, x)不为 0 时，我们可以通过使用赋值来轻松适应这种变化

其中函数被适当地向量化以使这个公式有效。在解决问题的代码方面，我们只需要包括函数的定义，然后将解决方案的循环更改如下：

for i in range(steps):
    u[i+1, :] = A @ u[i, :] + f(t[i], x)

从物理上讲，这个函数代表了杆上每个点的外部热源（或者热池）。这可能随时间变化，这就是为什么一般来说，这个函数应该有t和x作为参数（尽管它们不一定都被使用）。

我们在这个例子中给出的边界条件代表了杆的两端保持在恒定温度为 0。这种边界条件有时被称为Dirichlet边界条件。还有Neumann边界条件，其中函数u的导数在边界处给定。例如，我们可能已经给出了边界条件

这在物理上可以被解释为杆的两端被绝缘，因此热量不能通过端点逃逸。对于这种边界条件，我们需要稍微修改矩阵A，但方法本质上保持不变。实际上，在边界的左侧插入一个虚拟的x值，并在左边界（x = 0）使用向后有限差分，我们得到

使用这个二阶有限差分近似，我们得到

这意味着我们矩阵的第一行应包含1-r，然后是r，然后是 0。对右手极限进行类似的计算得到矩阵的最后一行。

diag = [1-r, *(1-2*r for _ in range(N-1)), 1-r]
abv_diag = [*(r for _ in range(N))]
blw_diag = [*(r for _ in range(N))]
A = sparse.diags([blw_diag, diag, abv_diag], (-1, 0, 1), shape=(N+1, N+1), dtype=np.float64, format="csr")

对于涉及偏微分方程的更复杂问题，可能更适合使用有限元求解器。有限元方法使用比我们在本篇中看到的有限差分方法更复杂的方法来计算解决方案，通常比偏微分方程更灵活。然而，这需要更多依赖更高级数学理论的设置。另一方面，有一个用于使用有限元方法解决偏微分方程的 Python 包，如 FEniCS (fenicsproject.org)。使用 FEniCS 等包的优势在于它们通常针对性能进行了调整，这在解决复杂问题时对高精度至关重要。

另请参阅

FEniCS 文档介绍了有限元方法以及使用该包解决各种经典偏微分方程的示例。有关该方法和理论的更全面介绍可在以下书籍中找到：

• Johnson, C. (2009).Numerical solution of partial differential equations by the finite element method. Mineola, N.Y.: Dover Publications.

有关如何使用 Matplotlib 生成三维曲面图的更多详细信息，请参阅第二章中的曲面和等高线图食谱。

使用离散傅立叶变换进行信号处理

来自微积分中最有用的工具之一是傅立叶变换。粗略地说，傅立叶变换以可逆的方式改变了某些函数的表示。这种表示的改变在处理作为时间函数的信号时特别有用。在这种情况下，傅立叶变换将信号表示为频率函数；我们可以将其描述为从信号空间到频率空间的转换。这可以用于识别信号中存在的频率以进行识别和其他处理。在实践中，我们通常会有信号的离散样本，因此我们必须使用离散傅立叶变换来执行这种分析。幸运的是，有一种计算效率高的算法，称为快速傅立叶变换（FFT），用于对样本应用离散傅立叶变换。

*我们将遵循使用 FFT 处理嘈杂信号的常见过程。第一步是应用 FFT 并使用数据计算信号的功率谱密度。然后我们识别峰值并滤除对信号贡献不足够大的频率。然后我们应用逆 FFT 来获得滤波后的信号。

在本篇中，我们使用 FFT 来分析信号的样本，并识别存在的频率并清除信号中的噪音。

准备工作

对于本篇，我们只需要导入 NumPy 和 Matplotlib 包，如往常一样命名为np和plt。

如何做…

按照以下说明使用 FFT 处理嘈杂信号：

1. 我们定义一个函数来生成我们的基础信号：

def signal(t, freq_1=4.0, freq_2=7.0):
    return np.sin(freq_1 * 2 * np.pi * t) + np.sin(freq_2 * 2 *
        np.pi * t)

1. 接下来，我们通过向基础信号添加一些高斯噪声来创建我们的样本信号。我们还创建一个数组，以便在以后方便时保存样本t值处的真实信号：

state = np.random.RandomState(12345)
sample_size = 2**7 # 128
sample_t = np.linspace(0, 4, sample_size)
sample_y = signal(sample_t) + state.standard_normal(sample_size)
sample_d = 4./(sample_size - 1) # Spacing for linspace array
true_signal = signal(sample_t)

1. 我们使用 NumPy 的fft模块来计算离散傅立叶变换。在开始分析之前，我们从 NumPy 中导入这个模块：

from numpy import fft

1. 要查看嘈杂信号的样子，我们可以绘制样本信号点并叠加真实信号：

fig1, ax1 = plt.subplots()
ax1.plot(sample_t, sample_y, "k.", label="Noisy signal")
ax1.plot(sample_t, signal(sample_t), "k--", label="True signal")
ax1.set_title("Sample signal with noise")
ax1.set_xlabel("Time")
ax1.set_ylabel("Amplitude")
ax1.legend()

此处创建的图表显示在图 3.6中。正如我们所看到的，嘈杂信号与真实信号几乎没有相似之处（用虚线表示）：

图 3.6：带真实信号的噪声信号样本

1. 现在，我们将使用离散傅立叶变换来提取样本信号中存在的频率。fft模块中的fft例程执行 FFT（离散傅立叶变换）：

spectrum = fft.fft(sample_y)

1. fft模块提供了一个用于构建适当频率值的例程，称为fftfreq。为了方便起见，我们还生成一个包含正频率出现的整数的数组：

freq = fft.fftfreq(sample_size, sample_d)
pos_freq_i = np.arange(1, sample_size//2, dtype=int)

1. 接下来，计算信号的功率谱密度（PSD）如下：

psd = np.abs(spectrum[pos_freq_i])**2 + np.abs(spectrum[-
        pos_freq_i])**2

1. 现在，我们可以绘制信号的正频率的 PSD，并使用这个图来识别频率：

fig2, ax2 = plt.subplots()
ax2.plot(freq[pos_freq_i], psd)
ax2.set_title("PSD of the noisy signal")
ax2.set_xlabel("Frequency")
ax2.set_ylabel("Density")

结果可以在图 3.7中看到。我们可以在这个图表中看到大约在 4 和 7 左右有尖峰，这些是我们之前定义的信号的频率：

图 3.7：使用 FFT 生成的信号的功率谱密度

1. 我们可以识别这两个频率，尝试从噪声样本中重建真实信号。所有出现的次要峰值都不大于 10,000，所以我们可以将其作为滤波器的截止值。现在，我们从所有正频率索引的列表中提取（希望是 2 个）对应于 PSD 中大于 10,000 的峰值的索引：

filtered = pos_freq_i[psd > 1e4]

1. 接下来，我们创建一个新的干净频谱，其中只包含我们从噪声信号中提取的频率。我们通过创建一个只包含 0 的数组，然后将spectrum的值从对应于滤波频率及其负值的索引复制过来来实现这一点：

new_spec = np.zeros_like(spectrum)
new_spec[filtered] = spectrum[filtered]
new_spec[-filtered] = spectrum[-filtered]

1. 现在，我们使用逆 FFT（使用ifft例程）将这个干净的频谱转换回原始样本的时间域。我们使用 NumPy 的real例程取实部以消除错误的虚部：

new_sample = np.real(fft.ifft(new_spec))

1. 最后，我们绘制这个滤波信号和真实信号进行比较：

fig3, ax3 = plt.subplots()
ax3.plot(sample_t, true_signal, color="#8c8c8c", linewidth=1.5, label="True signal")
ax3.plot(sample_t, new_sample, "k--", label="Filtered signal")
ax3.legend()
ax3.set_title("Plot comparing filtered signal and true signal")
ax3.set_xlabel("Time")
ax3.set_ylabel("Amplitude")

步骤 12的结果显示在图 3.8中。我们可以看到，滤波信号与真实信号非常接近，除了一些小的差异：

图 3.8：比较使用 FFT 和滤波生成的滤波信号与真实信号的图

工作原理…

函数f(t)的傅立叶变换由积分给出

离散傅立叶变换如下：

这里，f[k]值是复数形式的样本值。可以使用前面的公式计算离散傅立叶变换，但实际上这并不高效。使用这个公式计算的复杂度是O(N²)。FFT 算法将复杂度提高到O(N log N)，这显然更好。书籍Numerical Recipes（在进一步阅读部分给出完整的参考文献细节）对 FFT 算法和离散傅立叶变换有很好的描述。

我们将对从已知信号（具有已知频率模式）生成的样本应用离散傅立叶变换，以便我们可以看到我们获得的结果与原始信号之间的联系。为了保持这个信号简单，我们创建了一个只有两个频率分量（值为 4 和 7）的信号。从这个信号，我们生成了我们分析的样本。由于 FFT 的工作方式，最好是样本的大小是 2 的幂；如果不是这种情况，我们可以用零元素填充样本使其成为这种情况。我们向样本信号添加了一些高斯噪声，这是一个正态分布的随机数。

fft例程返回的数组包含N+1个元素，其中N是样本大小。索引为 0 的元素对应于 0 频率，或者直流偏移。接下来的N/2个元素是对应于正频率的值，最后的N/2个元素是对应于负频率的值。频率的实际值由采样点数N和采样间距确定，在这个例子中，采样间距存储在sample_d中。

频率ω处的功率谱密度由以下公式给出

其中H(ω)代表信号在频率ω处的傅里叶变换。功率谱密度测量了每个频率对整体信号的贡献，这就是为什么我们在大约 4 和 7 处看到峰值。由于 Python 索引允许我们对从序列末尾开始的元素使用负索引，我们可以使用正索引数组从spectrum中获取正频率和负频率元素。

在步骤 9中，我们提取了在图表上峰值超过 10,000 的两个频率的索引。对应于这些索引的频率分别是 3.984375 和 6.97265625，它们并不完全等于 4 和 7，但非常接近。这种差异的原因是我们使用有限数量的点对连续信号进行了采样。（使用更多点当然会得到更好的近似。）

在步骤 11中，我们提取了从逆 FFT 返回的数据的实部。这是因为从技术上讲，FFT 处理复杂数据。由于我们的数据只包含实数据，我们期望这个新信号也只包含实数据。然而，会有一些小错误产生，这意味着结果并非完全是实数。我们可以通过取逆 FFT 的实部来纠正这一点。这是合适的，因为我们可以看到虚部非常小。

我们可以在图 3.8中看到，滤波信号非常接近真实信号，但并非完全相同。这是因为，如前所述，我们正在用相对较小的样本来逼近连续信号。

还有更多…

在生产环境中的信号处理可能会使用专门的软件包，比如scipy中的signal模块，或者一些更低级的代码或硬件来执行信号的滤波或清理。这个示例更多地应该被看作是 FFT 作为处理从某种基础周期结构（信号）采样的数据的工具的演示。FFT 对于解决偏微分方程非常有用，比如在数值解偏微分方程食谱中看到的热方程。

另请参阅

有关随机数和正态分布（高斯）的更多信息可以在第四章中找到，处理随机性和概率。

进一步阅读

微积分是每门本科数学课程中非常重要的一部分。有许多关于微积分的优秀教科书，包括 Spivak 的经典教科书和 Adams 和 Essex 的更全面的课程：

• Spivak, M. (2006). Calculus. 3rd ed. Cambridge: Cambridge University Press
  
• Adams, R. and Essex, C. (2018). Calculus: A Complete Course. 9th ed. Don Mills, Ont: Pearson.Guassian
  

关于数值微分和积分的良好来源是经典的数值方法书，其中详细描述了如何在 C++中解决许多计算问题的理论概述：

• Press, W., Teukolsky, S., Vetterling, W. and Flannery, B. (2007). Numerical recipes: The Art of Scientific Computing. 3rd ed. Cambridge: Cambridge University Press****