|||
%% 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,[]);
得到的结果如下:未能正确恢复出相位,不知为何?已经搞了十天了。谢谢。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 08:43
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社