角光谱方法在菲涅尔衍射中的高效实现与验证

1. 从“硬算”到“巧算”：为什么我们需要角光谱方法？

如果你之前尝试过直接计算菲涅尔衍射积分，那你一定体会过什么叫“算力焦虑”。那个双重积分公式，看起来就让人头疼，更别说用程序去实现了。它需要对整个源平面上的每一个点进行遍历计算，计算量随着采样点数的增加呈爆炸式增长。我刚开始做光学仿真的时候，就曾试图用最笨的循环嵌套去实现它，结果一个512x512的仿真算了快十分钟，这还怎么玩？

所以，光学模拟领域的前辈们早就开始寻找更聪明的办法。角光谱方法（Angular Spectrum Method, ASM）就是这样一个“巧算”的典范。它的核心思想非常巧妙：把复杂的空域卷积，转换到频域里去做简单的乘法。这就像你要计算两个超大数字的卷积，直接算很麻烦，但如果你先把它们都转换成对数，相加之后再转换回来，就简单多了。角光谱方法利用的正是傅里叶变换这个“数学显微镜”，让我们能透过现象看本质。

简单来说，光波的传播过程，在频域视角下，可以看作是不同的空间频率分量（也就是不同方向的平面波，即“角光谱”）各自经历了一个相位延迟，然后再叠加起来。这个相位延迟的规律，就由一个叫做“传递函数”的东西来描述。这样一来，整个传播过程就变成了三步：1）把源平面的光场变到频域；2）在频域乘上传递函数；3）再变回空域。计算复杂度从O(N⁴)直接降到了O(N² log N)，这效率提升可不是一点半点。

这个方法特别适合我们这些需要快速验证光学设计、分析衍射图样的工程师和研究者。无论是设计一个新型衍射光学元件，还是模拟激光光束经过光学系统后的变化，角光谱方法都能提供既快速又足够精确的结果。接下来，我就带你亲手实现它，并验证它到底有多靠谱。

2. 角光谱方法的两种“武功秘籍”：脉冲响应与传递函数

角光谱方法在具体实现时，有两种在数学上等效但代码实现上略有差别的路径，我习惯把它们比作“武功秘籍”的两种不同练法。理解这两种方式，能让你更深刻地把握方法的本质，并在不同场景下做出最佳选择。

2.1 秘籍一：基于脉冲响应函数（Impulse Response）

脉冲响应函数 h(x, y)，物理意义非常直观：它描述了一个位于源平面原点的点光源，传播距离 z 后，在观察平面上形成的复振幅分布。你可以把它想象成光学系统的“单位脉冲响应”。整个源平面上的复杂光场，可以看作是无数个点光源的叠加。根据线性系统的叠加原理，观察平面的光场就是所有这些点光源产生的脉冲响应的加权（权重就是该点的复振幅）叠加之和。数学上，这就是一个卷积运算。

公式就是咱们熟悉的样子： U2(x, y) = U1(x, y) * h(x, y) 这里的 * 代表二维卷积。而卷积运算，在计算机里如果直接按定义去算，又回到了老路，效率低下。幸运的是，我们有卷积定理这个“法宝”：时域（空域）卷积等于频域乘积。所以，我们的计算流程就变成了：

1. 分别计算 U1 和 h 的傅里叶变换，得到 F_U1 和 H。
2. 在频域将两者相乘：F_U2 = F_U1 .* H。
3. 对 F_U2 做逆傅里叶变换，就得到了观察平面的光场 U2。

这种方法的思路很直接，但有一个关键的实现细节需要注意：在计算脉冲响应函数 h 的傅里叶变换 H 时，我们理论上得到的就是传递函数。但在数值计算中，我们需要对 h 进行离散采样和FFT。h 函数本身是一个在空域振荡的复指数函数，其“频率”与坐标的平方成正比。因此，我们必须确保我们的采样间隔 delta1 和采样范围足够大，以捕捉到 h 函数的所有重要振荡细节，否则就会产生混叠误差，导致仿真结果失真。这是使用脉冲响应法时需要特别留意的“坑”。

2.2 秘籍二：基于传递函数（Transfer Function）

传递函数 H(fx, fy) 是脉冲响应函数的傅里叶变换。它直接定义在频域，其物理意义更加清晰：它告诉我们在空间频率 (fx, fy) 上的平面波分量，在传播了距离 z 后，振幅不变，但相位会增加一个与频率平方相关的项 exp(-jπλz(fx²+fy²))。

基于传递函数的实现流程更为简洁：

1. 计算源平面光场 U1 的傅里叶变换 F_U1。
2. 直接在频域生成传递函数 H，并与 F_U1 相乘：F_U2 = F_U1 .* H。
3. 对 F_U2 做逆傅里叶变换，得到 U2。

