MATLAB角谱法仿真中光斑物理尺寸的精确标定与验证

1. 从“算出来”到“量出来”：一个光学仿真中的经典困惑

如果你用MATLAB做过光学衍射仿真，尤其是用角谱传播（Angular Spectrum Method）配合fft2计算光场传播，大概率会遇到一个让人挠头的问题：屏幕上那个亮斑，它到底有多大？我说的不是像素数，而是它对应的真实物理尺寸，比如毫米或者微米。你按照教科书上的公式，输入波长、传播距离、采样数，一通计算，imagesc一画，一个漂亮的光斑图案出来了。但当你试图回答“这个光斑直径是多少微米”时，往往会卡壳——程序给的是像素坐标，而像素和真实尺寸的映射关系，似乎总有点模糊。

这个问题看似基础，却是连接“数字仿真”与“物理现实”的关键桥梁。我最初接触角谱法时，也在这个问题上栽过跟头。当时我需要模拟一个激光光束经过透镜聚焦后的光斑，用来评估系统的分辨率。代码跑通了，光斑图案也看起来“很物理”，但当我试图将仿真结果与理论公式（比如艾里斑尺寸）对比时，却发现对不上，误差能达到百分之几十。问题就出在fft2变换后，那个横纵坐标轴到底代表什么，以及如何将它校准到真实的物理空间。

网上很多教程和代码示例，重点都放在推导角谱公式和实现fft2计算上，对于“尺寸校准”这一步往往一笔带过，或者默认读者已经清楚。但实际上，这正是新手（甚至一些有经验者）最容易混淆的地方。fft2带来的频率域采样、fftshift导致的坐标重排、以及物理尺寸与像素网格的对应关系，这几个环节环环相扣，任何一个理解偏差都会导致最终结果的尺寸失真。

所以，这篇内容我们就彻底把这个问题掰开揉碎。不谈高深的理论推导，就聚焦在实操层面：当你用MATLAB写完角谱传播的代码后，如何准确无误地确定并标定出仿真光斑的实际物理尺寸。我们会从角谱法的核心计算步骤出发，一步步追踪空间坐标的变换过程，并给出可直接“抄作业”的校准公式和验证方法。

2. 角谱传播与FFT2：核心计算流程与坐标陷阱

要理解尺寸问题，必须先搞清楚角谱传播法中fft2究竟扮演了什么角色，以及数据是如何流动的。角谱法的基本思想是：任何一个光场，都可以分解为不同方向传播的平面波（即角谱）。传播过程，就是让这些平面波在空间频率域乘以一个与传播距离相关的相位因子。

2.1 角谱法的标准计算步骤

假设我们有一个初始光场复振幅U0(x, y)，采样在一个Nx × Ny的网格上，网格的物理采样间隔是dx和dy（单位通常是米或微米）。我们希望计算它传播距离z后的光场Uz(x, y)。标准的角谱法步骤如下：

1. 二维傅里叶变换：对初始光场进行FFT，转换到空间频率域。

   U0_ft = fft2(U0);

   这一步之后，U0_ft是一个复数矩阵，其行和列索引对应的是空间频率分量。

2. 构建传递函数：角谱法的传递函数（或称为频域相位因子）为：H(fx, fy) = exp(i * 2*pi * z * sqrt(1/lambda^2 - fx^2 - fy^2))， 当fx^2 + fy^2 <= 1/lambda^2（传播波）；H(fx, fy) = 0， 当fx^2 + fy^2 > 1/lambda^2（倏逝波，通常忽略）。 其中lambda是波长，fx和fy是空间频率，单位是1/长度。

3. 频域相乘：

   Uz_ft = U0_ft .* H;

4. 逆傅里叶变换：

   Uz = ifft2(Uz_ft);

   得到的Uz就是传播后的光场复振幅。

