关于matlab的问题

Question

有没有能看懂这串代码的，学长给了串代码我看不懂
每一步具体在干什么，公式是在算什么
具体是关于高斯光束的光束整形

%%%%%%%%%%%%%%%%%%%%
%内能量超过90%的衍射圆
% 衍射距离。程序中的uint是毫米
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lamd=0.193*10^-3;%% 入射波长
w=1.5;%% 入射光束腰
R=1.5;%% 元件的极限孔径
D=700;%% ;衍射距离
r1=1.5;r2=3;%% 圆的内外半径
L0=10;%%mm
k=2*pi/lamd;%%% 波数
N=256;
%%%%%%%%%%%%% 判断（分析变换）
Judging=(sqrt(N*lamd*D)<=L0);%% 计算一个非负数的平方根
if Judging==0
disp('????ERROR');
disp('......Fresnel Analytic Transformation is not satisfied');%%菲涅尔变换不满足
% break;
end
%%%%%%%%%%%%%
x11=linspace(-L0/2,L0/2,N);y11=linspace(-L0/2,L0/2,N);%%y = linspace(x1,x2,n) 生成 n 个点。这些点的间距为 (x2-x1)/(n-1)。
[x1,y1]=meshgrid(x11,y11);
J1=zeros(N);
for m=1:N
for n=1:N
if x1(m,n)^2+y1(m,n)^2<=R^2
J1(m,n)=1;
end
end
end
A=exp(-(x1.^2+y1.^2)/w^2).*J1;
%%%%%
fx=1/L0*(-N/2:N/2-1);fy=1/L0*(-N/2:N/2-1);
[fx,fy]=meshgrid(fx,fy);
%%%%%
JJ=zeros(N);
for m=1:N
for n=1:N
if x1(m,n)^2+y1(m,n)^2>=r1^2&&x1(m,n)^2+y1(m,n)^2<=r2^2
JJ(m,n)=1;
end
end
end
a=sum(sum(A.^2))/sum(sum(JJ.^2));
J2=JJ*sqrt(a);
%imagesc(J2);axis square;colormap(gray)图像SC（J2）；轴线正方形；颜色图（灰色）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% G-S
%pha0=load('C:\Documents and Settings\wjq\My Documents\MATLAB\pha90.txt');
pha0=2*pi*rand(N)-pi;M=0;CC=[];eta=0;
while eta<0.90
M=M+1;
U1=A.*exp(i*pha0);
temp1=(fft2(U1));temp2=exp(i*k*D*(1-lamd^2/2*(fx.^2+fy.^2)));
U2=ifft2(temp1.*temp2);pha2=angle(U2);A2=abs(U2);I2=abs(U2).^2;
Err(M)=sum(sum((I2-J2.^2 ).^2));eta(M)=sum(sum(I2.*JJ))/sum(sum(I2));
subplot(221);plotyy(1:M,Err,1:M,eta)
title('Err (left) and eta (right)');xlabel('M');
subplot(222);imagesc(x11,y11,J2);axis square;title('Object function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude')
subplot(223);imagesc(x11,y11,A2);axis square;title('Output function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude')
%%%%%%%%
U2=J2.*exp(i*pha2);
tmep3=fft2(U2);temp4=exp(-i*k*D*(1-lamd^2/2*(fx.^2+fy.^2)));
U11=ifft2(tmep3.*temp4);pha0=angle(U11);
subplot(224);imagesc(x11,y11,pha0);axis square; title('phase of DOE

xlabel('x / mm');ylabel('y / mm')
pause(0.01)
end

Answer 1

我引用ChatGPT作答：这段代码实现的是高斯光束的光束整形，具体步骤如下：

1 定义输入光束的一些参数，包括入射波长，光束腰大小，元件的极限孔径等。

2 判断菲涅尔变换是否满足条件。

3 生成平面内网格点的坐标，根据孔径大小定义一个矩形区域，并将矩形区域内的像素值赋为1，其他像素为0。

4 计算高斯光束的振幅分布，即在平面内每个像素点上的振幅值，由高斯分布函数和孔径控制函数共同决定。

5 计算傅里叶变换后的频率分布，即在傅里叶空间内每个像素点上的振幅值。

6 定义一个圆形区域，并将圆形区域内的像素值赋为1，其他像素为0。

7 根据圆形区域内的能量百分比，对输入光束的振幅进行缩放，使其满足设定的能量约束条件。

8 通过迭代的方式，将输入光束经过傅里叶变换和逆傅里叶变换，得到目标光束的振幅和相位分布。在迭代的过程中，根据能量收敛的情况，判断是否达到了目标要求。

9 最后输出目标光束的相位分布。
下面是添加了注释的代码，您可以参考一下：

%内能量超过90%的衍射圆
% 衍射距离。程序中的uint是毫米

% 定义变量
lamd=0.193*10^-3; % 入射波长
w=1.5; % 入射光束腰
R=1.5; % 元件的极限孔径
D=700; % 衍射距离
r1=1.5; % 圆的内半径
r2=3; % 圆的外半径
L0=10; % 水平长度（以毫米为单位）
k=2*pi/lamd; % 波数
N=256; % 采样点数

%%%%%%%%%%%%% 判断（分析变换）
% 判断是否满足菲涅尔变换条件
Judging=(sqrt(N*lamd*D)<=L0); % 计算一个非负数的平方根
if Judging==0
    disp('????ERROR');
    disp('......Fresnel Analytic Transformation is not satisfied'); % 菲涅尔变换不满足
    % break;
end
%%%%%%%%%%%%%

% 生成网格点
x11=linspace(-L0/2,L0/2,N); % 在[-L0/2, L0/2]范围内生成N个点
y11=linspace(-L0/2,L0/2,N); % 在[-L0/2, L0/2]范围内生成N个点
[x1,y1]=meshgrid(x11,y11); % 生成网格点

% 生成光强函数J1
J1=zeros(N);
for m=1:N
    for n=1:N
        if x1(m,n)^2+y1(m,n)^2<=R^2 % 如果点在圆内部，则取值为1，否则为0
            J1(m,n)=1;
        end
    end
end

% 生成干涉图的振幅函数A
A=exp(-(x1.^2+y1.^2)/w^2).*J1;

%%%%%
% 生成频率域上的网格点
fx=1/L0*(-N/2:N/2-1);
fy=1/L0*(-N/2:N/2-1);
[fx,fy]=meshgrid(fx,fy);

%%%%%%
% 生成干涉图的掩模函数JJ
JJ=zeros(N);
for m=1:N
    for n=1:N
        if x1(m,n)^2+y1(m,n)^2>=r1^2 && x1(m,n)^2+y1(m,n)^2<=r2^2 % 如果点在圆环内，则取值为1，否则为0
            JJ(m,n)=1;
        end
    end
end

a=sum(sum(A.^2))/sum(sum(JJ.^2)); % 计算放缩系数
J2=JJ*sqrt(a); % 对JJ进行缩放，使其平均光强等于A

%imagesc(J2);axis square;colormap(gray)图像SC（J2）；轴线正方形；颜色图（灰色）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% G-S
%pha0=load('C:\Documents and Settings\wjq\My Documents\MATLAB\pha90.txt');
pha0=2pirand(N)-pi; % 生成一个N×N大小的矩阵，矩阵中的每个元素都是在[-pi, pi]范围内随机生成的一个角度，作为初相。
M=0;CC=[];eta=0;
while eta<0.90 % 当eta小于0.9时，进行循环计算
M=M+1; % 循环次数自加1
U1=A.exp(ipha0); % 通过矩阵乘法生成初始的光强分布
temp1=(fft2(U1)); % 对初始光强分布进行二维傅里叶变换
temp2=exp(ikD*(1-lamd^2/2*(fx.^2+fy.^2))); % 计算平面波前传播的相位调制因子
U2=ifft2(temp1.temp2); % 对相位调制后的光强分布进行二维傅里叶逆变换，得到通过透镜后的光强分布
pha2=angle(U2); % 对通过透镜后的光强分布求相位，得到相位分布
A2=abs(U2); % 对通过透镜后的光强分布求幅值，得到幅值分布
I2=abs(U2).^2; % 对通过透镜后的光强分布求强度，得到强度分布
Err(M)=sum(sum((I2-J2.^2 ).^2)); % 计算输出光强与目标光强之间的误差
eta(M)=sum(sum(I2.JJ))/sum(sum(I2)); % 计算输出光强与目标光强之间的相似度
subplot(221);plotyy(1:M,Err,1:M,eta) % 绘制误差和相似度的变化曲线
title('Err (left) and eta (right)');xlabel('M');
subplot(222);imagesc(x11,y11,J2);axis square;title('Object function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') % 绘制目标光强分布的图像
subplot(223);imagesc(x11,y11,A2);axis square;title('Output function')
xlabel('x / mm');ylabel('y / mm');zlabel('Amplitude') % 绘制透镜后的光强分布的图像
%%%%%%%%
U2=J2.exp(ipha2); % 将通过透镜后的光强分布的相位分布与目标光强分布相乘，得到下一次迭代的初始光强分布
tmep3=fft2(U2); % 对下一次迭代的初始光强分布进行二维傅里叶变换
temp4=exp(-ikD*(1-lamd^2/2*(fx.^2+fy.^2))); % 计算平面波后传播的相位调制因子
U11=ifft2(tmep3.*temp4); % 利用傅里叶变换将初始光强分布与相位调制因子相乘得到通过光阑后的光强分布，然后再进行傅里叶反变换得到下一次迭代的光强分布
subplot(224);imagesc(x11,y11,pha0);axis square; title('phase of DOE xlabel('x / mm');ylabel('y / mm') %绘图
pause(0.01)
end