你会发现，这种方法跳过了在空域生成 h 并对其做FFT的步骤。它直接在频域操作，概念上更贴近角光谱方法的物理本质——即对各个角谱分量进行相位调制。在代码实现上，它通常也更简洁，且避免了因对 h 采样不足而引入的混叠误差。因为 H 是一个光滑的复指数函数，我们只要在频域坐标上正确计算它即可，无需担心空域采样的限制。

那么问题来了：既然两种方法数学上等价，为什么结果会有细微差别？哪种更好？这正是我们接下来要通过代码实战来验证的核心问题。差别就源于我们刚才提到的离散化与数值采样。脉冲响应法受限于空域采样，传递函数法则受限于频域采样（即计算窗口的大小）。在不同的仿真参数下，这两种限制对结果的影响程度不同。

3. 手把手代码实战：从方形孔衍射看两种实现

光说不练假把式，咱们直接上MATLAB代码（Python的实现逻辑几乎完全一致）。我会用一个最经典的例子——方形孔的菲涅尔衍射，来演示这两种方法，并对比它们的结果。我们假设一束波长为532nm的绿光，照射在一个边长为1mm的方形孔上，观察距离为0.1米处的衍射图样。

3.1 准备工作：定义方形孔和辅助函数

首先，我们需要一个生成方形孔光场的函数，这里用到了一个自定义的 rect 函数，它生成一个在给定范围内值为1，其他为0的矩阵，用来模拟方形孔。

function y = rect(x, D)
    % 矩形函数，当 |x| <= D/2 时输出1，否则为0
    y = abs(x) <= D/2;
end

然后，是仿真参数的设置。这部分参数直接影响仿真精度和计算量，需要根据实际情况仔细选择。

clear; close all; clc

%% 仿真参数设置
a = 1e-3;          % 方形孔边长，单位：米
N = 1024;           % 采样点数（每边），必须是2的幂次以便FFT高效计算
L = 1e-2;           % 源平面的物理尺寸，单位：米。要大于孔径，并留出足够缓冲区域。
delta1 = L / N;     % 源平面的采样间隔（假设x和y方向相同）
lambda = 532e-9;    % 光波长，单位：米
z = 0.1;            % 传播距离，单位：米

%% 生成方形孔光场
[x1, y1] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-N/2 : 1: N/2 -1) * delta1);
U1 = rect(x1, a) .* rect(y1, a); % 二维方形孔

% 显示源平面光场（即方形孔）
figure(1);
imagesc(x1(1,:), y1(:,1), U1);
colormap('gray'); axis square; axis xy;
xlabel('x [m]'); ylabel('y [m]');
title('源平面：方形孔');

运行这段代码，你会看到一个位于中央的明亮方形，这就是我们的光源。接下来，好戏开场。

3.2 实现秘籍一：脉冲响应函数法

根据2.1节的原理，我们首先实现脉冲响应函数 h，然后利用卷积定理进行计算。

function U2 = angular_spectrum_impulse(U1, lambda, delta1, z)
    % 基于脉冲响应函数的角光谱方法
    % 输入：
    %   U1: 源平面复振幅分布
    %   lambda: 波长
    %   delta1: 源平面采样间隔
    %   z: 传播距离
    % 输出：
    %   U2: 观察平面复振幅分布

    [M, N] = size(U1);
    k = 2 * pi / lambda; % 波数

    % 生成源平面坐标网格
    [x, y] = meshgrid((-N/2 : 1: N/2 -1) * delta1, (-M/2 : 1: M/2 -1) * delta1);

    % 计算脉冲响应函数 h(x,y)
    h = exp(1j * k * z) / (1j * lambda * z) * exp((1j * k) / (2 * z) * (x.^2 + y.^2));

    % 利用卷积定理计算衍射
    F_U1 = fftshift(fft2(fftshift(U1))); % 对U1做FFT
    H = fftshift(fft2(fftshift(h))) * delta1^2; % 对h做FFT，并乘以采样间隔的平方以近似积分
    U2 = fftshift(ifft2(fftshift(F_U1 .* H))); % 频域相乘，再逆变换回空域
end

调用这个函数并显示结果：

U2_impulse = angular_spectrum_impulse(U1, lambda, delta1, z);

figure(2);
imagesc(x1(1,:), y1(:,1), abs(U2_impulse));
colormap('gray'); axis square; axis xy;
xlabel('x [m]'); ylabel('y [m]');
title('脉冲响应法：衍射光斑（强度）');

你会看到，原来清晰的方形光斑，经过一段距离传播后，变成了一个带有复杂条纹的衍射图样，中心最亮，边缘有明暗相间的衍射条纹。这就是菲涅尔衍射的典型特征。

