如何在 COMSOL Multiphysics 中实现傅里叶变换

在之前的博客中，我们讨论了如何模拟用于全息数据存储的聚焦激光束。在一个具体的示例中，通过对透镜入口处的电磁场振幅进行傅里叶变换得到由傅里叶透镜聚焦的电磁波。这篇博客，我们将通过一个夫琅禾费衍射示例，演示如何在 COMSOL Multiphysics 中实现带有积分运算的前处理和后处理。

通过夫琅禾费衍射示例理解傅里叶变换

傅里叶变换在许多仿真应用场景中都发挥着重要作用。除了傅里叶光学外，我们还再夫琅禾费衍射理论，信号处理中的频谱信息提取，图像处理中的降噪与过滤等方面使用了傅里叶变换。

在本文的示例中，我们模拟了一个交通灯光穿过网格状纱帘后的成像，如下图所示。为了简化模型，我们假设光是电场强度均匀的平面波，例如强度为 1V/m。假设纱帘是垂直于光传播方向的平面，将纱帘几何进行网格离散化，并用局部坐标 x 和 y 进行度量；在靠近眼睛的位置取一平行于纱帘的平面，通过局部坐标 u 和 v 测量成像图案。

用作纱帘中方形孔径的傅里叶变换的夫琅禾费衍射图。

根据夫琅禾费衍射理论，我们可以简单地通过传递函数的傅里叶变换来计算上述图像；如果网格是正方形的，则该传递函数是一个周期性的矩形函数。接下来，我们考虑一种简化情况，即传递函数为矩形函数的单一网格。我们将在后文中讨论周期性传递函数的情况。

当光线照射到一个正方形网格，并在网格中心透射时将被织物的锐边衍射，我们用二维矩形函数来描述传递函数。通过在 COMSOL Multiphysics 仿真中引入傅里叶变换，可以更全面地了解这个过程。

利用 COMSOL Multiphysics 数据集

为了学习如何实现傅里叶变换，我们需要先讨论数据集 或多维矩阵数据存储的概念。COMSOL Multiphysics 中有两种可能的数据集类型：解 和栅格。对于任何计算，COMSOL 软件都会创建一个位于结果 >数据集 节点下的数据集。

解数据集包含非结构化栅格，用于存储解数据。为了利用此数据集，指定每一列和每一行所对应的数据。如果我们指定解1（sol1），则矩阵维数对应于研究1 中模型的维数。例如，如果是瞬态问题，数据集有三维数组，可以写成 T(i,j,k) ，i=1,\cdots, N_t, \ j=1, \cdots, N_n, \ k = 1, \cdots, N_s   。其中, N_t 是存储的时步， N_n 是节点数，N_s 是空间维度。类似地，瞬态参数化研究的数据集由四维数组组成。同样，请注意，空间数据（时间和参数数据除外）与网格上的节点位置上存储的节点数据（不一定在规则栅格上）链接。

另一方面，栅格数据集具有规则栅格，用于函数和所有其他通用用途。存储在栅格数据集中的所有数值都链接到设置窗口定义的栅格中。在定义 节点中定义一个函数并单击 创建绘图 时，会自动创建此数据集。此操作将在数据集 节点中创建一维栅格数据集。

另外，我们还需要指定自变量的范围和分辨率。默认情况下，一维栅格数据集的分辨率设置为 1000。如果自变量（即 x）的范围从 0 到 1，则栅格数据集将准备 0, 0.001, 0.002，…，0.999 和 1 的数据序列。 对于二维和三维数据集，默认的分辨率分别为 100 和 30。对于傅里叶变换，使用栅格数据集。我们也可以将此数据集用作计算的独里工具，因为它没有强制指向解数据集。

实现傅里叶变换

在开始仿真前，先定义内置的一维矩形函数，如下图所示。

定义内置的一维矩形函数。

然后，在设置窗口中单击创建绘图 按钮，在结果 节点中创建一个单独的一维绘图组。

内置的一维矩形函数绘图。

下图为设置窗口的绘图。扩展一维绘图组 1 节点，然后单击线图 1 查看指向一维栅格 的数据集。在一维栅格节点设置中，可以看到数据集与 rect1 函数相关联。

设置内置的一维矩形函数。

设置一维栅格数据集。

