GS算法流程
GS（Gerchberg-Saxton）相位恢复算法是一种基于傅里叶变换的最速下降算法，可以通过输出平面和输入平面上光束的光强分布计算出光束的相位分布。图1是基于GS算法的涡旋光束畸变波前校正系统框图，在该框图中，已知输入平面涡旋光束振幅和相位以及输出平面的涡旋光束振幅，经过算法多次迭代运算并满足误差需求后可以得到需要的校正相位信息，进而实现校正效果。
图1：GS算法流程[1]
GS算法恢复涡旋光束畸变波前的具体流程如下[2]：
- 将LG涡旋光束作为输入平面初始参考光，计算出参考光的振幅和相位。
- 将初始参考光束进行快速傅里叶变换，得到其幅度谱A0(kx,ky)和相位谱Φ0(kx,ky)。
- 使用畸变光束的幅度谱替换A0(kx,ky)，将替换后的幅度谱与相位谱Φ0(kx,ky)组合成新的函数，再通过傅里叶逆变换即可得到光束的空域幅度a0(x,y)和相位H0(x,y)。
- 用参考光的振幅代替a0(x,y)，并与H0(x,y)组合成新函数，将新的函数做快速傅里叶变换后得到新的幅度谱A1(kx,ky)和相位谱Φ1(kx,ky)。将A1(kx,ky)与畸变涡旋光束的幅度谱做差，若满足误差则输出H(x,y)，否则就使用新的幅度谱和相位谱进行下一轮的计算。
- 得出预估相位与理想相位差距，对畸变相位进行校正。
仿真结果展示
文章中所用的光束为拉盖尔高斯涡旋光束，光束经大气湍流MATLAB仿真详见之前的文章：MATLAB仿真：涡旋光束经大气湍流后的光强及相位分布基于蒙特卡洛相位屏法
仿真参数：波长lambda=632.8e-9m；波失k = 2 * pi / lambda；束腰w=5e-3；拓扑荷l=3；外尺度L0=1.0m；内尺度l0=1.0e-3；传播距离z=1000m；湍流屏张数n=10
如图2展示了不同大气湍流下迭代次数为500轮的GS算法校正结果。由图可以看出GS算法可以有效校正大气湍流带来的光束畸变。图3展示了在不同迭代次数较强湍流中涡旋光束的校正结果。可以看到随着迭代次数增加，校正效果越好，能量越集中。
图2不同大气湍流下迭代次数为500轮的GS算法校正结果。
图3：不同迭代次数较强湍流中涡旋光束的GS校正结果。
参考文献：
- 《大气湍流下涡旋光束波前校正及模式检测研究》
- 《液晶相控阵波前相位校正算法研究》
MATLAB部分仿真代码：
喜欢的话给小编点个关注点个赞呀
有任何问题可以在评论区留言
更多完整代码进T店博士生牛马赚点小钱
clc;clear;close;
%% 坐标
N=512;
delta1=0.6/N;
[x,y]=meshgrid((-N/2:1:N/2-1).*delta1);
[phi,r1]=cart2pol(x,y);
deltaf1=1/(N*delta1);
%% 光束参数
lambda=632.8e-9;%m
k = 2 * pi / lambda;
w=5e-3;
l=3;
%% 模拟测量面经过湍流的LG光束光强
Cn2=1.0e-13;
L0=1.0;
l0=1.0e-3;
z=1000;
n=10;
dz=z/n;
r0=(0.423 * k^2 * Cn2 *dz) .^(-3/5);
E1= (sqrt(r1.^2)./w).^abs(l).*exp(-(r1.^2)./(w.^2)).*exp(1i.*l.*phi);
%% 初始光场假定未经过湍流的自由空间传播光强
E0= (sqrt(r1.^2)./w).^abs(l).*exp(-(r1.^2)./(w.^2)).*exp(1i.*l.*phi);
[xn, yn, G0]=ang_spec_prop(E0,lambda,delta1,delta1,z);
%gs循环
max_iter = 1000;%迭代次数
iter = 0;
epsilon = 1e-3;
ER = inf;figure(1)
subplot(2,3,1)
imagesc(I_vac);axis square;title('未畸变光强','FontSize',15);axis off
subplot(2,3,2 )
imagesc(I_tur);axis square;title('畸变光强','FontSize',15);axis off
subplot(2,3,3 )
imagesc(I_correct);axis square;title('校正后光强','FontSize',15);axis off
subplot(2,3,4 )
imagesc(p0);axis square ; title('未畸变相位','FontSize',15);axis off
subplot(2,3,5 )
imagesc(p12);axis square;title('畸变相位','FontSize',15);axis off
subplot(2,3,6 )
imagesc(pcor);axis square;title('校正后的相位','FontSize',15);axis off