3.3 实现秘籍二：传递函数法

现在，我们用更直接的传递函数法来实现。

function U2 = angular_spectrum_transfer(U1, lambda, delta1, z)
    % 基于传递函数的角光谱方法
    % 输入输出参数同上

    [M, N] = size(U1);
    k = 2 * pi / lambda;

    % 生成频域坐标网格 (空间频率)
    % 注意：频域坐标间隔 = 1 / (物理尺寸) = 1 / (N * delta1)
    [fx, fy] = meshgrid((-N/2 : 1: N/2 -1) * (1 / (N * delta1)), ...
                        (-M/2 : 1: M/2 -1) * (1 / (M * delta1)));

    % 计算传递函数 H(fx, fy)
    H = exp(1j * k * z) * exp(-1j * pi * lambda * z * (fx.^2 + fy.^2));

    % 计算衍射
    F_U1 = fftshift(fft2(fftshift(U1)));
    U2 = fftshift(ifft2(fftshift(F_U1 .* H)));
end

同样地，调用并显示结果：

U2_transfer = angular_spectrum_transfer(U1, lambda, delta1, z);

figure(3);
imagesc(x1(1,:), y1(:,1), abs(U2_transfer));
colormap('gray'); axis square; axis xy;
xlabel('x [m]'); ylabel('y [m]');
title('传递函数法：衍射光斑（强度）');

乍一看，图3的结果和图2几乎一模一样。这说明两种方法在宏观上都正确模拟了衍射过程。但“几乎一样”就意味着有细微差别，真相需要更严格的检验。

4. 终极验证：与“标准答案”对比，谁更精确？

仿真的结果再漂亮，我们也得问问：这跟物理现实符合吗？对于方形孔菲涅尔衍射，我们是幸运的，因为它存在解析解（Analytical Solution）。这个解析解可以用菲涅尔正弦和余弦积分来表示，是理论上精确的结果，可以作为我们数值模拟的“金标准”。

这里我们不深入解析解的复杂公式，直接使用一个实现好的函数 analytical_square_fresnel_propagation（其内部基于菲涅尔积分函数 fresnelc 和 fresnels 计算）。我们将数值模拟的结果与这个解析解进行对比。

% 计算方形孔菲涅尔衍射的解析解
U2_analytic = analytical_square_fresnel_propagation(x1, y1, a, lambda, z);

% 显示解析解的结果
figure(4);
imagesc(x1(1,:), y1(:,1), abs(U2_analytic));
colormap('gray'); axis square; axis xy;
xlabel('x [m]'); ylabel('y [m]');
title('解析解：衍射光斑（强度）');

% 一维切片对比：取穿过图像中心的一行数据（第513行，因为N=1024）
row_index = N/2 + 1; % 即513行
sim_impulse = abs(U2_impulse(row_index, :));
sim_transfer = abs(U2_transfer(row_index, :));
ana = abs(U2_analytic(row_index, :));
x_line = x1(1, :);

figure(5);
plot(x_line, sim_impulse, '--r', 'LineWidth', 1.5); hold on;
plot(x_line, sim_transfer, '--g', 'LineWidth', 1.5);
plot(x_line, ana, 'b', 'LineWidth', 1);
xlim([-0.005, 0.005]); % 聚焦中心区域
xlabel('x [m]'); ylabel('光场强度');
legend('脉冲响应法模拟', '传递函数法模拟', '解析解');
title('一维切片对比（中心行）');
grid on;

这张对比图会告诉你一切。你会看到三条曲线：蓝色的解析解曲线、红色的脉冲响应法曲线和绿色的传递函数法曲线。

在我的多次测试中，通常会发现：绿色的传递函数法曲线与蓝色的解析解曲线贴合得几乎严丝合缝，尤其是在主要的衍射瓣区域。而红色的脉冲响应法曲线，在整体趋势一致的前提下，在高频振荡的细节处（比如衍射条纹的峰谷值）有时会出现微小的偏差或衰减。这个偏差的根源，就是我们之前提到的混叠。

因为脉冲响应 h 是一个空域振荡函数，其“频率”随着离中心距离的平方增加而增加。在我们有限的采样窗口 L 和采样间隔 delta1 下，对于 h 函数边缘的高频振荡成分，可能无法满足奈奎斯特采样定理，导致这部分信息在离散化时失真。这个失真被带入了后续的FFT计算中，最终影响了仿真精度。

而传递函数 H 是频域的一个平滑二次相位因子，只要我们的频域坐标 fx, fy 计算正确，就不会有采样不足的问题。因此，在大多数参数设置下，基于传递函数的实现方式数值稳定性更好，精度更高。

