%% D-FFT  角谱理论衍射.GS（Gerchberg_Saxton）相位复原算法

%% 衍射光场的宽度为L=lamda*d*N/L0. (L0为物尺寸)

close all;

clear all;clc;

%% 衍射参数

d=100e3;%赋值观察屏到衍射面的距离,单位:米

row=512; column=row;        %给出衍射面上的抽样数 row行。column 列。

N=row;%采样点数

lamda=632e-9;        

k0=2*pi/lamda;   %赋值波长、波数

L0=4e-3;     %赋值衍射面尺寸，单位:米

widthObj=2e-3;%孔的尺寸

dxObj=L0/row;%每一个像素代表的实际尺寸

%% 生成  振幅分布

object=zeros(row,column);

object(row/2-round(widthObj/dxObj/2):row/2+round(widthObj/dxObj/2),...

  row/2-round(widthObj/dxObj/2):row/2+round(widthObj/dxObj/2))=1;  %生成衍射孔

% object= imresize((imread('cameraman.tif')),[N,N]);%以图像作为振幅

%% 加载相位

objectPhase = double(imread('westconcordorthophoto.png'));

objectPhase = 2*pi*imresize(objectPhase,[row,column])./max(max(objectPhase));

object=object.*exp(1i*(objectPhase));

Amplitude_object=abs(object);%测量到的物体的空间域中的振幅

%% 是否满足采样条件

Judge=(N>=2*L0^2/lamda/sqrt(d^2+2*L0^2));

%% 生成坐标

fmax=N/(2*L0);%由物理意义得出最高频率

fx1=linspace(-fmax/2,fmax/2,N); %给出频域坐标

fy1=fx1;

[fx,fy]=meshgrid(fx1,fy1);

%% 计算衍射屏上的光振幅分布

Hfxfy= (exp(1i.*k0.*d.*sqrt(1-(lamda.*fx).^2-(lamda.*fy).^2)));

FTobject=fftshift (fft2(object));

FH=FTobject.*Hfxfy;%物频谱与传递函数相乘

U1_amp=abs(ifft2((FH)));%得到衍射屏上应该有的振幅。

%figure,imshow(U1_amp,[]);%

%% GS算法主体

%% 产生随机相位

Phase_retrival=2*pi*rand(row , column);%产生随机相位

for LOOP=1:30               %迭代次数

       object_cal = Amplitude_object.*exp(1i*Phase_retrival);%振幅用实际值

       %% 下面开始用 角谱传播理论 计算衍射

       Hfxfy=exp(1i.*k0.*d.*sqrt(1-(lamda.*fx).^2-(lamda.*fy).^2));

       FTobject= fftshift(fft2(object_cal));

       FH=FTobject.*Hfxfy;%物频谱与传递函数相乘

       U1_cal=U1_amp.*exp(1i*angle(ifft2(FH)));%仅保留衍射面相位。振幅用目标值替换

       I1_cal=U1_cal.*conj(U1_cal);%观察屏上的光强分布

      %% 下面开始用角谱理论逆衍射。

      Hfxfy_reverse= (exp(-1i.*k0.*d.*sqrt(1-(lamda.*fx).^2-(lamda.*fy).^2)));%d取负值

       FTobject_reverse= fftshift(fft2(U1_cal));

       FH_

reverse=FTobject_reverse.*Hfxfy_reverse;%物频谱与传递函数相乘

       U0_cal=ifft2(FH_reverse);%逆衍射的光场

       Phase_retrival=angle(U0_cal);%只保留相位

       I0_cal=U0_cal.*conj(U0_cal);%逆观察屏上的光强分布

end

%% 显示

figure,

subplot(2,2,1),imshow( (Phase_retrival),[]),colormap(gray),title('恢复的相位');

subplot(2,2,2),imshow((abs(object)),[]),colormap(gray),title('物振幅');

subplot(2,2,3),imshow(angle(object),[]),colormap(gray),title('原图相位');

subplot(2,2,4),imshow(((U1_amp)),[]),colormap(gray),title('屏目标振幅');

% TEST = abs(ifft2(FH)) ;

% figure,imshow(TEST,[]);

得到的结果如下：未能正确恢复出相位，不知为何？已经搞了十天了。谢谢。

T_FFT_GS1.txt