MATLAB角谱传播仿真：从像素到物理尺寸的精确校准

1. 项目概述：从“算出来”到“量出来”的鸿沟

最近在光学仿真和图像处理社区里，一个看似基础但频繁引发讨论的问题又浮出水面了：用MATLAB的fft2函数配合角谱传播理论计算出的光场分布，屏幕上那个漂亮的光斑，它的尺寸到底对应现实世界中的多少毫米？这个问题，几乎每个从理论迈入实际光学系统设计或实验验证的工程师和研究者都会遇到。你可能会觉得，代码跑通了，衍射图样也画出来了，任务不就完成了吗？但当你试图把仿真结果和实验台上的CCD相机拍摄到的光斑进行对比时，才会发现两者之间隔着一道“单位换算”的鸿沟。

这个问题的核心，远不止是MATLAB里一个简单的缩放因子。它触及了离散傅里叶变换（DFT）的本质、角谱传播的物理假设，以及如何将抽象的“像素”世界与真实的“物理尺寸”世界进行校准。网上能找到的代码片段很多，但往往只给出了计算步骤，却鲜少深入解释每一步参数背后的物理意义和量纲转换逻辑。这就导致很多人“知其然，不知其所以然”，仿真结果看起来合理，但一和实验数据对不上，就陷入迷茫。

本文将从一个资深光学仿真工程师的视角，彻底拆解这个问题。我们不只给出“怎么做”的代码，更要深挖“为什么这么做”。我会带你一步步构建从光源到观测面的完整角谱传播模型，重点剖析fft2变换后，空间频率坐标与物理空间坐标的映射关系，并最终回答那个终极问题：如何从仿真图像中的一个像素值，准确无误地推算出它代表的实际物理尺寸。无论你是正在做课程设计的学生，还是需要进行光学系统前期仿真的工程师，理解这套校准流程，都能让你的仿真工作从“看起来对”提升到“经得起实验检验”的可靠水平。

2. 理论基础与核心概念澄清

在动手写代码之前，我们必须把几个关键概念和它们之间的关系彻底理清。角谱传播方法本身是标量衍射理论的一种数值实现方式，而fft2是我们进行快速计算的工具。混淆往往发生在这三者交汇的地方。

2.1 角谱传播法：在频率域中“行走”的光

角谱传播的基本思想非常优美：它不直接在空间域处理光波的传播，而是将其转换到空间频率域（即角谱域）。任何一个复杂的光场分布，都可以看作是由无数个不同方向传播的平面波（即角谱分量）叠加而成。传播过程，在频率域里就变成了每个平面波分量乘以一个简单的相位因子exp(i*k_z*z)，其中k_z是该平面波在传播方向上的波矢分量。

用数学公式表达，假设初始平面（z=0）的光场复振幅为U0(x, y)，其二维傅里叶变换（角谱）为A0(fx, fy)。那么传播到距离为z的平面后，其角谱A_z(fx, fy)为：A_z(fx, fy) = A0(fx, fy) * exp(i * 2π * z * sqrt(1/λ² - fx² - fy²))这里(fx, fy)就是空间频率坐标，单位是长度分之一（如 mm⁻¹）。之后，再对A_z做一次逆傅里叶变换，就得到了传播后的空间光场U_z(x, y)。

为什么选择角谱法？相比另一种常用的菲涅尔衍射积分（FFT实现），角谱法在理论上是一种精确的标量衍射计算方法，没有近场近似。只要采样满足奈奎斯特条件，它可以准确计算从很近到很远的传播距离。这对于需要高精度仿真，或者传播距离与孔径尺寸可比拟的情况尤为重要。

2.2 FFT2的“双刃剑”：离散化与周期性

MATLAB中的fft2是对离散二维序列进行快速傅里叶变换。这里隐藏了两个关键陷阱：

1. 离散化与采样定理：fft2处理的是定义在离散网格点上的数据。我们的物理光场U(x, y)被采样成了一个M×N的矩阵。根据采样定理，要想无混叠地表示这个光场，我们的采样间隔dx,dy（即一个像素代表的物理尺寸）必须小于等于光场中最高空间频率分量的一半。这在光学中对应着：dx ≤ λ/(2*NA)，其中NA是数值孔径。不满足这个条件，仿真结果就会出现虚假的条纹和失真。