问题就隐藏在第一步和第二步的衔接处。fft2输出的U0_ft，其矩阵元素的索引(m, n)对应的空间频率(fx, fy)是多少？如果我们不能正确建立(m, n)到(fx, fy)的映射，那么构建的传递函数H就是错的，最终结果自然也是错的。

2.2 FFT2输出的频率坐标约定

这是最关键也最易错的一点。MATLAB的fft2默认的输出频率顺序是“标准顺序”（0到正频率，再到负频率）。对于一个长度为N的序列，经过fft后，其频率分量的顺序是：[0, 1, 2, ..., N/2-1, -N/2, -N/2+1, ..., -1] * (采样频率Fs / N)， 当N为偶数时。

为了在可视化时将零频分量移到频谱中心，我们通常会使用fftshift。但请注意：在角谱法的计算中，我们通常不对U0_ft使用fftshift，而是直接在这种“标准顺序”的频率坐标下构建传递函数H。因为fft2和ifft2是成对出现的，它们内部遵循相同的坐标约定。如果我们对U0_ft做了fftshift，那么在乘以H之前，H也必须按同样的方式移位，并且在ifft2之前还需要用ifftshift将顺序恢复，否则会引入错误。为了简化且避免混淆，最稳妥的做法是：在整个计算流程中，都不使用fftshift，直到最后需要显示光场强度分布时，再对abs(Uz).^2使用fftshift以便将零频（即光轴附近）移到图像中心显示。

那么，如何为“标准顺序”的fft2输出生成正确的空间频率坐标fx和fy呢？

假设 x 方向有Nx个采样点，物理宽度为Lx = Nx * dx。那么 x 方向的采样频率Fs_x = 1/dx。根据 DFT 理论，fft2输出的第m个元素（MATLAB索引从1开始）对应的空间频率fx为：fx(m) = (m-1 - floor(Nx/2)) * (1/Lx)？错！这是fftshift之后的频率坐标。

对于“标准顺序”（即fft2直接输出，未经过fftshift），频率坐标应为：fx(m) = (m-1) / Lx， 当m-1 <= Nx/2时；fx(m) = (m-1 - Nx) / Lx， 当m-1 > Nx/2时。

更简洁的生成方法是使用fftfreq函数的思想，在MATLAB中可以这样实现：

% 生成标准顺序的空间频率向量 fx 和 fy fx = ((0:Nx-1) - floor(Nx/2)) / Lx; % 注意：这样生成的是fftshift后的顺序，不直接适用于H fx = ifftshift(((0:Nx-1) - floor(Nx/2)) / Lx); % 先生成中心化顺序，再ifftshift回标准顺序

但实际上，更直观且不易错的做法是：先构建一个网格，其坐标范围是关于零点对称的（即包含正负频率），然后对这个网格使用ifftshift，使其顺序与fft2输出的标准顺序匹配。

% 正确的构建传递函数H的频率网格方法 Nx = size(U0, 2); % 列数对应x方向 Ny = size(U0, 1); % 行数对应y方向 dx = Lx / Nx; % x方向采样间隔 dy = Ly / Ny; % y方向采样间隔 % 1. 先构建“视觉上中心化”的频率坐标（即零频在中间） fx_centered = (-floor(Nx/2):ceil(Nx/2)-1) / Lx; % 单位: 1/长度 fy_centered = (-floor(Ny/2):ceil(Ny/2)-1) / Ly; % 2. 生成网格 [Fx_centered, Fy_centered] = meshgrid(fx_centered, fy_centered); % 3. 使用 ifftshift 将网格顺序调整为 fft2 的标准输出顺序 Fx = ifftshift(Fx_centered); Fy = ifftshift(Fy_centered); % 4. 此时，Fx和Fy的每个元素，与fft2(U0)输出矩阵的对应位置元素，其空间频率是严格对应的。 % 用它们来构建传递函数H rho2 = Fx.^2 + Fy.^2; % 空间频率半径的平方 k = 2*pi / lambda; H = exp(1i * z * sqrt(k^2 - (2*pi)^2 * rho2)); H(rho2 >= (1/lambda)^2) = 0; % 滤除倏逝波

