相干光衍射原理&衍射积分代码

衍射计算

摘要

介绍相干光衍射积分的方法和原理，并附加角谱衍射积分代码、泊松亮斑仿真代码，四步相移共轴全息术代码。

衍射积分简介

衍射计算常被用于验证一些相干光现像。比如全息仿真，相位共轭仿真、全息重建等。衍射计算方法有诸多种，但是最准确的还是角谱衍射计算。角谱衍射的传输函数是麦克斯韦方程组(Maxwell方程组)的解析解，所以角谱衍射计算的准确性更高。因为Maxwell方程组有很多解释文章，所以一笔掠过。

衍射计算目前有： 基尔霍夫衍射积分、菲涅尔衍射积分、柯林斯积分、角谱衍射积分等等。目前最常用的是菲涅尔衍射积分和角谱衍射积分。而菲涅尔衍射积分要满足傍轴近似条件的，而且有些情况下，光学实验并不满足傍轴近似（经常遇到）。

衍射积分基础：Maxwell方程组

通过Maxwell(麦克斯韦)方程组可以计算光场分布。Maxwell方程组由四个部分组成：高斯定律、高斯磁定律、法拉第磁感应定律和全电流定律。

高斯定律：在静电场中，穿过任一封闭曲面的电场强度通量只与封闭曲面内的电荷的代数和有关，且等于封闭曲面的电荷的代数和除以真空中的电容率。
公式：
(1)
高斯磁定律：磁场的散度等于零。
公式：
（2）
电磁感应定律：电磁感应定律也叫法拉第电磁感应定律，电磁感应现象是指因磁通量变化产生感应电动势的现象，例如，闭合电路的一部分导体在磁场里做切割磁感线的运动时，导体中就会产生电流，产生的电流称为感应电流，产生的电动势（电压）称为感应电动势。
公式：
（3）
全电流定律：任意一个闭合回线上的总磁压等于被这个闭合回路所包围的面内穿过的全部电流的代数和。
公式：
（4）

根据公式(1)和(3)，以及矢量计算中的关系：(5)
可以得到如下公式：
(6)
为了得到麦克斯韦方程组的解析解，可以从频域的角度入手。定义物体的复振幅分布为U0，像平面的复振幅分布为U。相应的频谱用A和A0表示。
（7）
（8）
将公式(7)和(8)代入公式(6)得到
(9)
所以
自由空间相干传递函数：(10)
其中：

表示空间频率，高频信号沿着小角度传输，低频信号沿着大角度传输。角度α、β角度如图：

代码对图片要求

图片必须是方阵，即行列数量相等。如果不满足，可以自己用0元素扩充为方阵。

仿真四步相移全息

原图
blog.csdnimg.cn/2020052321535435.png)
衍射图：

图像传感器接收图

四步相移后重建图：

仿真泊松亮斑

仿真泊松亮斑所用的圆形屏障

仿真的泊松亮斑（热图）

衍射圆形屏障半径1mm，衍射距离100mm。遮挡物后的中心有明显的光斑。

MATLAB代码部分

衍射函数本体

function [img]=angulardiffraction(obj,sizex,sizey,z,lambda)  
    [pixely,pixelx]=size(obj);                                              %obj 的行和列 
    if pixelx~=pixely
        error('Pixel width and length must be equal');
    end
    k=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加速+采样条件检查），Matlab（2019a）的GPU加速目前只支持NVIDIA，如果没有硬件不符合，请将MaximunPixel赋值为NaN
所有涉及到的坐标系，x正方向为从左至右，y正方向为从上至下

function [img]=angulardiffraction(obj,sizex,sizey,z,lambda)  
    [pixely,pixelx]=size(obj);                                              %obj 的行和列 
    if pixelx~=pixely
        error('Pixel width and length must be equal');
    end
    k=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)，即[