菲涅尔转换的数值计算

菲涅尔转换的数值计算

在这篇论文中，我们通过均匀空间采样所获得的有限数值来处理计算菲涅尔转换积分这一问题。通常简单的数值采样规则能够让使用者对于任何散射距离都能计算出他的分布。这些规则是如何被扩展到利用基于快速傅里叶变换这一算法来提高计算效率的, 在这篇论文中也有所解释。这篇论文也把本方法与其他理论方法做了比较。

1、 引言

光的传播和衍射是个经典问题，在当下仍十分重要。麦斯威尔等式为检测光的行为提供了一个严格的框架，他把光描述成一种电磁波。电场和磁场组成部分的行为都应该特别的考虑，例如，当分析数值孔径高的镜头，就如在显微镜中或者是为平面等离子体匹配模式时发现的。然而在许多场合，对于一个给定的物理问题进行模式匹配时，电磁场的全矢量描述是不需要的，一个简单地标量模型可以用在这里。用一个标量模型去表示一个光学场是一个有用的近似表示。在光学制度中，它为光通过孔径（在这里孔径的大小远大于光波长）的衍射提供了一个高度准确的预测，可参见【1】的实验结果。其他相关理论和数值计算在【2,3】中被提出来了。在这篇论文中我们假设一个标量模型对于我们的目的是足够精确的。

在光传输的标量描述中，区分以下两个定义是非常有用的，即傍轴和非傍轴衍射积分。两类衍射积分来自于麦斯威尔等式，参见【4,5】。然而在傍轴处理过程中，在推倒时应做更多的近似处理。因此，傍轴模式更加脱离于真实的衍射基础物理过程，伴随着这个结论还有非傍轴衍射积分被认为是更加的准确。我们参见在【4】的3、4章中写的衍射处理，并且据此把基尔霍夫定律或瑞利-索末菲第一和第二定律归类为非傍轴衍射积分。菲涅尔转换是一个傍轴衍射积分。

在这篇论文中我们用有限的采样数值来研究菲涅尔转换的数值计算问题。这是因为数值模拟技术对模拟光传输和为没有解析解的普通衍射问题计算解析解是相当重要的。菲涅尔转换在数字全息术【6-12】、迭代相位检索【13-15】和估计相位分布的强度方式的传输中扮演着不可缺少的角色，估计相位分布的强度方式的传输在生物图像应用【16,17】中起着重要的作用。菲涅尔转换也广泛应用于针对波束成形【18,19】的衍射光学元件的设计中，在计量应用中衍射光斑区域的去相关也被应用于判定变形，例如，在测试时光学粗糙表面的压力、应变以及裂缝，【20,22】，或者是位移和倾斜的检测【23,24】。在这所有的应用当中，菲涅尔转换的数值计算是必要的，因此完全明白如何执行运算具有显著的理论和实践重要性。

或许在理论计算过程中遇到的第一个问题会是如何选择一个合适的均匀分布的抽样值来代替输入平面的复杂标量波场。在文献中有许多通过检测需要多少抽样值来处理这个问题的手稿，因此在检测内核或者是菲涅尔啁啾函数时尼奎斯特采样率（或者是一个相关的标准）被采用了，参见例子【25-31】。然而这个处理可以等效于准确的理论计算，它往往可以是不必要的数值密集型【28】。用于计算输入输出信号的范围大小被发现也影响光学采样率，然而在文献中也出现了一些关于理想采样率的意见分歧，例如【29-31】。通常也不是很清楚当违反采样规则时数值计算是为何变得不准确的。因为这些因素，很有必要来重新审视这个问题，即为衍射积分运算提供一个合适的采样率。尽管我们的分析是针对于菲涅

尔转换，但是这一结果可以概括到在【29-31】中讨论的傍轴积分中。

在此，我们以一种完全不同的方式来处理采样问题，即依靠一种对数值衍射过程的直观的物理描述。这一描述对分析普遍的衍射问题和准确的确定通过采样所引进的理想误差都是有帮助的。和文献中的其他手稿不同的是我们在决定合适的采样率时不考虑菲涅尔内核和啁啾函数在计算中的作用。相反地，我们依靠输入域的空间和空间频率范围并依靠衍射平面中重现波的位置来定义采样率。这些重现波的出现是由于采样操作（平面1中的复杂信号的离散化，看图1和第3部分）。这些重现波和他们特性之间的离散化取决于菲涅尔转换是如何运算的。鉴于我们正在计算的信号的输出范围小于相邻重现波之间的离散这一事实，我们认为我们的信号是被合适采样的。