关键提示：很多错误源于直接用fx = (0:Nx-1)/Lx这种“从0开始”的频率向量去构建网格。这会导致传递函数H的相位梯度方向错误，相当于给光场附加了一个不该有的倾斜相位，使得传播后的光斑位置发生不可预测的偏移，但光斑形状可能看起来“还行”，极具迷惑性。

3. 从像素坐标到物理尺寸：光斑尺寸的标定方法

经过上述步骤，我们得到了传播后的光场Uz。它的强度分布Iz = abs(Uz).^2。现在，Iz是一个Ny×Nx的矩阵。我们想画图，并给坐标轴标上真实的物理尺寸。

3.1 构建物理空间坐标网格

这比频率网格简单。我们的采样网格在实空间是等间距的。假设仿真区域的左下角（或中心）对应物理坐标的原点。通常，我们更关心相对位置，所以常将网格中心设为坐标原点，这样对称性好。

对于 x 方向，物理坐标向量可以这样构建：

x = (-Nx/2 : (Nx/2-1)) * dx; % 当Nx为偶数时，这样设置能使坐标对称且不包含重复的边界点 % 或者更通用的写法： x = ((0:Nx-1) - floor(Nx/2)) * dx;

y 方向同理。然后生成网格：

[X, Y] = meshgrid(x, y);

现在，Iz(i,j)这个像素点的强度，对应的物理位置就是(X(i,j), Y(i,j))。

3.2 可视化与尺寸测量

使用imagesc或pcolor进行可视化时，需要正确传入坐标参数：

imagesc(x*1e6, y*1e6, fftshift(Iz)); % 假设dx,dy单位是米，乘以1e6转换为微米显示 % 使用fftshift是为了将光场中心（通常是能量最强处）显示在图像中心 xlabel('x (um)'); ylabel('y (um)'); axis image; % 保证横纵坐标比例相等，圆斑看起来才是圆的 colorbar;

现在，图像上的坐标轴就是真实的物理尺寸了。你可以直接用鼠标在图上读取光斑的直径（比如从强度降为中心最大值的1/e^2或一半的位置读取距离）。

3.3 程序化测量光斑尺寸

更精确的方法是程序化提取。对于常见的高斯光束，光斑尺寸通常指强度降为中心最大值的1/e^2（约13.5%）处的半径w（束腰半径）。对于衍射极限的艾里斑，则是指第一暗环的半径。

以高斯光束为例：

1. 找到强度矩阵Iz的最大值位置[y0_idx, x0_idx]。
2. 提取通过该中心点的水平线（row = Iz(y0_idx, :)）和垂直线（col = Iz(:, x0_idx)）。
3. 对这些一维强度分布进行高斯拟合，或者直接寻找强度下降到max(Iz)/e^2的位置。
4. 将找到的像素索引差值，乘以采样间隔dx或dy，就得到了物理尺寸。

% 示例：测量水平方向1/e^2半径 I_line = Iz(y0_idx, :); % 中心水平线强度分布 I_max = max(I_line); threshold = I_max / exp(2); % 找到强度大于阈值的区间 above_threshold = I_line > threshold; % 找到区间边界（可能需要处理多个区间，取包含中心点的那个） indices = find(diff(above_threshold)); % ... (具体边界索引处理逻辑) left_idx = indices(1); % 假设处理得到左边界索引 right_idx = indices(2); % 右边界索引 % 计算物理尺寸 spot_radius_x = (right_idx - left_idx) * dx / 2; % 粗略估计半径 spot_diameter_x = spot_radius_x * 2;

注意：这种方法得到的尺寸精度受限于采样间隔dx。如果光斑很小，只覆盖了几个像素，测量误差会很大。这就是为什么在仿真时，需要根据预期的光斑尺寸来合理设置dx和仿真区域大小Lx，确保光斑有足够的像素点来描绘其轮廓。

