相干光衍射原理衍射积分代码
相干光衍射计算
衍射计算 摘要
介绍相干光衍射积分的方法和原理，并附加角谱衍射积分代码、泊松亮斑仿真代码，四步相移共轴全息术代码。
衍射积分简介
衍射计算常被用于验证一些相干光现像。比如全息仿真，相位共轭仿真、全息重建等。衍射计算方法有诸多种，但是最准确的还是角谱衍射计算。角谱衍射的传输函数是麦克斯韦方程组(方程组)的解析解，所以角谱衍射计算的准确性更高。因为方程组有很多解释文章，所以一笔掠过。
衍射计算目前有： 基尔霍夫衍射积分、菲涅尔衍射积分、柯林斯积分、角谱衍射积分等等。目前最常用的是菲涅尔衍射积分和角谱衍射积分。而菲涅尔衍射积分要满足傍轴近似条件的，而且有些情况下，光学实验并不满足傍轴近似（经常遇到）。
衍射积分基础：方程组
通过(麦克斯韦)方程组可以计算光场分布。方程组由四个部分组成：高斯定律、高斯磁定律、法拉第磁感应定律和全电流定律。
高斯定律：在静电场中，穿过任一封闭曲面的电场强度通量只与封闭曲面内的电荷的代数和有关，且等于封闭曲面的电荷的代数和除以真空中的电容率。
公式：
(1)
高斯磁定律：磁场的散度等于零。
公式：
（2）
电磁感应定律：电磁感应定律也叫法拉第电磁感应定律，电磁感应现象是指因磁通量变化产生感应电动势的现象，例如，闭合电路的一部分导体在磁场里做切割磁感线的运动时，导体中就会产生电流，产生的电流称为感应电流，产生的电动势（电压）称为感应电动势。
公式：
（3）
全电流定律：任意一个闭合回线上的总磁压等于被这个闭合回路所包围的面内穿过的全部电流的代数和。
公式：
（4）
根据公式(1)和(3)，以及矢量计算中的关系：
(5)
可以得到如下公式：
(6)
为了得到麦克斯韦方程组的解析解，可以从频域的角度入手。定义物体的复振幅分布为U0，像平面的复振幅分布为U。相应的频谱用A和A0表示。
（7）
（8）
将公式(7)和(8)代入公式(6)得到
(9)
所以
自由空间相干传递函数：
(10)
其中：
表示空间频率，高频信号沿着小角度传输，低频信号沿着大角度传输。角度α、β角度如图：
代码对图片要求
图片必须是方阵，即行列数量相等。如果不满足，可以自己用0元素扩充为方阵。
仿真四步相移全息
原图
/25435.png)
衍射图：
图像传感器接收图
四步相移后重建图：
仿真泊松亮斑
仿真泊松亮斑所用的圆形屏障
仿真的泊松亮斑（热图）
衍射圆形屏障半径1mm，衍射距离100mm。遮挡物后的中心有明显的光斑。
代码部分
衍射函数本体
function [img]=angulardiffraction(obj,sizex,sizey,z,lambda) [pixely,pixelx]=size(obj); %obj 的行和列 if pixelx~=pixelyerror('Pixel width and length must be equal');endk=2.*pi./lambda; %定义波矢量%% 计算频域坐标和obj频谱originpoint=ceil((pixelx+1)/2); %寻找0频率坐标dfx=1/sizex;fx=[1-originpoint:-1,0:pixelx-originpoint].*dfx; fx=ones(pixelx,1)*fx;dfy=1/sizey;fy=[1-originpoint:-1,0:pixelx-originpoint]'.*dfy;%注意，fft2后，y方向频谱的正方向是向下的fy=fy*ones(1,pixely);fobj=fftshift(fft2(obj)); %傅里叶变换clear('obj');%% 计算像面频谱fimg=fobj.*exp(1j.*k.*z.*sqrt(1-(lambda.*fx).^2-(lambda.*fy).^2)); %原图与真空CTF作用clear('fobj');img=ifft2(ifftshift(fimg)); %逆傅里叶变换
end
衍射函数（GPU加速+采样条件检查），（2019a）的GPU加速目前只支持，如果没有硬件不符合，请将赋值为NaN
所有涉及到的坐标系，x正方向为从左至右，y正方向为从上至下
function [img]=angulardiffraction(obj,sizex,sizey,z,lambda) [pixely,pixelx]=size(obj); %obj 的行和列 if pixelx~=pixelyerror('Pixel width and length must be equal');endk=2.*pi./lambda; %定义波矢量MaximumPixel=5000; %显存支持最大图片宽度gpuWork=false; %判断GPU加速的指示变量%% 计算频域坐标和obj频谱originpoint=ceil((pixelx+1)/2); %寻找0频率坐标dfx=1/sizex;fx=[1-originpoint:-1,0:pixelx-originpoint].*dfx; %按照原理，采样率为fs的频谱，其频谱范围为[0,fs)，根据Nyquist采样定理，只有%[0,fs/2)有信息，[fs/2,fs)为对称谱。以256为例，则信息为[0,128)像素，即0至%127，共128个像素点，频率为[0:127*df]；对称一侧为[128,256)，同为128个像素点，%频率为[128*df,256*df)，即[128*df,255*df]，对称后为[-128*df,-df]。%matlab的fftshift函数会将0频率移动至频谱中心，在处理频谱中心时候，会根据像素%的奇偶不同，中心点选择不同。originpoint=ceil((pixelx+1)/2)便是寻找中心%点位置，中心点对应的频率为0频。[1:pixel]个像素移动后，其分布为%[1,originpoint)[originpoint,pixel]，右侧为正频率，左侧为负频率。%采样频率由单个像素点的空间周期决定，即fs=pixel/size。最小频率fm=1/size。%fft2后，每个像素的频率跨度为fm，所以正频谱为%[originpoint-originpoint,pixel-originpoint]*fm，%即[0,pixelx-originpoint]*fm。%负频谱为[1-originpoint,originpoint-originpoint)*fm，即%[1-originpoint,0)*fm，改写为[1-originpoint,-1]*fm。fx=ones(pixelx,1)*fx;dfy=1/sizey;fy=[1-originpoint:-1,0:pixelx-originpoint]'.*dfy;%注意，fft2后，y方向频谱的正方向是向下的fy=fy*ones(1,pixely);if pixelx<=MaximumPixel %将RAM移动至显存obj=gpuArray(obj);fx=gpuArray(fx);fy=gpuArray(fy);gpuWork=true;endfobj=fftshift(fft2(obj)); %傅里叶变换clear('obj');%% 计算频谱丢失情况(可用，根据需求使用)alpha=acot(sizex/2/z);frequencyx=cos(alpha)/lambda;beta=acot(sizey/2/z);frequencyy=cos(beta)/lambda;% pixelxn=power(frequencyx,-1)/(sizex/pixelx);% pixelyn=power(frequencyy,-1)/(sizey/pixely);% disp(['x y分辨最小像素:（',num2str(pixelxn),',',num2str(pixelyn),')']);%% 根据频谱丢失情况，判断是否满足采样条件decay=1;PositiveMaxFx=uint32(find(fx>frequencyx*decay));NegativeMaxFx=uint32(find(fx<-frequencyx*decay));PositiveMaxFy=uint32(find(fy>frequencyy*decay));NegativeMaxFy=uint32(find(fy<-frequencyy*decay));if ~(isempty(PositiveMaxFx)&&isempty(NegativeMaxFx)&&isempty...(PositiveMaxFy)&&isempty(NegativeMaxFy))warning('Function "angulardiffraction" failed to sample. Frequency lost.');end%% 计算像面频谱fimg=fobj.*exp(1j.*k.*z.*sqrt(1-(lambda.*fx).^2-(lambda.*fy).^2)); %原图与真空CTF作用clear('fobj');img=ifft2(ifftshift(fimg)); %逆傅里叶变换%%if gpuWork==true %如果启用GPU加速，则将显存移至RAMimg=gather(img);end
end
泊松亮斑仿真主函数
clear;
clc;
nano=1e-9;
meter=1;
milli=1e-3;
micron=1e-6;%定义单位pixelx=2000;
pixely=2000;%定义像素点数量%计算总尺寸
sizex=5*micron*pixelx;%每个像素点的大小为5*micron
sizey=5*micron*pixely;%总大小=像素点尺寸*像素点数量lambda=532*nano;
k=2*pi/lambda;
disp([num2str(sizex/milli),'×',num2str(sizey/milli),'mm']);%% 泊松亮斑用的圆形掩模
m=pixely; %矩阵的函数
n=pixelx; %矩阵的列数
r=200; %生成圆的半径
m1=-m/2:m/2-1; %把圆心变到矩阵的中间
n1=-n/2:n/2-1;
[x,y]=meshgrid(m1,n1);
circle=x.^2+y.^2; %计算出每一点到圆心的距离的平方circ_mask=zeros(m,n);
circ_mask(find(circle<=r*r))=0; %找到圆内的元素，并复制为0
circ_mask(find(circle>r*r))=1; %找到圆外的元素，并复制为1
imshow(circ_mask,[]);
disp(['半径',num2str(r*5*micron/milli),'mm']);%% 衍射计算
z=100*milli;%衍射距离
obj=circ_mask;%衍射物体为圆屏
[img]=angulardiffraction(obj,sizex,sizey,z,lambda);%衍射
imshow(abs(img).^2,[]),colormap('hot');
四步相移全息术主函数：
%---注意:本文利用角谱传播，其中包括了傅里叶变换，注意采样条件----%
% Notice:This promgram adopts the angular diffraction. %
% It includes the Fourier Transform. Pay attention to %
%---the sampling condition ---%
%% 单位定义 Define unit %%
clear all;clc;
cm=0.01;um=1e-6;mm=0.001;%% 载入图像 load the image %%
load('eword.mat');
B=eword;
clear('eword');%% 将图片变为方阵 Change the object into a square matrix %%
[m,n]=size(B);
ms=max([m,n]);
A=zeros(ms,ms);
A(1+floor((ms-m)/2):m+floor((ms-m)/2),1+floor((ms-n)/2):n+floor((ms-n)/2))=B;times=10;
[m,n]=size(A);
multiples=2;
M=2.*m.*multiples;
N=2.*n.*multiples;
%% 设定仿真参数 Parameter setting of simulation----------%%
%z: 传输距离 Diffraction distance %
%lambda: 波长 Wavelength %
%k: 波矢量 Wave vector %
%hx&hy: 原始图大小 Size of image% %
%dhx&dhy: 尺寸对像素微分 Differential of size to pixel %
%--------------------------------------------------------%
z=100.*cm; %传输距离 Diffraction distance
lambda=0.63*um; %波长 Wavelength
k=2*pi/lambda; %波矢量 Wave vector
hx=5*cm;hy=5*cm; %原始图大小 Size of image
dhx=hx/N;
dhy=hy/M;
x=dhx.*(ones(N,1)*[-M/2:M/2-1]);
y=(dhy.*(ones(N,1)*[-M/2:M/2-1]))';%%
figure;imshow(abs(A).^2); %原始图
[lm,ln]=size(A);%% 扩充图像 Padding the image %%
f0=zeros(N,1)*[1:M];
light0=f0;
light=ones(lm,ln);
f0(M/2-M/4/multiples+1:M/2+M/4/multiples,N/2-N/4/multiples+1:N/2+N/4/multiples)=A; %扩充图像%% 设定激光器发出光束 Setting of laser ------------------------------- %%
% light0 定义了激光束的形状, 在此程序中, light0为256x256的方形，变量f0是 %
% 扩充后的物体图 %
% light0 decide the shape of laser. In this program, it's square with %
% 256x256 pixel. The variable f0 is padding image of object %
% ------------------------------------------------------------------ %
light0(M/2-M/4/multiples+1:M/2+M/4/multiples,N/2-N/4/multiples+1:N/2+N/4/multiples)=light;%设置激光器光束形状和大小（方形，256x256）
light0=angulardiffraction(light0,hx,hy,10.*cm,lambda);%激光器发出光束传播至被拍照物体
clear('light');
f0=light0.*f0;%激光与目标作用
[image] = angulardiffraction(f0,hx,hx,z,lambda);%被拍照物体衍射至全息干板
figure('name','diffraction');imshow(mat2gray(abs(image).^2));%% 设定参考光 Setting of reference light ------ %%% alpha :angle between x axis and light vector %% beta :angle between y axis and light vector %% alpha :x轴与光矢量之间的夹角 %% beta :y轴与光矢量之间的夹角 %% ----------------------------------------------%
alpha=90/180*pi;beta=(90)/180*pi;
pic=zeros(M,N);
pic(M/2-M/4+1:M/2+M/4,M/2-M/4+1:M/2+M/4)=1;
ur=exp(1j*k*(cos(alpha)*x+cos(beta)*y)).*pic;%refernce light,parallel light%% 四步相移 Four step phase shift----------------------%%% 参考:四步相移; 功能:从 hologram0~3提取图像复振幅 %% Method reference:This formula refers to four step %% phase shift method %% Function:Extract the origin image's complex %% amplitude from hologram0~3 %%----------------------------------------------------%
%reference light phase shift 四步相移
ur0=ur;
ur1=ur.*exp(1j.*pi./2);
ur2=ur.*exp(1j.*pi);
ur3=ur.*exp(1j.*pi.*3./2);
%inteference between object and reference，物光、参考光干涉 %
hologram0=ur0+image;
hologram1=ur1+image;
hologram2=ur2+image;
hologram3=ur3+image;hologram0=abs(hologram0).^2;
hologram1=abs(hologram1).^2;
hologram2=abs(hologram2).^2;
hologram3=abs(hologram3).^2;hologram=((hologram0-hologram2)+1j.*(hologram3-hologram1))./4;
%% 显示衍射结果 Show the results of diffraction %%phi=angle(hologram);
figure('name','hologram');imshow(mat2gray(abs(hologram).^2));
figure('name','phase map');imshow(mat2gray(phi));%% -----重建目标 Reconstruction the object-------------------%%
% 重建目标的原理参考全息术技术 %
% Reconstruction the object refer to principle of hologram %
% ------------------------------------------------------ %
image=hologram.*conj(ur);
figure('name','phase map after illustrate');imshow(mat2gray(angle(image)));
[imageR] = angulardiffraction(image,hx,hy,z,lambda);
phiR=angle(imageR);
imageR=abs(imageR).^2;
figure('name','reconstruction');imshow(mat2gray(imageR));
figure('name','phase of reconstruction');imshow(mat2gray(phiR));
clear('hologram0');
clear('hologram1');
clear('hologram2');
clear('hologram3');
clear('image');
clear('image_cut');
clear('light0');参考文献：
[1]苏显渝, 吕乃光, 陈家璧. 信息光学原理[M]. 电子工业出版社, 2010.
