简介：面向光学工程师、科研人员及高校学生，这份MATLAB源码包聚焦DOE（衍射光学元件）设计中的GS算法（基于傅里叶变换的IFTA）与衍射计算。压缩包大小仅1KB，内含3个MATLAB脚本（.m文件），分别对应算法初始化、迭代优化与结果输出环节，可直接运行并观察相位变换过程。已有7959人学习下载，是光学计算领域较受关注的基础实现代码。脚本完整覆盖目标光场定义、傅里叶变换、相位迭代更新、反变换及误差判断等核心步骤，并借助fft/ifft函数完成频域与空间域转换。通过阅读和实践这套源码，读者可深入理解GS算法的计算流程与收敛条件，掌握DOE相位分布的设计方法，并能迁移应用到激光整形、自由空间通信、光学加密等实际场景，对从事信息光学与光束控制的研究者尤为实用。
1. 从零到能出图的DOE设计：这个GS算法的MATLAB包能干什么
拿一块纯相位透明片让照明光在焦面上拼出一个目标图案，这种听起来像光学玄学的需求，现在用一份MATLAB源码就能跑通。这套工程的核心是GS算法，也叫Gerchberg-Saxton迭代相位恢复，工程里最常见的落点就是DOE（衍射光学元件）设计：给它一张目标强度图，它迭代几十到上百轮，给你一个纯相位分布，把这个相位刻成台阶或浮雕就是DOE。适合两类人：一类是刚接手DOE设计任务的开发者，需要一个能跑、能出图、能改的基线工程；另一类是研究衍射计算、想要IFTA迭代框架做二次开发的人。标题里“基于matlab设计DOE元件”不是虚的，它把从参数设置、迭代计算到重建验证的一整条链都放进了src目录，属于拿来就能跑出第一版结果的范畴。
2. 算法与载体：GS与IFTA的关系，以及MATLAB现场为什么顺手
2.1 一个请求的链路：平面波从DOE面到目标面经历了什么
DOE设计要先建立物理模型，这件事不绕开。常见场景设为一个单色平面波垂直入射到DOE面上，DOE面是一个纯相位元件，只改入射光的相位、不改振幅；光继续传播一段距离z，在目标面上形成强度分布。这个传播过程在标量衍射条件下，可以抽象成一次二维傅里叶变换：DOE面上的复振幅做FFT，得到目标面的复振幅。GS算法要做的就是反着来——已知目标面想要的振幅分布，反推出DOE面的相位分布。
工程里对传播模型的处理分两派。一派用S-FFT（单次快速傅里叶变换），直接把DFT当衍射积分，代码最短，大部分GS教学代码都用这种；另一派用角谱传播（ASM），需要把传递函数显式写出来，代码长一点但更可控。你先跑通S-FFT版本就够了，后面二次开发再换角谱法。
这个链路里有个容易被忽略的细节：目标面的采样间距并不等于DOE面的采样间距。用S-FFT时，目标面的网格间距dx2由公式dx2 = lambda * z / (N * dx1)决定，其中N是网格边长，dx1是DOE面采样间距。也就是说，你设置N、dx1、z这组参数的同时，也把目标面的物理尺寸定死了。工程里跑出重建图，如果发现点阵间距和你预先想象的不一致，先回头检查这三者的匹配关系，而不是怀疑算法写错了。
2.2 GS算法与IFTA：同一条迭代主线的两个名字
IFTA全称是迭代傅里叶变换算法，GS算法是它最经典、最朴素的一个实现。两者不是两个东西，是一族方法。IFTA的核心套路就四步：随机初始相位起手；前向传播到目标面，把幅值替换成目标幅值、保留相位；逆传播回DOE面，把幅值约束成均匀（理想平面波），保留相位；循环直到重建像的误差不再明显下降。
GS能收敛的本质，是它在两个约束集合之间做交替投影：一个集合是“DOE面振幅等于1”的约束，另一个是“目标面振幅等于sqrt(I_target)”的约束。每一次正向和逆向传播，就是把当前相位分布往这两个集合上各投影一次。投影到交集时，解就出现了。这个视角对调试很有用——如果你看到误差曲线在某个值附近震荡下不去，说明两个约束集合的交集里，没有更优解了，要么放宽约束（改成加权GS），要么改约束本身（在目标面允许一个误差容限）。
工程源码里没有把IFTA单独拆成一个包，它就是把GS迭代循环封装成了一个函数。模块命名沿用了IFTA的思路，函数名带forward / inverse / amplitude_constrain这类词，读起来不费劲。
2.3 载体选择与运行环境：MATLAB加哪些工具箱
用MATLAB做这个事，优势在数据结构和调试体验。DOE面上的一个像素，在代码里就是一个复数矩阵元素，复振幅天然对应复数矩阵，fft2、ifft2、angle、abs都是内置的，不需要像C++那样先封一层复数数组再手写FFT。调试的时候，在循环里设断点，直接看phase矩阵的数值和分布，直观程度比Python加NumPy还要顺一点。出图也更省事：imshow看强度，surf看相位，plot看误差曲线。
这份工程跑起来只需要基础MATLAB环境，我建议用2019a之后的版本。工具箱方面，核心计算只用到了基础功能，Image Processing Toolbox在显示重建图时会用到imshow、mat2gray，实际上即使没有这个工具箱，用imagesc也能替代。所以安装MATLAB时勾上默认的基础工具箱就够跑通，不需要额外装Simulink或专用光学工具包。
运行准备只需要三步：确认路径下能访问到src目录里的所有.m文件；在工作区准备好目标强度图变量I_target（可以用mat文件里的现成数据）；直接运行入口脚本。脚本会把参数区、迭代循环、重建验证全部串起来，这个过程在下一章拆开讲。
3. 核心模块拆解：从目标强度图到重建光斑的完整链路
3.1 主入口脚本：先看参数区，再跑通一次设计
工程把主流程放在一个入口脚本里，结构很清楚：参数区、调用迭代核心、重建验证。参数区是最值得先改的地方，因为DOE设计里大多数“结果不对”都不是算法问题，是参数没配对。先看一段典型的参数区代码。
% main_design.m % 目标：10x10 点阵分束器的相位设计 clear; clc; % ---------- 物理参数 ---------- lambda = 532e-9; % 波长，单位米，这里取绿光532nm z = 0.1; % DOE面到目标面的传播距离，单位米 dx1 = 8e-6; % DOE面采样间距，单位米 N = 256; % 网格边长（方阵），单位像素 iter = 100; % GS迭代轮数 % ---------- 目标强度分布 ---------- % 用5个点构造一行弱目标，方便快速验证迭代是否收敛 I_target = zeros(N, N); step = 40; I_target( (0:4)*step + round(N/2) - round(4*step/2), ... round(N/2) ) = 1; I_target = I_target / sum(I_target(:)); % 能量归一化 % ---------- 调用GS迭代核心 ---------- phase = gs_core(I_target, iter); % ---------- 重建验证 ---------- E_out = fftshift(fft2(ifftshift(exp(1i * phase)))); I_rec = abs(E_out).^2; I_rec = I_rec / sum(I_rec(:)); figure; imshow(mat2gray(I_rec)); title('重建焦面强度');这段脚本里，lambda、z、dx1、N四个参数是绑定的，因为它们共同决定目标面的物理尺寸。按S-FFT的采样公式，目标面的网格间距是lambda * z / (N * dx1)，代入这组参数得到约0.104mm。这意味着每个点阵点之间的间距是0.104mm，整个目标面覆盖约26.6mm见方。你要把目标改成不同物理尺寸，不要只调I_target，要先算目标面间距，再反推N或dx1。这是DOE设计里最容易忽略的物理自洽问题。
入口脚本里I_target用“点阵”构造而不是直接用照片，是有意的。纯相位DOE对高细节目标的重建效果一般，但点阵这类稀疏目标，GS能在几十轮内收敛到不错的结果，适合先验证链路通不通。等你把GS跑顺了，再换复杂目标图不迟。
3.2 gs_core函数：两次FFT与两个幅值约束交替
核心迭代函数不长，但它就是一个纯相位DOE设计的中枢。读这个函数，关键是分清三个变量：phase是输出（我们要的DOE相位），F是目标面复振幅，E是DOE面复振幅。三个变量在迭代中轮流被更新，互相不混淆。
function phase = gs_core(I_target, iter) % GS算法核心：从目标强度分布反推DOE纯相位 % 输入： % I_target : 目标面光强分布，N*N的double矩阵 % iter : 迭代次数，建议50到200 % 输出： % phase : 设计的DOE相位分布，范围[-pi, pi] N = size(I_target, 1); A_target = sqrt(I_target); % 目标振幅约束 phase = 2 * pi * rand(N, N); % 随机相位起手，打破对称性 for k = 1:iter % ---- 前向传播：DOE面 -> 目标面 ---- E_in = ifftshift(exp(1i * phase)); % DOE面复振幅，幅值默认为1 F_out = fft2(E_in); F_out = fftshift(F_out); % ---- 目标面约束：替换振幅，保留相位 ---- F_new = A_target .* exp(1i * angle(F_out)); % ---- 逆传播：目标面 -> DOE面 ---- E_back = ifft2(ifftshift(F_new)); E_back = fftshift(E_back); % ---- DOE面约束：只取相位，振幅重置为1 ---- phase = angle(E_back); end end代码不长，但四个操作都踩在了要害上。前向传播里那句exp(1i * phase)，同时把DOE面振幅设为1，对应理想的平面波入射，这是纯相位元件的定义约束。目标面约束那句A_target .* exp(1i * angle(F_out))，是GS的灵魂：把传播结果的幅值用目标幅值替换，但把相位原样保留，等价于往“目标振幅集合”上做投影。逆传播后，DOE面的振幅约束同理——直接用angle取相位，把振幅全部丢掉，等价于往“纯相位集合”上投影。两处约束交替来，就是GS的全部秘密。
两个shift函数值得单独说。MATLAB里fft2默认把零频放在矩阵左上角，但物理上我们习惯零频在中心。ifftshift把中心移到左上角再给fft2，处理完用fftshift把零频挪回中心。这四行必须成对出现，少一个都不行。如果不做shift，重建图案在方位上可能旋转180度或出现中心对称翻转，但强度数值本身不受影响——初看似乎没错位，细看排列规律就发现问题。工程里常见的排查项就在这里。
3.3 重建评估：不只是肉眼看图，用误差曲线判断收敛
GS迭代不是跑完就算了，工程里留了误差评估的钩子。建议你在gs_core循环里维护一条误差记录，每轮计算重建面强度与目标强度的归一化误差，跑完直接plot看下降趋势。
% 在迭代循环内追加，每轮记录归一化均方误差 I_out = abs(F_out).^2; I_out = I_out / sum(I_out(:)); err(k) = sum(abs(I_out(:) - I_target(:)).^2) / sum(I_target(:).^2);迭代结束后在命令行敲一句plot(err)，就能看到GS的典型收敛形态：前20轮误差快速下降，后面几十轮几乎是一条缓坡，甚至偶尔反弹。这不是实现错了，是GS交替投影的特性——越接近约束交集，收敛越慢。判断是否继续跑，看最后一段曲线斜率是否趋平，如果趋势已经不动了，再增加迭代次数意义不大。此时应当改算法策略，比如切到加权GS，而不是干等。
我这里还有一个自己的习惯：误差曲线要看，但更要看重建图上“不该亮的地方”。均匀背景上如果出现少量散斑亮点，这是纯相位重建的正常现象，不是bug；如果背景整体发亮且呈条纹状，那多半是传播模型或采样参数有问题，优先查dx1与z的匹配关系。
4. 使用避坑：跑这套DOE设计时最常见的四个翻车现场
4.1 现象到原因：直接抄代码时最容易踩的坑
第一坑：重建图是旋转180度后的镜像。现象是对应的点阵位置先后顺序反了，目标是左三右二，出来是右三左二。原因是IFTA实现里fftshift和ifftshift用乱了，正向传播用fft2(fftshift(x))，逆传播又用ifft2(ifftshift(x))，前后没有保持一致的坐标系定义。解决方法是固定一套：输入用ifftshift，输出用fftshift，正反向都按这个规则来，不要混。最简单的自查办法是构造一个不对称目标图，比如“L”形点阵，跑一次就能看出方位是否正确。
第二坑：目标面强度一直收敛不动，误差曲线卡在0.8以上。现象是迭代几十轮，重建图还是一个模糊光斑，点阵没有成形。原因多半是目标图分辨率远高于DOE可表达的物理极限，比如直接用512x512的照片做目标，纯相位元件重建不了这么高频的信息。解决方法是先降目标图的分辨率，或者改用稀疏目标验证——这也是我在3.1里坚持用点阵的原因。验证GS实现本身没问题之后，再上复杂目标不迟。
第三坑：重建图能量分布不对，点阵整体发暗，背景大面积亮斑。现象是全局能量守恒了，但能量明显往中心聚拢。原因一方面是目标图本身包含了零频附近的直流分量，GS会优先重建低频。另一方面是I_target没有做能量归一化，目标振幅尺度与DOE面约束不匹配，导致约束投影的权重失衡。解决方法是先对I_target做sum归一化，让目标和输入的物理能量在同一量级；如果中心亮点还是太强，在目标图上把中心附近像素置零，也就是在目标强度里加一个“暗场”约束。
第四坑：512x512网格100轮迭代，跑出结果像蜗牛爬。现象是循环里每轮有三个矩阵运算，理论上不该这么慢，实际却要等很久。原因大概率是循环体内无意写了逐像素操作，比如用foor循环遍历矩阵给每个像素设初始值。MATLAB里这种写法慢得离谱。解决方法是把初始化改成phase = 2pirand(N,N)，目标图构造用矩阵切片一次性赋值，绝不用双for循环。要是还想再快，把gs_core里的fft2换成GPU版的gpuArray调用，个人场景下收益明显。
第五坑：主程序能跑，但figure窗口里什么都看不到，一片全白。现象是imshow(I_rec)显示全白，或者全是黑的。原因多数是矩阵类型问题：I_rec是double且数值范围在0到1e-6，imshow会把double图默认映射到0到1，数值太小导致全黑；反之没归一化直接显示，又容易全白。解决方法是显示前统一接一句mat2gray(I_rec)转成0~1灰度，不要直接接original未归一化的矩阵。这套工程里的脚本都带了mat2gray，你二次开发新脚本时别忘。
5. 二次开发思路：把标准GS改成你能用的DOE设计工具
5.1 加权GS：解决纯相位重建散斑问题
标准GS在目标面约束上做的事情是“无条件替换振幅”，太硬了，导致每次迭代都会留下散斑噪声。把约束换成带权重的软约束，就是WGS（加权GS）的核心。做法在约束那一段加一行权重：
F_new = (1 - w) * F_out + w * (A_target .* exp(1i * angle(F_out)));w取0.5到0.8之间，相当于每次只向目标振幅走半步，而不是一步到位。这会让收敛速度略微变慢，但散斑噪声明显减少，尤其适合目标图里同时有亮区和暗区的场景。
5.2 多平面与角谱法：摆脱S-FFT的距离限制
S-FFT把传播近似成一次FFT，条件是大菲涅尔数，但工程里很多时候DOE和目标面离得很近，此时采样间距那个公式会给出不合理的间距。换成角谱法就能把近场传播接进来，核心替换就是前向传播那两行，把传播算符写成乘积，具体传递函数是H = exp(1i * k * z)乘以sqrt(1 - (lambdafx)^2 - (lambdafy)^2)形式的项。改法不复杂，核心代码从fft2一行变成三个矩阵运算，多几句而已。
5.3 参数化封装：把脚本改成可批量扫描的函数
我在这个项目上的改造习惯，是把main_design脚本拆成参数入参、迭代过程、重建验证三块。参数入参变成函数参数，这样批量扫描z、dx1、N的组合就方便了。批量扫描对消费级应用场景很实用，比如你要对比不同传播距离下DOE的容差，脚本里写个嵌套循环，每组参数跑一次gs_core，记录均匀性指标，最后画一条曲线，选最优解。
这个工程跑顺了以后，我再看类似的迭代相位设计项目，都会先确认约束写法是硬约束还是软约束，再看采样参数是否物理自洽。GS本身不复杂，翻车几乎都翻在约束写法和坐标shift上。从那以后我每次设计DOE都强制走一遍流程：先跑稀疏目标验证链路，再画误差曲线确认收敛，最后用归一化灰度图检查重建结果。希望帮到你。
本文还有配套的精品资源，点击获取
无人机三维路径规划，输出做短路径图和适应度曲线.)