这篇论文的论据基于以下考虑：

（1） 在空间和空间频率域有一个有效的有限的支持的信号被定义为使用来自在【27,28,32-35】中描述的空间-带宽产品的概念。这种信号包含定义的（或使用者选择的）最大的空间频率，f x ，和一个有限的输入范围。衍射的真正的物理过程使得这个信号在衍射平面中有一个无限的输出空间范围，我们后边会定义为SE 。

（2） 从一个有限的采样值而不是从1中考虑的连续情形来计算衍射分布从根本上改变了这一问题的性质。这一理论计算的信号将重新产生一个来自于能够产生范围为SE 的输出平面的点1的分析结果。然而除此之外，重现波的一个无限集产生于衍射平面。这些重现波的特性取决于这一运算是如何运行的。

（3） 对于数值计算衍射积分考虑了两种方法。直接计算衍射积分会产生彼λ此间隔为z 的复制平面，在此λ是波长，z 是传输距离，δX 是在输入空间域

X

的采样间隔。通过运算在输入域中信号的空间频率没有发生相位移动，并且不是通过尼奎斯特标准采样获得的空间频率也将会传播并会加强在衍射域中的场。基于光谱的运算在衍射平面中产生了一个彼此间隔的重现波的无限集，在此v

δv 是空间频率域的采样间隔。在光谱运算过程中衍射信号的空间频率成分经常被限制在空间频率域。

图1，衍射过程的图示。一个无线平面波入射进平面1中的孔径中，该平面以空间变量X 表示。考虑到如何计算P(x,z)中的复振幅，仅仅是开放孔径的光对P(x,z)中的复振幅有贡献，我们把开放孔径内的光设想

成由无穷多个点光源发光组成的，每一个点光源都带有强度值和相位值（图中只画了一部分点光源）。因为这个孔径是被理想平面波照明的，每一个第二级点光源都带有相同的强度值和相位值，尽管当考虑到距离

l1和l4时，由于光的路径不同每一点光源将积累不同的相位值。

（4）我们建议的抽样规则遵从以上1-3点：SE 比v 和λz X 都小。当这一

点确定后，这个衍射信号被认为是被很好抽样的。我们再一次标注，以上所采取的方法不同于文献中所提出来的方法，这是因为我们在思考时没有用到菲涅尔啁啾函数的特性。

在2-5节中，我们得出了以上的点1-4。在第六节中，我们把注意力转移到计算菲涅尔转换的高效计算算法，对于这一点我们用到了快速傅里叶变换。两种算法被明确的命名为：直接算法和光谱算法，其中光谱算法是我们在第4和5节中讨论过的直接和光谱计算的快速实现（基于快速傅里叶变换）。我们的理论结果与文献中的其他分析方法所得的结论做了对比。最后，在第七节中，我们考虑一个重要的案例在这个案例中我们在第4和5节中提出的抽样规则是不合适的：也就是说，在衍射区被抽样时我们希望对这个结果做菲涅尔逆变换，这后一种情况对数字全息照相术有重要的寓意，然后我们总结这篇论文的主要发现。

2、 一个分析解决方案

在随后的介绍中为了表述简单，我们执行一维分析。在图1中，我们简绘了问题的特性。平面1中的场我们以大写字母U （x ）表示，并假设它是已知的，我们的任务是推导出衍射场的表达式，即平面2的表示式u z (x ) 。通过【4】中的分析，我们现在写出菲涅尔转换的定义式：

u z (x ) =FST z {U (X )}(x )

=exp(j 2π) j z ⎰∞

-∞U (X ) j π(x -X ) 2]d X (1) λz

在这里λ是单色相干光（在时间和空间上的）的波长，z 是平面1和平面2的轴距离（参见图1），j=-1，变量X 和x 分别代表在平面1和平面2上的空间坐标。为了余下的分析我们将降低开始的相位项，即exp(j 2πz ) 。 在有些衍射问题中有可能能解出等式1的解析解。由于一些原因这些解决办法非常重要：他们洞察了这个衍射过程，他们也可以作为一种评估数值方法准确度的方法。这个可以通过将数值计算结果和通过对的分析方法而获得的结果作比较而获得。一些与菲涅尔转换有关的特性在【27,36-41】中讨论过了。我们用下面这个等式第一次检验等式1的分析方法：

