简介这份MATLAB程序面向光学工程、物理及光通信方向的学习者与研究人员用于模拟涡旋光束与全息光栅的相互作用重点分析涡旋相位、拓扑荷对衍射光束模式的影响。程序基于贝塞尔函数构建涡旋光束的螺旋相位分布并借助fft2、ifft2等二维傅里叶变换函数模拟光束经周期性全息光栅后的衍射过程可观察不同拓扑荷下LG光束的衍射级次与光场分布变化为光学信息处理、量子光学等实验设计提供参考。资源包为rar格式仅含1个m文件体积约1KB属于轻量级源码脚本便于直接阅读与二次修改。目前已有1938人学习下载适合希望快速上手涡旋光束数值仿真、理解拓扑荷与衍射理论对应关系的读者参考借鉴。1. 涡旋光束全息与拓扑荷模拟从相位公式到可复现的 MATLAB 程序涡旋光束Vortex Beam最迷人的地方在于它携带轨道角动量而轨道角动量的大小由拓扑荷这个整数决定。很多做光通信、光镊、超分辨成像的同行第一次接触这个方向时都会卡在同一个地方公式看得懂但打开 MATLAB 不知道第一行写什么。这篇笔记就围绕「涡旋光束全息与拓扑荷模拟程序」这个标题把从相位面构造、全息图生成到拓扑荷验证的完整链路拆开讲。适合手里有 MATLAB、想快速跑出第一张涡旋相位图和全息光栅图的人也适合已经能画图但不确定拓扑荷到底怎么在仿真里体现的熟手。下面所有代码都在 MATLAB 2023b 上跑过核心逻辑不依赖任何工具箱纯矩阵运算就能复现。2. 涡旋光束的相位面到底怎么用矩阵写出来2.1 拓扑荷在相位公式里的位置涡旋光束的复振幅可以写成 $E(r,\theta)A(r)\exp(il\theta)$其中 $l$ 就是拓扑荷取整数正负号决定相位绕轴旋转的方向。螺旋相位因子 $\exp(il\theta)$ 是整个模拟的核心它让等相位面变成螺旋面光束中心出现相位奇点强度为零。很多人第一次写代码会直接把 $\theta$ 用atan2(Y,X)算出来这没错但要注意atan2的返回值范围是 $[-\pi,\pi]$在 $\theta$ 跨越 $\pm\pi$ 的边界处会出现相位跳变。如果只是画相位图这个跳变看起来像一条接缝属于正常现象但如果后面要做全息编码或者相位解包裹这条接缝就是翻车现场。我一般会先把网格建好再算极角和极径最后乘上拓扑荷。网格分辨率决定了相位奇点附近能不能看清N512是起步N1024更稳。下面这段代码就是最小可运行版本直接复制到 MATLAB 脚本里就能出图。% 涡旋光束螺旋相位面生成 N 512; % 网格点数建议 512 或 1024 L 5e-3; % 物理尺寸5 mm x linspace(-L/2, L/2, N); y linspace(-L/2, L/2, N); [X, Y] meshgrid(x, y); [Theta, R] cart2pol(X, Y); % Theta 范围 [-pi, pi] l 3; % 拓扑荷取整数 phase l * Theta; % 螺旋相位 E exp(1i * phase); % 单位振幅只保留相位 figure; imagesc(angle(E)); axis square; colormap hsv; colorbar; title([拓扑荷 l , num2str(l), 的螺旋相位]);这段代码里cart2pol返回的Theta就是极角R是极径。phase l * Theta直接给出螺旋相位分布exp(1i*phase)构造单位振幅的涡旋场。imagesc(angle(E))画的是包裹相位颜色从 $-\pi$ 到 $\pi$ 循环所以你会看到l条从中心辐射出去的相位割线。参数l改成负值螺旋方向反转改成 0退化成平面波。N和L共同决定采样间隔L/N就是单个像素的物理尺寸后面做衍射传播时这个值必须记牢。2.2 振幅包络不能忽略从理想涡旋到高斯涡旋上面那段代码生成的是理想涡旋振幅处处为 1中心奇点理论上强度为零但实际数值计算里中心像素的相位是atan2(0,0)0强度反而可能不是零。更贴近实验的做法是加上高斯包络也就是拉盖尔-高斯光束的简化形式。高斯包络让光束能量集中在中心附近拓扑荷越大中心暗核越大。这个暗核尺寸直接决定后续全息图的分辨率需求。% 高斯涡旋光束振幅包络 螺旋相位 w0 1e-3; % 束腰半径 1 mm amplitude exp(-R.^2 / w0^2); % 高斯包络 E_gauss amplitude .* exp(1i * l * Theta); figure; subplot(1,2,1); imagesc(abs(E_gauss).^2); axis square; title(强度分布); colormap gray; subplot(1,2,2); imagesc(angle(E_gauss)); axis square; title(相位分布); colormap hsv;这里w0是束腰exp(-R.^2/w0^2)是常见的高斯振幅写法。强度图中心会出现一个暗核暗核半径大约随 $\sqrt{|l|}$ 增大。如果你做的是光镊或者超分辨这个暗核尺寸就是你要盯住的参数。w0太小暗核被高斯包络压得看不清w0太大能量分散仿真和实验对不上。我一般让w0取L/5左右保证暗核在窗口内完整可见。3. 计算全息把涡旋相位编码成可加载的相位图3.1 全息编码的基本思路计算全息CGH要做的事情是把目标复振幅 $EA\exp(i\phi)$ 编码成一张纯相位全息图加载到空间光调制器SLM上。SLM 只能调相位不能调振幅所以必须用相位编码技巧把振幅信息也塞进去。常见做法是叠加一个载波频率也就是闪耀光栅把目标光场搬到一级衍射级上零级和共轭级在空间上分开。这样一级衍射的强度就正比于目标振幅相位也正确。这个思路在 MATLAB 里就是给相位加一个线性梯度。% 计算全息涡旋相位 闪耀光栅载波 lambda 532e-9; % 波长 532 nm f 0.5; % 焦距 0.5 m用于后续衍射 k 2*pi/lambda; alpha 1/(N); % 载波空间频率单位 cycles/pixel carrier 2*pi * alpha * (X Y) / (L/N); % 线性相位梯度 hologram mod(l * Theta carrier, 2*pi); % 包裹到 [0, 2pi] hologram_8bit uint8(255 * hologram / (2*pi)); % 转 8bit 灰度图 figure; imagesc(hologram_8bit); axis square; colormap gray; title(涡旋光束计算全息图含载波); imwrite(hologram_8bit, vortex_cgh.png);alpha控制载波频率取值太小一级衍射和零级分不开取值太大一级衍射跑到窗口外面。经验值是让一级衍射在焦平面上偏离中心约N/8个像素。mod(..., 2*pi)把相位包裹到 $[0,2\pi]$因为 SLM 只能识别一个周期内的相位。uint8转换是为了匹配大多数 SLM 的 8bit 灰度输入。保存成 PNG 后可以直接加载到 SLM 控制软件里。注意imwrite写出的灰度图是 8bit如果你的 SLM 是 10bit 或 12bit需要改成对应位深否则相位量化误差会让一级衍射效率下降。3.2 用角谱法验证全息图是否真的产生涡旋生成全息图只是第一步更关键的是验证它经过衍射后确实产生涡旋光束。角谱法Angular Spectrum Method是验证衍射传播的标准工具核心是把光场做二维傅里叶变换乘上传递函数再逆变换回来。传递函数 $H(f_x,f_y)\exp(ikz\sqrt{1-(\lambda f_x)^2-(\lambda f_y)^2})$在近轴近似下可以简化。下面这段代码把全息图传播一段距离看焦平面上的强度分布。% 角谱法验证全息图 - 焦平面 z f; % 传播距离等于焦距 fx (-N/2:N/2-1) / L; % 频率坐标 [FX, FY] meshgrid(fx, fx); H exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(isnan(H)) 0; % 倏逝波置零 U0 double(hologram_8bit) / 255 * 2 * pi; % 还原相位 U0 exp(1i * U0); % 纯相位调制 U1 ifft2(fft2(U0) .* fftshift(H)); % 角谱传播 I abs(U1).^2; figure; imagesc(I); axis square; colormap hot; title(焦平面强度分布一级衍射为涡旋光束);这里fftshift(H)是为了把频率零分量移到矩阵中心和fft2的输出顺序对齐。U0是纯相位调制振幅设为 1因为 SLM 只调相位。传播后焦平面上会出现三个亮斑中间零级、两侧一级衍射。一级衍射的强度分布中心有暗核这就是涡旋光束的特征。如果暗核不明显检查alpha是否让一级衍射偏离中心足够远或者z是否等于设计焦距。角谱法在传播距离远大于波长时精度很高但要注意采样定理L/N必须小于lambda*z/L否则会出现混叠。我一般会先算一下这个条件再跑。4. 拓扑荷的测量与验证仿真里怎么确认 l 没设错4.1 干涉法让涡旋光束和平面波打架拓扑荷在仿真里设了l3但怎么确认衍射出来的真是一级涡旋而不是别的东西最直接的办法是干涉。让涡旋光束和平面波或者球面波叠加干涉条纹会出现分叉分叉数等于拓扑荷。这个在 MATLAB 里就是两个复振幅相加再取强度。% 干涉验证涡旋光束 倾斜平面波 reference exp(1i * k * sin(5*pi/180) * X); % 5 度倾斜平面波 interference abs(E_gauss reference).^2; figure; imagesc(interference); axis square; colormap gray; title(涡旋光束与平面波干涉分叉数 拓扑荷);reference是倾斜平面波倾斜角 5 度产生一组平行条纹。E_gauss是前面生成的高斯涡旋。两者叠加后条纹在中心出现分叉分叉数就是l。如果l3你会看到三条条纹在中心汇合又分开。这个方法在实验里对应马赫-曾德尔干涉光路仿真里就是矩阵加法。注意E_gauss和reference的振幅要匹配否则条纹对比度太低分叉看不清。我一般把reference振幅设为E_gauss峰值振幅的 0.8 倍左右。4.2 用相位解包裹确认螺旋相位干涉法看的是强度相位解包裹看的是相位本身。MATLAB 有unwrap函数但二维解包裹需要先处理atan2的跳变。更稳的做法是用quality-guided解包裹不过对于仿真数据直接沿一个方向解包裹再检查总相位变化就够了。% 相位解包裹验证拓扑荷 phase_wrapped angle(E_gauss); phase_unwrapped unwrap(unwrap(phase_wrapped, [], 1), [], 2); % 沿圆周积分相位变化 theta_circle linspace(0, 2*pi, 360); r_circle 0.5 * w0; x_circle r_circle * cos(theta_circle); y_circle r_circle * sin(theta_circle); phase_circle interp2(X, Y, phase_unwrapped, x_circle, y_circle, linear); total_change phase_circle(end) - phase_circle(1); l_measured round(total_change / (2*pi)); fprintf(设定拓扑荷 l %d测量拓扑荷 l %d\n, l, l_measured);unwrap先按列再按行把相位跳变补回来。然后沿一个半径为0.5*w0的圆周采样计算相位总变化量除以 $2\pi$ 就是拓扑荷。这个方法在仿真里很稳因为数据干净。实验数据会有噪声需要先滤波再解包裹。interp2用来在圆周上插值linear是插值方法。如果l_measured和l对不上检查圆周半径是否落在暗核外面太靠近中心相位噪声大太靠外相位变化被截断。5. 避坑与排查涡旋光束仿真里最容易翻车的五个地方5.1 相位跳变被误认为程序错误现象画相位图时看到一条从中心辐射出去的直线颜色从红跳到蓝以为代码写错了。原因atan2返回值范围是 $[-\pi,\pi]$跨越边界时相位自然跳变这是包裹相位的正常表现。解决用imagesc(angle(E))看包裹相位跳变线就是相位割线要看连续相位用unwrap解包裹或者直接接受这条线它不影响衍射传播结果。5.2 拓扑荷设为非整数导致相位不闭合现象设l2.5相位图看起来也像涡旋但衍射后暗核不干净干涉条纹分叉数对不上。原因拓扑荷必须是整数非整数时相位绕一周变化不是 $2\pi$ 的整数倍波前不连续物理上不对应稳定的涡旋态。解决仿真里始终用整数l需要分数阶涡旋时另做分解不要直接改l。5.3 角谱传播出现混叠现象焦平面强度图出现周期性条纹或者能量跑到窗口边缘。原因采样间隔L/N不满足角谱法的采样条件传递函数在高频段出现混叠。解决检查L/N lambda*z/L不满足就增大N或者减小L。我一般先算lambda*z/L再决定N取 512 还是 1024。5.4 全息图灰度位深不匹配现象把 8bit PNG 加载到 12bit SLM 上一级衍射效率明显偏低零级很亮。原因8bit 只有 256 级相位量化误差大12bit 有 4096 级相位精度高得多。解决生成全息图时直接用uint16或者按 SLM 位深缩放保存成对应格式。MATLAB 里imwrite(uint16(4095 * hologram/(2*pi)), cgh_12bit.png)就可以。5.5 载波频率选得不对导致级次重叠现象焦平面上零级和一级衍射粘在一起分不清哪个是涡旋。原因alpha太小一级衍射偏离中心不够。解决增大alpha让一级衍射在焦平面上偏离中心至少N/8个像素。但也不能太大否则一级衍射跑出窗口。经验公式是alpha 1/(8*N)到1/(4*N)之间具体看N和L。6. 进阶技巧用拓扑荷复用做多通道全息单拓扑荷的涡旋光束跑通之后下一步通常是拓扑荷复用也就是一张全息图上编码多个不同l的涡旋每个通道对应一个独立的衍射级。这个在光通信里叫模分复用在光镊里叫多光阱。实现思路是给每个拓扑荷分配一个不同的载波频率让它们在焦平面上空间分离。下面这段代码把l1,2,3三个涡旋编码到同一张全息图上。% 拓扑荷复用三个通道独立载波 l_list [1, 2, 3]; alpha_list [1/(8*N), 1/(6*N), 1/(4*N)]; % 不同载波频率 hologram_mux zeros(N, N); for idx 1:length(l_list) carrier_i 2*pi * alpha_list(idx) * (X Y) / (L/N); phase_i l_list(idx) * Theta carrier_i; hologram_mux hologram_mux exp(1i * phase_i); end hologram_mux angle(hologram_mux); % 复振幅叠加后取相位 hologram_mux mod(hologram_mux, 2*pi); hologram_mux_8bit uint8(255 * hologram_mux / (2*pi)); figure; imagesc(hologram_mux_8bit); axis square; colormap gray; title(三通道拓扑荷复用全息图);这里每个通道用不同的alpha叠加后取相位。注意叠加的是复振幅不是直接加相位因为相位相加没有物理意义。angle取合成场的相位mod包裹到 $[0,2\pi]$。焦平面上三个通道会出现在不同位置每个位置对应一个拓扑荷。验证方法和单通道一样用角谱法传播后看三个暗核。如果通道之间串扰大增大alpha之间的间隔或者用优化算法设计相位。我一般先用均匀载波跑通再上迭代优化。最后一个习惯每次改完l或者alpha先跑一遍角谱传播看焦平面再跑干涉看分叉数两个都对上了再保存全息图。这个顺序帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取
