1. 项目概述从“算出来”到“量出来”的鸿沟最近在光学仿真和图像处理社区里一个看似基础但频繁引发讨论的问题又浮出水面了用MATLAB的fft2函数配合角谱传播理论计算出的光场分布屏幕上那个漂亮的光斑它的尺寸到底对应现实世界中的多少毫米这个问题几乎每个从理论迈入实际光学系统设计或实验验证的工程师和研究者都会遇到。你可能会觉得代码跑通了衍射图样也画出来了任务不就完成了吗但当你试图把仿真结果和实验台上的CCD相机拍摄到的光斑进行对比时才会发现两者之间隔着一道“单位换算”的鸿沟。这个问题的核心远不止是MATLAB里一个简单的缩放因子。它触及了离散傅里叶变换DFT的本质、角谱传播的物理假设以及如何将抽象的“像素”世界与真实的“物理尺寸”世界进行校准。网上能找到的代码片段很多但往往只给出了计算步骤却鲜少深入解释每一步参数背后的物理意义和量纲转换逻辑。这就导致很多人“知其然不知其所以然”仿真结果看起来合理但一和实验数据对不上就陷入迷茫。本文将从一个资深光学仿真工程师的视角彻底拆解这个问题。我们不只给出“怎么做”的代码更要深挖“为什么这么做”。我会带你一步步构建从光源到观测面的完整角谱传播模型重点剖析fft2变换后空间频率坐标与物理空间坐标的映射关系并最终回答那个终极问题如何从仿真图像中的一个像素值准确无误地推算出它代表的实际物理尺寸。无论你是正在做课程设计的学生还是需要进行光学系统前期仿真的工程师理解这套校准流程都能让你的仿真工作从“看起来对”提升到“经得起实验检验”的可靠水平。2. 理论基础与核心概念澄清在动手写代码之前我们必须把几个关键概念和它们之间的关系彻底理清。角谱传播方法本身是标量衍射理论的一种数值实现方式而fft2是我们进行快速计算的工具。混淆往往发生在这三者交汇的地方。2.1 角谱传播法在频率域中“行走”的光角谱传播的基本思想非常优美它不直接在空间域处理光波的传播而是将其转换到空间频率域即角谱域。任何一个复杂的光场分布都可以看作是由无数个不同方向传播的平面波即角谱分量叠加而成。传播过程在频率域里就变成了每个平面波分量乘以一个简单的相位因子exp(i*k_z*z)其中k_z是该平面波在传播方向上的波矢分量。用数学公式表达假设初始平面z0的光场复振幅为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是对离散二维序列进行快速傅里叶变换。这里隐藏了两个关键陷阱离散化与采样定理fft2处理的是定义在离散网格点上的数据。我们的物理光场U(x, y)被采样成了一个M×N的矩阵。根据采样定理要想无混叠地表示这个光场我们的采样间隔dx,dy即一个像素代表的物理尺寸必须小于等于光场中最高空间频率分量的一半。这在光学中对应着dx ≤ λ/(2*NA)其中NA是数值孔径。不满足这个条件仿真结果就会出现虚假的条纹和失真。隐含的周期性假设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 (z0平面) % 示例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/21处在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^2FY^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); % 注意对于倏逝波区域H0意味着这些高频分量被滤除。 % 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_idx1), profile_x(left_idx:left_idx1), 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;尺寸校准的核心逻辑定义标尺在仿真开始时我们通过Lx, Ly, M, N定义了物理世界到数字网格的映射关系dx Lx/M。固定网格整个仿真过程中空间采样网格(X, Y)是固定的。无论光场如何传播每个像素点(i,j)始终对应着物理坐标(x(i), y(j))。结果解读仿真输出的光强矩阵Iz其第i行第j列的值就代表了物理位置(x(j), y(i))处的光强注意MATLAB矩阵索引是先行后列。因此当我们测量光斑在矩阵中占据了N_pixel个像素时其对应的物理尺寸就是N_pixel * dx。实操心得永远不要在代码里直接写“像素数”作为尺寸。任何关于距离、宽度、半径的计算都必须通过dx或dy转换为物理单位。养成这个习惯你的仿真结果才能和实验数据直接对话。5. 影响仿真精度与尺寸准确性的关键因素即使校准逻辑正确仿真结果也可能不准。以下是几个需要仔细排查的方面。5.1 采样不足与混叠这是最常见的问题。如果初始光场包含的高频成分例如硬边光阑产生的剧烈变化超出了采样能力就会发生混叠。现象仿真出的衍射图样出现非物理的高频条纹、毛刺或者光斑边缘出现锯齿状。诊断检查dx是否满足dx ≤ λ/2。对于包含硬边界的场景这个条件可能还不够严格需要更小的dx。解决增加采样点数M, N这是最直接的方法但会增加计算负担。增大仿真窗口Lx, Ly在M, N不变的情况下增大窗口也会减小dx但同时会降低对感兴趣区域的空间分辨率因为总像素数不变窗口大了单位距离像素变少。需要权衡。使用零填充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的边缘剖面图查看。解决确保仿真窗口足够大让初始光场及其传播一段距离后的场在到达计算边界前都已衰减到可忽略的程度。通常窗口边长取为初始光斑主要部分的4-8倍以上。使用切趾函数Apodization在窗口边缘乘以一个缓慢衰减到零的函数如 Tukey 窗强制边界值为零。但这会人为修改光场需谨慎使用。采用更高级的边界处理方法如吸收边界条件PML但这在角谱法中不常用更多用于FDTD等时域方法。5.3 相位因子的数值计算问题计算H exp(1i * kz * z)时对于空间频率满足fx^2fy^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 从仿真到实验的标定验证如何确认你的仿真尺寸是准确的一个很好的方法是进行“自洽验证”或与简单解析解对比。夫琅禾费衍射验证当传播距离z很大满足夫琅禾费近似条件时衍射图样是孔径函数的傅里叶变换的模平方。对于一个半径为a的圆孔其夫琅禾费衍射的爱里斑第一暗环半径约为1.22 * λ * z / (2a)。你可以用角谱法仿真一个大z的情况测量光斑第一暗环的半径与解析解对比。如果一致说明你的尺寸校准是正确的。高斯光束传播验证高斯光束在自由空间传播有严格的解析解。仿真一个高斯光束的传播在不同z位置测量其束宽与w(z) w0 * sqrt(1 (z/zr)^2)对比。6.3 常见问题速查表问题现象可能原因排查步骤与解决方法光斑尺寸与预期严重不符1.dx,dy计算或使用错误。2. 频率坐标fx, fy定义错误。3. 传播距离z单位错误。1. 重新检查dx Lx/M的计算和单位。2. 检查dfx 1/Lx的计算。3. 打印出x向量的前几个值和fx向量的前几个值看量级是否合理。4. 确认z是以米为单位输入。衍射图样出现明显锯齿或高频噪声1. 采样不足 (dx太大)。2. 初始光场包含高于奈奎斯特频率的成分如理想硬边。1. 检查dx ≤ λ/2是否满足。2. 尝试增加采样点数M, N。3. 对硬边光阑进行轻微软化如用误差函数过渡。光场在边界处出现“重影”周期性边界条件导致。初始光场在边界处值不为零。1. 增大仿真窗口Lx, Ly使光场在边界处衰减至零。2. 使用零填充。3. 慎用乘以切趾窗函数。仿真结果能量不守恒严重泄露1. 传播掩膜propagation_mask设置过窄滤除了过多有效高频分量。2. 倏逝波处理不当导致能量异常。1. 检查propagation_mask的条件(FX.^2 FY.^2) (1/lambda^2)。2. 尝试不应用掩膜直接计算H观察能量变化。对于远场掩膜是合理的。计算速度慢M, N过大。1. 对于探索性计算可先用较小的M, N如256x256。2. 利用MATLAB的GPU计算功能 (gpuArray) 加速FFT。3. 考虑是否需要全分辨率或许可以降低采样率。光斑位置偏移1. 坐标网格X, Y零点未对齐。2. 使用了fftshift和ifftshift但顺序错误。1. 确保空间域和频率域在使用fftshift前坐标原点定义一致。2. 遵循标准流程fftshift(fft2(ifftshift(U0)))和fftshift(ifft2(ifftshift(Az)))可以确保物理中心对齐。本文为简化未用ifftshift因坐标网格已以零点为中心定义。6.4 性能优化与代码封装建议对于需要反复运行仿真的情况可以考虑以下优化预计算相位因子如果波长λ和传播距离z不变只是改变初始光场U0那么传播相位因子H可以预先计算并存储避免重复计算sqrt和exp。使用单精度如果精度要求可接受使用single精度浮点数可以减半内存占用并提升计算速度。U0 single(U0);函数封装将角谱传播过程封装成一个函数输入(U0, lambda, dx, dy, z)输出Uz。这使代码更清晰易于复用和调试。并行计算如果需要对多个波长或多个距离进行扫描计算可以使用parfor循环进行并行处理。回到最初的问题“FFT2光斑实际尺寸是多少” 答案现在已经很清晰光斑的实际尺寸等于它在仿真图像中占据的像素数量乘以每个像素所代表的物理尺寸dx或dy。而dx和dy是由你定义的仿真窗口物理大小Lx, Ly和采样点数M, N共同决定的 (dx Lx/M)。整个仿真的可信度就建立在你对Lx, M, dx这一系列参数物理意义的深刻理解以及对其取值合理性的反复校验之上。我个人在多年的光学仿真中有一个深刻的体会仿真代码跑通只是第一步让仿真结果具有“物理可度量性”才是从玩具走向工具的关键。每次开始一个新的仿真项目我都会花相当的时间来设计和验证这个“像素-物理尺寸”的映射关系甚至会专门写一个小脚本用已知解析解的特例如圆孔夫琅禾费衍射来验证我的整个参数化流程是否正确。这个习惯帮我避免了许多后续与实验对比时的头疼时刻。最后一个小技巧是在可视化时养成使用imagesc(x, y, I)并正确设置xlabel和ylabel的习惯让图中的坐标轴直接显示物理单位这样任何看到图的人包括未来的你自己都能立刻理解图中的尺度而不需要再去代码里翻找换算关系。