简介:本资源面向光学仿真与光束整形方向的学习者与研究人员,围绕GS算法实现高斯光束向高阶高斯光、一阶空心高斯光及贝塞尔高斯光的转换,适合具备Matlab基础、希望深入理解相位恢复与光束调控的读者参考。压缩包共10个文件,包含3个m源码文件、1个mat数据文件,以及5张jpg和1张bmp运行结果图,整体约4.07MB,源码与效果图配套,便于对照验证算法输出。资源中主函数与调用函数分工明确,涵盖正向与逆向变换流程,可帮助读者掌握GS迭代的核心思路与实现细节。目前已有143人学习,适合用于光学课程设计、科研仿真复现及光束整形实验的入门与进阶练习,读者可据此快速搭建可运行的光学仿真环境并观察不同光束模式的转换效果。
1. 从一束高斯光到三种“变形光”:GS 算法到底在算什么
手里有一台 632.8nm 的氦氖激光器,出射的是一束标准基模高斯光,光斑圆、相位平、能量集中。但很多实验要的不是这束光——冷原子捕获想要中空的一阶空心光做暗势阱,光镊和超分辨想要贝塞尔高斯光的无衍射长焦深,光通信模式复用想要高阶高斯光凑模式基。买空间光调制器(SLM)贵,定制相位板又只能对应一种模式。这时候 GS 算法(Gerchberg-Saxton)就是那根救命稻草:它让你用一块普通的相位型 SLM,通过迭代计算出一张相位图,把入射高斯光整形成你想要的任意目标光场。
这套方案解决的核心问题是相位恢复:已知入射光振幅和目标光振幅,求 SLM 上该加载什么相位。它适合做光束整形、光镊、量子光学模式制备的从业者,也适合正在做光学仿真大作业的学生。Matlab 在这里的角色是纯计算引擎——迭代、FFT、误差评估全在矩阵里跑,跑完把相位图导出成 8bit 灰度图喂给 SLM 就行。下面从原理到代码,把三种目标光的生成路径拆开讲清楚。
2. GS 迭代的数学骨架与 Matlab 实现要点
2.1 为什么是“来回傅里叶变换”而不是直接反解
GS 算法的本质是一个约束满足问题。入射面(SLM 平面)和输出面(目标平面)之间是傅里叶变换关系,两个平面各有已知约束:入射面振幅已知(高斯分布),输出面振幅已知(你想要的目标光振幅)。但相位两边都未知,直接反解无解,只能迭代逼近。
一次迭代的流程是:入射面给一个初始相位猜测,做 FFT 到输出面,保留计算出的相位但把振幅替换成目标振幅,再 IFFT 回入射面,同样保留相位替换回入射振幅。如此往复,输出面的振幅分布会逐渐逼近目标。收敛判据通常用输出面实际振幅与目标振幅的均方误差(MSE)或衍射效率。
这里有个容易翻车的地方:FFT 默认把零频放在角落,而物理上光轴在中心。不写fftshift的话,你的目标光会跑到画面四个角拼起来,看起来像鬼影。我一般会在每次变换前后都套一层fftshift,虽然多两次内存搬运,但省掉大量调试时间。
% 基础 GS 迭代核心(单目标光场) N = 512; % 计算窗口分辨率,SLM 常见 512/1024 lambda = 632.8e-9; % 氦氖激光波长 pitch = 8e-6; % SLM 像素间距,按实际器件改 k = 2*pi/lambda; % 入射高斯光振幅(束腰 w0 覆盖约 1/4 窗口) x = (-N/2:N/2-1)*pitch; [X, Y] = meshgrid(x, x); w0 = N*pitch/4; A_in = exp(-(X.^2 + Y.^2)/w0^2); % 目标振幅:这里先用一个空心光占位,后面章节替换 A_target = double(sqrt(X.^2 + Y.^2) > N*pitch/8); phase = 2*pi*rand(N); % 随机初始相位,避免对称停滞 for iter = 1:200 U_in = A_in .* exp(1i*phase); U_out = fftshift(fft2(ifftshift(U_in))); % 到输出面 U_out = A_target .* exp(1i*angle(U_out)); % 替换振幅 U_back = fftshift(ifft2(ifftshift(U_out)));% 回入射面 phase = angle(U_back); % 保留相位 % 误差评估 err(iter) = sum(sum((abs(U_out)-A_target).^2)) / sum(sum(A_target.^2)); end phase_wrapped = mod(phase, 2*pi); % 归一化到 [0,2pi] imwrite(uint8(phase_wrapped/(2*pi)*255), 'phase_gs.png');逻辑说明:ifftshift/fftshift配对保证零频居中;A_target替换时只取相位不取振幅,这是 GS 的硬约束;err用归一化 MSE,一般 50 到 200 次迭代能降到 1% 以下。参数上,N越大光场越细腻但单次 FFT 越慢,512 是 SLM 和算力的平衡点;w0太大会让入射光溢出窗口,太小则高频信息不足,取窗口的四分之一是经验值。pitch必须和你的 SLM 实际像素尺寸一致,填错会导致目标光尺寸整体缩放。
2.2 迭代不收敛时先查这三处
第一,初始相位别用全零。全零相位在对称目标下会让迭代卡在对称解里出不来,随机相位或二次相位(透镜相位)都是好起点。第二,检查目标振幅是不是超出了系统的衍射极限——目标光的最小特征尺寸如果小于lambda*f/(N*pitch),物理上就传不过去,迭代只会震荡。第三,迭代次数不是越多越好,超过 300 次后误差可能因为数值噪声反而回升,我一般设 200 次并在误差曲线拐点处截断。
3. 三种目标光的振幅构造:高阶、一阶空心、贝塞尔高斯
3.1 高阶高斯光:用厄米多项式搭振幅模板
高阶高斯光(HG 模)的振幅是厄米多项式乘以高斯包络。Matlab 没有内置厄米多项式,但递推公式很简单:H0=1,H1=2x,Hn=2xH(n-1)-2(n-1)H(n-2)。构造时把归一化坐标代入即可。阶数 (m,n) 决定光斑的瓣数和节点数,比如 HG01 是两瓣,HG11 是四瓣。
function A = hg_amplitude(X, Y, w0, m, n) % 归一化坐标 xi = sqrt(2)*X/w0; eta = sqrt(2)*Y/w0; Hm = hermite_poly(m, xi); Hn = hermite_poly(n, eta); A = Hm .* Hn .* exp(-(X.^2+Y.^2)/w0^2); A = A / max(abs(A(:))); % 归一化,方便替换 end function H = hermite_poly(n, x) if n == 0, H = ones(size(x)); return; end if n == 1, H = 2*x; return; end Hp = ones(size(x)); Hc = 2*x; for k = 2:n Hn = 2*x.*Hc - 2*(k-1)*Hp; Hp = Hc; Hc = Hn; end H = Hc; end逻辑说明:xi、eta用sqrt(2)/w0归一化是为了让厄米多项式和高斯包络的尺度匹配,否则高阶模会撑得过大。A最后除以最大值做归一化,因为 GS 里目标振幅的绝对幅度不影响相位解,只影响误差的数值尺度。参数上,m、n越大光斑越分散,能量利用率越低,一般不超过 4 阶,再高衍射效率会掉到 50% 以下。w0要和入射高斯光的束腰匹配,差太多会导致模式串扰。
3.2 一阶空心高斯光:拉盖尔多项式加螺旋相位
一阶空心光(LG01 模)的特征是中心光强为零、相位绕轴一圈变化 2π。振幅用拉盖尔多项式 L_p^l,一阶对应 p=0、l=1,此时 L=1,振幅简化为径向坐标乘以高斯包络。相位部分直接叠加一个exp(1i*angle(X+1i*Y))的螺旋相位。
% 一阶空心光振幅(LG01) r = sqrt(X.^2 + Y.^2); A_hollow = (sqrt(2)*r/w0) .* exp(-r.^2/w0^2); A_hollow = A_hollow / max(A_hollow(:)); % 螺旋相位(如果目标本身带轨道角动量) spiral = exp(1i * atan2(Y, X)); % 注意:GS 迭代时目标相位若固定,则替换步骤只换振幅 % 若要让 GS 同时求螺旋相位,则 A_target 只给振幅,相位自由逻辑说明:A_hollow在 r=0 处为零,这正是空心光的物理特征。spiral是可选叠加项——如果你的目标光需要携带轨道角动量(比如光通信模式复用),就在目标场里固定这个相位;如果只是要空心振幅分布,让 GS 自由求相位即可。参数上,w0决定空心环的半径,环半径约等于w0/sqrt(2)。这里有个坑:空心光中心零点的相位是奇点,FFT 在奇点附近数值不稳定,迭代误差会比普通光场高一个量级,属于正常现象,别以为是代码写错了。
3.3 贝塞尔高斯光:用贝塞尔函数截断加高斯窗
贝塞尔高斯光(BG 光)的核心是无衍射——中心主瓣在传播一段距离内不扩散。它的振幅是零阶贝塞尔函数 J0 乘以高斯截断窗。Matlab 的besselj(0, x)直接可用,但要注意 J0 在 x 很大时振荡剧烈,必须用高斯窗截断,否则目标振幅会超出 SLM 的衍射能力。
% 贝塞尔高斯光振幅 kr = 2*pi*sin(0.5*pi/180)/lambda; % 锥角 0.5 度对应的径向波矢 A_bg = abs(besselj(0, kr*r)) .* exp(-r.^2/(2*w0^2)); A_bg = A_bg / max(A_bg(:)); % 检查目标特征尺寸是否可分辨 min_feature = lambda / (2*sin(0.5*pi/180)); fprintf('最小特征尺寸 %.2f um, 像素间距 %.2f um\n', min_feature*1e6, pitch*1e6);逻辑说明:kr由锥角决定,锥角越大主瓣越窄但无衍射距离越短,0.5 度是常见折中。besselj(0, kr*r)产生振荡环,exp(-r.^2/(2*w0^2))把外围环压掉,只留中心几个环。参数上,w0控制截断范围,太小则贝塞尔特性不明显,太大则外围环超出 SLM 调制范围。min_feature那行是自检——如果最小特征尺寸小于两倍像素间距,说明目标光的高频成分 SLM 根本表达不了,迭代再久也是白费,这时候要么加大锥角对应的特征尺寸,要么换更高分辨率的 SLM。
4. 把三种光场跑通:参数配置与结果验证
4.1 一套可复用的参数配置表
三种目标光共用同一套 GS 框架,差异只在A_target的构造函数和少数参数。下面这张表是我实际跑通后整理的配置,直接抄改即可。
| 参数 | 高阶高斯光 | 一阶空心光 | 贝塞尔高斯光 |
|---|---|---|---|
| 目标振幅函数 | hg_amplitude | A_hollow | A_bg |
| 关键参数 | m, n 阶数 | w0 环半径 | 锥角、w0 截断 |
| 典型迭代次数 | 150 | 200 | 250 |
| 收敛误差量级 | 0.5% | 2% | 1.5% |
| 衍射效率 | 70%~85% | 60%~75% | 55%~70% |
| 易翻车点 | 高阶能量分散 | 中心奇点数值噪声 | 外围环超分辨 |
衍射效率定义为目标区域能量除以总入射能量,高阶光因为瓣数多、能量集中度好,效率偏高;空心光和贝塞尔光有大量能量分布在环上,效率天然低一些,这是物理决定的,不是算法问题。
4.2 验证相位图对不对:三步自检
第一步,把算出的相位图重新加载回 GS 正向传播,看输出面振幅和目标振幅的相关系数,大于 0.95 才算合格。第二步,把相位图imwrite成 8bit 灰度 PNG,用图像查看器确认没有大面积纯黑或纯白——那意味着相位被截断,实际 SLM 上会出现相位跳变。第三步,如果条件允许,在光路里放一个 CCD,实际采集的光斑和仿真结果对比,中心位置和环半径应该一致,偏差超过 10% 就要回头查pitch和lambda是否填对。
% 正向验证:用最终相位重建输出面 U_final = A_in .* exp(1i*phase); U_check = fftshift(fft2(ifftshift(U_final))); I_check = abs(U_check).^2; corr_coef = corr2(I_check/max(I_check(:)), A_target/max(A_target(:))); fprintf('振幅相关系数: %.4f\n', corr_coef); imagesc(I_check); axis image; colormap hot; colorbar;逻辑说明:corr2算的是二维相关系数,比 MSE 更直观,0.95 以上说明相位解可用。imagesc画出的光斑应该和目标形状一致,如果出现十字条纹或同心圆伪影,多半是fftshift漏了或者目标振幅有硬边(硬边会产生振铃)。参数上,I_check归一化只是为了显示,不影响相关系数计算。
5. 避坑与排查:GS 光束整形最常见的五个翻车现场
现象一:迭代误差降到某个值就下不去,光斑形状始终差一点。原因通常是目标振幅超出了系统的数值孔径限制,高频成分无法传递。解决方法是把目标振幅做一次低通滤波,或者降低目标光的空间频率(比如减小贝塞尔锥角、降低高阶模阶数)。
现象二:仿真光斑很好,SLM 上实际出射光是散的。原因多半是相位图量化或像素匹配出了问题。SLM 是 8bit 相位调制,mod(phase,2*pi)后量化到 256 级,量化误差会引入额外衍射。解决方法是加一步误差扩散抖动(dithering),或者确认pitch和 SLM 实际像素一致——我见过有人把 8um 的 SLM 按 20um 算,光斑直接大了两倍多。
现象三:空心光中心不是暗的,有个亮点。原因是螺旋相位的奇点被数值噪声填平了。解决方法是构造目标振幅时在中心加一个极小的零区(半径 1~2 像素),强制中心振幅为零,迭代时就不会在奇点处产生伪相位。
现象四:高阶光出现旋转或镜像。这是 FFT 坐标约定问题。Matlab 的fft2默认零频在左上角,fftshift后到中心,但如果ifftshift和fftshift用反了,输出面会整体翻转。检查方法很简单:用一个偏心点源做测试,看输出位置对不对。
现象五:换一台电脑跑,结果不一样。随机初始相位没设种子。rand默认按时间播种,每次结果都不同。解决方法是迭代前加rng(42)固定种子,保证结果可复现。这个坑在写论文或做实验记录时特别致命,血泪经验。
6. 进阶技巧:用加权 GS 和相位平滑把效率再提一档
基础 GS 跑通后,如果想把衍射效率从 70% 推到 85% 以上,有两个实操技巧值得试。第一个是加权 GS:在误差评估时给目标光的关键区域(比如空心环、贝塞尔主瓣)加权重,让迭代优先保证这些区域的振幅精度,牺牲边缘区域。实现上只需在替换振幅那一步乘一个权重矩阵W,误差计算也乘W。第二个是相位平滑:GS 算出的相位图往往有孤立跳变点,SLM 上表现为散射中心。可以在每次迭代后用一次 3x3 中值滤波平滑相位,再mod回 2π,代价是收敛慢一点,但实际光斑质量明显更干净。
% 加权 GS 的替换步骤(其余不变) W = double(sqrt(X.^2+Y.^2) > N*pitch/8); % 空心环区域权重 1 W = 0.3 + 0.7*W; % 边缘保底权重 0.3 U_out = (W.*A_target + (1-W).*abs(U_out)) .* exp(1i*angle(U_out)); % 相位平滑 phase = angle(U_back); phase = medfilt2(phase, [3 3]);逻辑说明:W让目标区域权重为 1、背景为 0.3,迭代会优先满足环上的振幅。medfilt2平滑相位时要注意别在 2π 跳变处直接滤波,否则会把跳变抹成斜坡,正确做法是先解包裹再滤波再重新包裹,但简单场景下直接中值滤波也能用。参数上,权重比例 0.3 是经验值,太低会让背景区域完全失控,太高则加权效果不明显。
我自己做这套东西最大的教训是:别一上来就追求完美光斑,先把基础 GS 跑通、用相关系数验证相位解正确,再逐步加加权和平滑。很多翻车都是因为同时改了太多参数,最后不知道是哪个起了作用。每次只动一个变量,记录误差曲线,这是光学仿真里最笨也最靠谱的习惯。希望帮到你。
本文还有配套的精品资源,点击获取