基于S-FFT的相息图
基于S-FFT的相息图**
光波是一种电磁波,它满足麦克斯韦方程组,将光波表示成球面波的形式,根据推导,可以获得标量场下光波的菲涅尔衍射积分的表达式如下
(1)
上述式子中分别为衍射面和观察面的平面坐标系,为物平面函数,为波数,为光波波长。将上式改写为离散的傅里叶变换的形式为
(2)
其中表示傅里叶变换,设衍射平面尺寸为,观察平面为,采样数为,则两者的关系可以由下式表示
(3)
根据迭代的相关公式及知识,编写基于S-FFT的反复迭代相息图(MATLAB)程序具体如下:
% 功能: 调入一图像,用菲涅耳衍射的S-FFT计算形成
% 能够由空间光调制器LCOS显示图像的相息图
% 主要变量:
% h ——波长(mm); Ih ——相息图;
% L ——SLM宽度(mm); z0——记录相息
图的距离(mm);
% 取衍射距离1200mm,显示比例1,迭代次数15次
%-------------------------------------------------------------
clear;
close all;
chemin=‘D:\我的资料库\Pictures\3D\figs’;%图片路径
[nom,chemin]=uigetfile([chemin,’.’],[‘调入模拟物体图像’],100,100);
[XRGB,MAP]=imread([chemin,nom]);%获取图像
figure,imshow(XRGB);
X0=rgb2gray(XRGB); %彩色图像转换为灰度图像
figure,imshow(X0,[]);
[M0,N0]=size(X0); %获取灰度图像的像素数大小
N1=min(M0,N0);
N=1080; %相息图取样数, 可按需要修改
m0=0.5 %图像在重建周期中的显示比例,
m0=input(‘图像在重建周期中的显示比例(0->1)’);
X1=imresize(X0,N/N1*m0);
[M1,N1]=size(X1);
X=zeros(N,N);
X(N/2-M1/2+1:N/2+M1/2,N/2-N1/2+1:N/2+N1/2)=X1(1:M1,1:N1);
h=0.532e-3; %波长(mm), 可按需要修改
k=2pi/h;
pix=0.0064; %SLM像素宽度(mm), 可按需要修改
L=Npix; %SLM宽度(mm)
z0=1200 %----衍射距离(mm),
z0=input(‘衍射距离(mm)’);
L0=h*z0/pix %重建像平面宽度(mm)
Y=double(X);
a=ones(N,N);
b=rand(N,N)2pi;
U0=Y.exp(i.b); %叠加随机相位噪声,形成振幅正比于图像的初始场复振幅
X0=abs(U0); %初始场振幅,后面叠代运算用
figstr=strcat(‘SLM平面宽度=’,num2str(L),‘mm’);
figstr0=strcat(‘初始物平面宽度=’,num2str(L0),‘mm’);
figure(1),imshow(X,[]),colormap(gray); xlabel(figstr);title(‘物平面图像’);
np=input(‘叠代次数’);
for p=1:np+1 %叠代次数
%---------------菲涅耳衍射的S-FFT计算开始
n=1:N;
x=-L0/2+L0/N(n-1);
y=x;
[yy,xx] = meshgrid(y,x);
Fresnel=exp(ik/2/z0*(xx.2+yy.2));
f2=U0.Fresnel;
Uf=fft2(f2,N,N);
Uf=fftshift(Uf);
x=-L/2+L/N(n-1);%SLM宽度取样(mm)
y=x;
[yy,xx] = meshgrid(y,x);
phase=exp(ikz0)/(ihz0)exp(ik/2/z0*(xx.2+yy.2));
Uf=Uf.phase;
%---------------菲涅耳衍射的S-FFT计算结束
figstr=strcat(‘SLM宽度=’,num2str(L),‘mm’);
figure(2),imshow(abs(Uf),[]),colormap(gray); xlabel(figstr);title(‘到达SLM平面的物光振幅分布’);
Phase=angle(Uf)+pi;
Ih=uint8(Phase/2/pi255);%形成0-255灰度级的相息图
figure(3),imshow(Phase,[]),colormap(gray); xlabel(figstr);title(‘相息图’);
%---------------菲涅耳衍射的S-IFFT计算开始
U0=cos(Phase-pi)+isin(Phase-pi);
n=1:N;
x=-L/2+L/N(n-1);
y=x;
[yy,xx] = meshgrid(y,x);
Fresnel=exp(-ik/2/z0(xx.2+yy.2));
f2=U0.Fresnel;
Uf=ifft2(f2,N,N);
x=-L0/2+L0/N(n-1);%SLM宽度取样(mm)
y=x;
[yy,xx] = meshgrid(y,x);
phase=exp(-ikz0)/(-ihz0)exp(-ik/2/z0*(xx.2+yy.2));
Uf=Uf.phase;
figure(4),imshow(abs(Uf),[]),colormap(gray); xlabel(figstr0);title(‘逆运算重建的物平面振幅分布’);
%---------------保持相位不变,引用原图振幅,重新开始新一轮计算
Phase=angle(Uf);
U0=X0.(cos(Phase)+i*sin(Phase));
end;
figure(5),imshow(Ih,[]),colormap(gray); xlabel(figstr);title(‘相息图’);
参考书目
[1] 古德曼 ,傅里叶光学导论[M],北京:电子工业出版社,2006,第3版
[2]李俊昌,衍射计算及数字全息[M],北京:科学出版社,2014
[3] 钱晓凡,信息光学数字实验室[M],北京:科学出版社,2014
[4]陈家壁,苏显渝.光学信息技术原理及应用[M].北京:高等教育出版社,2001
上一篇: PyQt5 - Lable中显示gif动画(+作为包导入不显示问题)
下一篇: PCL点云显示