简介：这份MATLAB程序面向光学、物理与光学工程方向的学习者和研究者，用于模拟涡旋光束与全息光栅的相互作用，帮助理解螺旋波前、涡旋相位及拓扑荷对光束传播特性的影响。程序以贝塞尔函数构建涡旋光束相位分布，结合傅里叶光学方法模拟光束经全息光栅后的衍射过程，并通过fft2、ifft2等函数完成二维傅里叶变换与逆变换，直观展示不同衍射级次的光束分布。调整拓扑荷参数还可观察LG光束等不同模式的旋涡性质变化，为光学通信、信息处理与量子光学等场景提供实验设计参考。资源包为rar压缩格式，仅含1个m源码文件，体积约1KB，结构精简，便于直接阅读与二次修改。目前已有1938人学习下载，适合希望快速上手涡旋光束全息模拟、掌握拓扑荷调控思路的读者参考使用。

1. 涡旋光束全息与拓扑荷模拟：一份能直接跑通的 MATLAB 程序包

做光场调控的同行大概率都遇到过这种场景：想验证一个螺旋相位板产生的涡旋光到底是不是目标拓扑荷，实验室里没有 SLM 可以实时调，或者想先看看不同拓扑荷下相位分布和强度环的对应关系，结果卡在“从哪写起”上。这份涡旋光束全息与拓扑荷模拟程序就是冲着这个痛点来的——它把涡旋相位生成、全息图计算、拓扑荷扫描、光场传播这几件事用 MATLAB 串成了一条完整链路，输入拓扑荷和光束参数就能出相位图、全息图和衍射光斑。适合光学工程、光通信、超分辨成像方向的研究生和工程师，尤其是需要快速出图验证想法、又不想从零推公式的人。程序本身不依赖额外工具箱，纯 MATLAB 基础函数就能跑，对版本要求不高，2023b 到 2026b 实测都能用。

2. 涡旋相位与拓扑荷的数学建模：从螺旋相位板到全息图

2.1 涡旋光束的相位结构为什么用 exp(i·l·θ) 描述

涡旋光束的核心特征是相位波前绕光轴旋转，每绕一圈相位变化 2πl，l 就是拓扑荷。数学上写成 exp(i·l·θ)，其中 θ 是方位角，l 取整数时相位连续，取非整数时会出现相位突变——这个突变在模拟里会直接反映成衍射图样的不对称。程序里用atan2(Y, X)算 θ，比atan(Y./X)稳，因为后者在第二、三象限会翻车，导致相位图出现假的分支切割。我一般会在生成相位前先把网格中心对齐到矩阵的中间像素，否则涡旋核会偏，拓扑荷越大偏得越明显。

% 生成极坐标网格，中心对齐到矩阵中心 N = 512; % 网格尺寸，建议 256/512/1024 L = 5e-3; % 物理尺寸 5mm 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 为拓扑荷 l = 3; vortex_phase = exp(1i * l * theta); % 高斯振幅包络，w0 为束腰 w0 = 1e-3; amplitude = exp(-(r.^2) / w0^2); field = amplitude .* vortex_phase;

这段代码里cart2pol返回的 theta 是 [-π, π]，直接乘 l 没问题，因为 exp 对 2π 周期不敏感。但如果你后面要做相位解包裹或者算相位梯度，就得注意 theta 的跳变位置。w0控制光束尺寸，一般取物理尺寸的 1/5 到 1/3，太小会截断，太大会浪费网格。N取 512 是速度和精度的折中，1024 能看清高阶拓扑荷的细节，但内存翻四倍。

2.2 计算全息图：把相位编码成振幅或相位光栅

有了涡旋相位，下一步是把它变成 SLM 能加载的全息图。常见做法是加一个载波频率，把相位信息搬到空间频率域，形成叉形光栅。程序里用cos(l*theta + 2*pi*fx*X)生成全息图，fx 是载波频率，一般取 1/(10*像素尺寸) 左右，太小衍射级分不开，太大超出 SLM 带宽。

% 计算全息图，载波频率 fx lambda = 632.8e-9; % 波长 pixel_size = L / N; % 像素尺寸 fx = 1 / (20 * pixel_size); % 载波频率，约 20 像素一个周期 hologram = 0.5 + 0.5 * cos(l * theta + 2*pi*fx*X); % 显示 figure; subplot(1,3,1); imagesc(angle(vortex_phase)); axis image; colormap gray; title('涡旋相位'); subplot(1,3,2); imagesc(hologram); axis image; colormap gray; title('全息图'); subplot(1,3,3); imagesc(abs(field).^2); axis image; colormap hot; title('强度');

hologram的值在 0 到 1 之间，直接对应 SLM 的灰度。载波频率fx决定了衍射级的位置，重建时用傅里叶变换取一级衍射就行。这里有个参数坑：fx如果和l的符号不匹配，一级衍射会跑到另一边，重建时取错区域就得到共轭相位。我一般会在代码里加个判断，根据 l 的正负自动翻转 fx 符号。