2. 隐含的周期性假设：fft2默认其输入输出序列是周期性的。这意味着，你计算的U0矩阵，在fft2看来，其左右、上下都是首尾相连的无限重复。如果你的光场在边界处不是平滑衰减到零（例如一个硬边光阑），那么这种周期性就会在边界引入虚假的衍射效应，相当于光从一边衍射出去，又从对边绕回来。这是仿真中产生“鬼影”和误差的一个重要来源。通常的解决办法是使用“零填充”（Zero-padding），即在计算阵列周围补上一圈零，给衍射场足够的扩展空间，避免周期性折叠的影响。

2.3 空间坐标与频率坐标的映射关系

这是连接“像素索引”与“物理世界”的桥梁，也是回答光斑尺寸问题的核心。假设我们有一个M×N的采样网格，物理尺寸为Lx = M * dx,Ly = N * dy。

• 空间域坐标：经过fft2和ifft2变换后，MATLAB输出的矩阵，其行、列索引直接对应着空间位置。通常，我们会用fftshift将零频分量移到中心，这样坐标范围是[-floor(M/2) : floor((M-1)/2)] * dx和[-floor(N/2) : floor((N-1)/2)] * dy。这里的dx和dy就是每个像素对应的物理尺寸。

• 频率域坐标：fft2输出的矩阵，经过fftshift后，其行、列索引对应着空间频率。频率坐标的间隔（即频率分辨率）为dfx = 1 / Lx,dfy = 1 / Ly。频率坐标的范围是[-floor(M/2) : floor((M-1)/2)] * dfx和[-floor(N/2) : floor((N-1)/2)] * dfy。

关键推导：角谱传播公式中的exp(i * 2π * z * sqrt(1/λ² - fx² - fy²))要求fx² + fy² ≤ 1/λ²，否则根号内为负，对应倏逝波，在远场可忽略。因此，有效的仿真范围受限于|fx| ≤ 1/λ。这反过来对我们的采样提出了要求：为了能表征最高空间频率1/λ，根据采样定理，我们需要dx ≤ λ/2。这是一个非常重要的自检条件。

注意：很多初学者直接套用网上代码，却从不检查自己的dx是否满足dx ≤ λ/2。如果不满足，那么你的仿真从采样层面就已经失真了，后续计算再精确也无意义。这是第一个必须养成的习惯。

3. 从零构建角谱传播模型：参数设计与初始化

理论清晰后，我们开始动手构建模型。一个好的仿真，始于严谨的参数定义和初始化。

3.1 定义物理参数与仿真窗口

我们以一个典型的激光高斯光束通过圆形光阑衍射为例。首先，定义所有物理参数，并附上注释说明其单位和意义。

%% 1. 物理参数定义 lambda = 632.8e-9; % 波长，单位：米 (He-Ne激光) k = 2 * pi / lambda; % 波数 % 光源参数（假设为高斯光束） w0 = 1e-3; % 光束束腰半径，单位：米 z_prop = 0.5; % 传播距离，单位：米 % 仿真窗口参数 Lx = 10e-3; % 仿真窗口物理宽度，单位：米 (10mm) Ly = 10e-3; % 仿真窗口物理高度，单位：米 M = 1024; % x方向采样点数（建议为2的幂次，FFT效率高） N = 1024; % y方向采样点数 % 计算采样间隔 dx = Lx / M; % x方向单个像素的物理尺寸 dy = Ly / N; % y方向单个像素的物理尺寸 % ！！！关键检查：采样间隔是否满足奈奎斯特条件？ if dx > lambda/2 || dy > lambda/2 warning('采样间隔(%.2e m)大于 lambda/2(%.2e m)。可能发生混叠，建议增加采样点数M/N以减小dx/dy。', dx, lambda/2); end

参数设计心得：

• 仿真窗口大小Lx, Ly：它必须大于你关心的光场分布区域，并预留足够空间给衍射场扩展，避免周期性边界效应。通常取为初始光斑尺寸的4-8倍。可以先估算，后调整。
• 采样点数M, N：在内存允许范围内，越多越好。增加点数可以减小dx，提高空间分辨率，同时也能提供更大的频率域范围。但点数翻倍，计算量和内存消耗会显著增加。1024×1024是一个兼顾精度和效率的常用起点。
• dx, dy的计算：这是整个仿真的标尺。后续所有从“像素数”到“物理尺寸”的转换，都依赖于这两个值。务必准确理解dx = Lx / M的含义：它表示在生成的图像中，相邻两个像素点所代表的实际空间距离。

