如何在光学应用中使用新的空间 FFT 特征

2022年 8月 19日

快速傅里叶变换 (FFT) 是一种有用并且强大的数值方法。COMSOL Multiphysics ® 软件的最新 6.0 版本增加了与此方法相关的新特征:空间 FFT 特征。在这篇博客中,我们将讨论如何将这一新功能用于光学应用并展示了一些应用案例。

术语和定义

首先,我们来明确一些术语和定义。需要区分三个术语:傅里叶变换 (FT)、离散傅里叶变换 (DFT) 和 FFT。一个函数的 FT定义为

\hat{u}(\xi) = \int_{-\infty}^{\infty} u(x)e^{-2\pi i \xi x} dx,

式中,x\xi 分别是物理空间和傅里叶空间中的变量。当物理空间变量为时间 t 时,变量 \xi 称为频率。在光学中,\xi 被称为空间频率,通常与波长和焦距成比例(我们将在后面讨论),而 x 是用于描述感兴趣的光学结构附近位置的物理空间坐标。

在之前的博客 “如何在 COMSOL Multiphysics 中实现傅里叶变换”和“如何由计算解实现傅里叶变换”中,我们讨论了如何在 COMSOL ®中进行傅里叶变换。我们可以使用一种数值技术,通过基于辛普森法则的直接数值积分来执行傅里叶变换公式。在本文后面的部分,我们可以将其称为“数值积分的 FT”。

DFT 是 FT 的离散版本,它是对一组离散的点进行运算。在 COMSOL ®中,它被定义为

\hat{u}(\xi_k) = \sum_{j=0}^{N-1}u(x_j) e^{-2\pi ijk/N}, \ \ \ k=0, \cdots, N-1

FFT 是一种计算 DFT 的有效算法。

请注意,这些 FT 和 DFT 的定义是最通用的定义,但符号约定与 COMSOL 的波动方程符号约定不一致,即 \exp(i(\omega t-{\mathbf k}\cdot{\mathbf x}))。当使用这些符号定义弗劳恩霍夫和菲涅耳衍射公式时,请注意不要弄错。符号不一致性不影响稳态解。

如何使用空间 FFT 特征

接下来,我们将演示如何在 COMSOL® 中将新的空间 FFT 特征 用于光学。可以分别使用步骤 1 和 2 设置和执行 FFT 特征:

  • 步骤1:准备数据集
    • 通过右键单击 数据集 → 更多数据集 添加空间 FFT 数据(这定义了傅立叶空间)
    • 选择合适的源数据集作为物理空间,然后进行变换
    • 空间分辨率 设置为手动
    • 采样分辨率 设置为适当的数字
    • 空间布局 中选择使用补零 并将 x 填充 设置为适当的数字
    • 傅里叶空间变量 中选择频率
    • 取消选中屏蔽 DC
  • 步骤2:使用 fft() 算子 绘图
    • 在绘图设置中调整 x 轴数据的空间频率缩放

矩形函数示例

矩形函数是光学中最常用的函数之一,因为它代表了一个硬边光阑。当有硬边光阑时,总是涉及矩形函数的 FT。矩形函数的 FT 可以很容易地手动计算出来,如下所示:

{\rm rect}(x/a)=
\begin{cases}
0 & |x|>a/2 \\
1 & |x|\le a/2
\end{cases}

 

{\mathcal F}[ {\rm rect(x/a)}](f) = a\:{\rm sinc}(\pi f a),

式中,\mathcal F 代表傅里叶变换算子,a 是一个常数,{\rm sinc}(x) = \sin x/x 是正弦函数。

让我们看看如何使用 COMSOL® 新的空间 FFT 特征 计算该 FT 。

左侧是模型构建器的屏幕截图,其中选择了 Grid 1D 节点和相应的 Settings 窗口,其中 Data、Parameter Bounds 和 Grid 部分展开。 右侧是模型构建器的屏幕截图,其中选择了 Spatial FFT 节点和相应的 Settings 窗口,其中展开了 Data 和 Transformation 部分。
矩形函数(左)及其 FT (右)的数据集设置示例。

矩形函数是定义 > 函数 下的内置函数。如果点击创建绘图 按钮,就会在结果下的数据库 节点下自动为该函数创建一个新数据集。 默认情况下,范围和分辨率也是自动设置的。执行 FFT 时,自己控制这些参数很重要。傅里叶空间分辨率由物理空间范围的倒数和物理空间数据的零填充确定。傅里叶空间范围由物理空间范围和傅里叶空间采样数决定。FFT 结果的大小因物理空间范围和傅里叶空间采样数而异。下表是 FFT 参数表达式的汇总,包括与上图中显示的 FFT 设置对应的参数值。