我们可以通过在 定义 节点下定义一个解析函数如 rect1(x)*rect1(y)，来创建二维矩形函数。为了便于学习，我们将创建和定义一个二维网格数据集，并通过手动绘制而不是自动绘制。结果如下文中的图像所示。

在二维栅格设置中，选择 全部函数 作为数据源，因为二维矩形函数还引用了另一个函数 rect1。将 x 和 y （之前定义的纱帘的局部坐标）设置为第一、第二参数，并将分辨率设置为 64 以加速测试。设置完毕后将其重命名为二维栅格（源空间），并在绘图组设置窗口，选择该数据集绘制结果。 在二维栅格设置中定义函数。

在二维栅格设置中定义函数。

创建和定义一个二维数据集。

为二维矩形函数设置二维绘图组。

二维矩形函数的二维绘图。

现在，我们通过计算下式来实现该函数的傅里叶变换：

其中，如我们之前讨论的 u 和 v 代表目标空间（傅里叶/频率空间）的独立变量。

由于已经为 x 和 y 创建了二维数据集，如下图所示，现在我们可以为 u 和 v 创建一个二维栅格数据集 [被重命名为二维栅格(目标空间)]。从 源 中选择 函数，并从 函数 中选择 全部，因为函数 rect 也会调用函数 rect1。为了更快的计算，我们可以像处理二维数据集一样，在此处将分辨率更改为 64。

为傅里叶空间设置二维栅格数据集。

现在，我们处于仿真阶段，可以使用 integrate 算子写入方程。

输入二维矩形函数的傅里叶变换方程。

如下图所示，我们最终得到了傅里叶变换后的结果，将其（更准确地说是其正方形）与纱帘照片中的每个闪烁彩色光进行比较。实际上，我们还没有真正看到透射后的图像。要计算出眼睛可以看到的图像（最终目标），我们需要再进行一次傅里叶变换。

二维矩形函数的傅里叶变换。

结语

在 COMSOL Multiphysics 中，我们可以在主计算之前或之后将数据集功能和 integrate 算子，用作方便快捷的独立计算工具以及前处理和后处理工具。请注意，此处讨论的傅里叶变换不是离散傅里叶变换（FFT）。我们仍然使用离散数学，但是使用辛普森法则（Simpson’s rule）进行数值积分。该函数在 COMSOL Multiphysics 的积分算子中使用，而离散傅里叶变换是通过对数字序列进行运算而形成的。因此，不必担心混叠问题，傅里叶空间分辨率问题或傅里叶空间偏移问题。

关于这个主题还有更多的问题需要讨论，先分析一下之前简化的两种情况。我们进行了单个网格计算，在实际应用中，纱帘由有限数量的周期性正方形开口组成。这看似必须针对周期性情况重新进行计算，但幸运的是，最终结果仅与周期性相关的包络函数而有所不同。有关详细信息，Hecht 在光学 模块对该主题进行了很好地概述。

第二种简化是假定网格传递函数为非平滑矩形函数。在 COMSOL Multiphysics 中，出于数值稳定性和准确性的原因，除用户定义的函数外，所有函数都会被做某种程度的平滑处理。您可能已经注意到，我们的矩形函数的斜率很小。这是一种复杂而不是简化，因为最简单的情况是一个斜率无穷大的矩形函数，并且我们对矩形函数进行了平滑处理。

我们知道，有两种极端情况的傅里叶变换。即，将斜率无穷大的矩形函数转换为辛格（sinc）函数（sin(x)/x），以及将一个高斯函数转换为另一个高斯函数。辛格函数在中心周围有波动，代表衍射效应，而高斯函数在衰减时没有任何波动。我们的平滑矩形函数处于这两个极端之间，因此它的傅里叶变换也在辛格函数和高斯函数之间。正如我们前文所提到的，纱帘织物不能有锐边，因此对于此示例，我们的结果可能更准确。

延伸阅读

• 查看下列有关全息数据存储系统模拟的博客：
  • 在 COMSOL Multiphysics 中模拟全息数据存储
  • 如何模拟全息页面数据存储系统
• 阅读下列光学入门书籍，了解更多信息：
  • J.W. Goodman, Introduction to Fourier Optics, W. H. Freeman, 2004.
  • E. Hecht, Optics, Pearson Education Limited, 2014.