X 2

U(x)=exp(-2) cos(2πf x X), (2) a i

它也可以表述为下式：

1X 2

U (X ) =exp(-2)[exp(j 2πf x X ) +exp(-j 2πf x X )] (3) 2a i

在上式中，a i 和U(x)的空间扩展有关。然而严格的说，U （x ）扩展到一个无穷的区域，信号的超过99.9%的能量存在于以下的范围中：-2a i ≤X ≤2a i 【33】。U(x)也有空间频率成分f x 。将等式3带入等式1就得到了如下结果：

u z (x ) =

并且 1+[u z (x ) +u z -(x )] (4) 2

(x -λzf x ) 2

l q 2 u (x ) =K z exp[-]exp[j (φx +φx )] （5） z z 2a z +z

a z =(πa i

2a 4i +z 2λ2) -1 （6）

jz λ-1-j πa 4i π2f 2x z λK z =(+2) 24) （7） πa i πa i +z 2λ2

2π3f x a i 4 φ=24 （8） πa i +z 2λ2l z

φ=q

z πz 2λ2πλza +z λ24

i 33 （9）

我们注意输入的调制高斯函数已经被转换成带有4a z 的空间扩展的高斯函数。高斯函数U z (x ) +的中心位置在x=λzf x [U Z (x ) _存在于x=-λzf x ]。因此在菲涅尔转换之后我们会看到U(x)在平面2将有一个大小为SE 的空间扩展。

SE=2(λzf x +2a z ). (10) 这代表了我们的第一个重要结论。它表示一个带有有效有限支持并知道它的空间频率成分f x 的信号由于他的两个独立分量：f x 和a z ，所以它将随着它的传播程度而增强。由于结果相对简单我们选择了这个分析形势，它使得人可以看到和特定空间频率有关的能量是如何随着传播而远离光学轴的，然而仍然有部分处于高斯外壳以内，参见【42】和【43】的18.6节。空间频率越高它随着传播偏离空间轴越远。尽管如此，对于相对简单的衍射问题这也是一个解决办法，它的结论可以被扩展到其他更多的复杂案例中。

几乎所有的我们希望去数值估计的衍射问题都在输入平面内有一个有限的扩展，和等式2很相似。总的来说这些信号在输入平面内都有一些和他们相似的空间结构。当我们试图去推广等式10的结果到更多的普通信号时出现了一个问题。这是因为等式10的结果来自于等式2，等式2有一个明确的空间频率成分。对于一般信号，我们如何去选择一个合适的空间频率值f x 呢？这个问题更详细的

在第五节被检验了，在那儿我们展示了通过研究信号在它的傅里叶区的能量是如何分布的去为f x 选择一个合适的值，也可以参看【27,34,42,44】。我们假设如果信号的傅里叶转换能量具有多余一个具体数值存在于某个空间频率下，这上一个空间频率值就变为f x 。通过计算输入波在整个输入面的卷积可以计算出这个信号在空间区的总能量。

3、 数值计算的输入

在接下来的一节中，我们介绍可以用来计算菲涅尔转换的两种不同的方法。我们以直接计算和光谱计算来提出它们。这些方法不同于接下来将要在第6节中介绍的快速数值计算，因为它们没有利用快速傅里叶转换算法并且输出空间的坐标轴x 可以被任意选取。这一节的目的是清楚的描述我们数值计算方法的输入，这个输入将采用两个坐标矢量的形式：一个空间矢量和一个代表我们输入信号的复杂值。这些值一旦被选取，我们就可以比较每个计算方法的相关表现。

为了开始我们的数值分析，我们假设我们有一个复杂值矢量U ，它对应于函数U(x)在均匀空间采样点上的值（采样点以δX 间隔分开）。并用空间矢量定义。

U =[U (X 1), U (X 2),... U (X N ), ],

=[U 1, U 2,... U N , ] （11）

在此

X =[X 1, X 2,... X N , ] （12）

并且在此N 是采样总数。在这篇论文的其他地方，我们用黑体字符来表示一个矢量。这两个矢量还有已经给出的λ和z 我们定义为是数值计算的输入参数。

4、 直接计算