3.2 构建空间坐标网格与初始光场

基于上面定义的参数，我们创建空间坐标网格，并生成初始光场（这里以高斯光束为例）。

%% 2. 创建空间坐标网格 % 生成一维坐标向量 x = (-M/2 : M/2-1) * dx; % 坐标原点在中心 y = (-N/2 : N/2-1) * dy; % 创建二维网格 [X, Y] = meshgrid(x, y); %% 3. 构建初始光场 U0 (z=0平面) % 示例1：理想高斯光束 U0_gauss = exp(-(X.^2 + Y.^2) / w0^2); % 示例2：高斯光束通过圆形光阑 aperture_radius = 2e-3; % 光阑半径，2mm aperture = (X.^2 + Y.^2) <= aperture_radius^2; U0_apertured = U0_gauss .* aperture; % 选择初始光场 U0 = U0_apertured; % 这里使用加光阑的例子 % 可视化初始光场和光阑 figure; subplot(1,2,1); imagesc(x*1e3, y*1e3, abs(U0).^2); % 显示光强，坐标转换为mm axis image; xlabel('x (mm)'); ylabel('y (mm)'); title('初始光场光强 |U0|^2'); colorbar; subplot(1,2,2); imagesc(x*1e3, y*1e3, aperture); axis image; xlabel('x (mm)'); ylabel('y (mm)'); title('光阑形状'); colorbar; colormap('gray');

坐标网格构建的细节： 注意我们使用(-M/2 : M/2-1)来生成坐标。这是因为当M为偶数时（通常如此），fft2输出的零频分量位于索引M/2+1处（在MATLAB中索引从1开始）。fftshift操作后，零频会移到矩阵中心。我们这样定义x和y，使得网格坐标X和Y的中心就是(0,0)，与fftshift后的频谱中心对齐，这在物理上更直观，也便于后续分析。

4. 角谱传播的核心计算与尺寸校准

这是最核心的部分，我们将实现角谱传播，并重点解决尺寸校准问题。

4.1 计算角谱与传播相位因子

%% 4. 角谱传播核心计算 % 4.1 对初始光场进行二维傅里叶变换，得到角谱 A0 = fft2(U0); % 4.2 创建空间频率坐标网格 (fx, fy) % 频率坐标间隔（频率分辨率） dfx = 1 / Lx; % fx方向频率间隔 dfy = 1 / Ly; % fy方向频率间隔 % 生成一维频率向量（同样以零频为中心） fx = (-M/2 : M/2-1) * dfx; fy = (-N/2 : N/2-1) * dfy; % 创建二维频率网格 [FX, FY] = meshgrid(fx, fy); % 4.3 计算传播相位因子 H % 角谱传播的传递函数 kz = 2*pi * sqrt((1/lambda)^2 - FX.^2 - FY.^2); % 计算 k_z % 对于满足 (FX^2+FY^2) > 1/lambda^2 的点，kz为虚数，对应倏逝波。 % 在大多数远场或非超分辨情况下，这些分量迅速衰减，可以强制置零或保留（但传播后会衰减）。 H = exp(1i * kz * z_prop); % 更稳健的做法：只对传播波分量应用相位因子 % 创建一个掩膜，标识传播波分量（空间频率在“光锥”内） propagation_mask = (FX.^2 + FY.^2) < (1/lambda^2); H = propagation_mask .* exp(1i * kz .* propagation_mask * z_prop); % 注意：对于倏逝波区域，H=0，意味着这些高频分量被滤除。 % 4.4 在频率域应用传播因子，并变换回空间域 Az = A0 .* H; Uz = ifft2(Az); % 注意：经过fft2和ifft2后，结果会有一个整体的缩放因子 (M*N)。 % 对于能量守恒的检查很重要，但通常我们更关心相对分布。 % 如果需要严格的能量校准，需要除以 sqrt(M*N)。这里先不处理，专注于尺寸。

4.2 关键步骤：从像素到物理尺寸的校准