2.3 拓扑荷扫描：批量生成不同 l 的相位与全息图

单次模拟只能看一个 l，实际分析往往要扫 l = -5 到 5，看衍射环半径和拓扑荷的关系。程序里用 for 循环批量生成，存成 cell 数组或者三维矩阵。注意每次循环要重新算 theta 和 r 吗？不用，网格不变，只变 l 就行。

l_list = -5:5; num_l = length(l_list); phase_stack = zeros(N, N, num_l); hologram_stack = zeros(N, N, num_l); for idx = 1:num_l l = l_list(idx); phase_stack(:,:,idx) = angle(exp(1i * l * theta)); hologram_stack(:,:,idx) = 0.5 + 0.5 * cos(l * theta + 2*pi*fx*X); end % 保存为 mat 文件，方便后续加载 save('vortex_stack.mat', 'phase_stack', 'hologram_stack', 'l_list', 'x', 'y');

phase_stack存的是相位角，范围 [-π, π]，用angle提取。hologram_stack存全息图灰度。保存成 mat 文件后，下次可以直接 load，不用重算。如果 l 范围更大，比如 -20 到 20，建议只存全息图，相位可以现算，省空间。这里l_list用整数，非整数拓扑荷也能算，但衍射图样会分裂，程序里没做特殊处理，需要自己加判断。

3. 光场传播与衍射计算：从全息图到远场光斑

3.1 角谱传播法：为什么比菲涅尔更稳

模拟涡旋光传播，常见方法有菲涅尔衍射积分和角谱传播。菲涅尔在近场精度可以，但远场或者大数值孔径下会失真。角谱法直接在频域乘传递函数，近场远场都适用，程序里默认用角谱。传递函数写成 exp(i·kz·sqrt(1 - (λfx)^2 - (λfy)^2))，其中 fx、fy 是频率坐标。

% 角谱传播 k = 2*pi / lambda; z = 0.5; % 传播距离 0.5m fx_vec = (-N/2:N/2-1) / L; [FX, FY] = meshgrid(fx_vec, fx_vec); H = exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); H(abs(lambda*FX) > 1 | abs(lambda*FY) > 1) = 0; % 倏逝波截断 % 对全息图传播 field_holo = hologram_stack(:,:,6); % 取 l=0 的全息图 field_far = ifft2(fft2(field_hlo) .* H); intensity_far = abs(field_far).^2;

H是传递函数，sqrt里出现负数时对应倏逝波，程序里直接置零，避免复数溢出。z取 0.5m 是典型远场距离，实际按实验光路改。field_hlo取的是全息图，不是涡旋场本身，因为 SLM 加载的是全息图，传播后一级衍射才是涡旋光。这里有个细节：全息图是实振幅，传播后要取一级衍射，程序里用fft2后取中心偏移的区域，偏移量由载波频率决定。

3.2 一级衍射提取：叉形光栅的衍射级分离

全息图传播后会出现多个衍射级，零级在中心，一级在 ±fx 位置。提取一级衍射就是取对应区域的复振幅。程序里用circshift把一级移到中心，再截取。

% 提取一级衍射 shift_pix = round(fx * L); % 载波对应的像素偏移 field_1st = circshift(field_far, [0, -shift_pix]); field_1st = field_1st(round(N/4):round(3*N/4), round(N/4):round(3*N/4)); % 看相位，应该是涡旋 figure; imagesc(angle(field_1st)); axis image; colormap gray; title('一级衍射相位');

shift_pix是载波频率对应的像素数，circshift把一级移到中心。截取范围取中间一半，避免边缘噪声。如果shift_pix算错，一级没对齐，相位图就是乱的。我一般会先跑一个 l=0 的看零级和一级位置，确认shift_pix再批量跑。

3.3 拓扑荷与衍射环半径的定量关系

涡旋光远场强度是环状，环半径和拓扑荷近似成正比：r ≈ l·λ·z/(2π·w0)。程序里可以自动测环半径，用regionprops或者找强度峰值。

% 测环半径 intensity_1st = abs(field_1st).^2; [~, idx] = max(intensity_1st(:)); [cy, cx] = ind2sub(size(intensity_1st), idx); r_ring = sqrt((cx - size(intensity_1st,2)/2)^2 + (cy - size(intensity_1st,1)/2)^2) * pixel_size; % 理论值 r_theory = abs(l) * lambda * z / (2*pi*w0); fprintf('l=%d, 实测半径=%.4f mm, 理论=%.4f mm\n', l, r_ring*1e3, r_theory*1e3);

