简介一份围绕高斯光束Matlab仿真的完整技术文档适合激光原理、光学工程等课程的学生与科研人员参考。文档从高斯光束数学模型出发推导强度分布公式并给出谐振腔中位置处的归一化强度分布仿真方法覆盖二维、三维强度分布图绘制。同时结合CCD采集的真实激光光斑图像通过imread读取强度数据与理论高斯曲线进行对比展示数学仿真与实际测量的差异。进一步利用角谱衍射法模拟高斯光束在自由空间传播过程中光强、光斑有效截面半径及等相位面的变化包含核心Matlab代码与参数输入示例方便读者复现。资源为单份docx文档体积432KB内容结构清晰既有公式推导也有程序实现可直接用于实验报告撰写或课程设计参考。已有461人学习下载适合需要快速完成高斯光束仿真作业或理解光场传播原理的读者。1. 高斯光束的Matlab仿真先画强度分布再谈传播用Matlab画高斯光束十行代码就能出一张三维图但把CCD拍到的光斑和理论曲线叠在一起时峰值位置、半宽、背景噪声都会对不上三维图还会有一边整体竖起来。问题往往不在公式而在矩阵读取方向、灰度通道和白边处理。下面这套流程来自激光原理课程里的高斯光束仿真作业解决两个具体问题一是用imread读取实际光斑图像生成二维、三维归一化强度分布并与理论高斯曲线对比二是用角谱传播模型仿真光束传播不同距离后光斑半径和峰值强度的变化。适合正在做激光实验、光电课程设计或者想搞懂Matlab图像处理与傅里叶光学仿真边界的读者。2. 高斯光束数学模型与二维强度矩阵的生成2.1 基模高斯光束的强度表达式与参数约定谐振腔中常见的基模场直角坐标下振幅分布写作E(x,y)E0 exp(-(x^2y^2)/w0^2)对应的强度分布是 |E|^2写成I(x,y)I0 exp(-2(x^2y^2)/w0^2)这里的 w0 是束腰半径在振幅定义里是振幅下降到中心值 1/e 的位置在强度定义里1/e^2 处对应同样半径。很多Matlab示例直接用 exp(-r^2/w0^2) 画“强度分布”严格说那画的是振幅场。如果只做归一化展示两者形状相近但和实测光斑做半宽对比时强度曲线的半高位置会在 sqrt(2) 倍处出现偏差。做这一题之前先把要画的是振幅还是强度定下来后面所有公式、阈值和束宽计算才不会乱。这份课程作业文档资料里的 M 文件把理论曲线写成a2exp(-x2.^2/10)形式上也是 exp 函数但它没有按严格的 1/e^2 强度定义写更像是一个拟合用高斯函数。作业里可以沿用但在程序注释里最好写明w_fitsqrt(10)而不是w0避免后续传播仿真和二维图对比时概念混淆。2.2 用meshgrid生成二维高斯强度场二维高斯分布最简单的写法是用 meshgrid 构建网格再把径向距离平方代入 exp。下面的代码生成 224×224 的归一化强度场和后面 CCD 图像尺寸保持一致。N 224; % 与CCD图像行数保持一致方便对比 L 10e-3; % 仿真区域边长单位m wo 1e-3; % 束腰半径单位m dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; % 以0为中心的横坐标 [X, Y] meshgrid(x, x); % 若按强度分布指数项为 -2*(X.^2Y.^2)/wo^2 Gau exp(-(X.^2 Y.^2) / wo^2); % 振幅场 Gau_norm Gau / max(Gau(:)); % 归一化 figure; surf(X*1e3, Y*1e3, Gau_norm, EdgeColor, none); xlabel(x (mm)); ylabel(y (mm)); zlabel(归一化强度);这段代码里N224不是必须的后面做 FFT 传播时也会用 500。关键是(-N/2 : N/2-1)这种写法让光斑中心落在网格中央而不是数组的(1,1)位置如果直接用1:N后续fft2和surf都要处理中心偏移麻烦得多。dxL/N决定空间分辨率L取 10mm 时单个像素相当于约 44.6µm这个量级和实验中 CCD 的像素尺寸接近适合做毫米级光斑仿真。绘图时乘 1e3 是为了把米换成毫米否则横轴刻度会显示成0.001这类小数值。2.3 参数表波长、束腰、瑞利长度和采样点传播仿真中需要提前确定的参数如下。参数符号本资源中的值作用波长lambda0.568 µm决定波数 k0 和瑞利长度束腰半径wo1 mm初始光斑尺寸二维图半宽来源瑞利长度z_R约 5530 mm近场/远场分界传播 z 的标尺采样点数N100500影响 FFT 精度太小会混叠仿真区域L10 mm限制频域分辨率为 1/L影响传播相位精度瑞利长度由zR k0*wo^2/2计算代入k02*pi/lambda后等于pi*wo^2/lambda。用上面的参数算出来约 5.53m所以文档中z100000mm即 100m 时已经远大于瑞利长度光斑半径接近线性扩展峰值强度明显下降。运行结果里峰值从 1 降到 0.21是这个参数下的正常物理趋势不是程序错误。2.4 为什么要做归一化CCD 读出的灰度值不是光功率绝对值单位是“数”和理论公式中的 I0 没有可比性。所以画二维强度分布图前一定要把数据除以自身最大值。更稳妥的做法是先减背景再归一化I_norm(I-min(I))/(max(I)-min(I))。原始光斑照片如果有暗电流或环境光底部不是 0直接除以最大值会把背景也纳入比较导致理论曲线和实验曲线的边缘不贴合。这一行处理在后面的 imread 实验数据中同样适用。3. 用imread读取CCD光斑二维/三维强度分布与白边处理3.1 读取图片前先看矩阵维度imread是 Matlab 读取位图最常用的命令但读到的是什么尺寸取决于图片本身是灰度图还是 RGB 图。A imread(D:\documents\作业\激光原理与应用\高斯.bmp); disp(size(A)); imshow(A); axis off; title(CCD采集的高斯光束光强分布);如果size(A)返回224 244说明是单通道灰度图如果返回224 244 3说明是 RGB 三通道后面所有取行、取列的操作都必须指定第三维下标否则A(:,122)拿到的是 244×3 的矩阵plot会画出三条线。原文档写的是 224×244 矩阵但代码里又出现了[high, width, color] size(A)说明实际文件很可能是 RGB 格式。常见做法是先用size确认然后取第一通道或者用rgb2gray转成灰度。imshow对uint8数据会自动按灰度范围映射直接看即可不要先用double转换再 imshow否则显示结果可能变黑。3.2 提取直径方向数据绘制二维强度分布原始程序的读取核心是两行A1 A(:, 122); % 取第122列长度224 x1 1:1:224; figure; plot(x1, A1, LineWidth, 1, Color, r); title(理论高斯曲线);这里要特别注意“中间一行”和“中间一列”的差别。题目文字说的是“取中间一行122行”但代码A(:,122)取的是第 122 列。如果矩阵是 224 行 × 244 列那么第 122 列对应的剖面长度是 224和第 122 行对应的剖面长度 244 并不相同。两个方向理论上都可以反映高斯剖面但横轴刻度必须和取的方向一致。更好的写法是if size(A,3) 3 A A(:,:,1); % 取单通道后续都是二维矩阵 end line_data double(A(122, :)); % 第122行长度244 line_data line_data / max(line_data); % 归一化到[0,1] Nx length(line_data); xpix (1:Nx) - (Nx1)/2; % 以图像中心为原点 figure; plot(xpix, line_data, r, LineWidth, 1);double转换是必须的因为imread返回的是uint8直接做除法时 Matlab 会按整数运算取整。归一化后纵轴可以统一到 01方便和理论曲线放在同一张图里。3.3 理论高斯曲线与实验曲线的叠加文档里的理论曲线是x2-100:1:100; a2exp(-x2.^2/10)。它的 x 范围是 [-100,100]和实验曲线的像素范围 [1,224] 不在同一坐标基准上所以原图把两条曲线分开画看起来各自都像高斯叠在一起却对不上。叠加时要做两件事把实验曲线中心移到 0再让理论曲线的横轴范围和实验曲线匹配。x_theory -100:0.1:100; I_theory exp(-x_theory.^2 / 10); figure; plot(xpix, line_data, r, LineWidth, 1); hold on; plot(x_theory, I_theory, b--, LineWidth, 1); legend(实验数据, 理论曲线); xlabel(相对中心像素); ylabel(归一化强度);如果发现实验峰值不在 x0 附近说明光斑中心没有落在图像中心可以用[~, idx] max(line_data)把峰值位置平移到 0 后再比较。像素坐标到毫米坐标需要标定像素间距原题没有给出该参数所以这一步只能做“形”的比较。严格的做法是用光斑的实际直径和像素数算出pixel_size_mm再把xpix乘上它。3.4 用mesh命令画三维强度分布并去白边二维剖面只能看到一条直径上的分布完整光斑要用整幅矩阵画三维图。原文档里的做法是[high, width, color] size(A); x 1:width; y 1:high-1; mesh(x, y, double(A(2:224,:,1))); grid on; xlabel(x); ylabel(y); zlabel(z); title(三维强度分布);这段代码已经处理了顶部白边A(2:224,:,1)从第 2 行开始取跳过了第 1 行的亮线。但它写死了 224如果图片尺寸变化会越界。我一般会写成A_gray double(A(:,:,1)); A_crop A_gray(2:end-1, :); % 去掉顶部和底部各一行 [Ny, Nx] size(A_crop); x 1:Nx; y 1:Ny; figure; mesh(x, y, A_crop); grid on; xlabel(x (pixel)); ylabel(y (pixel)); zlabel(强度); title(三维强度分布);这里2:end-1比2:224更通用。白边在 mesh 图里的表现是 z 轴一整边被拉高因为白边灰度接近 255而光斑中心灰度可能不到 200如果不剔除三维图会有一条边整体竖起来光斑本身的形状反而看不清楚。如果白边在侧边而不是顶部可以通过imshow(A(1:5,:))或size(A)确认再相应调整列方向的范围。注意mesh 的输入必须是 double 类型uint8 矩阵在三维坐标下可能不按数值直接渲染会出现颜色被当作索引的问题。4. 角谱法传播仿真不同z位置的光斑半径与强度变化4.1 为什么是角谱法而不是直接积分高斯光束在自由空间中传播理论上可以用菲涅尔衍射积分但离散化之后每个观察点都要做一次二维求和计算量很大。角谱法的思路是把初始光场分解成平面波乘上自由空间传递函数再逆变换回空间域。由于只涉及两次 FFT 和一次矩阵乘法扫描多个 z 值时速度很快非常适合画“传播 10m、20m、50m”这组图。近轴条件下传递函数可以写成H(fx,fy)exp(jz(kx^2ky^2)/(2*k0))其中 kx2pifxky2pify。这个公式是原始 M 文件的核心也是高斯光束从束腰传播到远场时振幅和束宽变化的数值基础。适用前提是传播角足够小、介质均匀、没有增益或损耗激光器谐振腔外检测正好满足这些条件。4.2 角谱传播代码与fftshift配合下面这段脚本可以直接替换原 M 文件中的传播部分解决中心偏移和 ifftshift 缺失的问题。clear; close all; N 500; % 采样点数100~500 L 10e-3; % 仿真区域宽度单位m dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; y x; [X, Y] meshgrid(x, y); lambda 0.568e-6; % 波长单位m k0 2*pi / lambda; wo 1e-3; % 束腰半径单位m zR k0 * wo^2 / 2; % 瑞利长度 fprintf(Rayleigh range: %.2f m\n, zR); Gau exp(-(X.^2 Y.^2) / wo^2); % 初始振幅分布 fx (-N/2 : N/2-1) / (N*dx); % 空间频率单位1/m [FX, FY] meshgrid(fx, fx); z input(Propagation distance z (m): ); H exp(1j * z / (2*k0) * ((2*pi*FX).^2 (2*pi*FY).^2)); FGau fftshift(fft2(Gau)); % 零频移到矩阵中心 Gau_pro ifft2(ifftshift(FGau .* H)); % 乘H后逆shift再ifft2 figure; subplot(1,2,1); mesh(x*1e3, y*1e3, abs(Gau)); title(Initial Gaussian Beam); xlabel(x (mm)); ylabel(y (mm)); zlabel(Amplitude); subplot(1,2,2); mesh(x*1e3, y*1e3, abs(Gau_pro)); title([z, num2str(z), m]); xlabel(x (mm)); ylabel(y (mm)); zlabel(Amplitude);代码的要点是频域坐标必须和fft2的输出顺序一致。fft2的零频在矩阵的(1,1)位置fftshift后零频在中心同时H是按零频在中心构造的所以两者可以直接相乘。乘完之后要再ifftshift把零频挪回(1,1)再送进ifft2。原始 M 文件最容易忽略这一步结果是光斑在空域发生循环平移峰值下降数值束宽也可能从理论上的十几毫米变成一两毫米。fx (-N/2 : N/2-1) / (N*dx)是空间频率轴单位为 1/m。频域间隔等于 1/(N*dx)1/L也就是仿真区域越宽频域越密角谱法的分辨率越高。H中z/(2*k0)这个系数决定了相位累积的快慢当 z 很大时H 在频域边缘振荡很快N 不够就会发生仿真发散。4.3 扫描z值并比较数值束宽与理论束宽要观察光斑随传播距离的变化不能只跑一次可以把 z 写成一个数组循环。理论束宽用wz wo*sqrt(1(z/zR)^2)。z_list [0.1, 1, 5, 10, 20, 50, 100]; w_num zeros(size(z_list)); w_theory zeros(size(z_list)); FGau fftshift(fft2(Gau)); for k 1:length(z_list) zk z_list(k); Hk exp(1j * zk / (2*k0) * ((2*pi*FX).^2 (2*pi*FY).^2)); Gk ifft2(ifftshift(FGau .* Hk)); Ik abs(Gk).^2; Ik Ik / max(Ik(:)); line_prof Ik(N/21, :); pos1 find(line_prof exp(-2), 1, first); pos2 find(line_prof exp(-2), 1, last); w_num(k) (x(pos2) - x(pos)) / 2; w_theory(k) wo * sqrt(1 (zk / zR)^2); end figure; semilogy(z_list, w_num*1e3, ro-, z_list, w_theory*1e3, b^-); legend(数值束宽, 理论束宽, Location, northwest); xlabel(z (m)); ylabel(束宽 (mm));line_prof取的是通过中心的水平强度剖面阈值exp(-2)对应强度降到峰值 1/e^2 的位置。find找到左右两个越过阈值的点相减再除以 2 就是光斑半径。如果pos1为空说明光斑已经扩展出仿真区域需要增大L或N。原文档运行结果中100m 处的数值束宽是 1.94mm理论值 18.107mm二者差一个数量级最可能的原因就是fftshift/ifftshift不配对修正后 100m 处会明显接近理论束宽。不同传播距离下的束宽变化趋势可以总结成下表供报告排版参考。z (m)z/z_R理论 w(z) (mm)主要特征0.10.018约 1.00接近束腰强度近高斯10.18约 1.02变化很小50.90约 1.35开始明显扩展101.81约 2.03远场区近似直线扩展509.04约 9.09峰值降低明显10018.08约 18.11扩展约18倍表中的理论值按wzwo*sqrt(1(z/zR)^2)计算z_R5.53m。数值仿真的峰值强度也会随 z 增大而下降但要注意abs(Gau_pro)是振幅画强度图时要取平方再归一化。5. 调试细节白边剔除、采样点与束宽计算的一致性5.1 白边不一定要手动数行数三维图画出来一边高先不要急着改数据。用下面几行定位白边位置imshow(A(1:5, :, 1)); % 查看顶部5行 figure; imagesc(double(A(:,:,1))); colorbar; % 完整显示强度范围白边在imagesc下会是一条接近 colorbar 顶部的亮色条比直接看 mesh 图更直观。确认白边只有顶部一行时用A(2:end-1, :)去掉首尾如果白边有十行就把前面范围改成A(11:end-10, :)。不要写死A(2:224,:,1)换一张图就容易越界。5.2 仿真发散时的三个参数调整高频相位因子在频域边缘振荡时会出现边缘条纹或能量不衰减这就是常见的“仿真发散”。按以下顺序调整提高采样点数N从 100 加到 500频域采样更密缩小仿真区域L太大时光斑只占很少像素相位变化集中在少数网格上建议L为光斑直径的 35 倍分段传播z 超过几十米时把 z 拆成若干段循环传播每段不超过瑞利长度可以显著减少混叠。每次调整后检查中心行剖面如果边缘仍有周期性起伏说明还没收敛。5.3 束宽定义统一才能对比数值和理论定义阈值位置实际用途振幅半径 w(z)振幅降到 1/e理论公式 wzw0*sqrt(1(z/zR)^2)强度半径强度降到 1/e^2CCD 光斑实测常用FWHM强度降到半高工程测量常见数值比 w(z) 小理论束宽公式里的 w(z) 对应振幅降到 1/e 的位置换算成强度就是 1/e^2。所以代码里数值束宽也应该按强度阈值exp(-2)找边缘才能和wzwo*sqrt(1(z/zR)^2)对齐。如果按半高 FWHM 找结果会差约 1.18 倍不要混用。把像素坐标换成空间坐标时先标定像素间距。CCD 图像里的一个像素对应多少毫米应通过光斑实际直径和像素数算出来而不是对矩阵下标直接取abs。否则二维图、三维图和传播仿真四张图之间的束宽永远对不上。本文还有配套的精品资源点击获取 高斯光束Matlab仿真：强度分布、图像处理与角谱传播
简介一份围绕高斯光束Matlab仿真的完整技术文档适合激光原理、光学工程等课程的学生与科研人员参考。文档从高斯光束数学模型出发推导强度分布公式并给出谐振腔中位置处的归一化强度分布仿真方法覆盖二维、三维强度分布图绘制。同时结合CCD采集的真实激光光斑图像通过imread读取强度数据与理论高斯曲线进行对比展示数学仿真与实际测量的差异。进一步利用角谱衍射法模拟高斯光束在自由空间传播过程中光强、光斑有效截面半径及等相位面的变化包含核心Matlab代码与参数输入示例方便读者复现。资源为单份docx文档体积432KB内容结构清晰既有公式推导也有程序实现可直接用于实验报告撰写或课程设计参考。已有461人学习下载适合需要快速完成高斯光束仿真作业或理解光场传播原理的读者。1. 高斯光束的Matlab仿真先画强度分布再谈传播用Matlab画高斯光束十行代码就能出一张三维图但把CCD拍到的光斑和理论曲线叠在一起时峰值位置、半宽、背景噪声都会对不上三维图还会有一边整体竖起来。问题往往不在公式而在矩阵读取方向、灰度通道和白边处理。下面这套流程来自激光原理课程里的高斯光束仿真作业解决两个具体问题一是用imread读取实际光斑图像生成二维、三维归一化强度分布并与理论高斯曲线对比二是用角谱传播模型仿真光束传播不同距离后光斑半径和峰值强度的变化。适合正在做激光实验、光电课程设计或者想搞懂Matlab图像处理与傅里叶光学仿真边界的读者。2. 高斯光束数学模型与二维强度矩阵的生成2.1 基模高斯光束的强度表达式与参数约定谐振腔中常见的基模场直角坐标下振幅分布写作E(x,y)E0 exp(-(x^2y^2)/w0^2)对应的强度分布是 |E|^2写成I(x,y)I0 exp(-2(x^2y^2)/w0^2)这里的 w0 是束腰半径在振幅定义里是振幅下降到中心值 1/e 的位置在强度定义里1/e^2 处对应同样半径。很多Matlab示例直接用 exp(-r^2/w0^2) 画“强度分布”严格说那画的是振幅场。如果只做归一化展示两者形状相近但和实测光斑做半宽对比时强度曲线的半高位置会在 sqrt(2) 倍处出现偏差。做这一题之前先把要画的是振幅还是强度定下来后面所有公式、阈值和束宽计算才不会乱。这份课程作业文档资料里的 M 文件把理论曲线写成a2exp(-x2.^2/10)形式上也是 exp 函数但它没有按严格的 1/e^2 强度定义写更像是一个拟合用高斯函数。作业里可以沿用但在程序注释里最好写明w_fitsqrt(10)而不是w0避免后续传播仿真和二维图对比时概念混淆。2.2 用meshgrid生成二维高斯强度场二维高斯分布最简单的写法是用 meshgrid 构建网格再把径向距离平方代入 exp。下面的代码生成 224×224 的归一化强度场和后面 CCD 图像尺寸保持一致。N 224; % 与CCD图像行数保持一致方便对比 L 10e-3; % 仿真区域边长单位m wo 1e-3; % 束腰半径单位m dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; % 以0为中心的横坐标 [X, Y] meshgrid(x, x); % 若按强度分布指数项为 -2*(X.^2Y.^2)/wo^2 Gau exp(-(X.^2 Y.^2) / wo^2); % 振幅场 Gau_norm Gau / max(Gau(:)); % 归一化 figure; surf(X*1e3, Y*1e3, Gau_norm, EdgeColor, none); xlabel(x (mm)); ylabel(y (mm)); zlabel(归一化强度);这段代码里N224不是必须的后面做 FFT 传播时也会用 500。关键是(-N/2 : N/2-1)这种写法让光斑中心落在网格中央而不是数组的(1,1)位置如果直接用1:N后续fft2和surf都要处理中心偏移麻烦得多。dxL/N决定空间分辨率L取 10mm 时单个像素相当于约 44.6µm这个量级和实验中 CCD 的像素尺寸接近适合做毫米级光斑仿真。绘图时乘 1e3 是为了把米换成毫米否则横轴刻度会显示成0.001这类小数值。2.3 参数表波长、束腰、瑞利长度和采样点传播仿真中需要提前确定的参数如下。参数符号本资源中的值作用波长lambda0.568 µm决定波数 k0 和瑞利长度束腰半径wo1 mm初始光斑尺寸二维图半宽来源瑞利长度z_R约 5530 mm近场/远场分界传播 z 的标尺采样点数N100500影响 FFT 精度太小会混叠仿真区域L10 mm限制频域分辨率为 1/L影响传播相位精度瑞利长度由zR k0*wo^2/2计算代入k02*pi/lambda后等于pi*wo^2/lambda。用上面的参数算出来约 5.53m所以文档中z100000mm即 100m 时已经远大于瑞利长度光斑半径接近线性扩展峰值强度明显下降。运行结果里峰值从 1 降到 0.21是这个参数下的正常物理趋势不是程序错误。2.4 为什么要做归一化CCD 读出的灰度值不是光功率绝对值单位是“数”和理论公式中的 I0 没有可比性。所以画二维强度分布图前一定要把数据除以自身最大值。更稳妥的做法是先减背景再归一化I_norm(I-min(I))/(max(I)-min(I))。原始光斑照片如果有暗电流或环境光底部不是 0直接除以最大值会把背景也纳入比较导致理论曲线和实验曲线的边缘不贴合。这一行处理在后面的 imread 实验数据中同样适用。3. 用imread读取CCD光斑二维/三维强度分布与白边处理3.1 读取图片前先看矩阵维度imread是 Matlab 读取位图最常用的命令但读到的是什么尺寸取决于图片本身是灰度图还是 RGB 图。A imread(D:\documents\作业\激光原理与应用\高斯.bmp); disp(size(A)); imshow(A); axis off; title(CCD采集的高斯光束光强分布);如果size(A)返回224 244说明是单通道灰度图如果返回224 244 3说明是 RGB 三通道后面所有取行、取列的操作都必须指定第三维下标否则A(:,122)拿到的是 244×3 的矩阵plot会画出三条线。原文档写的是 224×244 矩阵但代码里又出现了[high, width, color] size(A)说明实际文件很可能是 RGB 格式。常见做法是先用size确认然后取第一通道或者用rgb2gray转成灰度。imshow对uint8数据会自动按灰度范围映射直接看即可不要先用double转换再 imshow否则显示结果可能变黑。3.2 提取直径方向数据绘制二维强度分布原始程序的读取核心是两行A1 A(:, 122); % 取第122列长度224 x1 1:1:224; figure; plot(x1, A1, LineWidth, 1, Color, r); title(理论高斯曲线);这里要特别注意“中间一行”和“中间一列”的差别。题目文字说的是“取中间一行122行”但代码A(:,122)取的是第 122 列。如果矩阵是 224 行 × 244 列那么第 122 列对应的剖面长度是 224和第 122 行对应的剖面长度 244 并不相同。两个方向理论上都可以反映高斯剖面但横轴刻度必须和取的方向一致。更好的写法是if size(A,3) 3 A A(:,:,1); % 取单通道后续都是二维矩阵 end line_data double(A(122, :)); % 第122行长度244 line_data line_data / max(line_data); % 归一化到[0,1] Nx length(line_data); xpix (1:Nx) - (Nx1)/2; % 以图像中心为原点 figure; plot(xpix, line_data, r, LineWidth, 1);double转换是必须的因为imread返回的是uint8直接做除法时 Matlab 会按整数运算取整。归一化后纵轴可以统一到 01方便和理论曲线放在同一张图里。3.3 理论高斯曲线与实验曲线的叠加文档里的理论曲线是x2-100:1:100; a2exp(-x2.^2/10)。它的 x 范围是 [-100,100]和实验曲线的像素范围 [1,224] 不在同一坐标基准上所以原图把两条曲线分开画看起来各自都像高斯叠在一起却对不上。叠加时要做两件事把实验曲线中心移到 0再让理论曲线的横轴范围和实验曲线匹配。x_theory -100:0.1:100; I_theory exp(-x_theory.^2 / 10); figure; plot(xpix, line_data, r, LineWidth, 1); hold on; plot(x_theory, I_theory, b--, LineWidth, 1); legend(实验数据, 理论曲线); xlabel(相对中心像素); ylabel(归一化强度);如果发现实验峰值不在 x0 附近说明光斑中心没有落在图像中心可以用[~, idx] max(line_data)把峰值位置平移到 0 后再比较。像素坐标到毫米坐标需要标定像素间距原题没有给出该参数所以这一步只能做“形”的比较。严格的做法是用光斑的实际直径和像素数算出pixel_size_mm再把xpix乘上它。3.4 用mesh命令画三维强度分布并去白边二维剖面只能看到一条直径上的分布完整光斑要用整幅矩阵画三维图。原文档里的做法是[high, width, color] size(A); x 1:width; y 1:high-1; mesh(x, y, double(A(2:224,:,1))); grid on; xlabel(x); ylabel(y); zlabel(z); title(三维强度分布);这段代码已经处理了顶部白边A(2:224,:,1)从第 2 行开始取跳过了第 1 行的亮线。但它写死了 224如果图片尺寸变化会越界。我一般会写成A_gray double(A(:,:,1)); A_crop A_gray(2:end-1, :); % 去掉顶部和底部各一行 [Ny, Nx] size(A_crop); x 1:Nx; y 1:Ny; figure; mesh(x, y, A_crop); grid on; xlabel(x (pixel)); ylabel(y (pixel)); zlabel(强度); title(三维强度分布);这里2:end-1比2:224更通用。白边在 mesh 图里的表现是 z 轴一整边被拉高因为白边灰度接近 255而光斑中心灰度可能不到 200如果不剔除三维图会有一条边整体竖起来光斑本身的形状反而看不清楚。如果白边在侧边而不是顶部可以通过imshow(A(1:5,:))或size(A)确认再相应调整列方向的范围。注意mesh 的输入必须是 double 类型uint8 矩阵在三维坐标下可能不按数值直接渲染会出现颜色被当作索引的问题。4. 角谱法传播仿真不同z位置的光斑半径与强度变化4.1 为什么是角谱法而不是直接积分高斯光束在自由空间中传播理论上可以用菲涅尔衍射积分但离散化之后每个观察点都要做一次二维求和计算量很大。角谱法的思路是把初始光场分解成平面波乘上自由空间传递函数再逆变换回空间域。由于只涉及两次 FFT 和一次矩阵乘法扫描多个 z 值时速度很快非常适合画“传播 10m、20m、50m”这组图。近轴条件下传递函数可以写成H(fx,fy)exp(jz(kx^2ky^2)/(2*k0))其中 kx2pifxky2pify。这个公式是原始 M 文件的核心也是高斯光束从束腰传播到远场时振幅和束宽变化的数值基础。适用前提是传播角足够小、介质均匀、没有增益或损耗激光器谐振腔外检测正好满足这些条件。4.2 角谱传播代码与fftshift配合下面这段脚本可以直接替换原 M 文件中的传播部分解决中心偏移和 ifftshift 缺失的问题。clear; close all; N 500; % 采样点数100~500 L 10e-3; % 仿真区域宽度单位m dx L / N; % 空间采样间隔 x (-N/2 : N/2-1) * dx; y x; [X, Y] meshgrid(x, y); lambda 0.568e-6; % 波长单位m k0 2*pi / lambda; wo 1e-3; % 束腰半径单位m zR k0 * wo^2 / 2; % 瑞利长度 fprintf(Rayleigh range: %.2f m\n, zR); Gau exp(-(X.^2 Y.^2) / wo^2); % 初始振幅分布 fx (-N/2 : N/2-1) / (N*dx); % 空间频率单位1/m [FX, FY] meshgrid(fx, fx); z input(Propagation distance z (m): ); H exp(1j * z / (2*k0) * ((2*pi*FX).^2 (2*pi*FY).^2)); FGau fftshift(fft2(Gau)); % 零频移到矩阵中心 Gau_pro ifft2(ifftshift(FGau .* H)); % 乘H后逆shift再ifft2 figure; subplot(1,2,1); mesh(x*1e3, y*1e3, abs(Gau)); title(Initial Gaussian Beam); xlabel(x (mm)); ylabel(y (mm)); zlabel(Amplitude); subplot(1,2,2); mesh(x*1e3, y*1e3, abs(Gau_pro)); title([z, num2str(z), m]); xlabel(x (mm)); ylabel(y (mm)); zlabel(Amplitude);代码的要点是频域坐标必须和fft2的输出顺序一致。fft2的零频在矩阵的(1,1)位置fftshift后零频在中心同时H是按零频在中心构造的所以两者可以直接相乘。乘完之后要再ifftshift把零频挪回(1,1)再送进ifft2。原始 M 文件最容易忽略这一步结果是光斑在空域发生循环平移峰值下降数值束宽也可能从理论上的十几毫米变成一两毫米。fx (-N/2 : N/2-1) / (N*dx)是空间频率轴单位为 1/m。频域间隔等于 1/(N*dx)1/L也就是仿真区域越宽频域越密角谱法的分辨率越高。H中z/(2*k0)这个系数决定了相位累积的快慢当 z 很大时H 在频域边缘振荡很快N 不够就会发生仿真发散。4.3 扫描z值并比较数值束宽与理论束宽要观察光斑随传播距离的变化不能只跑一次可以把 z 写成一个数组循环。理论束宽用wz wo*sqrt(1(z/zR)^2)。z_list [0.1, 1, 5, 10, 20, 50, 100]; w_num zeros(size(z_list)); w_theory zeros(size(z_list)); FGau fftshift(fft2(Gau)); for k 1:length(z_list) zk z_list(k); Hk exp(1j * zk / (2*k0) * ((2*pi*FX).^2 (2*pi*FY).^2)); Gk ifft2(ifftshift(FGau .* Hk)); Ik abs(Gk).^2; Ik Ik / max(Ik(:)); line_prof Ik(N/21, :); pos1 find(line_prof exp(-2), 1, first); pos2 find(line_prof exp(-2), 1, last); w_num(k) (x(pos2) - x(pos)) / 2; w_theory(k) wo * sqrt(1 (zk / zR)^2); end figure; semilogy(z_list, w_num*1e3, ro-, z_list, w_theory*1e3, b^-); legend(数值束宽, 理论束宽, Location, northwest); xlabel(z (m)); ylabel(束宽 (mm));line_prof取的是通过中心的水平强度剖面阈值exp(-2)对应强度降到峰值 1/e^2 的位置。find找到左右两个越过阈值的点相减再除以 2 就是光斑半径。如果pos1为空说明光斑已经扩展出仿真区域需要增大L或N。原文档运行结果中100m 处的数值束宽是 1.94mm理论值 18.107mm二者差一个数量级最可能的原因就是fftshift/ifftshift不配对修正后 100m 处会明显接近理论束宽。不同传播距离下的束宽变化趋势可以总结成下表供报告排版参考。z (m)z/z_R理论 w(z) (mm)主要特征0.10.018约 1.00接近束腰强度近高斯10.18约 1.02变化很小50.90约 1.35开始明显扩展101.81约 2.03远场区近似直线扩展509.04约 9.09峰值降低明显10018.08约 18.11扩展约18倍表中的理论值按wzwo*sqrt(1(z/zR)^2)计算z_R5.53m。数值仿真的峰值强度也会随 z 增大而下降但要注意abs(Gau_pro)是振幅画强度图时要取平方再归一化。5. 调试细节白边剔除、采样点与束宽计算的一致性5.1 白边不一定要手动数行数三维图画出来一边高先不要急着改数据。用下面几行定位白边位置imshow(A(1:5, :, 1)); % 查看顶部5行 figure; imagesc(double(A(:,:,1))); colorbar; % 完整显示强度范围白边在imagesc下会是一条接近 colorbar 顶部的亮色条比直接看 mesh 图更直观。确认白边只有顶部一行时用A(2:end-1, :)去掉首尾如果白边有十行就把前面范围改成A(11:end-10, :)。不要写死A(2:224,:,1)换一张图就容易越界。5.2 仿真发散时的三个参数调整高频相位因子在频域边缘振荡时会出现边缘条纹或能量不衰减这就是常见的“仿真发散”。按以下顺序调整提高采样点数N从 100 加到 500频域采样更密缩小仿真区域L太大时光斑只占很少像素相位变化集中在少数网格上建议L为光斑直径的 35 倍分段传播z 超过几十米时把 z 拆成若干段循环传播每段不超过瑞利长度可以显著减少混叠。每次调整后检查中心行剖面如果边缘仍有周期性起伏说明还没收敛。5.3 束宽定义统一才能对比数值和理论定义阈值位置实际用途振幅半径 w(z)振幅降到 1/e理论公式 wzw0*sqrt(1(z/zR)^2)强度半径强度降到 1/e^2CCD 光斑实测常用FWHM强度降到半高工程测量常见数值比 w(z) 小理论束宽公式里的 w(z) 对应振幅降到 1/e 的位置换算成强度就是 1/e^2。所以代码里数值束宽也应该按强度阈值exp(-2)找边缘才能和wzwo*sqrt(1(z/zR)^2)对齐。如果按半高 FWHM 找结果会差约 1.18 倍不要混用。把像素坐标换成空间坐标时先标定像素间距。CCD 图像里的一个像素对应多少毫米应通过光斑实际直径和像素数算出来而不是对矩阵下标直接取abs。否则二维图、三维图和传播仿真四张图之间的束宽永远对不上。本文还有配套的精品资源点击获取