现在，我们有了传播后的光场复振幅Uz。它是一个M×N的矩阵。我们如何知道第(i,j)个像素的强度值，对应着实际空间中哪个位置的光强呢？答案就在我们最初定义的x和y坐标向量里。

%% 5. 结果可视化与物理尺寸标注 % 计算观测面光强分布 Iz = abs(Uz).^2; % 找到光斑中心（用于后续分析） [~, max_index] = max(Iz(:)); [center_y_idx, center_x_idx] = ind2sub(size(Iz), max_index); % 注意MATLAB是行(y)列(x) % 创建与Uz对应的物理坐标网格（与初始坐标相同，因为采样网格没变） % X_grid 和 Y_grid 已经定义，直接使用。 % 绘制光强分布图，并标注物理坐标轴 figure; imagesc(x*1e3, y*1e3, Iz); % x, y 单位转换为毫米进行显示 axis image; xlabel('x (mm)'); % 坐标轴标签直接显示物理单位 ylabel('y (mm)'); title(sprintf('传播距离 z = %.2f m 处的光强分布', z_prop)); colorbar; colormap('hot'); % 在图上叠加一条通过光斑中心的水平线，用于提取一维剖面 hold on; plot(x*1e3, y(center_y_idx)*ones(size(x))*1e3, 'w--', 'LineWidth', 1); hold off; %% 6. 光斑尺寸计算（以束腰或半高全宽FWHM为例） % 提取通过光斑中心的一维光强剖面（水平方向） horizontal_profile = Iz(center_y_idx, :); % 对应的物理坐标就是 x 向量 profile_x = x * 1e3; % 转换为毫米 % 方法1：高斯拟合求束腰半径（适用于高斯或近高斯光斑） % 使用曲线拟合工具箱或自定义拟合函数 % 这里简单演示使用 findpeaks 和半高宽计算 [peak_intensity, peak_idx] = max(horizontal_profile); half_max = peak_intensity / 2; % 找到强度降至一半的点（更稳健的方法是插值） left_idx = find(horizontal_profile(1:peak_idx) <= half_max, 1, 'last'); right_idx = peak_idx + find(horizontal_profile(peak_idx:end) <= half_max, 1, 'first') - 1; if ~isempty(left_idx) && ~isempty(right_idx) % 线性插值以提高精度 x_left = interp1(horizontal_profile(left_idx:left_idx+1), profile_x(left_idx:left_idx+1), half_max, 'linear'); x_right = interp1(horizontal_profile(right_idx-1:right_idx), profile_x(right_idx-1:right_idx), half_max, 'linear'); fwhm_x = abs(x_right - x_left); % 半高全宽 (FWHM)，单位：mm % 对于高斯光束，束腰半径 w = FWHM / (2*sqrt(ln2)) ≈ FWHM / 1.665 w_x = fwhm_x / (2*sqrt(log(2))); fprintf('计算得到的光斑尺寸（水平方向）：\n'); fprintf(' 半高全宽 FWHM = %.3f mm\n', fwhm_x); fprintf(' 等效高斯束腰半径 w = %.3f mm\n', w_x); fprintf(' （注意：此结果基于当前采样，精度受dx=%.3f mm限制）\n', dx*1e3); end % 将关键尺寸标注在图上 figure(gcf); % 回到当前图窗 hold on; plot([x_left, x_right], [y(center_y_idx)*1e3, y(center_y_idx)*1e3], 'b-', 'LineWidth', 2, 'DisplayName', sprintf('FWHM=%.2fmm', fwhm_x)); legend('Location', 'best'); hold off;

尺寸校准的核心逻辑：

1. 定义标尺：在仿真开始时，我们通过Lx, Ly, M, N定义了物理世界到数字网格的映射关系dx = Lx/M。
2. 固定网格：整个仿真过程中，空间采样网格(X, Y)是固定的。无论光场如何传播，每个像素点(i,j)始终对应着物理坐标(x(i), y(j))。
3. 结果解读：仿真输出的光强矩阵Iz，其第i行第j列的值，就代表了物理位置(x(j), y(i))处的光强（注意MATLAB矩阵索引是先行后列）。因此，当我们测量光斑在矩阵中占据了N_pixel个像素时，其对应的物理尺寸就是N_pixel * dx。