4. 关键参数的影响与仿真设置经验谈

光斑尺寸的仿真精度，不仅仅取决于坐标标定是否正确，更取决于整个仿真参数的设置是否合理。这里分享几个从踩坑中总结出的经验。

4.1 采样间隔dx与仿真区域Lx：如何平衡精度与计算量？

这是一对矛盾。根据奈奎斯特采样定理，要无混叠地表示一个光场，采样频率必须大于其最高空间频率的两倍。在角谱法中，传递函数H有一个截止频率1/lambda（对应传播波的最大空间频率）。因此，采样间隔必须满足：dx < lambda / 2这是一个非常严格的下限。例如，对于波长633nm的氦氖激光，dx必须小于316.5nm。但这只是保证频率域不混叠。

在实空间，我们需要仿真区域Lx = Nx * dx足够大，以容纳传播后的光场，特别是当光束有较大发散角或传播距离较远时。如果Lx太小，光场能量会“溢出”到计算窗口边界，由于FFT的周期性假设，这会导致“卷绕误差”（aliasing in space），即从一边溢出的光会从另一边绕回来，严重干扰结果。

实操建议：

1. 初始设置：dx可以设为lambda/4到lambda/10之间，这是一个比较安全的范围，既能保证采样足够，又不会让计算量爆炸。
2. 区域大小：Lx至少应大于初始光束尺寸的4-5倍，并且要预估传播后光束的展宽。一个粗略的估计是Lx > 初始光斑尺寸 + 2 * z * NA，其中NA是光束的数值孔径。
3. 动态调整：对于长距离传播，一种稳健的方法是使用“分步传播”，或者采用“吸收边界”（在计算窗口边缘乘以一个渐变的衰减函数）来抑制卷绕误差。但这会引入额外复杂度。对于大多数初步分析和设计，确保窗口足够大是最简单有效的方法。

4.2 传播距离z与菲涅耳数：何时角谱法依然有效？

角谱法是一种严格的标量衍射计算方法，理论上适用于任何传播距离。但是，当传播距离非常小（近场）或非常大（远场/夫琅禾费区）时，可能会遇到数值问题。

• 近场问题：当z非常小，接近零时，传递函数H的相位变化过于剧烈，空间频率高的分量对应的相位因子exp(i*2*pi*z/lambda * sqrt(1-(lambda*fx)^2))中，根号部分接近1，相位随z线性变化。此时对采样要求极高，否则相位采样不足会导致严重误差。实际上，当z小到与波长可比拟时，可能需要考虑矢量衍射效应，标量理论本身可能已不适用。
• 远场简化：当传播距离z足够大，满足夫琅禾费近似条件时，角谱法的传递函数可以简化为一个二次相位因子。更重要的是，此时输出平面的光场分布正比于初始光场傅里叶变换的模平方。这就是标题中“夫琅禾费衍射”的关联点。在MATLAB中，如果你用角谱法计算远场，你会发现，除了一个相位弯曲因子，强度分布Iz确实与abs(fft2(U0)).^2的形状一致。此时，光斑尺寸（如艾里斑半径）有明确的理论公式：半径 = 1.22 * lambda * z / D，其中D是孔径直径。你可以用这个公式来验证你的仿真标定是否正确。

验证步骤：

1. 设置一个圆形孔径的初始光场U0。
2. 选择一个足够大的传播距离z，使其满足夫琅禾费条件（例如，菲涅耳数F = D^2/(4*lambda*z) << 1）。
3. 用你的角谱法程序计算Iz。
4. 测量仿真结果中艾里斑第一暗环的半径r_sim。
5. 计算理论值r_theory = 1.22 * lambda * z / D。
6. 对比r_sim和r_theory。如果两者吻合得很好（比如误差<2%），那么恭喜你，你的坐标标定和仿真设置很可能是正确的。如果不吻合，就需要回头检查频率网格构建、fftshift使用等环节。