参数 表达式 示例值
实际上的总范围 X_0 2
傅里叶空间采样数 N_x 16
补零 N_z 8
傅里叶空间总范围 N_x/X_0 8
傅里叶空间分辨率 \frac{N_x}{X_0} \cdot \frac{1}{N_x+2N_z} 1/4
FT 归一化因子 X_0/N_x 1/8

通过上述设置,下图是矩形函数 rect1(x) ,其 FT的 绝对值 abs(fft(rect1(x)) 由 FFT 特征计算。傅里叶空间总范围是 N_x/X_0 = 16/2 = 8,即从 -4 到 4。可以看到傅里叶空间的采样点总数为 N_x+2N_z = 32。

为什么是 2N_z? 那是因为在补零中,N_z 零被添加到物理空间数据的两侧。傅里叶空间分辨率为 8/32 = 0.25。在没有归一化的情况下,FFT 运算结果因子为 N_x/X_0。所以,我们需要将结果乘以得到一个单位峰值。稍后,我们将对各种公式使用 FFT,每个公式都有不同的乘法常数。因此,我们必须将 FFT 结果归一化。

显示矩形函数的折线图,其中 a = 1(蓝线)及其 FT 的绝对值(红线)
矩形函数 a = 1 及其 FT 的绝对值,由 FFT 和前面显示的设置确定。

在这个示例中,采样数被有意设置为较低的数字,以便可以参考前面的公式。不过,我们可以看到 FT,{\rm sinc}(\pi \xi) 是用一个很好的近似值计算出来的。使用更合适的参数(例如 X_0 = 3,N_x = 128,N_z = 512),我们得到以下结果,这是完美的。将数值积分的 FT 结果叠加以进行比较。当然,这两种方法应该匹配!

将矩形函数与 a = 1(蓝线)、其 FT 的绝对值(红线)和数值积分的 FT(绿线)进行比较的折线图。
a = 1时,由高分辨率的 FFT 确定的矩形函数的 FT 绝对值和由数值积分确定的 FT 绝对值的对比

在光学中应用FT

现在,我们已经学习了如何为矩形函数(一维分析函数)设置和使用空间 FFT 特征。接下来,让我们将它用于光学中的一些实际应用示例。

在光学中,将光电场的时间信号与其光谱(频率或波长)相关联的时频傅里叶变换可能更为人所知。空间 FT 用于电场从一个平面传播到另一个平面的各种传播(变换)方法。在这个例子中,空间傅里叶变换将一个平面中电场的空间形状与另一个平面中的形状(称为空间频率)联系起来。让我们考虑一个标量电场或矢量电场的一个分量,它入射到一个平面中的扰动上,例如光圈或镜头,然后到达另一个平面,例如焦平面或像平面,如下图所示:

光学中传播方法的坐标系示意图。
光学中传播(变换)方法的坐标系。

让我们来表征扰动后平面内的电场 E(x,y,0)。然后,根据不同的目标,使用四种传播方法之一计算了另一个平面的电场 \hat{E}(x’,y’,L)。下表总结了四种方法。这些公式由 FT 的简单相位函数符号 \mathcal{F}[\cdot] 表示。

理论 公式(简单符号) 应用
1. 夫琅禾费衍射理论 \hat{E}(x’) = {\mathcal F}[E(x)] 夫琅禾费衍射条件下的标量远场——观察者距离衍射物体*很远,用于孔径、光栅和傅里叶光学等应用。
2. 菲涅耳衍射理论 \hat{E}(x’) = {\mathcal F}[E(x)\exp(-ikx^2/(2L))] 菲涅耳衍射条件**下的标量近场至远场,适用于低数值孔径 (NA) 透镜系统等应用。
3. 角谱法*** \hat{E}(x’) ={\mathcal F}^{-1}[{\mathcal F}[E(x)]\exp(ik\sqrt{1-\alpha^2})] 适用于任何系统(例如高数值孔径透镜系统)的严格单向标量场解决方案(不考虑反射)。
4. 部分相干理论(Schell 模型)**** I_{pc}(x’,L) = {\mathcal F}^{-1}[{\mathcal F}[I_c(x)] {\mathcal F}[\gamma(x)]^\ast ] 非干扰或低干扰光源,例如 LED 和太阳光,使用在在夫琅禾费或菲涅耳衍射近似下的互相干函数的 Schell 模型假设。

脚注:

夫琅禾费衍射条件
** 菲涅耳衍射条件
*** \alpha 是方向余弦
**** I_{pc}是部分相干强度,I_c 是相干强度,并且 \gamma 是互相干函数

夫琅禾费衍射

夫琅禾费衍射公式用于计算满足夫琅禾费条件时,从物体衍射的远场。

以下是完整的公式:

\hat{E}
(x’,y’,L) = \frac{e^{i k L}}{i \lambda L} \iint_{-\infty}^{\infty}E(x,y,0) e^{-i 2\pi (x’ x+y’ y) / (\lambda L)}dxdy

该公式用于计算孔径、光栅的远场和傅里叶光学焦平面内的场(参考文献 1)。物体是一个具有均匀照明的方形孔径。孔径出口平面的电场是一个二维矩形函数,远场由 FFT 计算。这形成了一种衍射图案,可能是一个熟悉的景象,类似于从网状窗帘后面观察路灯时的样子。请注意,需要将图中的 x 轴数据缩放为 \lambda L,因为空间频率被缩放为 f=x’/\lambda L。 通过数值积分使用 FT 计算二维 FT 需要很长时间,但 FFT 可以非常快速地完成这项工作。

左侧为方形孔径的图像,右侧为其衍射图案。
方形孔径(左)及其衍射图案(右)。

模拟方形孔径的设置非常简单:

 

菲涅耳衍射

第二个应用,菲涅耳衍射公式,可用于计算远场以及近场干扰。这个近似值的完整公式是:

\hat{E}
(x’,y’,L) = \frac{e^{ikL}}{i\lambda L}e^{ik(x’^2+y’^2)/(2L)}\iint_{-\infty}^{\infty}E(x,y,0)e^{-ik(x^2+y^2)/(2L)} e^{-i 2\pi (x’ x +y’ y)/ (\lambda L)}dxdy

请注意,x 轴数据需要按因子 \lambda L 进行缩放。这个应用在菲涅耳透镜模型中实现,其中使用了 波动光学,频域 接口计算透镜内部的电场。使用菲涅耳衍射公式通过数值积分由傅里叶变换计算焦平面中的场。如下图所示,可以在该模型中使用 FFT 特征,并通过数值积分得到与 FT 相同的结果。

显示菲涅耳模型焦平面中的电场模的折线图,宝蓝色线代表积分菲涅耳近似,绿线代表亥姆霍兹方程(ewfd),红线代表亥姆霍兹方程(ewbe)和代表 FFT 菲涅耳近似的水线。
菲涅耳透镜模型焦平面中的电场模。FFT 特征用于计算菲涅耳衍射公式,并与其他方法进行比较。

A screenshot of the COMSOL Multiphysics UI showing the Line Graph settings window with the Label of FFT Fresnel Approximation. The Data, y-Axis Data, and x-Axis Data sections are all expanded.
菲涅耳衍射公式 FFT 特征的后处理设置。请注意,y 轴数据是标准化的,x 轴数据是按比例缩放的。

角谱法

第三个应用,角谱法,实现起来有点麻烦,因为它需要使用两次 FT,就像它的完整公式:

\hat{E}(x’,y’,L) = \iint_{-\infty}^{\infty}
A \left(\frac{\alpha}{\lambda},\frac{\beta}{\lambda},0\right) e^{ikL\sqrt{1-\alpha^2-\beta^2}} e^{i2\pi (\alpha x+\beta y) /\lambda} d\frac{\alpha} {\lambda} d\frac{\beta}{\lambda}
,

式中,

A \left(\frac{\alpha} {\lambda},\frac{\beta}{\lambda}
,0\right) = \iint_{-\infty}^{\infty} E(x,y,0) e^{-i2\pi (\alpha x + \beta y)/\lambda} dxdy,

\alpha\beta 是方向余弦。

在前面提到的博客文章中,我们介绍了如何模拟大型光学器件;即可以使用波动光学、频域 接口来计算光学元件周围的小域,然后使用弗劳恩霍夫或菲涅耳衍射公式,或者使用光束包络 接口来模拟整个域。然而,这两种方法仅适用于慢速(低 NA)镜头,因为快速(高 NA)镜头需要大量网格单元,而夫琅禾费和菲涅耳公式无法给出很好的近似值。波矢太陡,光束包络 接口无法对计算域进行网格划分。

模拟大型高数值孔径透镜的唯一方法是使用角谱法(ASM)。这是一种传播方法,与夫琅禾费和菲涅耳衍射公式属于同一类型的数值方法。一旦知道一个平面中的场,就可以计算另一个平面中的场。角谱法是严格的,因为它满足亥姆霍兹方程。该方法可以结合波动光学模块计算某个域中的场,然后使用角谱法将场传播到更远的平面。

下面是一个高 NA 镜头 (NA = 0.66) 的示例,它比 DVD 光读取快得多。镜头半径为 16μm,后焦距(镜头第二面与焦平面之间的距离)为 10μm。使用 几何光学 接口与优化模块相结合,对该透镜头进行了优化,使其在 0.66 μm 的波长下具有衍射极限。透镜被特意设计得很小,以便波动光学,频域 接口可以计算出严格的解进行比较。我们将演示如何使用角谱法将场从该透镜的出射面传播到焦平面。

左侧是使用射线光学模块和优化模块设计的高数值孔径镜头模型。 右边是镜头的全波模拟。
使用射线光学模块和优化模块(左)设计的 NA = 0.66 镜头。使用 波动光学,频域接口(右)模拟的透镜的全波模拟。注意代表透镜出射平面的线,场从该平面传播到最右边的边缘,即焦平面。

比较菲涅耳衍射公式(蓝线)、亥姆霍兹溶液(粉线)和 ASM(绿线)的高 NA 透镜光斑轮廓线图。
NA = 0.66 透镜模型的光斑轮廓比较菲涅耳衍射公式;波动光学,频域接口的严格解;和角谱法。请注意,对于这个透镜,菲涅耳衍射公式不再准确。(为了更好的比较,斑点轮廓显示为 11 μm 而不是 10 μm。)

为了执行两次 FT,我们需要将第一次 FT 存储在数据集中。这是因为 fft() 算子只是一个后处理算子,不是通用算子,比如可以在物理设置中使用的算子 integrate。目前,在当前版本的 COMSOL Multiphysics 中(在未来版本中,fft() 算子将被提升为通用算子),我们仍然需要在第一次 FT 的物理场设置中通过数值积分来使用FT,然后将 fft() 算子用于后处理中的第二次 FT。边界常微分和微分代数方程分布式常微分方程 节点的接口被定义在透镜出射平面上,通过数值积分 FT 执行第一次 FT,并将结果存储为函数 u,如下图所示:

模型开发器的屏幕截图,其中选择了分布式 ODE 节点,并展开了相应的设置窗口,其中边界选择、源项、阻尼或质量系数和质量系数部分展开。
ASM 中第一个 FT 的 边界常微分和微分代数方程设置的屏幕截图。请注意,我们通过透镜半径 D/2 对傅里叶空间进行了归一化,进行适当的缩放。

COMSOL Multiphysics UI 的屏幕截图,显示了带有 Line Graph 1 标签的 Line Graph 设置窗口。数据、y轴数据 和 x轴数据 部分均已展开。
后处理中 ASM 中第二次 FT 设置的屏幕截图。对于第二次 FT,注意方向余弦 ay 轴数据中由归一化的 y 坐标表示,\alpha = y/(D/2),x 轴数据中的归一化因子 1/wl 来自变量的微分 d(\alpha/\lambda)。另请注意,空间频率名称 y2 是在 空间 FFT 数据集中任意选择的。

值得一提的是,第二次 FT 其实就是逆 FT,但是 FT 的绝对值和逆 FT 到常数之间没有区别。看到 ASM 作为亥姆霍兹解给出了准确的结果,我们可以将这种技术用于其他高 NA 透镜系统,例如大型高 NA 菲涅尔透镜。

结束语

在这篇博客中,我们了解了如何设置新的空间 FFT 特征的基础知识以及如何在光学的一些重要应用中使用它的示例。我们将在本系列的下一篇博客中讨论第四个应用,即部分相干光束计算(公式与此处用于第三个应用的公式相同)。

参考文献

  1. J. W. Goodman, Introduction to Fourier Optics, 3rd ed., Roberts and Company Publishers, 2005.
  2. M. Born and E. Wolf, Principles of Optics, 7th ed., McGraw-Hill, 1968.
  3. A. C. Schell, “A technique for the determination of the radiation pattern of a partially coherent aperture,” IEEE Trans. Antennas Propag., vol. 15, no. 1, pp. 187–188, 1967.

评论 (0)

正在加载...
浏览 COMSOL 博客