实操心得：永远不要在代码里直接写“像素数”作为尺寸。任何关于距离、宽度、半径的计算，都必须通过dx或dy转换为物理单位。养成这个习惯，你的仿真结果才能和实验数据直接对话。

5. 影响仿真精度与尺寸准确性的关键因素

即使校准逻辑正确，仿真结果也可能不准。以下是几个需要仔细排查的方面。

5.1 采样不足与混叠

这是最常见的问题。如果初始光场包含的高频成分（例如硬边光阑产生的剧烈变化）超出了采样能力，就会发生混叠。

• 现象：仿真出的衍射图样出现非物理的高频条纹、毛刺，或者光斑边缘出现锯齿状。
• 诊断：检查dx是否满足dx ≤ λ/2。对于包含硬边界的场景，这个条件可能还不够严格，需要更小的dx。
• 解决：
  1. 增加采样点数M, N：这是最直接的方法，但会增加计算负担。
  2. 增大仿真窗口Lx, Ly：在M, N不变的情况下，增大窗口也会减小dx，但同时会降低对感兴趣区域的空间分辨率（因为总像素数不变，窗口大了，单位距离像素变少）。需要权衡。
  3. 使用零填充（Zero-padding）：在计算FFT前，将初始光场矩阵用零包围。这并不能增加原始数据的频率信息，但可以增加频率域的采样点数，使频谱看起来更平滑，有时能缓解因周期性边界条件引起的混叠。在角谱传播中，零填充也为衍射场提供了扩展空间。

% 零填充示例：将U0放置在更大的零矩阵中心 M_pad = 2*M; % 填充后的大小 N_pad = 2*N; U0_padded = padarray(U0, [(M_pad-M)/2, (N_pad-N)/2], 0, 'both'); % 注意：填充后，dx, dy不变，但总网格点数增加，频率分辨率dfx, dfy会变好。 % 后续计算需要使用新的 M_pad, N_pad, Lx, Ly (Lx, Ly 是否变？这里概念容易混淆) % 更严谨的做法是：保持物理窗口Lx, Ly不变，通过增加M,N来减小dx。零填充通常用于处理卷积/相关时的边界效应。

5.2 周期性边界条件导致的误差

如前所述，fft2的周期性假设会带来问题。

• 现象：如果初始光场在边界处强度不为零，你会看到光好像从一边“绕”到了另一边，在对面边界处出现虚假的光场。
• 诊断：观察初始光场U0在仿真窗口边缘的值是否已经衰减到接近零。绘制U0的边缘剖面图查看。
• 解决：
  1. 确保仿真窗口足够大：让初始光场及其传播一段距离后的场，在到达计算边界前都已衰减到可忽略的程度。通常窗口边长取为初始光斑主要部分的4-8倍以上。
  2. 使用切趾函数（Apodization）：在窗口边缘乘以一个缓慢衰减到零的函数（如 Tukey 窗），强制边界值为零。但这会人为修改光场，需谨慎使用。
  3. 采用更高级的边界处理方法，如吸收边界条件（PML），但这在角谱法中不常用，更多用于FDTD等时域方法。

5.3 相位因子的数值计算问题

计算H = exp(1i * kz * z)时，对于空间频率满足fx^2+fy^2 > 1/λ^2的点，kz是虚数。直接计算会导致这些高频分量的幅度指数增长或衰减。

• 处理策略：
  • 方法A（滤波法）：如前面代码所示，创建一个传播掩膜propagation_mask，只对传播波分量应用相位因子，对倏逝波分量直接置零。这是最常用且物理上合理的方法，相当于一个理想低通滤波器。
  • 方法B（保留衰减）：即使对于倏逝波，也计算exp(1i * kz * z)，其中kz为虚数，结果是一个实指数衰减项。这对于模拟近场超分辨效应可能是必要的，但通常衰减极快，在远场可忽略。
  • 注意数值稳定性：当z很大时，exp(1i*kz*z)中的相位变化非常剧烈，可能导致数值精度问题。但角谱法本身对远场计算是稳定的。

5.4 能量守恒检查

一个快速的 sanity check 是计算传播前后的总光强（能量）是否大致守恒。注意，由于离散采样和数值误差，完全守恒很难，但不应有数量级上的差异。