因为我们现在仅使用抽样中的一个有限值来表示平面1中的输入场，因此等式1的衍射积分推导得出一个有限和式

δX u (x ) =j z S

z ∑U

n =1N n j π(x -X n ) 2] λz

2N j πX n -j 2πxX n δX j πx 2

=) ⨯{∑U n ) )} （13） λz λz λz j z n =1

在这里n 是一个整数。注意到u S z (x ) 是一个关于x 的连续函数，并且它是通过一个有限的数据点集计算出来的，我们用上标S 来表示它。因此等式13可以等同于以下描述：乘以一个啁啾函数，经过一个傅里叶的尺度变换，再乘以另一个啁啾函数。

我们现在用一个具体例子来检查一下u S z (x )

U (X ) =p L (X ) cos(2πΓX ) （14）

在这里

, p L (x ) ={1

0, when x

otherwise （15）

我们现在来抽样U(x)，取N=100，变化区间是-L ≤X ≤L 在此L=0.1mm，已知δX 为0.2μm 。其它数据取做z=1cm,λ=505.7nm,Γ=30条。在图2和3中，我们分别表示出了在区间-≤X ≤λX +0. 5 mm 和区间0≤X ≤δX mm 上的计算结果。在图2中我们表示震级分布，在图3中表示相位分布。

图2. 衍射计算的结果。重现波的出现在x=2.5mm处是清晰可见的。红色虚线表示衍射区域的范围。以黑线画出的分布已被空间滤波，不含有高于35条的位置，解释的空间频率。注意信号的有限范围和蓝线作比较。（报告，以蓝色绘制的中央级重现波所画的范围不会超过区间- ≤x ≤）

图3. 衍射运算的相位分布结果。两条线都超出了范围δv 。左图是零级重现波的相位。右图是一级重现波的，在右图中有一个大小为π的恒定偏移量和一个斜率为X 的线性偏移量。

对于这一数值例子计算SE 是有指导意义的。对于这个计算，我们最初选择f x =Γ条, 随后我们将会看到，等式14在空间频率域包含的能量高于Γ，这是由于矩形函数它限制了信号有限的空间扩展。为了计算出a z 的值，我们需要将a i 和等式15确定的矩形函数的有限扩展联系起来，例如我们设定4a i =2L。结

果画在了图2中。在图2中我们看到在x =λX =2. 5 mm 处有一个波形的再现。比较一阶重现和中心处的相位，我们可以发现一阶相位分布有一个额外的线性相位和一个恒定相位。我们注意到由于运算在平面2中的波形重现中有一个无穷量，它可以通过在一个很广的区间x 上绘制Uz (x ) S 的图形而验证。接下来的关系能够保持并能通过附录A 中的方法概论推导出来：

m λz j πλzm 2j 2πmx u (x ) =∑u z (x -) ⨯) ) （16） 2δX δX δX m =-∞S z ∞

从这些分析中我们可以发现下面的采样规则证明了它本身，例如：

X （17）

为了保证相邻的重现波不发生相互重叠，我们在数值计算衍射分布时引入了误差。

当DC 被用作计算菲涅尔积分时，所有的空间频率都通过媒介传播并且对平面2中的复杂分布有贡献，甚至是那些通过尼奎斯特采样率获得的空间频率也这样。从图2中发现空间频率成分在空间坐标大于的地方出现是很明显的，在这里一个振荡信号在x >时向外扩展。这个震荡是由于输入信号的空间频率高于Γ。对于f x =Γ这一不理想的取值，导致了在处没能反映出衍射信号准确的输出扩展。选择一个合适的f x 值应该基于输入信号的能量在傅里叶区的分布。在下一节中，我们将验证更多的基于信号产生的空间-带宽概念来为f x 选择一个合适的值。

然而，让信号在傅里叶区通过一个空间滤波系统来限制u S z (x ) 的空间频率扩展并因此控制SE 是很有可能的。在这里我们移动所有的频率到高于35条 SE

【连续的Uz (x ) S 的采样率非常高，这些样本是使用了FFT 操作的傅里叶转换，并且这个结果的空间频率分布被乘了一个矩形孔径；这一滤波空间频率分布后又经过了快速傅里叶逆变换。】。在图2中，经过这样一个滤波系统的结果可以通过检测黑点来获得。这一分布有一个非常清晰有限的空间扩展并且它与等式10以

