尧图网站建设 
简介基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台面向光学工程、信息光学及图像处理方向的学习者与研究者以光强数据反演相位分布为切入点系统展示了从相位恢复到波前重构的完整过程。包内共30个文件主要包括10个Matlab脚本.m、5个图形文件.fig、5张结果图片.jpg另有备份文件、数据文件及txt/md说明文档压缩包整体约53.55MB目录结构清晰便于按模块查阅。目前已有62人学习了该资源适合个人自学、课程实践或科研前期参考。通过修改输入参数或调整算法结构用户可深入分析GS算法在不同条件下的系统响应配合说明文档与示例图片能够快速掌握算法原理与实现细节并借助样例数据进行结果验证进一步提升光学仿真和算法开发的实践能力。1. 相位恢复问题的病态性Gerchberg-Saxton 迭代为何有效CCD 只记录强度波前相位被丢弃这是光学检测和成像里最难绕过的一道坎。Gerchberg-Saxton 迭代不依赖额外测量设备只靠物面振幅和频谱振幅这两个约束来回傅里叶变换就能把相位逼出来。手上这套 Matlab 仿真工程正好把这条链路拆开Gerchberg_Saxton_Algorithm.m是主循环FT2Dc.m、IFT2Dc.m负责正逆变换target.jpg、DOE.jpg提供目标图和衍射元件singleslit.m、difraction_grating.m则把单缝、光栅这类经典物面约束做成可切换场景。适合刚开始接触相位恢复的人跑通原理也适合做过 GS 但想研究收敛边界的人对照源码调参数。2. 傅里叶光学基础与频谱约束FT2Dc 的 fftshift 陷阱GS 的核心是空域和频域之间的交替投影但 Matlab 里的fft2和光学里的连续傅里叶变换并不一致尤其是坐标原点。这个平台把FT2Dc.m和IFT2Dc.m独立拆出来就是提示不要忽略这一步写错一个fftshift迭代出来的相位图会整体偏一个常数梯度看起来像倾斜波前实际上只是坐标绕乱了。2.1 正/逆傅里叶变换怎么和采样坐标对齐夫琅禾费衍射对应的连续变换是[ U(f_x, f_y)\iint u(x,y)\exp[-j2\pi(f_x x f_y y)]dxdy ]但离散 FFT 默认把零频放在矩阵的(1,1)位置而光学仿真中零频对应光轴中心。于是多数仿真代码会先对输入矩阵做ifftshift让阵列中心的光学原点移动到 FFT 约定的原点算完再fftshift把零频挪回中心。我一般把这类变换封装成独立函数避免在 GS 主循环里反复写错function F FT2Dc(u) % 二维光学傅里叶变换输入 u 的矩阵中心对应光轴原点 % 输出 F 的矩阵中心对应零频 F fftshift(fft2(ifftshift(u))); end function u IFT2Dc(F) % 对应的逆变换FT2Dc 与 IFT2Dc 可逆 u fftshift(ifft2(ifftshift(F))); end这段代码里ifftshift的作用是把中心原点挪到 FFT 的左上角原点fftshift负责把算完的零频移到矩阵中心。ifftshift和fftshift在偶数尺寸时行为相同在奇数尺寸时是严格互逆的所以不要用同一个函数替代另一个。真正验证这对函数是否写对的方法很简单对一个全 1 矩阵做FT2Dc理想频谱应该只有一个中心主峰直流分量集中在矩阵中央而不是四角。这部分在工程里容易被忽略因为单独画频谱时即使坐标偏了abs的图形看起来也只是整体平移。可一旦进入 GS 迭代物面约束和频谱约束会交替使用如果正逆变换的坐标约定不自洽相位恢复会收敛到一组看起来平滑但物理上完全错误的解。平台里每个.m文件都调用FT2Dc.m说明作者希望所有模块共享同一套坐标约定这也是光学仿真代码最值得先看的部分。函数在 GS 中的作用常见误用fft2/ifft2正逆傅里叶变换忘记先处理原点位置fftshift零频移到矩阵中心奇数尺寸时当成ifftshift用ifftshift中心原点移回 FFT 原点偶尺寸时不报错奇尺寸时错位abs取振幅约束对复数场直接取实部angle取相位约束相位跳变处需要包裹处理2.2 振幅约束与相位自由度的关系相位恢复问题的困难在于只有两个振幅约束可用来反推相位物面振幅 (A(x,y)) 和频谱振幅 (B(f_x,f_y))。它们之间的映射由傅里叶变换连接但相位本身不直接进入测量值。由于相位有大量自由度直接求逆是不适定的GS 算法的思路是先随机给一个初始相位然后反复执行两个投影在空域保持已知振幅、替换相位在频域同样保持已知振幅、替换相位。这里要注意物面的振幅约束并不一定等于目标图本身。这个工程把target.jpg当作已知振幅约束是一种有监督演示方便你比对重建结果和原图真正用于散射成像相位恢复时物面约束往往是照明光斑形状或一个支撑域掩模比如圆孔、狭缝、矩形区域。支撑域给得越大未知自由度越多GS 越容易陷入局部极小给得太小又会截断有效信号。工程里的singleslit.m、difraction_grating.m就是在构造不同类型的支撑域我建议你把这些文件里的 mask 改成圆形或环形再跑一轮观察收敛速度的变化。2.3 为什么选 GS 而不是直接求逆直接求逆需要同时知道复振幅而实验只能拿到强度。替代方案里有基于传输方程的 TIE、基于菲涅耳迭代的混合输入输出算法但 GS 的优势在于实现简单、每步只有两次 FFT特别适合二维相位恢复的快速验证。它的代价是没有严格的全局收敛保证初始相位离真实解太远时迭代可能停在某一组局部解上表现为重建图像边缘出现高频波纹。这个平台上把GS.zip和b_GS.m单独拆出来说明当时是边对照标准 GS 流程边修改的。标准 GS 流程里不设松弛因子只用硬替换工程里如果要改善收敛可以在替换相位时引入部分更新比如 (u^{(k1)}(1-\alpha)u^{(k)}\alpha u_{\text{proj}})其中 (\alpha) 通常取 0.5 到 0.8。这个改动不要直接写进主循环的初版先跑通无松弛版本再用误差曲线判断是否需要加。3. 在 Matlab 中落地 GS 主循环从 target.jpg 到重建相位这一层才是真正能动手改代码的地方。平台里的Untitled.m和b_GS.m应该是不同时期的主程序文件结构里既有.m又有.mat明显是实验过程中不断备份留下的。我复现时不会直接打开最大的脚本从头看到尾而是先拆出三块读数据、初始相位、迭代循环。3.1 读入振幅与生成初始相位target.jpg是二维灰度图可以直接当作物面振幅约束。读取时要注意图像可能不是方形GS 对非方形矩阵也能运行但后续傅里叶变换尺寸不变频谱坐标会受非对称采样影响相位误差统计也容易混入边缘效应。常见做法是先变成灰度再裁剪成和采样网格一致的尺寸并归一化到 0 到 1% 读取目标图作为物面振幅约束并统一尺寸 A im2double(imread(target.jpg)); if size(A, 3) 1 A rgb2gray(A); % 转为单通道强度 end A imresize(A, [256 256]); % 保证与频谱网格一致 A A / max(A(:)); % 归一化振幅保持能量量纲一致 % 初始相位均匀随机分布比常数更能避开对称解 rng(7); phi0 pi * (2 * rand(size(A)) - 1); u A .* exp(1i * phi0);这里的im2double把 uint8 数据转到双精度浮点避免后续exp和 FFT 时出现整数溢出。初始相位范围取 (-\pi) 到 (\pi)是因为相位本质上是周期量如果把范围缩小到 (-\pi/2) 到 (\pi/2)部分真实相位超出范围时迭代会明显变慢。rng(7)是为了固定随机种子方便对比不同参数下的误差曲线。3.2 GS 迭代主循环接下来是标准的 GS 交替投影主循环。物面振幅约束是A频谱振幅约束是B在有监督演示中B可以由A的频谱取模得到也可以读取a_sample_amplitude.bin文件获得实验测量振幅。主循环如下maxIter 200; err zeros(maxIter, 1); for k 1:maxIter U FT2Dc(u); % 从物面变换到频谱面 U B .* exp(1i * angle(U)); % 保持频谱相位替换频谱振幅 u IFT2Dc(U); % 变换回物面 u A .* exp(1i * angle(u)); % 保持物面相位替换物面振幅 err(k) sum(abs(abs(FT2Dc(u)) - B).^2, all) / sum(B(:).^2); end这段代码中angle(U)提取的是当前迭代下的频谱相位exp(1i * angle(U))把相位恢复成单位复振幅再乘上已知频谱振幅B。物面侧的替换同样只保留相位振幅被强制重置为A。误差指标用的是相对能量误差分母中的sum(B(:).^2)使得误差量纲与频谱总能量解耦便于比较不同尺寸和不同亮度图像下的收敛情况。这里有一个参数容易被忽略B和A是否满足能量一致性。GS 要求物面振幅和频谱振幅的能量一致也就是帕塞瓦尔定理成立如果A来源于target.jpg而B直接采用实验衍射图两边的总能量很可能差几个数量级这时迭代会稳定在噪声上。常见做法是在循环外B B / sqrt(sum(B(:).^2)) * sqrt(sum(A(:).^2))把能量对齐后再进入迭代。参数推荐值影响maxIter100500小于 50 通常欠收敛大于 500 提升有限初始随机相位范围([-π, π])范围过小容易停在对称解物面振幅归一化最大值归一化避免频谱动态范围过大频谱振幅能量匹配帕塞瓦尔归一化不匹配时误差曲线下不去支撑域 mask圆孔或矩形约束越强收敛越快但会引入截断伪影3.3 误差监测和相位可视化迭代结束时除了看err(end)我通常还会把重建相位和原始目标图的相位放在一起对比。target.jpg本身是强度图所以它并没有一个真正意义上的“真实相位”可直接比对工程里真正存放相位真值的是a_sample_phase.bin。读取这类二进制文件时要确认它保存的是单精度还是双精度浮点以及是按行还是按列写入。否则相位图会直接变成噪声fid fopen(a_sample_phase.bin, rb); phi_truth fread(fid, [256 256], single); fclose(fid); figure; subplot(1,2,1); imagesc(angle(u)); title(recovered phase); axis image; colormap(jet); colorbar; subplot(1,2,2); imagesc(phi_truth); title(sample phase); axis image; colormap(jet); colorbar;这里fread的single对应 32 位浮点如果实际文件是 64 位就会读取长度减半矩阵形状错乱。判断方法是先读取文件字节数除以 4观察是否能被 256×256 整除不能整除时改为double再试。这套判别步骤是排错里最常用的手段比直接看图像更可靠。4. 单缝、光栅与随机相位不同物面约束下的收敛行为GS 不是一个黑箱它的收敛行为直接由物面约束决定。平台里同时出现了singleslit.m、difraction_grating.m、RP.m和Gussianlight.m说明作者想把常规 GS、衍射光栅、随机相位屏和高斯光照明放在同一个框架里对比。这些文件之间共享FT2Dc.m和IFT2Dc.m差异只在于物面振幅A的构造。4.1 singleslit.m 与 difraction_grating.m把物面约束换成可调掩膜单缝衍射的物面约束是典型的矩形支撑域。若缝宽太小频谱会拉得很宽GS 在频域替换后容易出现能量泄漏到支撑域外缝宽太宽约束变弱相位恢复又多解。常见写法是在一维网格上构造布尔向量再扩展到二维例如N 512; L 5e-3; % 物面尺寸 5 mm x linspace(-L/2, L/2, N); w 0.3e-3; % 缝宽 0.3 mm slitMask abs(x) w/2; A double(slitMask); % 单缝在竖直方向无限的简化模型这段代码把x方向做成单缝y方向保持全 1得到的是竖直缝的物面振幅。注意实际仿真中缝宽和网格像素宽度需要匹配采样定理缝宽至少占到 4 个像素以上否则衍射图案的高频分量混叠GS 的频谱振幅约束本身就不可信。difraction_grating.m则会构造周期结构通常用余弦或方波函数生成周期光栅。光栅周期越小频谱上各级次之间的距离越远GS 需要的迭代次数也越多若周期接近采样间隔频谱级次会直接落在 Nyquist 边界之外相位恢复结果出现类似摩尔纹的重影。场景物面振幅特征频谱特征收敛行为单缝稀疏支撑域连续 sin 型包络收敛快但缝宽过窄时能量泄漏光栅周期调制分离的各级衍射峰级次间隔大时收敛慢随机相位屏均匀振幅随机相位散斑状频谱容易陷入局部解需增加迭代DOE连续浮雕相位非对称高光谱对初始相位敏感建议多随机重启4.2 散射成像相位恢复中的随机相位和迭代稳定性RP.m应该是生成随机相位屏的脚本这类约束在散射成像相位恢复里非常常见。随机相位屏的高频成分会把能量打得非常散GS 每次循环都在替换振幅、保留相位但物面侧的相位不断被随机初始值扰动导致误差曲线一开始快速下降随后进入一段长平台。这个平台期不是算错了而是 GS 正在慢慢调整中低频相位。处理这种场景我一般会做三次随机重启每次用不同rng种子跑 100 次迭代保留误差最小的结果。这个做法能显著降低散斑图案带来的局部极小问题。另一个技巧是在频域振幅约束上做一个低频增强把频谱根号强度作为替换模值而不是直接用功率谱强度这能削弱高频噪声对相位更新的干扰在散射成像相位恢复时效果尤其明显。4.3 a_simulate_DP.m 中常用的传播仿真链a_simulate_DP.m这个命名 DP 可能对应 Digital Propagation 或衍射传播它的角色是把物面复振幅传播到探测器平面。有了它就能在无实验数据的情况下生成合成衍射强度图再用这套强度图反向做 GS 相位恢复。工程里的a_sample_amplitude.bin和a_sample_phase.bin很可能就是这段链路生成的样本。如果要在自己的数据上复现这一链路通常流程是先把a_sample_phase.bin读取成相位矩阵加上高斯光振幅Gussianlight.jpg得到初始复振幅 (u_0)然后用角谱传播函数传播到探测器平面。角谱传播的传递函数为 (H(f_x,f_y)\exp(j2\pi z\sqrt{1-\lambda^2 f_x^2-\lambda^2 f_y^2}))。这一步在实际工程里经常用repmat或meshgrid构建频域坐标再与FT2Dc配合使用但要注意传递函数在 evanescent 波区域需要置零否则数值会发散。5. 收敛判据与病态配置GS 仿真的四个边界检查GS 跑完之后不能只盯着一张相位图。最后这一层是工程里最容易忽视的边界检查。我在拿到这套代码后会在输出前检查四项内容采样间隔是否满足 Nyquist、物面约束是否具备支撑域、频谱振幅是否做了能量对齐、最终相位是否出现棋盘格伪影。棋盘格伪影是 GS 最常见的失败模式表现是相位图中出现一个像素间隔的高频跳变。出现它说明频谱振幅约束里混入了混叠成分通常发生在L固定而N偏小时。检查方法是把采样间隔 (dxL/N) 和最高空间频率 (f_{\max}1/(2dx)) 打印出来再把频谱振幅B的有效截止频率和它对比。如果B在高频处还有显著能量说明采样不足。第二个边界检查是物面约束的支撑域。GS 需要足够强的约束才能收敛到物理正确的解。一个简单测试是把物面振幅A的像素总和除以全零数组的面积得到支撑域占比占比低于 5% 时收敛快但容易出现截断伪影高于 50% 时约束变弱误差曲线会有明显平台。对于随机相位屏这类振幅均匀的场景我一般会在物面侧额外加一个圆域掩模把支撑域收紧。第三个检查放在误差曲线末端。好的收敛曲线应当单调下降并进入平台如果末端出现周期震荡说明频谱能量和物面能量不匹配或者FT2Dc和IFT2Dc的坐标约定互相不对应。此时先跑 10 次迭代把每次的err打印出来震荡通常在偶数次和奇数次之间交替这是替换振幅时没有保持能量归一化造成的。最后一个技巧是把重建对象当成系统响应的验证实验只用一张已知的target.jpg手动给相位加上一个离焦二次项再看 GS 能否恢复出这个二次相位。这样做的优势是即便你的实验数据里没有真值相位也能验证传播链和 GS 主循环的物理正确性。检查时改动一个参数就够了改变波长或光栅周期前先确认max(fx) 1/(2*dx)否则 GS 会把混叠当真值收敛。本文还有配套的精品资源点击获取 Gerchberg-Saxton相位恢复算法：Matlab仿真与fftshift陷阱解析
简介基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台面向光学工程、信息光学及图像处理方向的学习者与研究者以光强数据反演相位分布为切入点系统展示了从相位恢复到波前重构的完整过程。包内共30个文件主要包括10个Matlab脚本.m、5个图形文件.fig、5张结果图片.jpg另有备份文件、数据文件及txt/md说明文档压缩包整体约53.55MB目录结构清晰便于按模块查阅。目前已有62人学习了该资源适合个人自学、课程实践或科研前期参考。通过修改输入参数或调整算法结构用户可深入分析GS算法在不同条件下的系统响应配合说明文档与示例图片能够快速掌握算法原理与实现细节并借助样例数据进行结果验证进一步提升光学仿真和算法开发的实践能力。1. 相位恢复问题的病态性Gerchberg-Saxton 迭代为何有效CCD 只记录强度波前相位被丢弃这是光学检测和成像里最难绕过的一道坎。Gerchberg-Saxton 迭代不依赖额外测量设备只靠物面振幅和频谱振幅这两个约束来回傅里叶变换就能把相位逼出来。手上这套 Matlab 仿真工程正好把这条链路拆开Gerchberg_Saxton_Algorithm.m是主循环FT2Dc.m、IFT2Dc.m负责正逆变换target.jpg、DOE.jpg提供目标图和衍射元件singleslit.m、difraction_grating.m则把单缝、光栅这类经典物面约束做成可切换场景。适合刚开始接触相位恢复的人跑通原理也适合做过 GS 但想研究收敛边界的人对照源码调参数。2. 傅里叶光学基础与频谱约束FT2Dc 的 fftshift 陷阱GS 的核心是空域和频域之间的交替投影但 Matlab 里的fft2和光学里的连续傅里叶变换并不一致尤其是坐标原点。这个平台把FT2Dc.m和IFT2Dc.m独立拆出来就是提示不要忽略这一步写错一个fftshift迭代出来的相位图会整体偏一个常数梯度看起来像倾斜波前实际上只是坐标绕乱了。2.1 正/逆傅里叶变换怎么和采样坐标对齐夫琅禾费衍射对应的连续变换是[ U(f_x, f_y)\iint u(x,y)\exp[-j2\pi(f_x x f_y y)]dxdy ]但离散 FFT 默认把零频放在矩阵的(1,1)位置而光学仿真中零频对应光轴中心。于是多数仿真代码会先对输入矩阵做ifftshift让阵列中心的光学原点移动到 FFT 约定的原点算完再fftshift把零频挪回中心。我一般把这类变换封装成独立函数避免在 GS 主循环里反复写错function F FT2Dc(u) % 二维光学傅里叶变换输入 u 的矩阵中心对应光轴原点 % 输出 F 的矩阵中心对应零频 F fftshift(fft2(ifftshift(u))); end function u IFT2Dc(F) % 对应的逆变换FT2Dc 与 IFT2Dc 可逆 u fftshift(ifft2(ifftshift(F))); end这段代码里ifftshift的作用是把中心原点挪到 FFT 的左上角原点fftshift负责把算完的零频移到矩阵中心。ifftshift和fftshift在偶数尺寸时行为相同在奇数尺寸时是严格互逆的所以不要用同一个函数替代另一个。真正验证这对函数是否写对的方法很简单对一个全 1 矩阵做FT2Dc理想频谱应该只有一个中心主峰直流分量集中在矩阵中央而不是四角。这部分在工程里容易被忽略因为单独画频谱时即使坐标偏了abs的图形看起来也只是整体平移。可一旦进入 GS 迭代物面约束和频谱约束会交替使用如果正逆变换的坐标约定不自洽相位恢复会收敛到一组看起来平滑但物理上完全错误的解。平台里每个.m文件都调用FT2Dc.m说明作者希望所有模块共享同一套坐标约定这也是光学仿真代码最值得先看的部分。函数在 GS 中的作用常见误用fft2/ifft2正逆傅里叶变换忘记先处理原点位置fftshift零频移到矩阵中心奇数尺寸时当成ifftshift用ifftshift中心原点移回 FFT 原点偶尺寸时不报错奇尺寸时错位abs取振幅约束对复数场直接取实部angle取相位约束相位跳变处需要包裹处理2.2 振幅约束与相位自由度的关系相位恢复问题的困难在于只有两个振幅约束可用来反推相位物面振幅 (A(x,y)) 和频谱振幅 (B(f_x,f_y))。它们之间的映射由傅里叶变换连接但相位本身不直接进入测量值。由于相位有大量自由度直接求逆是不适定的GS 算法的思路是先随机给一个初始相位然后反复执行两个投影在空域保持已知振幅、替换相位在频域同样保持已知振幅、替换相位。这里要注意物面的振幅约束并不一定等于目标图本身。这个工程把target.jpg当作已知振幅约束是一种有监督演示方便你比对重建结果和原图真正用于散射成像相位恢复时物面约束往往是照明光斑形状或一个支撑域掩模比如圆孔、狭缝、矩形区域。支撑域给得越大未知自由度越多GS 越容易陷入局部极小给得太小又会截断有效信号。工程里的singleslit.m、difraction_grating.m就是在构造不同类型的支撑域我建议你把这些文件里的 mask 改成圆形或环形再跑一轮观察收敛速度的变化。2.3 为什么选 GS 而不是直接求逆直接求逆需要同时知道复振幅而实验只能拿到强度。替代方案里有基于传输方程的 TIE、基于菲涅耳迭代的混合输入输出算法但 GS 的优势在于实现简单、每步只有两次 FFT特别适合二维相位恢复的快速验证。它的代价是没有严格的全局收敛保证初始相位离真实解太远时迭代可能停在某一组局部解上表现为重建图像边缘出现高频波纹。这个平台上把GS.zip和b_GS.m单独拆出来说明当时是边对照标准 GS 流程边修改的。标准 GS 流程里不设松弛因子只用硬替换工程里如果要改善收敛可以在替换相位时引入部分更新比如 (u^{(k1)}(1-\alpha)u^{(k)}\alpha u_{\text{proj}})其中 (\alpha) 通常取 0.5 到 0.8。这个改动不要直接写进主循环的初版先跑通无松弛版本再用误差曲线判断是否需要加。3. 在 Matlab 中落地 GS 主循环从 target.jpg 到重建相位这一层才是真正能动手改代码的地方。平台里的Untitled.m和b_GS.m应该是不同时期的主程序文件结构里既有.m又有.mat明显是实验过程中不断备份留下的。我复现时不会直接打开最大的脚本从头看到尾而是先拆出三块读数据、初始相位、迭代循环。3.1 读入振幅与生成初始相位target.jpg是二维灰度图可以直接当作物面振幅约束。读取时要注意图像可能不是方形GS 对非方形矩阵也能运行但后续傅里叶变换尺寸不变频谱坐标会受非对称采样影响相位误差统计也容易混入边缘效应。常见做法是先变成灰度再裁剪成和采样网格一致的尺寸并归一化到 0 到 1% 读取目标图作为物面振幅约束并统一尺寸 A im2double(imread(target.jpg)); if size(A, 3) 1 A rgb2gray(A); % 转为单通道强度 end A imresize(A, [256 256]); % 保证与频谱网格一致 A A / max(A(:)); % 归一化振幅保持能量量纲一致 % 初始相位均匀随机分布比常数更能避开对称解 rng(7); phi0 pi * (2 * rand(size(A)) - 1); u A .* exp(1i * phi0);这里的im2double把 uint8 数据转到双精度浮点避免后续exp和 FFT 时出现整数溢出。初始相位范围取 (-\pi) 到 (\pi)是因为相位本质上是周期量如果把范围缩小到 (-\pi/2) 到 (\pi/2)部分真实相位超出范围时迭代会明显变慢。rng(7)是为了固定随机种子方便对比不同参数下的误差曲线。3.2 GS 迭代主循环接下来是标准的 GS 交替投影主循环。物面振幅约束是A频谱振幅约束是B在有监督演示中B可以由A的频谱取模得到也可以读取a_sample_amplitude.bin文件获得实验测量振幅。主循环如下maxIter 200; err zeros(maxIter, 1); for k 1:maxIter U FT2Dc(u); % 从物面变换到频谱面 U B .* exp(1i * angle(U)); % 保持频谱相位替换频谱振幅 u IFT2Dc(U); % 变换回物面 u A .* exp(1i * angle(u)); % 保持物面相位替换物面振幅 err(k) sum(abs(abs(FT2Dc(u)) - B).^2, all) / sum(B(:).^2); end这段代码中angle(U)提取的是当前迭代下的频谱相位exp(1i * angle(U))把相位恢复成单位复振幅再乘上已知频谱振幅B。物面侧的替换同样只保留相位振幅被强制重置为A。误差指标用的是相对能量误差分母中的sum(B(:).^2)使得误差量纲与频谱总能量解耦便于比较不同尺寸和不同亮度图像下的收敛情况。这里有一个参数容易被忽略B和A是否满足能量一致性。GS 要求物面振幅和频谱振幅的能量一致也就是帕塞瓦尔定理成立如果A来源于target.jpg而B直接采用实验衍射图两边的总能量很可能差几个数量级这时迭代会稳定在噪声上。常见做法是在循环外B B / sqrt(sum(B(:).^2)) * sqrt(sum(A(:).^2))把能量对齐后再进入迭代。参数推荐值影响maxIter100500小于 50 通常欠收敛大于 500 提升有限初始随机相位范围([-π, π])范围过小容易停在对称解物面振幅归一化最大值归一化避免频谱动态范围过大频谱振幅能量匹配帕塞瓦尔归一化不匹配时误差曲线下不去支撑域 mask圆孔或矩形约束越强收敛越快但会引入截断伪影3.3 误差监测和相位可视化迭代结束时除了看err(end)我通常还会把重建相位和原始目标图的相位放在一起对比。target.jpg本身是强度图所以它并没有一个真正意义上的“真实相位”可直接比对工程里真正存放相位真值的是a_sample_phase.bin。读取这类二进制文件时要确认它保存的是单精度还是双精度浮点以及是按行还是按列写入。否则相位图会直接变成噪声fid fopen(a_sample_phase.bin, rb); phi_truth fread(fid, [256 256], single); fclose(fid); figure; subplot(1,2,1); imagesc(angle(u)); title(recovered phase); axis image; colormap(jet); colorbar; subplot(1,2,2); imagesc(phi_truth); title(sample phase); axis image; colormap(jet); colorbar;这里fread的single对应 32 位浮点如果实际文件是 64 位就会读取长度减半矩阵形状错乱。判断方法是先读取文件字节数除以 4观察是否能被 256×256 整除不能整除时改为double再试。这套判别步骤是排错里最常用的手段比直接看图像更可靠。4. 单缝、光栅与随机相位不同物面约束下的收敛行为GS 不是一个黑箱它的收敛行为直接由物面约束决定。平台里同时出现了singleslit.m、difraction_grating.m、RP.m和Gussianlight.m说明作者想把常规 GS、衍射光栅、随机相位屏和高斯光照明放在同一个框架里对比。这些文件之间共享FT2Dc.m和IFT2Dc.m差异只在于物面振幅A的构造。4.1 singleslit.m 与 difraction_grating.m把物面约束换成可调掩膜单缝衍射的物面约束是典型的矩形支撑域。若缝宽太小频谱会拉得很宽GS 在频域替换后容易出现能量泄漏到支撑域外缝宽太宽约束变弱相位恢复又多解。常见写法是在一维网格上构造布尔向量再扩展到二维例如N 512; L 5e-3; % 物面尺寸 5 mm x linspace(-L/2, L/2, N); w 0.3e-3; % 缝宽 0.3 mm slitMask abs(x) w/2; A double(slitMask); % 单缝在竖直方向无限的简化模型这段代码把x方向做成单缝y方向保持全 1得到的是竖直缝的物面振幅。注意实际仿真中缝宽和网格像素宽度需要匹配采样定理缝宽至少占到 4 个像素以上否则衍射图案的高频分量混叠GS 的频谱振幅约束本身就不可信。difraction_grating.m则会构造周期结构通常用余弦或方波函数生成周期光栅。光栅周期越小频谱上各级次之间的距离越远GS 需要的迭代次数也越多若周期接近采样间隔频谱级次会直接落在 Nyquist 边界之外相位恢复结果出现类似摩尔纹的重影。场景物面振幅特征频谱特征收敛行为单缝稀疏支撑域连续 sin 型包络收敛快但缝宽过窄时能量泄漏光栅周期调制分离的各级衍射峰级次间隔大时收敛慢随机相位屏均匀振幅随机相位散斑状频谱容易陷入局部解需增加迭代DOE连续浮雕相位非对称高光谱对初始相位敏感建议多随机重启4.2 散射成像相位恢复中的随机相位和迭代稳定性RP.m应该是生成随机相位屏的脚本这类约束在散射成像相位恢复里非常常见。随机相位屏的高频成分会把能量打得非常散GS 每次循环都在替换振幅、保留相位但物面侧的相位不断被随机初始值扰动导致误差曲线一开始快速下降随后进入一段长平台。这个平台期不是算错了而是 GS 正在慢慢调整中低频相位。处理这种场景我一般会做三次随机重启每次用不同rng种子跑 100 次迭代保留误差最小的结果。这个做法能显著降低散斑图案带来的局部极小问题。另一个技巧是在频域振幅约束上做一个低频增强把频谱根号强度作为替换模值而不是直接用功率谱强度这能削弱高频噪声对相位更新的干扰在散射成像相位恢复时效果尤其明显。4.3 a_simulate_DP.m 中常用的传播仿真链a_simulate_DP.m这个命名 DP 可能对应 Digital Propagation 或衍射传播它的角色是把物面复振幅传播到探测器平面。有了它就能在无实验数据的情况下生成合成衍射强度图再用这套强度图反向做 GS 相位恢复。工程里的a_sample_amplitude.bin和a_sample_phase.bin很可能就是这段链路生成的样本。如果要在自己的数据上复现这一链路通常流程是先把a_sample_phase.bin读取成相位矩阵加上高斯光振幅Gussianlight.jpg得到初始复振幅 (u_0)然后用角谱传播函数传播到探测器平面。角谱传播的传递函数为 (H(f_x,f_y)\exp(j2\pi z\sqrt{1-\lambda^2 f_x^2-\lambda^2 f_y^2}))。这一步在实际工程里经常用repmat或meshgrid构建频域坐标再与FT2Dc配合使用但要注意传递函数在 evanescent 波区域需要置零否则数值会发散。5. 收敛判据与病态配置GS 仿真的四个边界检查GS 跑完之后不能只盯着一张相位图。最后这一层是工程里最容易忽视的边界检查。我在拿到这套代码后会在输出前检查四项内容采样间隔是否满足 Nyquist、物面约束是否具备支撑域、频谱振幅是否做了能量对齐、最终相位是否出现棋盘格伪影。棋盘格伪影是 GS 最常见的失败模式表现是相位图中出现一个像素间隔的高频跳变。出现它说明频谱振幅约束里混入了混叠成分通常发生在L固定而N偏小时。检查方法是把采样间隔 (dxL/N) 和最高空间频率 (f_{\max}1/(2dx)) 打印出来再把频谱振幅B的有效截止频率和它对比。如果B在高频处还有显著能量说明采样不足。第二个边界检查是物面约束的支撑域。GS 需要足够强的约束才能收敛到物理正确的解。一个简单测试是把物面振幅A的像素总和除以全零数组的面积得到支撑域占比占比低于 5% 时收敛快但容易出现截断伪影高于 50% 时约束变弱误差曲线会有明显平台。对于随机相位屏这类振幅均匀的场景我一般会在物面侧额外加一个圆域掩模把支撑域收紧。第三个检查放在误差曲线末端。好的收敛曲线应当单调下降并进入平台如果末端出现周期震荡说明频谱能量和物面能量不匹配或者FT2Dc和IFT2Dc的坐标约定互相不对应。此时先跑 10 次迭代把每次的err打印出来震荡通常在偶数次和奇数次之间交替这是替换振幅时没有保持能量归一化造成的。最后一个技巧是把重建对象当成系统响应的验证实验只用一张已知的target.jpg手动给相位加上一个离焦二次项再看 GS 能否恢复出这个二次相位。这样做的优势是即便你的实验数据里没有真值相位也能验证传播链和 GS 主循环的物理正确性。检查时改动一个参数就够了改变波长或光栅周期前先确认max(fx) 1/(2*dx)否则 GS 会把混叠当真值收敛。本文还有配套的精品资源点击获取