energy_before = sum(abs(U0(:)).^2) * dx * dy; % 近似积分 energy_after = sum(abs(Uz(:)).^2) * dx * dy; fprintf('传播前能量: %.6e\n', energy_before); fprintf('传播后能量: %.6e\n', energy_after); fprintf('能量比 (后/前): %.6f\n', energy_after/energy_before);

如果能量比严重偏离1（如<0.9或>1.1），就需要回头检查采样、边界条件或相位因子计算是否正确。

6. 进阶应用与常见问题排查

掌握了基础模型后，可以应对更复杂的场景，并快速定位问题。

6.1 处理倾斜入射与离轴情况

角谱法的一个巨大优势是易于处理倾斜入射的平面波。只需在初始光场U0上乘以一个倾斜相位因子即可。

% 假设平面波入射方向与z轴夹角为 theta_x, theta_y (弧度) theta_x = 5 * pi/180; % 5度 theta_y = 0; % 对应的初始相位倾斜 tilt_phase = exp(1i * k * (sin(theta_x)*X + sin(theta_y)*Y)); U0_tilted = U0 .* tilt_phase; % 倾斜入射的初始光场 % 后续角谱传播计算完全不变

这相当于在角谱（频率域）中将整个频谱平移(sin(theta_x)/lambda, sin(theta_y)/lambda)。计算后，光斑的中心位置会相应移动。

6.2 从仿真到实验的标定验证

如何确认你的仿真尺寸是准确的？一个很好的方法是进行“自洽验证”或与简单解析解对比。

1. 夫琅禾费衍射验证：当传播距离z很大，满足夫琅禾费近似条件时，衍射图样是孔径函数的傅里叶变换的模平方。对于一个半径为a的圆孔，其夫琅禾费衍射的爱里斑第一暗环半径约为1.22 * λ * z / (2a)。你可以用角谱法仿真一个大z的情况，测量光斑第一暗环的半径，与解析解对比。如果一致，说明你的尺寸校准是正确的。
2. 高斯光束传播验证：高斯光束在自由空间传播有严格的解析解。仿真一个高斯光束的传播，在不同z位置测量其束宽，与w(z) = w0 * sqrt(1 + (z/zr)^2)对比。

6.3 常见问题速查表

6.4 性能优化与代码封装建议

对于需要反复运行仿真的情况，可以考虑以下优化：

1. 预计算相位因子：如果波长λ和传播距离z不变，只是改变初始光场U0，那么传播相位因子H可以预先计算并存储，避免重复计算sqrt和exp。
2. 使用单精度：如果精度要求可接受，使用single精度浮点数可以减半内存占用并提升计算速度。U0 = single(U0);
3. 函数封装：将角谱传播过程封装成一个函数，输入(U0, lambda, dx, dy, z)，输出Uz。这使代码更清晰，易于复用和调试。
4. 并行计算：如果需要对多个波长或多个距离进行扫描计算，可以使用parfor循环进行并行处理。

回到最初的问题：“FFT2光斑实际尺寸是多少？” 答案现在已经很清晰：光斑的实际尺寸等于它在仿真图像中占据的像素数量，乘以每个像素所代表的物理尺寸dx或dy。而dx和dy是由你定义的仿真窗口物理大小Lx, Ly和采样点数M, N共同决定的 (dx = Lx/M)。整个仿真的可信度，就建立在你对Lx, M, dx这一系列参数物理意义的深刻理解，以及对其取值合理性的反复校验之上。

我个人在多年的光学仿真中有一个深刻的体会：仿真代码跑通只是第一步，让仿真结果具有“物理可度量性”才是从玩具走向工具的关键。每次开始一个新的仿真项目，我都会花相当的时间来设计和验证这个“像素-物理尺寸”的映射关系，甚至会专门写一个小脚本，用已知解析解的特例（如圆孔夫琅禾费衍射）来验证我的整个参数化流程是否正确。这个习惯帮我避免了许多后续与实验对比时的头疼时刻。最后一个小技巧是，在可视化时，养成使用imagesc(x, y, I)并正确设置xlabel和ylabel的习惯，让图中的坐标轴直接显示物理单位，这样任何看到图的人（包括未来的你自己）都能立刻理解图中的尺度，而不需要再去代码里翻找换算关系。