★ FEATURED ARTICLE 简介这份MATLAB程序面向光学、物理与光学工程方向的学习者和研究者用于模拟涡旋光束与全息光栅的相互作用帮助理解螺旋波前、涡旋相位及拓扑荷对光束传播特性的影响。程序以贝塞尔函数构建涡旋光束相位分布结合傅里叶光学方法模拟光束经全息光栅后的衍射过程并通过fft2、ifft2等函数完成二维傅里叶变换与逆变换直观展示不同衍射级次的光束分布。调整拓扑荷参数还可观察LG光束等不同模式的旋涡性质变化为光学通信、信息处理与量子光学等场景提供实验设计参考。资源包为rar压缩格式仅含1个m源码文件体积约1KB结构精简便于直接阅读与二次修改。目前已有1938人学习下载适合希望快速上手涡旋光束全息模拟、掌握拓扑荷调控思路的读者参考使用。1. 涡旋光束全息与拓扑荷模拟一份能直接跑通的 MATLAB 程序包做光场调控的同行大概率都遇到过这种场景想验证一个螺旋相位板产生的涡旋光到底是不是目标拓扑荷实验室里没有 SLM 可以实时调或者想先看看不同拓扑荷下相位分布和强度环的对应关系结果卡在“从哪写起”上。这份涡旋光束全息与拓扑荷模拟程序就是冲着这个痛点来的——它把涡旋相位生成、全息图计算、拓扑荷扫描、光场传播这几件事用 MATLAB 串成了一条完整链路输入拓扑荷和光束参数就能出相位图、全息图和衍射光斑。适合光学工程、光通信、超分辨成像方向的研究生和工程师尤其是需要快速出图验证想法、又不想从零推公式的人。程序本身不依赖额外工具箱纯 MATLAB 基础函数就能跑对版本要求不高2023b 到 2026b 实测都能用。2. 涡旋相位与拓扑荷的数学建模从螺旋相位板到全息图2.1 涡旋光束的相位结构为什么用 exp(i·l·θ) 描述涡旋光束的核心特征是相位波前绕光轴旋转每绕一圈相位变化 2πll 就是拓扑荷。数学上写成 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); % 取 l0 的全息图 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算错一级没对齐相位图就是乱的。我一般会先跑一个 l0 的看零级和一级位置确认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 扫到 50MATLAB 报内存不足。原因三维矩阵N×N×100在 N1024 时占 800MB 以上。解决只存全息图不存相位或者每算一个 l 就存成单独文件用matfile追加。4.4 传播后相位不是涡旋而是共轭现象一级衍射相位螺旋方向和预期相反。原因载波频率符号和 l 符号不匹配取了一级衍射的共轭级。解决根据 l 正负翻转fx符号或者在提取一级时取另一侧。4.5 环半径和理论对不上差一个数量级现象实测环半径比理论大 10 倍。原因pixel_size算错用了L/N但 L 单位是 mm 不是 m。解决统一用米L5e-3pixel_size L/N输出时再转 mm。5. 进阶技巧用 GPU 加速和相位解包裹验证拓扑荷5.1 GPU 加速批量传播拓扑荷扫描时每个 l 都要做一次 FFT 和角谱传播CPU 跑 100 个 l 要几分钟。有 GPU 的话把fft2换成gpuArray的fft2速度能快 5 到 10 倍。注意显存要够N1024 时单精度复数矩阵约 16MB100 个 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); endgpuArray把数据搬到显存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 参数表常用配置速查参数典型值说明N512网格尺寸1024 更精细但慢L5e-3 m物理尺寸按 SLM 靶面改lambda632.8e-9 m氦氖激光波长w01e-3 m束腰取 L/5 左右z0.5 m传播距离远场取 1m 以上fx1/(20*pixel_size)载波频率20 像素一周期l-5:5拓扑荷范围按需扩这张表是我跑了几十次后总结的大部分场景直接抄。如果换波长lambda和fx要同步改否则衍射角变一级提取位置偏。z小于 0.1m 时角谱法仍然准但菲涅尔也够用看个人习惯。从那以后我每次跑涡旋模拟都先拿 l1 跑一遍完整链路确认相位螺旋、全息图叉形、一级衍射环、解包裹拓扑荷四个环节都对再批量扫参数。这个习惯帮我省了至少两周的返工时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