4.3 一个完整的、可验证的仿真示例

下面提供一个从参数设置、坐标构建、计算到结果验证的完整代码框架，并附上关键注释。

%% 角谱传播仿真与光斑尺寸标定验证示例 clear; close all; clc; % ==================== 1. 物理参数设置 ==================== lambda = 633e-9; % 波长，单位：米 (He-Ne激光) z = 0.1; % 传播距离，单位：米 (10cm) D = 2e-3; % 圆形孔径直径，单位：米 (2mm) % ==================== 2. 仿真网格参数设置 ==================== % 采样间隔：需要满足奈奎斯特条件 dx < lambda/2，这里取更保守的值 dx = lambda / 8; % x方向采样间隔 dy = dx; % y方向采样间隔，通常设成一样 % 仿真区域大小：需要足够大以容纳远场衍射图样 % 预估远场衍射角 theta ~ lambda/D，光斑半径 ~ z * theta % 因此区域半径至少需要几倍的光斑半径 spot_radius_estimate = 1.22 * lambda * z / D; Lx = 20 * spot_radius_estimate; % 仿真区域宽度，取估计值的20倍以确保充足 Ly = Lx; % 计算采样点数（取为2的整数次幂，FFT效率高） Nx = 2^nextpow2(ceil(Lx / dx)); Ny = 2^nextpow2(ceil(Ly / dy)); % 根据调整后的Nx, Ny，重新计算精确的Lx, Ly和dx, dy Lx = Nx * dx; Ly = Ny * dy; fprintf('仿真参数：\n'); fprintf(' 采样点数 Nx x Ny = %d x %d\n', Nx, Ny); fprintf(' 物理尺寸 Lx x Ly = %.3f mm x %.3f mm\n', Lx*1e3, Ly*1e3); fprintf(' 采样间隔 dx = dy = %.3f um\n', dx*1e6); % ==================== 3. 构建坐标网格 ==================== % 3.1 实空间坐标（以网格中心为原点） x = ((-Nx/2):(Nx/2-1)) * dx; y = ((-Ny/2):(Ny/2-1)) * dy; [X, Y] = meshgrid(x, y); % 3.2 空间频率坐标（用于构建传递函数H） % 注意：这里构建的是与fft2输出顺序匹配的频率网格 fx = ((-Nx/2):(Nx/2-1)) / Lx; % 中心化顺序的频率向量 fy = ((-Ny/2):(Ny/2-1)) / Ly; [Fx_centered, Fy_centered] = meshgrid(fx, fy); % 使用ifftshift调整到fft2的标准输出顺序 Fx = ifftshift(Fx_centered); Fy = ifftshift(Fy_centered); % ==================== 4. 创建初始光场（圆形孔径） ==================== R = sqrt(X.^2 + Y.^2); U0 = double(R <= D/2); % 圆形孔径，内部振幅为1，外部为0 % U0 = U0 .* exp(1i * 0); % 可以添加初始相位，这里为零相位 % ==================== 5. 角谱传播计算 ==================== % 5.1 FFT U0_ft = fft2(U0); % 5.2 构建角谱传递函数 k = 2 * pi / lambda; % 计算频率半径的平方 rho2 = Fx.^2 + Fy.^2; % 传递函数 H = exp(1i * z * sqrt(k^2 - (2*pi)^2 * rho2)); % 滤除倏逝波（空间频率超过1/lambda的分量） H(rho2 >= (1/lambda)^2) = 0; % 5.3 频域相乘并逆变换 Uz_ft = U0_ft .* H; Uz = ifft2(Uz_ft); Iz = abs(Uz).^2; % 传播后的光强 % ==================== 6. 可视化与尺寸标定 ==================== % 6.1 显示结果 figure('Position', [100, 100, 1200, 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, U0); % 显示初始孔径 axis image; xlabel('x (mm)'); ylabel('y (mm)'); title('初始光场 (孔径)'); colormap('gray'); colorbar; subplot(1,3,2); imagesc(x*1e3, y*1e3, fftshift(log10(Iz+1e-10))); % 显示对数尺度光强，fftshift使中心在图中 axis image; xlabel('x (mm)'); ylabel('y (mm)'); title(sprintf('传播后光强 (z=%.2f m)', z)); colormap('jet'); colorbar; % 6.2 提取中心水平线光强分布，用于测量 Iz_shifted = fftshift(Iz); % 为了便于分析，将中心移到数组中心 center_row = round(Ny/2); center_line = Iz_shifted(center_row, :); subplot(1,3,3); plot(x*1e3, center_line / max(center_line), 'b-', 'LineWidth', 1.5); xlabel('x (mm)'); ylabel('归一化光强'); title('中心水平线光强分布'); grid on; % ==================== 7. 光斑尺寸测量与理论验证 ==================== % 7.1 从仿真结果测量艾里斑第一暗环半径（第一个极小值位置） [peaks, locs] = findpeaks(-center_line); % 找极小值（负峰值） if length(locs) >= 1 % 第一个极小值位置（像素索引） first_min_idx = locs(1); % 转换为物理坐标 first_min_x = x(first_min_idx); % 艾里斑半径通常指第一暗环半径，即第一个极小值对应的x坐标绝对值 r_sim = abs(first_min_x); fprintf('\n仿真测量结果：\n'); fprintf(' 第一暗环位置 x = %.4f mm\n', first_min_x*1e3); fprintf(' 艾里斑半径 r_sim = %.4f mm\n', r_sim*1e3); else fprintf('未找到明显的极小值，可能传播距离不够远或窗口太小。\n'); r_sim = NaN; end % 7.2 计算夫琅禾费衍射理论值 r_theory = 1.22 * lambda * z / D; % 艾里斑半径理论公式 fprintf('\n理论计算结果（夫琅禾费近似）：\n'); fprintf(' 艾里斑半径 r_theory = 1.22 * lambda * z / D = %.4f mm\n', r_theory*1e3); % 7.3 比较与误差分析 if ~isnan(r_sim) error_percent = abs(r_sim - r_theory) / r_theory * 100; fprintf('\n误差分析：\n'); fprintf(' 绝对误差 = %.6f mm\n', abs(r_sim - r_theory)*1e3); fprintf(' 相对误差 = %.2f%%\n', error_percent); if error_percent < 5 fprintf(' -> 仿真与理论吻合良好，坐标标定可能正确。\n'); else fprintf(' -> 误差较大，请检查频率网格构建、fftshift使用和采样参数。\n'); fprintf(' 常见原因：仿真区域Lx不足导致卷绕误差，或传播距离z不满足远场条件。\n'); % 检查菲涅耳数 F = D^2 / (4 * lambda * z); fprintf(' 当前菲涅耳数 F = D^2/(4*lambda*z) = %.3f\n', F); if F > 0.01 fprintf(' F > 0.01，可能不满足严格的夫琅禾费条件，会引入误差。\n'); end end end

运行这段代码，你会得到三个子图：初始孔径、传播后的光强分布（对数尺度）、以及中心线光强剖面。通过比较仿真测得的艾里斑半径与理论值，你可以直观地验证你的整个仿真链路（包括尺寸标定）是否正确。如果误差在可接受范围内（比如<2%），那么你对fft2光斑实际尺寸的把握就基本到位了。

这个过程中最容易出错的点，我最后再强调一次：频率网格Fx, Fy的顺序必须与fft2输出的数据顺序严格匹配。使用ifftshift对“视觉中心化”的网格进行处理，是避免这个坑的最有效方法。当你能够用夫琅禾费衍射这个经典案例验证通你的仿真时，再去模拟更复杂的光场（如高斯光束、涡旋光束、经过像差的光束），心里就会踏实很多。毕竟，仿真工作的第一步，永远是确保你的尺子是准的。