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)。标准的角谱法步骤如下二维傅里叶变换对初始光场进行FFT转换到空间频率域。U0_ft fft2(U0);这一步之后U0_ft是一个复数矩阵其行和列索引对应的是空间频率分量。构建传递函数角谱法的传递函数或称为频域相位因子为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/长度。频域相乘Uz_ft U0_ft .* H;逆傅里叶变换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/21, ..., -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束腰半径。对于衍射极限的艾里斑则是指第一暗环的半径。以高斯光束为例找到强度矩阵Iz的最大值位置[y0_idx, x0_idx]。提取通过该中心点的水平线row Iz(y0_idx, :)和垂直线col Iz(:, x0_idx)。对这些一维强度分布进行高斯拟合或者直接寻找强度下降到max(Iz)/e^2的位置。将找到的像素索引差值乘以采样间隔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即从一边溢出的光会从另一边绕回来严重干扰结果。实操建议初始设置dx可以设为lambda/4到lambda/10之间这是一个比较安全的范围既能保证采样足够又不会让计算量爆炸。区域大小Lx至少应大于初始光束尺寸的4-5倍并且要预估传播后光束的展宽。一个粗略的估计是Lx 初始光斑尺寸 2 * z * NA其中NA是光束的数值孔径。动态调整对于长距离传播一种稳健的方法是使用“分步传播”或者采用“吸收边界”在计算窗口边缘乘以一个渐变的衰减函数来抑制卷绕误差。但这会引入额外复杂度。对于大多数初步分析和设计确保窗口足够大是最简单有效的方法。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是孔径直径。你可以用这个公式来验证你的仿真标定是否正确。验证步骤设置一个圆形孔径的初始光场U0。选择一个足够大的传播距离z使其满足夫琅禾费条件例如菲涅耳数F D^2/(4*lambda*z) 1。用你的角谱法程序计算Iz。测量仿真结果中艾里斑第一暗环的半径r_sim。计算理论值r_theory 1.22 * lambda * z / D。对比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(Iz1e-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对“视觉中心化”的网格进行处理是避免这个坑的最有效方法。当你能够用夫琅禾费衍射这个经典案例验证通你的仿真时再去模拟更复杂的光场如高斯光束、涡旋光束、经过像差的光束心里就会踏实很多。毕竟仿真工作的第一步永远是确保你的尺子是准的。