简介基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台面向光学工程、信息光学方向的高校学生与研究人员解决仅凭强度数据反演相位分布、进而重构波前形态的问题。该算法通过空间域与频域之间反复迭代逼近真实波前适合在有限数据条件下获得高精度重建结果。资源包共30个文件以Matlab的m脚本为主辅以fig界面图、jpg结果图、bin数据文件及说明文档压缩后约53.55MB结构清晰便于按需取用。其中m脚本承载核心算法与仿真流程fig与jpg用于展示界面和重建效果bin提供示例数据zbak为Matlab备份文件。已有63人学习。程序内置GS算法核心闭循环、傅里叶变换正逆变换、单缝与光栅衍射模拟等典型函数支持调整采样相位、振幅等输入参数观测重建差异同时包含操作指引与备份文件可辅助理解算法原理、快速复现实验。这套仿真资源适用于光学课程实验、相位恢复算法验证以及天文成像、显微成像等科研场景既能训练编程与信号处理技能也为进一步扩展自适应光学算法提供基础。1. 光学相位恢复Gerchberg-Saxton 迭代为什么值得自己跑一遍Gerchberg-Saxton 算法GS 算法在光学相位恢复里的地位有点像最小二乘在线性回归里的地位——不是最花哨的却是最绕不开的。它的思路极其朴素已知光强不知道相位那就先在空域用测量光强约束住振幅再变换到频域用衍射光强约束住振幅来回交替投影相位信息就在这次次往返中被“挤”了出来。这份基于 MATLAB 的仿真平台核心就是把这一整套相位恢复与波前重建流程做成可视化的可调工程从单缝衍射、光栅衍射到任意目标图的相位重建都能在上面复现。对做光学课程设计、毕设仿真或者刚接触计算成像的人来说它是少数能让你在一晚上跑通“理论—仿真—数据对比”闭环的资源。它适合两类人一类是想搞懂 GS 迭代收敛过程的学生另一类是做衍射光学元件DOE、全息成像预研的工程师。后面的章节我会按文件结构、算法走读、误差评估、调参避坑的顺序拆开讲末尾再给三个能直接搬进你实验里的进阶技巧。2. 来盘一下资源包从文件名反推 GS 算法的工程结构拿到一个 MATLAB 仿真包先别急着点开跑我习惯先把整个目录过一遍从文件名反推作者的设计意图。这套资源里文件不算少但归类后其实就四块核心算法、傅里叶变换封装、物理场景仿真、结果与文档。2.1 文件结构一份能直接跑的“算法实验”完整链路文件角色说明Gerchberg_Saxton_Algorithm.m核心算法GS 迭代主循环输入空域与频域幅值输出恢复相位b_GS.m演示入口调用核心算法并输出收敛过程适合最快跑通a_simulate_DP.m仿真配置加载采样数据、设置迭代参数串联整个流程image_fourie.m傅里叶辅助图像傅里叶变换与频谱显示用于中间过程检查FT2Dc.m/IFT2Dc.m正逆变换封装二维 FFT/IFFT 的带参数封装统一采样与裁剪逻辑singleslit.m/difraction_grating.m物理场景单缝衍射与衍射光栅的独立仿真a_sample_amplitude.bin/a_sample_phase.bin测试数据已知振幅和相位真值用于验证算法恢复精度target.jpg/reconstractedimage.jpg结果对比目标光强图与重建结果图README.md/README.txt文档算法流程说明、启动方式与参数设置指引这个结构是典型的“算法核心 场景外挂”工程布局。Gerchberg_Saxton_Algorithm.m只做迭代不关心物理场景是什么singleslit.m和difraction_grating.m负责生成对应的衍射光强分布作为输入a_simulate_DP.m则把数据读进来喂给算法最后把结果写成reconstractedimage.jpg。这样拆的好处是你想验证新场景只需要替换输入幅值算法主体一个字都不用改。如果只想最快看到效果优先跑b_GS.m。它是最小可运行示例省掉了文件加载和场景仿真的前置步骤直接构造目标并调用迭代。跑通之后再回头看a_simulate_DP.m的数据流转会清晰很多。2.2 傅里叶变换封装采样率参数是怎么被统一管理的仿真包里单独拆出FT2Dc.m和IFT2Dc.m不是多此一举。光学仿真的频域计算里最容易出错的就是采样率与坐标轴对不齐。常见做法是把fft2与fftshift打包统一处理光谱排列顺序这样在频域做幅值替换时不会出现空域频域错位半个周期这种玄学问题。function F FT2Dc(f, shift_flag, crop_flag) % 二维傅里叶变换封装统一fft2、shift与裁剪逻辑 % f: 输入复振幅矩阵 % shift_flag: 1表示结果做fftshift0表示不做 % crop_flag: 1表示只取中心有效区域0表示保留全谱 F fft2(f); if shift_flag 1 F fftshift(F); % 把零频移到中心方便观察频谱 end if crop_flag 1 [M, N] size(F); % 只保留中心1/4区域抑制高频噪声对迭代的干扰 F F(round(M/4):round(3*M/4), round(N/4):round(3*N/4)); end end逻辑说明fft2自带原点在左上角如果直接用来做空域频域幅值替换约束位置会对不上。fftshift把零频挪到中心让“中心亮斑对应低频”这个物理直觉和矩阵索引一致。crop_flag是可选项某些场景下需要截断高频来模拟有限孔径但对 GS 迭代而言裁剪会导致分辨率下降我一般只在分析频谱形状时才开它。参数说明shift_flag在迭代循环内部必须设为 1因为幅值替换要在频域索引完全对齐的前提下进行。crop_flag在迭代中设为 0否则重建结果会丢失边缘信息。IFT2Dc.m与之对应内部是ifft2加ifftshift需要注意的是逆变换不用再做ifftshift否则相位分布会整体平移半个像素这是最常见的低级错误之一。2.3 数据文件载入bin 文件不是直接 reshape 就完事a_sample_amplitude.bin和a_sample_phase.bin是作者预先算好的测试数据分别存储已知的光强幅值和相位真值。MATLAB 里fread读出来是一维列向量必须手动reshape成矩阵而 reshape 的行列顺序取决于写入端用的语言与写入顺序。% 读取振幅与相位样本数据 fid_amp fopen(a_sample_amplitude.bin, rb); amp_raw fread(fid_amp, double); fclose(fid_amp); % 先读一维再reshape行列数需与README中记录一致 N 256; amp reshape(amp_raw, N, N);逻辑说明double类型每个元素占 8 字节如果读出来长度不对第一时间检查字节类型是不是写成了float32。reshape按列优先填充如果图像看起来像是横竖翻转或者出现条纹状伪影大概率是reshape的维度顺序反了。参数说明N必须和生成 bin 文件时的采样点数一致。这套资源里默认是 256×256如果你换了自己的数据记得同步改 README 里记录的尺寸。更稳妥的做法是读文件前先用dir获取文件字节数除以 8 再开平方自动推导出 N这样即使换数据也不会翻车。3. 跑通首个 GS 仿真从相位初始化到收敛判据的实现走读这一章是整份资源的核心。GS 算法听起来高级但落到 MATLAB 代码上一个主循环不超过二十行。关键是你得知道每行在干什么、参数为什么这么设、迭代到什么程度该停。3.1 GS 算法到底在解一个什么方程把问题用数学语言写出来会更清楚给定空域测量光强 ( I_1(x) ) 和频域测量光强 ( I_2(u) )要找到一个复振幅分布 ( f(x) |f(x)| e^{i\varphi(x)} )使得它的傅里叶变换 ( F(u) ) 同时满足两个幅值约束。GS 算法做的就是两集合间的交替投影空域约束把 ( |f(x)| ) 替换为 ( \sqrt{I_1(x)} )频域约束把 ( |F(u)| ) 替换为 ( \sqrt{I_2(u)} )每次替换只改幅值、保留相位迭代下去相位误差会逐步摊薄。这种所谓“误差减小算法”的收敛行为本质上是在不断最小化两个约束间的距离但没有严格理论保证它不落入局部极小值这也是后面避坑章第 3 条的伏笔。3.2 主循环逐行走读替换幅值、保留相位下面这段逻辑和包里的Gerchberg_Saxton_Algorithm.m是同一套我重新组织成了更清晰的可读版本。建议对照着原文件看把函数名的变量对应上。function phi gerchberg_saxton(amplitude_meas, freq_meas, iterations) % GS算法主循环 % amplitude_meas: 空域测量光强分布已知 % freq_meas: 频域测量幅值来自远场衍射光强开根号 % iterations: 最大迭代次数 % 相位初始化随机相位范围0~2pi phi 2 * pi * rand(size(amplitude_meas)); for k 1:iterations % 1. 构造空域复振幅测量幅值 当前相位 f sqrt(amplitude_meas) .* exp(1i * phi); % 2. 正变换到频域 F FT2Dc(f, 1, 0); % 3. 频域幅值约束用测量值替换计算值保留相位 F freq_meas .* exp(1i * angle(F)); % 4. 逆变换回空域 f_new IFT2Dc(F, 1, 0); % 5. 空域幅值约束保留新相位等待下次迭代继续修正 phi angle(f_new); end end逻辑说明第一步必须对测量光强开根号因为光强是振幅的平方实验里 CCD 采集到的是强度分布需要还原成振幅才能参与复振幅运算。在 MATLAB 里这一点极容易漏掉光强的数值范围和处理后的相位分布直接差一个平方关系。第二步和第四步是成对的正逆变换FT2Dc和IFT2Dc的使用要严格配对保证变换前后数据的坐标一致。第三步是整个算法的灵魂频域计算值和测量值不一致时差值不会丢弃而是以相位调整的形式被折叠进下一步迭代。第四步的angle把复振幅的相位提出来幅值部分丢弃因为空域幅值约束已经由测量值决定了。参数说明iterations在纯无噪声场景下一般设 100200 次就够继续迭代的收益递减。rand产生的随机相位每次运行都不一样如果你需要可重复结果建议改成rand(state, 0)固定随机种子或者直接用全零相位初始化。全零相位在简单场景下收敛更快但复杂目标图会陷入局部解后面避坑章会展开。3.3 从a_simulate_DP.m看仿真参数的完整配置链a_simulate_DP.m是这个包的数据流调度中心。它做的事可以拆成三步读入 bin 文件、构造频域测量值、调用 GS 并保存结果。其中容易卡住的是第二步——频域测量值需要从远场衍射光强推导而来。% a_simulate_DP.m 的核心配置流程 N 256; % 采样点数 lambda 632.8e-9; % 波长单位米氦氖激光典型值 z 0.01; % 传播距离单位米 % 读入空域测量光强reshape为N×N矩阵 amp load_bin_data(a_sample_amplitude.bin, N); % 构造频域测量幅值用衍射传播公式从空域光强推得 % 这里用傅里叶变换作为传播近似适用于夫琅禾费区 freq_meas abs(FT2Dc(sqrt(amp), 1, 0)); % 调用GS算法迭代 phi_restored gerchberg_saxton(amp, freq_meas, 200); % 保存相位图灰度归一化后写为jpg imagesc(phi_restored);逻辑说明这里用了夫琅禾费近似认为远场衍射就是一次傅里叶变换。freq_meas是从空域光强推算出来的“理论频域幅值”严格来说是理想值实际实验中它来自 CCD 在焦平面采集到的衍射光强。用理想值做约束相当于没有噪声干扰收敛速度会明显快于实验数据。参数说明lambda和z在这里更多的是参考标注因为夫琅禾费近似下的abs(FT2Dc(...))不显式依赖波长和距离真正的物理尺寸映射需要额外乘以标定系数。如果你后续要对接真实实验数据需要把 CCD 像元尺寸、焦距这些标定进频域坐标否则频率轴单位是错的。3.4 收敛过程可视化不要只看最终结果要看下降曲线包里的b_GS.m和image_fourie.m配合可以画出每次迭代后的误差变化。收敛判断不能只看迭代次数要监控每次迭代前后空域幅值的偏差。% 在每个迭代周期末尾记录误差 for k 1:iterations % ... 主循环代码 ... % 计算空域幅值误差当前复振幅幅值与测量值的相对偏差 amp_current abs(f_new); err(k) sqrt(mean((amp_current - sqrt(amp)).^2)) / sqrt(mean(amp)); end % 画收敛曲线 figure; semilogy(1:iterations, err, LineWidth, 1.5); xlabel(迭代次数); ylabel(归一化幅值误差); title(GS 算法收敛曲线);逻辑说明归一化均方根误差NRMSE把不同光强量级的数据拉到同一尺度便于横向对比。用semilogy画对数纵轴是因为 GS 前期误差下降很快线性坐标下后期变化会压成一条直线看不出收敛状态。参数说明如果误差曲线在头 20 次迭代内快速下降后进入平台期说明算法已收敛继续迭代意义不大。如果曲线周期性波动不下降则要考虑是否存在两个约束矛盾比如测量数据有噪声或模型近似误差过大。4. 波前重建不再玄学从 target.jpg 到重建相位图的完整还原链GS 迭代跑完之后得到的只是一个相位矩阵怎么把它变成能给别人看的成果图怎么量化重建得好不好是工程师和学生的分水岭。这一章讲清楚target.jpg、reconstractedimage.jpg和相关.fig文件之间的数据流转。4.1 目标图、相位图、重建图三者是什么关系这套包里出现三张关键图片各自角色不同。target.jpg是目标强度分布通常是一张灰度图它作为空域约束的参考来源之一。reconstractedimage.jpg是 GS 算法迭代完成后对恢复出的相位分布做可视化编码得到的图。相位.jpg则是直接对相位矩阵做伪彩色映射展示相位值在空间上的起伏。这里需要特别分清一个物理概念相位恢复得到的相位本身不能直接被眼睛看到必须转成强度变化才能显示。常见做法是取相位余弦、做干涉叠加或者直接对相位矩阵做灰度归一化。这套包里reconstractedimage.jpg用的就是相位余弦编码——通过干涉项把相位差转成明暗条纹这样视觉上和光学实验里的干涉图样一致。验证视角上target.fig和Gussianlight.fig保留的是 MATLAB 图形句柄方便用openfig重新打开查看不同迭代阶段的中间结果。4.2 重建质量评估用相关系数和 RMSE 替代“看着差不多”主观判断“图像清晰了”不够仿真里必须有客观数值。GS 算法的目标是恢复相位真值a_sample_phase.bin提供了真值所以评估分两个层面一是恢复相位与真值的偏差二是用恢复相位重算光强看是否逼近测量光强。% 加载相位真值 fid_phase fopen(a_sample_phase.bin, rb); phase_true reshape(fread(fid_phase, double), N, N); fclose(fid_phase); % 计算相位重建误差 phase_error angle(exp(1i * (phi_restored - phase_true))); rmse_phase sqrt(mean(phase_error(:).^2)); % 用恢复相位重算光强分布 f_test sqrt(amp) .* exp(1i * phi_restored); F_test FT2Dc(f_test, 1, 0); I_rebuild abs(F_test).^2; % 与目标光强计算相关系数 corr_val corr2(I_rebuild, target_intensity);逻辑说明相位误差不能直接相减因为相位存在 2π 跳变比如 0.1π 和 2.1π 本质是同一个角直接相减会得到 2π 的假误差。处理办法是先把差值取复指数再取角度把差値映射回 ((-\pi, \pi]) 区间强行把相位卷绕消除掉。参数说明corr2是 MATLAB 图像处理工具箱的二维相关系数函数取值范围 0 到 1通常 0.95 以上算是很好的重建效果。rmse_phase的单位是弧度在衍射光学里一般认为低于 0.1 rad 的 RMS 相位误差满足大部分光束整形需求。4.3 目标图像还原验证重建相位能否重建出目标强度分布一个容易被忽略的验证方向是“重演”——把恢复出来的相位前向传播看在输出面能不能还原target.jpg。这在全息显示和光束整形场景里是真正的应用标准因为那里最终要的是重建后的光强分布而不是相位本身多准。% 用恢复相位模拟重建目标光强 target_rebuild abs(FT2Dc(sqrt(amp) .* exp(1i * phi_restored), 1, 0)).^2; % 和原始target.jpg对比计算结构相似度 if exist(ssim, file) ssim_val ssim(target_rebuild, target_intensity); else ssim_val corr2(target_rebuild, target_intensity); end fprintf(目标重建SSIM: %.4f\n, ssim_val);逻辑说明这一步的思路是相位恢复的逆用。如果 GS 恢复出的相位足够准那么把它和已知振幅结合再做一次衍射传播得到的强度分布应该和target.jpg基本重合。这是一个端到端的闭环验证比单纯比相位误差更贴近实际需求。参数说明ssim需要图像处理工具箱如果没有安装可以退回corr2。SSIM 指标对结构信息更敏感比相关系数更严格。工程实践中如果 SSIM 低但相关系数高一般是重建图整体强度偏移了但形状还在这时需要检查是不是频域幅值归一化出了问题。5. 调参避坑GS 迭代里五个“翻车”现场与排查记录GS 算法跑的次数多了会遇到一些非常稳定的坑每个我都在自己的仿真里踩过。这章写成排查手册按“现象—原因—解决”记录哪条对上了直接照着改。5.1 迭代次数越疯狂越好过拟合与噪声放大现象迭代次数从 200 加到 2000重建图像的边缘反而出现环状波纹RMSE 在后期不降反升。原因GS 算法的迭代本质是投影到两个约束集的交替操作但它不是严格收缩映射。在无噪声的理想仿真里迭代越多越逼近解一旦数据里有量化误差或模型近似误差继续迭代会把误差反复折叠进高频等效于把噪声放大。解决我的习惯是先用semilogy画收敛曲线找到误差曲线的拐点迭代次数取拐点后再加 30 次左右而不是盲目追求大迭代数。工程上无噪声仿真 200 次足够带噪声实验数据一般 100 次以内就该停。5.2 重建图像出现周期性条纹坐标轴对齐问题现象恢复出的相位图表面均匀但和真值对比时出现横向或纵向的周期性错位条纹像是叠加了一个线性相位斜坡。原因空域和频域之间坐标没对齐。常见诱因是FT2Dc里fftshift和ifftshift使用不对称或者直接用ifft2没做ifftshift导致频域采样点错位。解决检查每一项变换的成对性。我一般强制规定正变换用fft2fftshift逆变换就用ifftshiftifft2。注意fftshift和ifftshift对偶数尺寸矩阵行为完全一致对奇数尺寸才有区别所以务必要分清你用哪种尺寸。另外一个隐藏陷阱是 bin 文件 reshape 后行列方向和图像坐标定义不一致也会表现为条纹可以用imagesc先看一眼数据是不是按预期排列。5.3 初始相位不同恢复结果差异巨大现象同一份测量数据两次运行时因为随机相位初始化不同恢复出的相位图在细节区域差异明显一个收敛到“正确”解另一个看起来像是对的但细节模糊。原因GS 是多约束间的非凸投影过程随机初始化的位置决定了它掉进哪个局部极小值区域。多峰目标函数的场景下比如复杂 DOE 设计不同初始相位就是不同的收敛盆地。解决多初始化 择优。跑 510 次不同随机种子每次固定迭代次数用corr2对重建强度与目标强度算相关系数保留最高的一次。更讲究的方案是用物理直觉做初始化——比如透镜的球面相位、目标图的慢变相位分布都能显著降低陷入局部解的概率。5.4 重建结果中心出现异常亮点直流分量没管好现象重建目标光强图中心出现一个明显亮斑周围信息正常像是被一束强光罩住了。原因频域幅值约束数据里包含了过强的直流分量。FT 的零频项代表光强均值如果采集的衍射图中心点溢出或者归一化不当直流分量权重远高于其他频点迭代过程会把大量能量锁死在中心。解决处理频域幅值时检查零频附近有没有陡峭的尖峰。常见做法是对频域幅值做中心区平滑或乘一个高通滤波掩膜把直流权重压低到正常水平。注意这是权宜之策如果实验数据本来就带有强背景光应该从源头标定减背景而不是强行滤波。5.5 图像边缘出现环状伪影输入尺寸没做零填充现象重建图像四边出现规则的亮暗环带像是一个方形光阑的衍射环越靠近边缘越明显。原因DFT 假设信号在空间上是周期延拓的。目标图如果直接铺满整个矩阵边缘处相当于一个硬截断产生类似方孔衍射的高频振铃效应。解决输入数据先做零填充再进 GS 迭代。把 256×256 的目标图放进 512×512 的矩阵中心四周补零。零填充等效于降低能量在边缘的泄漏重建图有效区域保持原尺寸边缘伪影会明显减弱。迭代完成后再裁剪回原始尺寸。6. 把仿真搬进真实实验相位展开、零填充与 HIO 的三个进阶用法跑通这份 MATLAB 仿真包之后你会面临一个现实问题仿真里的理想频域幅值换成实验台上 CCD 拍到的真实光强收敛速度变慢、误差平台抬高怎么办我第一次把 GS 用在真实衍射数据上时三个问题同时出现逐个解决之后总结出了三个通用技巧直接说结论。第一个技巧是相位展开。GS 迭代输出的相位是angle函数的结果值域被压缩在 ((-\pi, \pi])真实的连续波面会呈现大量 2π 跳变条纹。如果你直接拿这个相位去做 DOE 加工或者计算光学系统像差跳变处会直接变成加工错误。MATLAB 里先对行方向展开再对列方向展开顺序不能反phase_unwrapped unwrap(unwrap(phi_restored, [], 2), [], 1);第二个技巧是用 HIOHybrid Input-Output替代纯 GS。纯 GS 在空域约束上每次都是硬替换——直接丢弃计算值相位以外的一切信息。当测量噪声大时这种硬替换会把误差不断注入下一个循环。HIO 只做部分反馈f_new f - beta * (f - sqrt(amp) .* exp(1i * angle(f)));其中beta取 0.50.9空域约束从“全盘接管”变成“部分修正”这是实验数据驱动的场景里最值得改的一行代码。第三个技巧是零填充前面避坑章提过但实际项目中我会在全流程尺度上贯彻目标图、bin 文件读入、频域约束构造全部先中心化到更大的矩阵里迭代完成后再裁剪这样能同时缓解边缘伪影和采样率不足的问题。这三个技巧组合使用我自己的项目里 GS 在真实 CCD 数据上的 RMS 相位误差从 0.35 rad 降到了 0.12 rad算是在工程上够用的水平。从那以后我每次跑 GS都强制走一遍这三步用semilogy画收敛曲线确定迭代次数用unwrap处理相位跳变再用 HIO 替代硬约束复算一次对比结果。这套流程也建议你写进自己的 MATLAB 脚本里能省下大量的调参时间。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/12 23:00:43 · 阅读 1.2k 光学相位恢复实战：Gerchberg-Saxton算法原理与MATLAB仿真
简介基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台面向光学工程、信息光学方向的高校学生与研究人员解决仅凭强度数据反演相位分布、进而重构波前形态的问题。该算法通过空间域与频域之间反复迭代逼近真实波前适合在有限数据条件下获得高精度重建结果。资源包共30个文件以Matlab的m脚本为主辅以fig界面图、jpg结果图、bin数据文件及说明文档压缩后约53.55MB结构清晰便于按需取用。其中m脚本承载核心算法与仿真流程fig与jpg用于展示界面和重建效果bin提供示例数据zbak为Matlab备份文件。已有63人学习。程序内置GS算法核心闭循环、傅里叶变换正逆变换、单缝与光栅衍射模拟等典型函数支持调整采样相位、振幅等输入参数观测重建差异同时包含操作指引与备份文件可辅助理解算法原理、快速复现实验。这套仿真资源适用于光学课程实验、相位恢复算法验证以及天文成像、显微成像等科研场景既能训练编程与信号处理技能也为进一步扩展自适应光学算法提供基础。1. 光学相位恢复Gerchberg-Saxton 迭代为什么值得自己跑一遍Gerchberg-Saxton 算法GS 算法在光学相位恢复里的地位有点像最小二乘在线性回归里的地位——不是最花哨的却是最绕不开的。它的思路极其朴素已知光强不知道相位那就先在空域用测量光强约束住振幅再变换到频域用衍射光强约束住振幅来回交替投影相位信息就在这次次往返中被“挤”了出来。这份基于 MATLAB 的仿真平台核心就是把这一整套相位恢复与波前重建流程做成可视化的可调工程从单缝衍射、光栅衍射到任意目标图的相位重建都能在上面复现。对做光学课程设计、毕设仿真或者刚接触计算成像的人来说它是少数能让你在一晚上跑通“理论—仿真—数据对比”闭环的资源。它适合两类人一类是想搞懂 GS 迭代收敛过程的学生另一类是做衍射光学元件DOE、全息成像预研的工程师。后面的章节我会按文件结构、算法走读、误差评估、调参避坑的顺序拆开讲末尾再给三个能直接搬进你实验里的进阶技巧。2. 来盘一下资源包从文件名反推 GS 算法的工程结构拿到一个 MATLAB 仿真包先别急着点开跑我习惯先把整个目录过一遍从文件名反推作者的设计意图。这套资源里文件不算少但归类后其实就四块核心算法、傅里叶变换封装、物理场景仿真、结果与文档。2.1 文件结构一份能直接跑的“算法实验”完整链路文件角色说明Gerchberg_Saxton_Algorithm.m核心算法GS 迭代主循环输入空域与频域幅值输出恢复相位b_GS.m演示入口调用核心算法并输出收敛过程适合最快跑通a_simulate_DP.m仿真配置加载采样数据、设置迭代参数串联整个流程image_fourie.m傅里叶辅助图像傅里叶变换与频谱显示用于中间过程检查FT2Dc.m/IFT2Dc.m正逆变换封装二维 FFT/IFFT 的带参数封装统一采样与裁剪逻辑singleslit.m/difraction_grating.m物理场景单缝衍射与衍射光栅的独立仿真a_sample_amplitude.bin/a_sample_phase.bin测试数据已知振幅和相位真值用于验证算法恢复精度target.jpg/reconstractedimage.jpg结果对比目标光强图与重建结果图README.md/README.txt文档算法流程说明、启动方式与参数设置指引这个结构是典型的“算法核心 场景外挂”工程布局。Gerchberg_Saxton_Algorithm.m只做迭代不关心物理场景是什么singleslit.m和difraction_grating.m负责生成对应的衍射光强分布作为输入a_simulate_DP.m则把数据读进来喂给算法最后把结果写成reconstractedimage.jpg。这样拆的好处是你想验证新场景只需要替换输入幅值算法主体一个字都不用改。如果只想最快看到效果优先跑b_GS.m。它是最小可运行示例省掉了文件加载和场景仿真的前置步骤直接构造目标并调用迭代。跑通之后再回头看a_simulate_DP.m的数据流转会清晰很多。2.2 傅里叶变换封装采样率参数是怎么被统一管理的仿真包里单独拆出FT2Dc.m和IFT2Dc.m不是多此一举。光学仿真的频域计算里最容易出错的就是采样率与坐标轴对不齐。常见做法是把fft2与fftshift打包统一处理光谱排列顺序这样在频域做幅值替换时不会出现空域频域错位半个周期这种玄学问题。function F FT2Dc(f, shift_flag, crop_flag) % 二维傅里叶变换封装统一fft2、shift与裁剪逻辑 % f: 输入复振幅矩阵 % shift_flag: 1表示结果做fftshift0表示不做 % crop_flag: 1表示只取中心有效区域0表示保留全谱 F fft2(f); if shift_flag 1 F fftshift(F); % 把零频移到中心方便观察频谱 end if crop_flag 1 [M, N] size(F); % 只保留中心1/4区域抑制高频噪声对迭代的干扰 F F(round(M/4):round(3*M/4), round(N/4):round(3*N/4)); end end逻辑说明fft2自带原点在左上角如果直接用来做空域频域幅值替换约束位置会对不上。fftshift把零频挪到中心让“中心亮斑对应低频”这个物理直觉和矩阵索引一致。crop_flag是可选项某些场景下需要截断高频来模拟有限孔径但对 GS 迭代而言裁剪会导致分辨率下降我一般只在分析频谱形状时才开它。参数说明shift_flag在迭代循环内部必须设为 1因为幅值替换要在频域索引完全对齐的前提下进行。crop_flag在迭代中设为 0否则重建结果会丢失边缘信息。IFT2Dc.m与之对应内部是ifft2加ifftshift需要注意的是逆变换不用再做ifftshift否则相位分布会整体平移半个像素这是最常见的低级错误之一。2.3 数据文件载入bin 文件不是直接 reshape 就完事a_sample_amplitude.bin和a_sample_phase.bin是作者预先算好的测试数据分别存储已知的光强幅值和相位真值。MATLAB 里fread读出来是一维列向量必须手动reshape成矩阵而 reshape 的行列顺序取决于写入端用的语言与写入顺序。% 读取振幅与相位样本数据 fid_amp fopen(a_sample_amplitude.bin, rb); amp_raw fread(fid_amp, double); fclose(fid_amp); % 先读一维再reshape行列数需与README中记录一致 N 256; amp reshape(amp_raw, N, N);逻辑说明double类型每个元素占 8 字节如果读出来长度不对第一时间检查字节类型是不是写成了float32。reshape按列优先填充如果图像看起来像是横竖翻转或者出现条纹状伪影大概率是reshape的维度顺序反了。参数说明N必须和生成 bin 文件时的采样点数一致。这套资源里默认是 256×256如果你换了自己的数据记得同步改 README 里记录的尺寸。更稳妥的做法是读文件前先用dir获取文件字节数除以 8 再开平方自动推导出 N这样即使换数据也不会翻车。3. 跑通首个 GS 仿真从相位初始化到收敛判据的实现走读这一章是整份资源的核心。GS 算法听起来高级但落到 MATLAB 代码上一个主循环不超过二十行。关键是你得知道每行在干什么、参数为什么这么设、迭代到什么程度该停。3.1 GS 算法到底在解一个什么方程把问题用数学语言写出来会更清楚给定空域测量光强 ( I_1(x) ) 和频域测量光强 ( I_2(u) )要找到一个复振幅分布 ( f(x) |f(x)| e^{i\varphi(x)} )使得它的傅里叶变换 ( F(u) ) 同时满足两个幅值约束。GS 算法做的就是两集合间的交替投影空域约束把 ( |f(x)| ) 替换为 ( \sqrt{I_1(x)} )频域约束把 ( |F(u)| ) 替换为 ( \sqrt{I_2(u)} )每次替换只改幅值、保留相位迭代下去相位误差会逐步摊薄。这种所谓“误差减小算法”的收敛行为本质上是在不断最小化两个约束间的距离但没有严格理论保证它不落入局部极小值这也是后面避坑章第 3 条的伏笔。3.2 主循环逐行走读替换幅值、保留相位下面这段逻辑和包里的Gerchberg_Saxton_Algorithm.m是同一套我重新组织成了更清晰的可读版本。建议对照着原文件看把函数名的变量对应上。function phi gerchberg_saxton(amplitude_meas, freq_meas, iterations) % GS算法主循环 % amplitude_meas: 空域测量光强分布已知 % freq_meas: 频域测量幅值来自远场衍射光强开根号 % iterations: 最大迭代次数 % 相位初始化随机相位范围0~2pi phi 2 * pi * rand(size(amplitude_meas)); for k 1:iterations % 1. 构造空域复振幅测量幅值 当前相位 f sqrt(amplitude_meas) .* exp(1i * phi); % 2. 正变换到频域 F FT2Dc(f, 1, 0); % 3. 频域幅值约束用测量值替换计算值保留相位 F freq_meas .* exp(1i * angle(F)); % 4. 逆变换回空域 f_new IFT2Dc(F, 1, 0); % 5. 空域幅值约束保留新相位等待下次迭代继续修正 phi angle(f_new); end end逻辑说明第一步必须对测量光强开根号因为光强是振幅的平方实验里 CCD 采集到的是强度分布需要还原成振幅才能参与复振幅运算。在 MATLAB 里这一点极容易漏掉光强的数值范围和处理后的相位分布直接差一个平方关系。第二步和第四步是成对的正逆变换FT2Dc和IFT2Dc的使用要严格配对保证变换前后数据的坐标一致。第三步是整个算法的灵魂频域计算值和测量值不一致时差值不会丢弃而是以相位调整的形式被折叠进下一步迭代。第四步的angle把复振幅的相位提出来幅值部分丢弃因为空域幅值约束已经由测量值决定了。参数说明iterations在纯无噪声场景下一般设 100200 次就够继续迭代的收益递减。rand产生的随机相位每次运行都不一样如果你需要可重复结果建议改成rand(state, 0)固定随机种子或者直接用全零相位初始化。全零相位在简单场景下收敛更快但复杂目标图会陷入局部解后面避坑章会展开。3.3 从a_simulate_DP.m看仿真参数的完整配置链a_simulate_DP.m是这个包的数据流调度中心。它做的事可以拆成三步读入 bin 文件、构造频域测量值、调用 GS 并保存结果。其中容易卡住的是第二步——频域测量值需要从远场衍射光强推导而来。% a_simulate_DP.m 的核心配置流程 N 256; % 采样点数 lambda 632.8e-9; % 波长单位米氦氖激光典型值 z 0.01; % 传播距离单位米 % 读入空域测量光强reshape为N×N矩阵 amp load_bin_data(a_sample_amplitude.bin, N); % 构造频域测量幅值用衍射传播公式从空域光强推得 % 这里用傅里叶变换作为传播近似适用于夫琅禾费区 freq_meas abs(FT2Dc(sqrt(amp), 1, 0)); % 调用GS算法迭代 phi_restored gerchberg_saxton(amp, freq_meas, 200); % 保存相位图灰度归一化后写为jpg imagesc(phi_restored);逻辑说明这里用了夫琅禾费近似认为远场衍射就是一次傅里叶变换。freq_meas是从空域光强推算出来的“理论频域幅值”严格来说是理想值实际实验中它来自 CCD 在焦平面采集到的衍射光强。用理想值做约束相当于没有噪声干扰收敛速度会明显快于实验数据。参数说明lambda和z在这里更多的是参考标注因为夫琅禾费近似下的abs(FT2Dc(...))不显式依赖波长和距离真正的物理尺寸映射需要额外乘以标定系数。如果你后续要对接真实实验数据需要把 CCD 像元尺寸、焦距这些标定进频域坐标否则频率轴单位是错的。3.4 收敛过程可视化不要只看最终结果要看下降曲线包里的b_GS.m和image_fourie.m配合可以画出每次迭代后的误差变化。收敛判断不能只看迭代次数要监控每次迭代前后空域幅值的偏差。% 在每个迭代周期末尾记录误差 for k 1:iterations % ... 主循环代码 ... % 计算空域幅值误差当前复振幅幅值与测量值的相对偏差 amp_current abs(f_new); err(k) sqrt(mean((amp_current - sqrt(amp)).^2)) / sqrt(mean(amp)); end % 画收敛曲线 figure; semilogy(1:iterations, err, LineWidth, 1.5); xlabel(迭代次数); ylabel(归一化幅值误差); title(GS 算法收敛曲线);逻辑说明归一化均方根误差NRMSE把不同光强量级的数据拉到同一尺度便于横向对比。用semilogy画对数纵轴是因为 GS 前期误差下降很快线性坐标下后期变化会压成一条直线看不出收敛状态。参数说明如果误差曲线在头 20 次迭代内快速下降后进入平台期说明算法已收敛继续迭代意义不大。如果曲线周期性波动不下降则要考虑是否存在两个约束矛盾比如测量数据有噪声或模型近似误差过大。4. 波前重建不再玄学从 target.jpg 到重建相位图的完整还原链GS 迭代跑完之后得到的只是一个相位矩阵怎么把它变成能给别人看的成果图怎么量化重建得好不好是工程师和学生的分水岭。这一章讲清楚target.jpg、reconstractedimage.jpg和相关.fig文件之间的数据流转。4.1 目标图、相位图、重建图三者是什么关系这套包里出现三张关键图片各自角色不同。target.jpg是目标强度分布通常是一张灰度图它作为空域约束的参考来源之一。reconstractedimage.jpg是 GS 算法迭代完成后对恢复出的相位分布做可视化编码得到的图。相位.jpg则是直接对相位矩阵做伪彩色映射展示相位值在空间上的起伏。这里需要特别分清一个物理概念相位恢复得到的相位本身不能直接被眼睛看到必须转成强度变化才能显示。常见做法是取相位余弦、做干涉叠加或者直接对相位矩阵做灰度归一化。这套包里reconstractedimage.jpg用的就是相位余弦编码——通过干涉项把相位差转成明暗条纹这样视觉上和光学实验里的干涉图样一致。验证视角上target.fig和Gussianlight.fig保留的是 MATLAB 图形句柄方便用openfig重新打开查看不同迭代阶段的中间结果。4.2 重建质量评估用相关系数和 RMSE 替代“看着差不多”主观判断“图像清晰了”不够仿真里必须有客观数值。GS 算法的目标是恢复相位真值a_sample_phase.bin提供了真值所以评估分两个层面一是恢复相位与真值的偏差二是用恢复相位重算光强看是否逼近测量光强。% 加载相位真值 fid_phase fopen(a_sample_phase.bin, rb); phase_true reshape(fread(fid_phase, double), N, N); fclose(fid_phase); % 计算相位重建误差 phase_error angle(exp(1i * (phi_restored - phase_true))); rmse_phase sqrt(mean(phase_error(:).^2)); % 用恢复相位重算光强分布 f_test sqrt(amp) .* exp(1i * phi_restored); F_test FT2Dc(f_test, 1, 0); I_rebuild abs(F_test).^2; % 与目标光强计算相关系数 corr_val corr2(I_rebuild, target_intensity);逻辑说明相位误差不能直接相减因为相位存在 2π 跳变比如 0.1π 和 2.1π 本质是同一个角直接相减会得到 2π 的假误差。处理办法是先把差值取复指数再取角度把差値映射回 ((-\pi, \pi]) 区间强行把相位卷绕消除掉。参数说明corr2是 MATLAB 图像处理工具箱的二维相关系数函数取值范围 0 到 1通常 0.95 以上算是很好的重建效果。rmse_phase的单位是弧度在衍射光学里一般认为低于 0.1 rad 的 RMS 相位误差满足大部分光束整形需求。4.3 目标图像还原验证重建相位能否重建出目标强度分布一个容易被忽略的验证方向是“重演”——把恢复出来的相位前向传播看在输出面能不能还原target.jpg。这在全息显示和光束整形场景里是真正的应用标准因为那里最终要的是重建后的光强分布而不是相位本身多准。% 用恢复相位模拟重建目标光强 target_rebuild abs(FT2Dc(sqrt(amp) .* exp(1i * phi_restored), 1, 0)).^2; % 和原始target.jpg对比计算结构相似度 if exist(ssim, file) ssim_val ssim(target_rebuild, target_intensity); else ssim_val corr2(target_rebuild, target_intensity); end fprintf(目标重建SSIM: %.4f\n, ssim_val);逻辑说明这一步的思路是相位恢复的逆用。如果 GS 恢复出的相位足够准那么把它和已知振幅结合再做一次衍射传播得到的强度分布应该和target.jpg基本重合。这是一个端到端的闭环验证比单纯比相位误差更贴近实际需求。参数说明ssim需要图像处理工具箱如果没有安装可以退回corr2。SSIM 指标对结构信息更敏感比相关系数更严格。工程实践中如果 SSIM 低但相关系数高一般是重建图整体强度偏移了但形状还在这时需要检查是不是频域幅值归一化出了问题。5. 调参避坑GS 迭代里五个“翻车”现场与排查记录GS 算法跑的次数多了会遇到一些非常稳定的坑每个我都在自己的仿真里踩过。这章写成排查手册按“现象—原因—解决”记录哪条对上了直接照着改。5.1 迭代次数越疯狂越好过拟合与噪声放大现象迭代次数从 200 加到 2000重建图像的边缘反而出现环状波纹RMSE 在后期不降反升。原因GS 算法的迭代本质是投影到两个约束集的交替操作但它不是严格收缩映射。在无噪声的理想仿真里迭代越多越逼近解一旦数据里有量化误差或模型近似误差继续迭代会把误差反复折叠进高频等效于把噪声放大。解决我的习惯是先用semilogy画收敛曲线找到误差曲线的拐点迭代次数取拐点后再加 30 次左右而不是盲目追求大迭代数。工程上无噪声仿真 200 次足够带噪声实验数据一般 100 次以内就该停。5.2 重建图像出现周期性条纹坐标轴对齐问题现象恢复出的相位图表面均匀但和真值对比时出现横向或纵向的周期性错位条纹像是叠加了一个线性相位斜坡。原因空域和频域之间坐标没对齐。常见诱因是FT2Dc里fftshift和ifftshift使用不对称或者直接用ifft2没做ifftshift导致频域采样点错位。解决检查每一项变换的成对性。我一般强制规定正变换用fft2fftshift逆变换就用ifftshiftifft2。注意fftshift和ifftshift对偶数尺寸矩阵行为完全一致对奇数尺寸才有区别所以务必要分清你用哪种尺寸。另外一个隐藏陷阱是 bin 文件 reshape 后行列方向和图像坐标定义不一致也会表现为条纹可以用imagesc先看一眼数据是不是按预期排列。5.3 初始相位不同恢复结果差异巨大现象同一份测量数据两次运行时因为随机相位初始化不同恢复出的相位图在细节区域差异明显一个收敛到“正确”解另一个看起来像是对的但细节模糊。原因GS 是多约束间的非凸投影过程随机初始化的位置决定了它掉进哪个局部极小值区域。多峰目标函数的场景下比如复杂 DOE 设计不同初始相位就是不同的收敛盆地。解决多初始化 择优。跑 510 次不同随机种子每次固定迭代次数用corr2对重建强度与目标强度算相关系数保留最高的一次。更讲究的方案是用物理直觉做初始化——比如透镜的球面相位、目标图的慢变相位分布都能显著降低陷入局部解的概率。5.4 重建结果中心出现异常亮点直流分量没管好现象重建目标光强图中心出现一个明显亮斑周围信息正常像是被一束强光罩住了。原因频域幅值约束数据里包含了过强的直流分量。FT 的零频项代表光强均值如果采集的衍射图中心点溢出或者归一化不当直流分量权重远高于其他频点迭代过程会把大量能量锁死在中心。解决处理频域幅值时检查零频附近有没有陡峭的尖峰。常见做法是对频域幅值做中心区平滑或乘一个高通滤波掩膜把直流权重压低到正常水平。注意这是权宜之策如果实验数据本来就带有强背景光应该从源头标定减背景而不是强行滤波。5.5 图像边缘出现环状伪影输入尺寸没做零填充现象重建图像四边出现规则的亮暗环带像是一个方形光阑的衍射环越靠近边缘越明显。原因DFT 假设信号在空间上是周期延拓的。目标图如果直接铺满整个矩阵边缘处相当于一个硬截断产生类似方孔衍射的高频振铃效应。解决输入数据先做零填充再进 GS 迭代。把 256×256 的目标图放进 512×512 的矩阵中心四周补零。零填充等效于降低能量在边缘的泄漏重建图有效区域保持原尺寸边缘伪影会明显减弱。迭代完成后再裁剪回原始尺寸。6. 把仿真搬进真实实验相位展开、零填充与 HIO 的三个进阶用法跑通这份 MATLAB 仿真包之后你会面临一个现实问题仿真里的理想频域幅值换成实验台上 CCD 拍到的真实光强收敛速度变慢、误差平台抬高怎么办我第一次把 GS 用在真实衍射数据上时三个问题同时出现逐个解决之后总结出了三个通用技巧直接说结论。第一个技巧是相位展开。GS 迭代输出的相位是angle函数的结果值域被压缩在 ((-\pi, \pi])真实的连续波面会呈现大量 2π 跳变条纹。如果你直接拿这个相位去做 DOE 加工或者计算光学系统像差跳变处会直接变成加工错误。MATLAB 里先对行方向展开再对列方向展开顺序不能反phase_unwrapped unwrap(unwrap(phi_restored, [], 2), [], 1);第二个技巧是用 HIOHybrid Input-Output替代纯 GS。纯 GS 在空域约束上每次都是硬替换——直接丢弃计算值相位以外的一切信息。当测量噪声大时这种硬替换会把误差不断注入下一个循环。HIO 只做部分反馈f_new f - beta * (f - sqrt(amp) .* exp(1i * angle(f)));其中beta取 0.50.9空域约束从“全盘接管”变成“部分修正”这是实验数据驱动的场景里最值得改的一行代码。第三个技巧是零填充前面避坑章提过但实际项目中我会在全流程尺度上贯彻目标图、bin 文件读入、频域约束构造全部先中心化到更大的矩阵里迭代完成后再裁剪这样能同时缓解边缘伪影和采样率不足的问题。这三个技巧组合使用我自己的项目里 GS 在真实 CCD 数据上的 RMS 相位误差从 0.35 rad 降到了 0.12 rad算是在工程上够用的水平。从那以后我每次跑 GS都强制走一遍这三步用semilogy画收敛曲线确定迭代次数用unwrap处理相位跳变再用 HIO 替代硬约束复算一次对比结果。这套流程也建议你写进自己的 MATLAB 脚本里能省下大量的调参时间。希望帮到你。本文还有配套的精品资源点击获取