但这并不意味着脉冲响应法一无是处。在某些特定场景，比如需要考虑像差或非均匀介质时，脉冲响应的概念可能更容易与物理模型结合。而且，当你需要非常关注计算窗口边缘的效应时，两种方法可能呈现出不同的误差特性，需要具体分析。

5. 效率与精度权衡：让你的仿真又快又准

经过上面的对比，你可能已经倾向于使用传递函数法了。但在实际项目中，我们不仅要精度，还要速度。角光谱方法虽然快，但里面依然有可以优化的门道，也藏着一些需要避开的“坑”。

计算效率分析： 两种方法的计算核心都是三次FFT（一次正变换、一次逆变换，脉冲响应法多一次对 h 的变换）。对于N×N的网格，FFT的计算复杂度是O(N² log N)。因此，两种方法在时间复杂度上属于同一量级，传递函数法通常略快一点点，因为它少了一次FFT（但需要计算频域坐标和传递函数矩阵，这部分开销很小）。真正的效率差异微乎其微，精度才是关键区别。

关键参数选择与“避坑”指南：

1. 采样点数N：必须是2的幂次（如256, 512, 1024, 2048）。这不仅是FFT算法最高效的要求，也能保证频域坐标计算的对称性。点数越多，分辨率越高，但计算越慢。需要根据待模拟特征的尺寸来权衡。
2. 物理尺寸L：这是最容易出错的地方。L 必须足够大，要能完整包含源平面光场 U1，并且周围要留有充足的“零填充”区域。为什么？因为角光谱方法在数学上假设光场在计算窗口外是周期延拓的。如果你的光场（比如一个孔）紧贴着窗口边缘，那么周期延拓会导致虚假的“邻居”光场，引入严重的边界误差，表现为衍射图样中出现不该有的干涉条纹。一个经验法则是：L 至少是孔径特征尺寸的3-5倍。
3. 采样间隔delta1与频域范围：delta1 = L / N 决定了空域分辨率。它必须足够小，以分辨光场中最小的特征。同时，它决定了频域的最大范围 f_max = 1/(2*delta1)。这个 f_max 必须大于光场中包含的最高空间频率，否则也会发生混叠。对于传递函数法，你需要确保在 [-f_max, f_max] 的范围内，传递函数 H 被正确计算。
4. 传播距离z的限制：角光谱方法基于标量衍射理论中的角谱传播公式，该公式本身是精确的。但是，当我们用FFT进行数值计算时，实际上隐含了傍轴近似（即传播方向与光轴夹角较小）。如果 z 太小，或者孔径相对于 z 太大，导致衍射角很大，那么角谱方法的高频分量（对应大角度）的传播公式可能不再准确，或者会超出我们频域坐标所能表示的范围。通常，应满足 (孔径尺寸)^2 / (λz) 不太大。

为了更直观，我把这些关键参数的影响整理成了下表：

参数	作用	设置不当的后果	调整建议
采样点数 N	决定空域/频域网格分辨率	N太小：细节丢失，出现“像素化”；N太大：计算慢	从512或1024开始，根据需求倍增
物理尺寸 L	定义计算窗口大小	L太小：边界效应严重，出现周期性误差；L太大：计算资源浪费	确保孔径周围有足够宽的零背景区域（例如，孔径占窗口宽度1/3以内）
传播距离 z	模拟的衍射距离	z太小：不满足傍轴条件，精度下降；z太大：衍射图样扩散，需要更大的L来捕捉	根据实际光学系统设定，并检查 (a^2)/(λz) 的值

一个实用的调试流程： 当我拿到一个新的仿真任务时，通常会这样操作：首先，根据孔径大小预估一个 L（比如孔径宽度的4倍）和一个初始的 N（比如512）。用传递函数法跑一次仿真。然后，将N加倍（例如从512到1024），再跑一次。对比两次结果的中心区域衍射图样。如果主要特征（如主瓣宽度、前几个衍射条纹的位置和强度）没有明显变化，说明采样基本足够。如果变化很大，则需要继续增加N或调整L。这个“收敛性测试”是确保你仿真结果可靠的必要步骤。

最后，关于两种方法的选择，我的个人经验是：对于绝大多数标准的菲涅尔衍射模拟，直接使用基于传递函数的角光谱方法。它代码更简洁，数值上更稳健，精度也更有保障。你可以把脉冲响应法的实现作为理解原理的练习，但在生产代码或严肃的研究中，传递函数法是更可靠的选择。我踩过的坑告诉我，在参数空间的广袤范围内，传递函数法能提供更一致和可预测的精度，让你把精力更多地放在物理问题本身，而不是和数值误差作斗争。