计算出的SE 非常一致。这一分布（以黑线画出）有一个遵从于第2节中说过的结论的空间扩展，然而这一分布当与正确的分析方法做比较时并不是那么准确。滤波系统非常明显的移动了位于x =2. 5mm 处的重现波的结构细节的一部分。

尽管如此，如果图2所表示出的空间滤波结果（这一分布以黑线表示）被认为是满足于我们的要求的，那么我们要注意到可以通过减少抽样数量来得出一个f x =Γ=30 条

几乎相同的结论。对于一个给定的输入，减少抽样数量就降低了重现波之间的间隔；空间频率滤波系统控制SE 。在第5节中，我们发现对于f x 的一个合理的取值是60 条，它在确保重现波之间有少许交叉时来保持高度准确度提供了一个合理的平衡点。

我们建议当输入信号的空间和空间频率扩展已经知道时，等式17的关系式可以作为一个抽样规则来使用。对于等式14所讨论的具体案例，我们发现对于z 和Γ的一系列不同取值，信号至少有98%的能量存在于-≤x ≤上，当f x =60 条时，参见等式18。尽管SE 的标准在等式14中已经被具体验证了，我们发现对于一系列不同的信号类型它仍然满足，包括有限空间频率扩展的随机信号。

与第一抽样规则有关的重要一点是当z 变小时，重现波彼此之间距离减小。这一效应可以通过增加N （仍保持2L 的输入扩展）来补偿，因此降低了δX 。然而当过了某一确定点后，这一计算变得极其的数值密集，因此对于实际计算很难。

在这一节的最后我们用下式来计算输入信号的总能量，P i :

P i =⎰(X ) d X (18) -∞∞2

对于在这里出现的数值案例平面1的任意单元都达到了P i =0.0001。

5、 光谱运算

为了执行菲涅尔转换这个方法使用了傅里叶转换。对于一些复杂信号f (x ) 在这里定义正向和逆向傅里叶变换是很合适的：

~f (v ) =FT {f (x )}(v ) =⎰f (x ) exp(-j 2πxv ) d x （19） -∞∞

f (x ) =IFT (f (v ) )(v ) =⎰f (x ) exp(j 2πxv ) d v （20） -∞~∞

在这里v 空间频率坐标，函数f (x ) 和f (v ) 是傅里叶转换对儿。

从【4】中等式1也可以写成以下形式： ~

u z (x ) =IFT {exp(-j πλzv 2) FT {U (X )}(v )}(x ) （21）

等式21和1是数学上的等价命题。在下面的小节中我们研究当使用有限的抽样数值来对等式21进行数值计算时出现的不同的步骤。为此介绍离散形式的傅里叶变换是有帮助的，在这里再一次声明大写的S 表示这一连续函数来自于由复杂数据组成的有限集，在我们的案例中是来自于矢量U ，在第三节中明确定义：

U (v ) =δX ∑U n exp(-j 2πvX n ) （22）

n =1~S N

在这里n 是一个整数。

我们注意到等式22描述的分布也包含了一个无穷数量的重现波，他们在傅里叶区是彼此分离的，彼此间隔为∆B =X 。为了把它们和在第4节讨论的重现波区别开来我们这些重现波为傅里叶重现波。

A. 傅里叶逆变换操作

用等式22来替换等式21中的函数FT {U (X )}(v ) ，我们得到下式：

u (x ) =⎰exp(-j πλzv ) U (v ) exp(j 2πxv ) d v （23） -∞S z ∞2~S

在等式23中我们已经明确包含了傅里叶逆变换积分来强调积分是发生在整个傅里叶平面，因此这个积分就包含了来自这个平面的所有傅里叶重现波的贡献。这一沿着整个傅里叶平面的分布乘以一个有无限扩展和有效焦距为-(λ2z ) 的透镜或啁啾函数，更多细节参见【40】。

并且我们用这个矢量来抽样连续场

~

U =[U (v 1) , U (v 2) ,... U (v N ) , ]=[U 1, U 2,... U N , ]（26）

将这些矢量带入等式23中，并用有限求和来替换这些积分我们可以得到以下结果：

2 u z SSC (x ) =δv {∑U exp(-j πλzv n ) exp(+j 2πxv n )} （27）

n =1n N ~S S S ~S ~S ~~S ~S ~S