衍射函数（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);endend
泊松亮斑仿真主函数
clear;clc;nano=1e-9;meter=1;milli=1e-3;micron=1e-6;%定义单位pixelx=2000;pixely=2000;%定义像素点数量%计算总尺寸sizex=5*micron*pixelx;%每个像素点的大小为5*micronsizey=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;%找到圆内的元素，并复制为0circ_mask(find(circle>r*r))=1;%找到圆外的元素，并复制为1imshow(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');