如何在 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)，来创建二维矩形函数。为了便于学习，我们将创建和定义一个二维网格数据集，并通过手动绘制而不是自动绘制。结果如下列图像所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

结语

在 COMSOL Multiphysics 中，我们可以在主计算之前或之后将数据集功能和integrate运算符，用作便捷的独立计算工具以及前处理和后处理工具。请注意，此处讨论的傅立叶变换不是离散傅立叶变换（FFT）。我们仍然使用离散数学，但是使用辛普森法则进行数值积分。该函数在 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.