GS算法简介

起源

在光学领域，因为光波的频率过快，一般的探测器不能直接探测到其相位信息，只能获得强度信息，那么如何从强度信息中得到相位信息成为了长时间困扰光学研究者的一个问题。1972年（Gerchberg-Saxton)GS算法被提出，用于从两个已知的强度信息中恢复相位信息。

基本流程

1. 给出物体的初始估计$g_0(x,y)=f(x,y)e^{(j\ast \varphi (x,y))}$,其中$f(x,y)$是物体的振幅分布估计，$\varphi (x,y)$是相位估计，可以是[0,$2\pi$]之间的随机值。
2. 对其作傅里叶变换得到$G_k(f_x,f_y)=\Im (g_0(x,y))=F(f_x,f_y)e^{(j\ast \Phi (f_x,f_y))}$;
3. 在频域内作适当的约束，可以为幅值约束或者支持域约束。对于已知衍射光斑的强度分布$|G^{\prime }(f_x,f_y)|$情况下，可以将频域的模值$F(f_x,f_y)$替换为$G^{\prime }(f_x,f_y)$,保持相位不变；如果对于一个透镜系统而言，因为透镜衍射作用高频信息分量将丢失，相当于一个低通滤波器，可以定义一个光瞳函数
   $$P(f_x,f_y)=\{\begin{array}{rl}& 1，截止频率内\\ & 0，截至频率外\end{array}$$
   随后，光瞳函数乘上零频处于中心的频谱(MATLAB中可以用fftshift(fft2($g_0(x,y)$实现）即可作支持域约束。作完约束之后得到新的频谱分布$$G_k^{\prime }(f_x,f_y)=|G^{\prime }(f_x,f_y)|e^{(j\ast \Phi (f_x,f_y))}$$
   4.对新的频谱进行逆傅里叶变换得到新的空域分布
   $$g_k^{\prime }(x,y)={\Im }^{-1}(G^{\prime }(f_x,f_y))$$
   在空域内再次作恰当的约束，作为下一次迭代的初始分布。
   5.重复上述的4步，知道满足条件输出，条件可以是迭代次数、频域误差、空域误差等。
   

原始GS算法与Fienup算法

现有的大部分相位恢复算法，本质上都是对上述步骤的第4步的约束方法做了更改。
例如，最原始的GS算法，在空域内做约束时，直接用想获得的振幅$f(x,y)$做为振幅约束，但是这种方法容易陷入局部最优解，甚至难以收敛。
而Fienup算法则对约束条件做了修改，假定当前迭代得到的新的空域幅值为$g_k^{\prime }(x,y)$，其与想要得到的振幅的差值就表示为$\Delta g(x，y)=f(x,y)-g_k^{\prime }(x,y)$，约束条件变为$f(x,y)+\alpha \ast \Delta g(x,y)$，其中，$\alpha$是步长，取值范围可以是[0,1]，取值为0时，Fienup算法就变成了原始GS算法。这一操作为迭代过程引入了反馈调节，从直观上来说，如果更新后的振幅为$g_k^{\prime }(x,y)$比想要的振幅$f(x,y)$大，那么振幅的差值$\Delta g(x)$就是个负数，约束条件相当于将初始振幅$f(x,y)$加上一个负数，那么下一次得到的振幅就会变小。正是引入反馈调节这一个操作，能让迭代过程快速收敛。

MATLAB代码

%%% Copyright. first release on April 05, 2021, revised on Sep 21, 2023. 
%%% All rights reserved


clc;
clearvars;
close all;
gray = double(imread('test2.bmp'));               %读入初始图,作为初始迭代分布振幅 
Amplitude = imresize(gray,[512,512]);             %适应调制器分辨率，可注释掉
Amplitude = Amplitude./(max(max(Amplitude)));     %归一化
phase = 2*pi*rand(512,512);                  %产生随机相位
g0_Fie = Amplitude.*exp(1i*phase);           %Fienup算法初始复振幅分布
g0_GS = Amplitude.*exp(1i*phase);            %GS算法初始复振幅分布
  

%fienup算法,对GS算法加入反馈调节量ger*k;
step_size = 0.1;   %设置反馈参量，范围为[0,1]，step_size=0时为GS算法
itera  = 100;
RMS_GS = zeros(itera,1);                       %计算GS算法均方根误差
RMS_Fie = zeros(itera,1);                   %计算Fienup算法均方根误差     
for n = 1:itera    %设置最大迭代次数 
%    Fienup算法
   G0_Fie = fftshift(fft2(g0_Fie));            %傅立叶变换到频域
   G0_FieNew = 1*G0_Fie./abs(G0_Fie); %取相位值,频域作全1幅值约束，相位全息图
   g0_FieNew = ifft2(ifftshift(G0_FieNew));      %逆傅里叶变换返回空域
   g_er = abs(Amplitude) - abs(g0_FieNew)/max(abs(g0_FieNew(:)));  %计算误差矩阵，确定收敛方向
   RMS_Fie(n)=sqrt(mean2((g_er.^2)));        %计算均方根误差
   g0_Fie=(abs(Amplitude)+g_er*step_size).*(g0_FieNew./abs(g0_FieNew)); %引入反馈调节
  
%  GS算法
   G0_GS = fftshift(fft2(g0_GS));   %傅立叶变换到频域
   G0_GSNew = 1*G0_GS./abs(G0_GS); %取相位值,频域作全1幅值约束，相位全息图
   g0_GSNew = ifft2(ifftshift(G0_GSNew));      %逆傅里叶变换返回空域
   g_er = abs(Amplitude) - abs(g0_GSNew)/max(abs(g0_GSNew(:)));     %计算误差矩阵
   RMS_GS(n)=sqrt(mean2((g_er.^2)));        %计算均方根误差
   g0_GS=abs(Amplitude).*(g0_GSNew./abs(g0_GSNew)); %直接用初始振幅约束，不做改变
end

figure(1)
subplot(321);imshow(Amplitude,[]);title('原图');
subplot(323);imshow(angle(G0_FieNew),[]);title('衍射元件相位分布(Fie)');
subplot(325);imshow(abs(g0_FieNew),[]);title('模拟衍射输出(Fie)');
subplot(322);imshow(Amplitude,[]);title('原图');
subplot(324);imshow(angle(G0_GSNew),[]);title('相位原件分布(GS)');
subplot(326);imshow(abs(g0_GSNew),[]);title('模拟衍射输出(GS)');

figure(2)
subplot(121);plot(1:itera ,RMS_Fie);xlabel('循环次数');ylabel('RMS误差(Fie)');
subplot(122);plot(1:itera ,RMS_GS);xlabel('循环次数');ylabel('RMS误差(GS)');

运行结果

图1.Fienup算法和GS算法对比：衍射元件相位分布与模拟输出振幅