简介：这份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$ 的边界处会出现相位跳变。如果只是画相位图，这个跳变看起来像一条接缝，属于正常现象；但如果后面要做全息编码或者相位解包裹，这条接缝就是翻车现场。

我一般会先把网格建好，再算极角和极径，最后乘上拓扑荷。网格分辨率决定了相位奇点附近能不能看清，N=512是起步，N=1024更稳。下面这段代码就是最小可运行版本，直接复制到 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）要做的事情，是把目标复振幅 $E=A\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 干涉法：让涡旋光束和平面波打架

拓扑荷在仿真里设了l=3，但怎么确认衍射出来的真是一级涡旋而不是别的东西？最直接的办法是干涉。让涡旋光束和平面波或者球面波叠加，干涉条纹会出现分叉，分叉数等于拓扑荷。这个在 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。如果l=3，你会看到三条条纹在中心汇合又分开。这个方法在实验里对应马赫-曾德尔干涉光路，仿真里就是矩阵加法。注意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 拓扑荷设为非整数导致相位不闭合

现象：设l=2.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的涡旋，每个通道对应一个独立的衍射级。这个在光通信里叫模分复用，在光镊里叫多光阱。实现思路是给每个拓扑荷分配一个不同的载波频率，让它们在焦平面上空间分离。下面这段代码把l=1,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，先跑一遍角谱传播看焦平面，再跑干涉看分叉数，两个都对上了再保存全息图。这个顺序帮我省了很多后悔药。希望帮到你。

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