regionprops更稳，但这里用峰值法够用。r_ring是像素转物理尺寸。理论公式里w0是束腰，z是传播距离。实测和理论差 10% 以内算正常，差太多检查w0和z是否对。程序里对每个 l 都跑一遍，最后画 r-l 曲线，线性度好的话说明模拟靠谱。

4. 避坑与排查：涡旋光模拟里最容易翻车的五个点

4.1 相位图出现径向条纹，不是涡旋

现象：生成的相位图不是螺旋，而是放射状条纹。原因：theta算错，用了atan(Y./X)而不是atan2(Y,X)，导致象限跳变。解决：统一用cart2pol或atan2，并在生成后imagesc(angle(field))确认螺旋方向。

4.2 全息图衍射级分不开，一级和零级重叠

现象：传播后光斑糊成一团，看不到分离的衍射级。原因：载波频率fx太小，一级衍射角小于衍射极限。解决：增大fx，一般取 1/(10~20 像素尺寸)，同时确认 SLM 像素尺寸和波长匹配。

4.3 拓扑荷扫描时内存爆掉

现象：l 从 -50 扫到 50，MATLAB 报内存不足。原因：三维矩阵N×N×100在 N=1024 时占 800MB 以上。解决：只存全息图不存相位，或者每算一个 l 就存成单独文件，用matfile追加。

4.4 传播后相位不是涡旋，而是共轭

现象：一级衍射相位螺旋方向和预期相反。原因：载波频率符号和 l 符号不匹配，取了一级衍射的共轭级。解决：根据 l 正负翻转fx符号，或者在提取一级时取另一侧。

4.5 环半径和理论对不上，差一个数量级

现象：实测环半径比理论大 10 倍。原因：pixel_size算错，用了L/N但 L 单位是 mm 不是 m。解决：统一用米，L=5e-3，pixel_size = L/N，输出时再转 mm。

5. 进阶技巧：用 GPU 加速和相位解包裹验证拓扑荷

5.1 GPU 加速批量传播

拓扑荷扫描时，每个 l 都要做一次 FFT 和角谱传播，CPU 跑 100 个 l 要几分钟。有 GPU 的话，把fft2换成gpuArray的fft2，速度能快 5 到 10 倍。注意显存要够，N=1024 时单精度复数矩阵约 16MB，100 个 l 就是 1.6GB，显存小的分批跑。

% GPU 加速版 if gpuDeviceCount > 0 field_gpu = gpuArray(field_hlo); H_gpu = gpuArray(H); field_far = ifft2(fft2(field_gpu) .* H_gpu); field_far = gather(field_far); else field_far = ifft2(fft2(field_hlo) .* H); end

gpuArray把数据搬到显存，gather搬回来。注意H也要转 GPU，否则每次传播都要传，反而慢。如果显存不够，用pagefun批量处理，但代码复杂些，新手先用循环。

5.2 相位解包裹验证拓扑荷

涡旋相位在 [-π, π] 之间有跳变，直接算梯度会出无穷大。用unwrap解包裹后，沿闭合路径积分相位梯度，除以 2π 就是拓扑荷。这个方法能验证模拟结果是不是真的 l 阶涡旋。

% 相位解包裹验证拓扑荷 phase_wrapped = angle(field_1st); phase_unwrapped = unwrap(unwrap(phase_wrapped, [], 1), [], 2); % 沿圆周积分 theta_circle = linspace(0, 2*pi, 360); r_circle = round(size(phase_unwrapped,1)/4); cx = size(phase_unwrapped,2)/2; cy = size(phase_unwrapped,1)/2; phase_circle = interp2(phase_unwrapped, cx + r_circle*cos(theta_circle), cy + r_circle*sin(theta_circle)); % 算相位梯度积分 dphi = diff(phase_circle); dphi(dphi > pi) = dphi(dphi > pi) - 2*pi; dphi(dphi < -pi) = dphi(dphi < -pi) + 2*pi; l_measured = round(sum(dphi) / (2*pi)); fprintf('设定 l=%d, 实测 l=%d\n', l, l_measured);

unwrap先按行再按列，避免二维跳变。interp2在圆周上采样，diff算梯度，处理跳变后求和除以 2π。l_measured应该等于设定值，差 1 以内算正常，差太多检查解包裹是否失败。我一般会在批量扫描后随机抽几个 l 做这个验证，确认程序没跑偏。

5.3 参数表：常用配置速查

这张表是我跑了几十次后总结的，大部分场景直接抄。如果换波长，lambda和fx要同步改，否则衍射角变，一级提取位置偏。z小于 0.1m 时角谱法仍然准，但菲涅尔也够用，看个人习惯。

从那以后我每次跑涡旋模拟，都先拿 l=1 跑一遍完整链路，确认相位螺旋、全息图叉形、一级衍射环、解包裹拓扑荷四个环节都对，再批量扫参数。这个习惯帮我省了至少两周的返工时间。希望帮到你。

本文还有配套的精品资源，点击获取