简介本资源是一套面向光学专业本科生与初学者的Matlab仿真教学工具包聚焦Gerchberg-SaxtonGS迭代算法在光学相位恢复与波前重建中的原理实现与可视化验证专为中国科学技术大学光学课程作业设计解决“仅知强度分布、如何反演相位信息”这一核心教学难点。压缩包共27个文件含10个核心Matlab脚本如Gerchberg_Saxton_Algorithm.m、FT2Dc.m、RP.m等、5幅算法中间过程及结果图.fig、5张关键输出图像.jpg含目标图、重建像、相位分布等以及说明文档.txt、理论补充.docx和二进制样例数据.bin总大小39.03MB。已有75人学习下载资源结构清晰主程序模块化、参数可调、每步迭代可视化附赠文档详述算法推导与计算全息关联配套fig与jpg直观呈现空间域与频域交替约束效果便于理解收敛行为与重建质量评估。1. 项目概述从课程作业到光学成像的核心技术看到这个项目标题很多光学、物理或者电子信息专业的同学应该会心一笑。这几乎是中国科学技术大学光学相关课程里一个经典的“大作业”或课程设计题目。它的核心是利用Gerchberg-SaxtonGS迭代算法在MATLAB环境中仿真实现光学相位恢复与波前重建。简单来说我们平时用相机拍照记录的是光波的强度信息也就是明暗但丢失了至关重要的相位信息。而相位决定了光波的波前形状包含了物体深度、折射率分布等关键信息。GS算法就是一种非常巧妙的方法它允许我们仅通过采集到的光强信息比如两张不同位置的光斑图通过数学迭代反推出丢失的相位从而重建出完整的光波场。这个项目绝不仅仅是为了完成作业。它在现代光学的前沿领域有着举足轻重的地位。计算全息、自适应光学用于校正大气湍流对天文望远镜的影响、相位显微成像用于观察透明生物样品、光学加密以及新型的光场调控等领域其底层核心都离不开相位恢复技术。通过亲手用MATLAB实现GS算法你不仅能深刻理解傅里叶光学、衍射理论的核心概念更能掌握一套解决实际逆向问题的计算工具。这个过程是从理论公式走向实际应用的关键一步你会遇到算法收敛性、初始猜测敏感性、噪声影响等一系列非常实际的问题而解决它们的过程正是能力提升的体现。2. 核心原理拆解GS算法如何“无中生有”出相位要理解GS算法我们得先搞清楚光波场是什么。一个单色光波在空间某一点的复振幅可以表示为 U(x, y) A(x, y) * exp(i * φ(x, y))。其中A是振幅amplitude直接对应我们探测到的光强 I |A|²φ就是相位phase它描述了光波在该点相对于某个参考点的振动状态。探测器如CCD相机是“相位盲”的只能记录光强I宝贵的相位φ在记录瞬间就丢失了。GS算法的天才之处在于它利用了两个域输入面和输出面通常是物面和频谱面之间的傅里叶变换关系以及我们对这两个面上光强分布的已知信息通过迭代来逼近真实的相位。其核心流程是一个闭合的迭代循环我们可以把它想象成一个不断自我修正的猜谜游戏。算法的基本迭代步骤一个循环如下初始猜测我们从物面输入面开始。已知物面的光强分布 I_object比如一个物体的透射率。但对于相位我们一无所知。因此我们随机生成一个初始相位猜测 φ_guess通常在0到2π之间均匀分布与已知的振幅 sqrt(I_object) 结合构造一个初始的复振幅场U1 sqrt(I_object) * exp(i * φ_guess)。正向传播约束变换将构造的物面场 U1 传播到输出面通常是焦平面或衍射场。这个传播过程在傍轴近似下通过一次傅里叶变换FFT来实现U2 FFT(U1)。此时U2 同时具有振幅 A2 和相位 φ2。施加输出面约束这是我们拥有的另一个已知信息。在输出面例如透镜的后焦面我们通过实验测量得到了该处的光强分布 I_measured。GS算法的关键操作来了我们保留计算得到的相位 φ2但用测量得到的振幅 sqrt(I_measured) 替换掉计算出的振幅 A2。即更新输出面场为U2‘ sqrt(I_measured) * exp(i * φ2)。这个步骤强制让计算出的波场满足实际的物理测量结果。反向传播逆变换将更新后的输出面场 U2‘ 传播回物面。通过逆傅里叶变换IFFT实现U1‘ IFFT(U2‘)。此时我们得到了一个新的物面复振幅场包含振幅 A1‘ 和相位 φ1‘。施加输入面约束同样我们保留新计算出的相位 φ1‘但用我们最初已知的物面振幅 sqrt(I_object) 替换掉计算出的振幅 A1‘。即更新物面场为U1_new sqrt(I_object) * exp(i * φ1‘)。这个步骤强制让波场满足物面的已知条件。至此一个迭代循环完成。U1_new 将作为下一个循环的起点重复步骤2到5。如此反复迭代每一次循环都让计算出的波场在输入面和输出面同时满足各自的振幅强度约束。理论上随着迭代次数增加计算得到的相位 φ1‘ 会收敛到真实的物面相位。注意这里描述的是最基本的“两平面”GS算法。在实际应用中根据不同的物理设置如多平面相位恢复、基于角谱的衍射传播等正向和反向传播算子可能不是简单的FFT/IFFT而是更复杂的衍射积分计算。但核心的“变换-施加振幅约束-反变换”的迭代思想是完全一致的。为什么这样迭代会收敛直观上理解这个算法在寻找一个解这个解在经过特定变换如傅里叶变换后其振幅分布能与我们在另一个平面上测量到的结果匹配。迭代过程就像在解一个巨大的方程组每次施加约束都是将解向可行解集合投影。GS算法本质是一种交替投影算法在两个约束集合满足输入面振幅的复函数集合和满足输出面振幅的复函数集合之间来回投影最终希望找到两个集合的交点即同时满足两个约束的解。3. MATLAB仿真系统设计与实现要点用MATLAB搭建这个仿真系统目的很明确在已知“真实”相位的情况下模拟光波的传播和探测过程生成仿真的“测量光强”然后假装我们不知道这个相位仅用“测量光强”和GS算法来恢复它最后与“真实”相位对比评估算法性能。这个过程让我们可以在完全可控的环境下研究算法的特性。3.1 系统模块分解一个结构清晰的仿真系统通常包含以下几个模块参数定义与初始化模块设定仿真环境的基本参数。这包括网格参数仿真区域的尺寸Lx,Ly单位米或毫米采样点数Nx,Ny。采样点数一般取2的整数次幂如2565121024以便利用FFT的高效性。网格间距dx Lx/Nx。光学参数波长lambda传播距离z如果是角谱传播透镜焦距f如果是4f系统等。目标物体定义输入面的振幅分布A_obj如一个圆形孔径、一个字母图案和真实的相位分布phi_true这是我们希望恢复的可以是一个简单的球面波相位、一个随机相位板或者一个复杂的像差分布。生成模拟的真实波场U_true A_obj .* exp(1i * phi_true)。正向传播模拟模块这个模块模拟真实物理过程即光从物面传播到探测面。根据实验光路选择合适的传播模型。对于最常见的4f滤波系统或焦面探测使用傅里叶变换U_prop fft2(U_true)或fftshift(fft2(fftshift(U_true)))注意零频分量移到中心。对于非焦面的自由空间衍射可能需要使用角谱传播或菲涅尔衍射积分。角谱传播更精确其核心是AS fft2(U_true); H exp(1i * 2*pi*z / lambda .* sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2)); U_prop ifft2(AS .* H);其中fx,fy是空间频率坐标。计算探测面的光强模拟测量结果I_measured abs(U_prop).^2。为了更真实可以在此处加入噪声如高斯噪声I_measured_noisy imnoise(I_measured, gaussian, 0, 0.01)添加方差为0.01的高斯噪声。GS算法迭代恢复模块这是系统的核心。初始化设定最大迭代次数max_iter收敛阈值threshold。生成随机的初始相位phi_guess 2*pi*rand(Nx, Ny)。构造初始估计场U_est A_obj .* exp(1i * phi_guess)。迭代循环for iter 1:max_iter正向传播至探测面U_prop_est Propagation_Model(U_est)使用与正向模拟模块相同的传播模型。施加探测面振幅约束A_prop_measured sqrt(I_measured)。U_prop_est_constrained A_prop_measured .* exp(1i * angle(U_prop_est))。angle()函数用于获取相位。反向传播回物面U_est_back Inverse_Propagation_Model(U_prop_est_constrained)。施加物面振幅约束U_est_new A_obj .* exp(1i * angle(U_est_back))。计算误差评估当前估计与约束的匹配程度。常用误差函数是探测面振幅误差error(iter) sum(sum(abs(abs(U_prop_est) - A_prop_measured).^2)) / sum(sum(A_prop_measured.^2))。判断收敛如果error(iter) threshold跳出循环。更新估计场U_est U_est_new。输出循环结束后从最终的U_est中提取恢复的相位phi_recovered angle(U_est)。结果可视化与评估模块这是验证算法效果的关键。图像显示并排显示真实相位phi_true、恢复相位phi_recovered以及两者的差值残差图。由于相位是周期性的模2π通常使用imagesc显示并注意相位包裹phase wrapping问题可能需要使用angle()函数直接显示范围 -π 到 π或使用wrapToPi等函数处理。误差分析绘制迭代误差error随迭代次数的下降曲线观察收敛速度和稳定性。定量指标计算恢复相位与真实相位之间的均方根误差RMSE、峰值信噪比PSNR或结构相似性SSIM。对于存在相位包裹的情况比较前需要先进行解包裹phase unwrapping或计算梯度。3.2 关键代码片段与实操心得下面给出一些核心步骤的MATLAB代码片段和编写时的注意事项1. 生成带有相位的测试物体% 参数定义 N 512; % 采样点数 L 10e-3; % 物理尺寸 10mm dx L/N; x linspace(-L/2, L/2, N); [X, Y] meshgrid(x, x); % 振幅物体一个圆形孔径 radius 1e-3; % 1mm半径 A_obj sqrt((X.^2 Y.^2) radius^2); % 振幅1代表透光0代表不透光 % 真实相位一个球面波相位模拟透镜叠加一些像差 lambda 632.8e-9; % He-Ne激光波长 k 2*pi/lambda; R 0.5; % 球面波曲率半径 phi_true mod(k/(2*R) * (X.^2 Y.^2) 0.5*sin(4*pi*X/L) .* cos(2*pi*Y/L), 2*pi); % 使用mod函数将其范围限制在[0, 2π)避免数值过大。实际仿真中相位值范围不重要重要的是其空间变化。 % 第二部分加入了一个简单的正弦型像差。 % 构造真实复振幅场 U_true A_obj .* exp(1i * phi_true);实操心得在定义相位phi_true时要特别注意其空间频率不能超过奈奎斯特频率即每个采样点间隔内相位变化不能超过π否则会产生混叠导致仿真失真。对于快速变化的相位需要确保采样足够密集。2. 角谱传播函数用于自由空间衍射模拟function Uout AngularSpectrumPropagation(Uin, lambda, dx, dy, z) % Uin: 输入场 % lambda: 波长 % dx, dy: 输入面采样间隔 % z: 传播距离 % Uout: 输出场 [Ny, Nx] size(Uin); % 生成空间频率坐标 fx linspace(-1/(2*dx), 1/(2*dx), Nx); fy linspace(-1/(2*dy), 1/(2*dy), Ny); [FX, FY] meshgrid(fx, fy); % 计算角谱传递函数 H exp(1i * 2*pi/lambda * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 对于频率分量满足 (lambda*fx)^2 (lambda*fy)^2 1 的部分sqrt结果为虚数对应倏逝波通常直接置零或做衰减处理。 H(sqrt((lambda*FX).^2 (lambda*FY).^2) 1) 0; % 进行角谱传播 Uin_f fft2(fftshift(Uin)); % 将输入场频谱零频移到角落便于与H相乘 Uout_f Uin_f .* H; Uout ifftshift(ifft2(Uout_f)); % 反变换并移回中心 end注意角谱传播是精确的标量衍射模型但计算时要注意传递函数H中根号内的项为负的情况对应高频倏逝波处理不当会引入数值不稳定。通常的做法是将其滤除置零如上面代码所示。这相当于一个低通滤波是物理合理的。3. GS算法迭代核心循环max_iter 200; threshold 1e-6; error_history zeros(max_iter, 1); % 初始猜测随机相位 phi_guess 2*pi * rand(N, N); U_est A_obj .* exp(1i * phi_guess); % 已知的测量面振幅模拟得到已加噪声 A_measured sqrt(I_measured_noisy); for iter 1:max_iter % 1. 正向传播 (使用角谱传播或FFT需与正向模拟一致) U_prop_est AngularSpectrumPropagation(U_est, lambda, dx, dx, z); % 2. 施加测量面振幅约束 A_prop_est abs(U_prop_est); U_prop_est_constrained A_measured .* exp(1i * angle(U_prop_est)); % 3. 反向传播 U_est_back AngularSpectrumPropagation(U_prop_est_constrained, lambda, dx, dx, -z); % 注意距离取负 % 4. 施加物面振幅约束 U_est_new A_obj .* exp(1i * angle(U_est_back)); % 5. 计算误差测量面振幅误差 error_history(iter) sum(sum((A_prop_est - A_measured).^2)) / sum(sum(A_measured.^2)); % 6. 判断收敛 if error_history(iter) threshold fprintf(在 %d 次迭代后收敛。\n, iter); break; end % 7. 更新估计 U_est U_est_new; end phi_recovered angle(U_est);实操心得反向传播函数理论上应该是正向传播的逆运算。对于角谱传播传递函数是exp(-i*...)即距离取负-z。对于FFT模型4f系统反向传播就是IFFT。务必保持正向模拟和迭代恢复中使用完全相同的传播模型和参数否则算法可能无法收敛到正确解。4. 影响算法性能的关键因素与调优策略GS算法原理简洁但实际应用中其性能受到多种因素影响。理解这些因素并学会调优是完成高质量相位恢复的关键。4.1 初始相位猜测GS算法本质上是一个非线性优化过程初始值会影响收敛速度和最终结果尤其是在解不唯一或存在噪声的情况下。随机相位最常用的方法。2*pi*rand()生成均匀分布的随机相位。优点是简单有助于算法跳出局部极小值。缺点是每次运行结果可能略有不同收敛所需的迭代次数不稳定。常数相位零相位zeros(N, N)。作为最简单的猜测在物体相位变化平缓且信噪比较高时可能收敛很快。但如果真实相位变化剧烈零相位猜测可能引导算法收敛到一个错误的解局部极小值。先验信息如果对物体相位有部分了解例如知道它是一个缓慢变化的像差可以将其作为初始猜测能显著加速收敛并提高准确性。策略在仿真中可以尝试多种初始猜测观察收敛结果的一致性。在实际实验中如果没有先验信息通常采用随机相位并多次独立运行算法选取误差最小的结果或对多次结果进行平均以增强鲁棒性。4.2 振幅约束的施加与支持域Support Constraint在基础GS算法中我们对物面的振幅约束是严格的A_obj。但在很多情况下我们只知道物体的大致范围支持域而不知道其内部精确的振幅分布。例如在晶体学中只知道分子的大致形状。紧支持域定义一个二值掩膜support在物体可能存在的区域为1其他区域为0。物面约束变为U_est_new (A_obj .* exp(1i*angle(U_est_back))) .* support (U_est_back .* (1-support)) * beta。这里beta是一个松弛因子0beta1用于在非支持域缓慢衰减场而不是直接置零这有助于提高收敛性这就是著名的混合输入输出HIO算法的核心思想。松弛振幅约束对于物面振幅如果不是精确已知可以施加一个范围约束而不是固定值。策略如果你的项目涉及非均匀振幅物体或支持域信息强烈建议实现HIO或其变种如OSP、RAAR算法它们比纯GS算法具有更强的抗噪声能力和更广的收敛域。4.3 噪声的影响与处理实验测量的光强I_measured必然包含噪声散粒噪声、读出噪声等。噪声会破坏振幅约束的有效性导致GS算法振荡或不收敛到正确解。现象在噪声较大时误差曲线error_history不会平稳下降到零而是在某个值附近波动。恢复的相位图中会出现明显的随机斑点或条纹。处理方法多次独立运行取平均用不同的随机初始相位多次运行GS算法将最终恢复的相位进行平均可以有效抑制随机噪声带来的波动。正则化在迭代过程中加入正则化项惩罚相位的高频变化即假设相位是平滑的。例如在每次迭代后对恢复的相位进行轻微的高斯滤波。使用更鲁棒的算法如前所述的HIO算法对噪声的容忍度比GS更高。数据预处理对测量的光强图像进行平滑滤波或降噪处理如小波去噪、非局部均值去噪但要注意不能过度平滑而损失细节。4.4 收敛性与停滞问题GS算法可能陷入停滞stagnation即误差不再下降但恢复的结果明显不正确。这通常是因为算法陷入了局部极小值。诊断观察误差曲线如果曲线在迭代几十次后早早地进入平台期且恢复的相位与真实相位相差甚远可能就是陷入了停滞。解决方案改变初始猜测换一个随机种子重新开始。引入扰动在迭代过程中每隔一定次数人为地对当前估计的相位加入一个小的随机扰动帮助算法跳出局部极小。这被称为“扰动GS”或“随机梯度”思想的一种简单实现。切换算法当GS停滞时切换到HIO算法迭代几十次然后再切换回GS往往能打破僵局。这种GS-HIO交替迭代的策略非常有效。多平面相位恢复如果条件允许采集多个不同传播距离或离焦量下的光强图像。对每个平面的光强同时施加约束大大增加了算法的约束条件能显著改善收敛性和唯一性。这就是著名的多平面迭代相位恢复技术。5. 仿真结果分析、可视化与报告撰写完成了算法实现和迭代如何呈现你的工作成果至关重要。一份优秀的课程报告或项目总结需要清晰、直观地展示仿真过程和结论。5.1 结果可视化技巧相位显示直接使用imagesc(phi_recovered)显示相位由于相位通常被包裹在[-π, π]或[0, 2π]区间会呈现出明显的2π跳变条纹。为了更直观地看到相位的连续变化需要进行相位解包裹。% 使用MATLAB的 unwrap函数适用于1D对2D效果有限或第三方解包裹算法 % 一个简单的2D解包裹思路质量图引导 % phi_wrapped angle(U_est); % 包裹相位 % 可以使用基于最小二乘的解包裹算法如 unwrapPhaseMesh (需要自己实现或找工具箱) % 或者直接显示相位梯度如x和y方向的差分梯度图不受2π跳变影响能反映相位的局部变化。 dx_phase diff(phi_recovered, 1, 2); % x方向差分 dy_phase diff(phi_recovered, 1, 1); % y方向差分对于课程作业如果真实相位本身变化范围不超过2π也可以直接比较包裹后的相位。但务必在报告中说明这一点。误差分析图迭代误差曲线semilogy(1:iter, error_history(1:iter))。使用对数坐标可以更清楚地看到误差前期的快速下降和后期的缓慢收敛。相位残差图imagesc(wrapToPi(phi_recovered - phi_true))。显示恢复相位与真实相位的差值包裹到[-π, π]。理想的恢复结果应该是一张均匀的灰度图差值接近0。剖面线对比在相位图上选取一条水平或垂直线绘制真实相位和恢复相位沿该线的变化曲线进行直接对比。line_index N/2; % 中间一行 figure; plot(x, phi_true(line_index, :), b-, LineWidth, 1.5, DisplayName, 真实相位); hold on; plot(x, phi_recovered(line_index, :), r--, LineWidth, 1.5, DisplayName, 恢复相位); xlabel(位置 (m)); ylabel(相位 (rad)); legend; title(相位剖面线对比);定量指标计算% 1. 均方根误差 (RMSE) - 需解包裹或处理跳变 phase_diff phi_recovered - phi_true; % 处理2π跳变将差值映射到[-π, π]区间 phase_diff_wrapped mod(phase_diff pi, 2*pi) - pi; rmse sqrt(mean(phase_diff_wrapped(:).^2)); fprintf(相位RMSE: %.4f rad\n, rmse); % 2. 相关系数 (Correlation Coefficient) C corrcoef(phi_true(:), phi_recovered(:)); correlation C(1,2); fprintf(相位相关系数: %.4f\n, correlation); % 3. 对于振幅恢复可以计算相对误差 amp_recovered abs(U_est); amp_error norm(amp_recovered(:) - A_obj(:)) / norm(A_obj(:)); fprintf(振幅相对误差: %.4f\n, amp_error);5.2 探索性仿真实验设计为了深入理解算法不要只满足于一个成功运行的例子。设计一系列对比实验能让你的报告脱颖而出噪声水平影响实验固定其他参数改变添加到I_measured上的高斯噪声方差例如从0到0.1观察恢复相位的RMSE和迭代收敛曲线如何变化。绘制“噪声水平-恢复误差”关系图。初始猜测敏感性实验对同一个物体和噪声水平分别用随机相位、零相位和一种有偏的初始猜测例如使用一个离焦的球面波相位运行GS算法。比较它们的最终误差、收敛速度和恢复结果的主观质量。支持域约束实验如果物体是一个简单形状如矩形尝试在物面约束中引入支持域。比较使用精确振幅约束、宽松支持域约束HIO算法以及无支持域信息时算法的恢复能力。可以尝试恢复一个在支持域外也有微弱信号的物体。欠采样与混叠实验故意减少采样点数N或者增加相位物体的空间频率使其超过奈奎斯特频率。观察混叠如何导致恢复失败并解释原因。5.3 项目报告与代码组织心得最后将你的工作整理成一份可交付的作业或项目报告。代码组织建议使用脚本.m脚本作为主程序调用多个函数文件.m函数。将核心功能模块化create_test_object.m,forward_propagation.m,gs_algorithm.m,hio_algorithm.m,plot_results.m。使用清晰的变量名添加充分的注释。特别是对于光学参数和坐标定义注释其物理意义和单位。在主脚本开头集中定义所有可调参数方便进行不同的实验。使用savefig或exportgraphics将生成的图表高质量保存便于插入报告。报告撰写要点引言简要说明相位恢复的意义、GS算法的背景和本项目目标。原理清晰阐述GS算法的数学原理和迭代步骤最好配以流程图。方法详细描述你的MATLAB仿真模型包括光学设置、参数、噪声模型以及你实现的GS算法细节是否包含HIO如何处理支持域。结果与分析这是核心部分。展示关键实验结果图原始物体、测量光强、恢复相位、误差曲线、剖面对比等。结合图表分析算法性能讨论前面提到的关键因素初始猜测、噪声、支持域的影响。展示你设计的探索性实验的结果和结论。讨论与结论总结GS算法的优缺点你在实现过程中遇到的主要挑战和解决方案以及对未来改进的设想如尝试更先进的算法Fienup算法、PIE、ePIE等。参考文献引用关键的教材、论文或网络资源。通过这样一个从理论到实践从实现到分析从基础到拓展的完整过程你完成的就不仅仅是一个课程作业而是一个扎实的光学计算研究小课题。这份经验对你理解更复杂的计算成像算法乃至从事相关领域的研究或开发工作都将是一块重要的基石。本文还有配套的精品资源点击获取 MATLAB实现GS算法：从相位恢复原理到光学成像仿真实践
简介本资源是一套面向光学专业本科生与初学者的Matlab仿真教学工具包聚焦Gerchberg-SaxtonGS迭代算法在光学相位恢复与波前重建中的原理实现与可视化验证专为中国科学技术大学光学课程作业设计解决“仅知强度分布、如何反演相位信息”这一核心教学难点。压缩包共27个文件含10个核心Matlab脚本如Gerchberg_Saxton_Algorithm.m、FT2Dc.m、RP.m等、5幅算法中间过程及结果图.fig、5张关键输出图像.jpg含目标图、重建像、相位分布等以及说明文档.txt、理论补充.docx和二进制样例数据.bin总大小39.03MB。已有75人学习下载资源结构清晰主程序模块化、参数可调、每步迭代可视化附赠文档详述算法推导与计算全息关联配套fig与jpg直观呈现空间域与频域交替约束效果便于理解收敛行为与重建质量评估。1. 项目概述从课程作业到光学成像的核心技术看到这个项目标题很多光学、物理或者电子信息专业的同学应该会心一笑。这几乎是中国科学技术大学光学相关课程里一个经典的“大作业”或课程设计题目。它的核心是利用Gerchberg-SaxtonGS迭代算法在MATLAB环境中仿真实现光学相位恢复与波前重建。简单来说我们平时用相机拍照记录的是光波的强度信息也就是明暗但丢失了至关重要的相位信息。而相位决定了光波的波前形状包含了物体深度、折射率分布等关键信息。GS算法就是一种非常巧妙的方法它允许我们仅通过采集到的光强信息比如两张不同位置的光斑图通过数学迭代反推出丢失的相位从而重建出完整的光波场。这个项目绝不仅仅是为了完成作业。它在现代光学的前沿领域有着举足轻重的地位。计算全息、自适应光学用于校正大气湍流对天文望远镜的影响、相位显微成像用于观察透明生物样品、光学加密以及新型的光场调控等领域其底层核心都离不开相位恢复技术。通过亲手用MATLAB实现GS算法你不仅能深刻理解傅里叶光学、衍射理论的核心概念更能掌握一套解决实际逆向问题的计算工具。这个过程是从理论公式走向实际应用的关键一步你会遇到算法收敛性、初始猜测敏感性、噪声影响等一系列非常实际的问题而解决它们的过程正是能力提升的体现。2. 核心原理拆解GS算法如何“无中生有”出相位要理解GS算法我们得先搞清楚光波场是什么。一个单色光波在空间某一点的复振幅可以表示为 U(x, y) A(x, y) * exp(i * φ(x, y))。其中A是振幅amplitude直接对应我们探测到的光强 I |A|²φ就是相位phase它描述了光波在该点相对于某个参考点的振动状态。探测器如CCD相机是“相位盲”的只能记录光强I宝贵的相位φ在记录瞬间就丢失了。GS算法的天才之处在于它利用了两个域输入面和输出面通常是物面和频谱面之间的傅里叶变换关系以及我们对这两个面上光强分布的已知信息通过迭代来逼近真实的相位。其核心流程是一个闭合的迭代循环我们可以把它想象成一个不断自我修正的猜谜游戏。算法的基本迭代步骤一个循环如下初始猜测我们从物面输入面开始。已知物面的光强分布 I_object比如一个物体的透射率。但对于相位我们一无所知。因此我们随机生成一个初始相位猜测 φ_guess通常在0到2π之间均匀分布与已知的振幅 sqrt(I_object) 结合构造一个初始的复振幅场U1 sqrt(I_object) * exp(i * φ_guess)。正向传播约束变换将构造的物面场 U1 传播到输出面通常是焦平面或衍射场。这个传播过程在傍轴近似下通过一次傅里叶变换FFT来实现U2 FFT(U1)。此时U2 同时具有振幅 A2 和相位 φ2。施加输出面约束这是我们拥有的另一个已知信息。在输出面例如透镜的后焦面我们通过实验测量得到了该处的光强分布 I_measured。GS算法的关键操作来了我们保留计算得到的相位 φ2但用测量得到的振幅 sqrt(I_measured) 替换掉计算出的振幅 A2。即更新输出面场为U2‘ sqrt(I_measured) * exp(i * φ2)。这个步骤强制让计算出的波场满足实际的物理测量结果。反向传播逆变换将更新后的输出面场 U2‘ 传播回物面。通过逆傅里叶变换IFFT实现U1‘ IFFT(U2‘)。此时我们得到了一个新的物面复振幅场包含振幅 A1‘ 和相位 φ1‘。施加输入面约束同样我们保留新计算出的相位 φ1‘但用我们最初已知的物面振幅 sqrt(I_object) 替换掉计算出的振幅 A1‘。即更新物面场为U1_new sqrt(I_object) * exp(i * φ1‘)。这个步骤强制让波场满足物面的已知条件。至此一个迭代循环完成。U1_new 将作为下一个循环的起点重复步骤2到5。如此反复迭代每一次循环都让计算出的波场在输入面和输出面同时满足各自的振幅强度约束。理论上随着迭代次数增加计算得到的相位 φ1‘ 会收敛到真实的物面相位。注意这里描述的是最基本的“两平面”GS算法。在实际应用中根据不同的物理设置如多平面相位恢复、基于角谱的衍射传播等正向和反向传播算子可能不是简单的FFT/IFFT而是更复杂的衍射积分计算。但核心的“变换-施加振幅约束-反变换”的迭代思想是完全一致的。为什么这样迭代会收敛直观上理解这个算法在寻找一个解这个解在经过特定变换如傅里叶变换后其振幅分布能与我们在另一个平面上测量到的结果匹配。迭代过程就像在解一个巨大的方程组每次施加约束都是将解向可行解集合投影。GS算法本质是一种交替投影算法在两个约束集合满足输入面振幅的复函数集合和满足输出面振幅的复函数集合之间来回投影最终希望找到两个集合的交点即同时满足两个约束的解。3. MATLAB仿真系统设计与实现要点用MATLAB搭建这个仿真系统目的很明确在已知“真实”相位的情况下模拟光波的传播和探测过程生成仿真的“测量光强”然后假装我们不知道这个相位仅用“测量光强”和GS算法来恢复它最后与“真实”相位对比评估算法性能。这个过程让我们可以在完全可控的环境下研究算法的特性。3.1 系统模块分解一个结构清晰的仿真系统通常包含以下几个模块参数定义与初始化模块设定仿真环境的基本参数。这包括网格参数仿真区域的尺寸Lx,Ly单位米或毫米采样点数Nx,Ny。采样点数一般取2的整数次幂如2565121024以便利用FFT的高效性。网格间距dx Lx/Nx。光学参数波长lambda传播距离z如果是角谱传播透镜焦距f如果是4f系统等。目标物体定义输入面的振幅分布A_obj如一个圆形孔径、一个字母图案和真实的相位分布phi_true这是我们希望恢复的可以是一个简单的球面波相位、一个随机相位板或者一个复杂的像差分布。生成模拟的真实波场U_true A_obj .* exp(1i * phi_true)。正向传播模拟模块这个模块模拟真实物理过程即光从物面传播到探测面。根据实验光路选择合适的传播模型。对于最常见的4f滤波系统或焦面探测使用傅里叶变换U_prop fft2(U_true)或fftshift(fft2(fftshift(U_true)))注意零频分量移到中心。对于非焦面的自由空间衍射可能需要使用角谱传播或菲涅尔衍射积分。角谱传播更精确其核心是AS fft2(U_true); H exp(1i * 2*pi*z / lambda .* sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2)); U_prop ifft2(AS .* H);其中fx,fy是空间频率坐标。计算探测面的光强模拟测量结果I_measured abs(U_prop).^2。为了更真实可以在此处加入噪声如高斯噪声I_measured_noisy imnoise(I_measured, gaussian, 0, 0.01)添加方差为0.01的高斯噪声。GS算法迭代恢复模块这是系统的核心。初始化设定最大迭代次数max_iter收敛阈值threshold。生成随机的初始相位phi_guess 2*pi*rand(Nx, Ny)。构造初始估计场U_est A_obj .* exp(1i * phi_guess)。迭代循环for iter 1:max_iter正向传播至探测面U_prop_est Propagation_Model(U_est)使用与正向模拟模块相同的传播模型。施加探测面振幅约束A_prop_measured sqrt(I_measured)。U_prop_est_constrained A_prop_measured .* exp(1i * angle(U_prop_est))。angle()函数用于获取相位。反向传播回物面U_est_back Inverse_Propagation_Model(U_prop_est_constrained)。施加物面振幅约束U_est_new A_obj .* exp(1i * angle(U_est_back))。计算误差评估当前估计与约束的匹配程度。常用误差函数是探测面振幅误差error(iter) sum(sum(abs(abs(U_prop_est) - A_prop_measured).^2)) / sum(sum(A_prop_measured.^2))。判断收敛如果error(iter) threshold跳出循环。更新估计场U_est U_est_new。输出循环结束后从最终的U_est中提取恢复的相位phi_recovered angle(U_est)。结果可视化与评估模块这是验证算法效果的关键。图像显示并排显示真实相位phi_true、恢复相位phi_recovered以及两者的差值残差图。由于相位是周期性的模2π通常使用imagesc显示并注意相位包裹phase wrapping问题可能需要使用angle()函数直接显示范围 -π 到 π或使用wrapToPi等函数处理。误差分析绘制迭代误差error随迭代次数的下降曲线观察收敛速度和稳定性。定量指标计算恢复相位与真实相位之间的均方根误差RMSE、峰值信噪比PSNR或结构相似性SSIM。对于存在相位包裹的情况比较前需要先进行解包裹phase unwrapping或计算梯度。3.2 关键代码片段与实操心得下面给出一些核心步骤的MATLAB代码片段和编写时的注意事项1. 生成带有相位的测试物体% 参数定义 N 512; % 采样点数 L 10e-3; % 物理尺寸 10mm dx L/N; x linspace(-L/2, L/2, N); [X, Y] meshgrid(x, x); % 振幅物体一个圆形孔径 radius 1e-3; % 1mm半径 A_obj sqrt((X.^2 Y.^2) radius^2); % 振幅1代表透光0代表不透光 % 真实相位一个球面波相位模拟透镜叠加一些像差 lambda 632.8e-9; % He-Ne激光波长 k 2*pi/lambda; R 0.5; % 球面波曲率半径 phi_true mod(k/(2*R) * (X.^2 Y.^2) 0.5*sin(4*pi*X/L) .* cos(2*pi*Y/L), 2*pi); % 使用mod函数将其范围限制在[0, 2π)避免数值过大。实际仿真中相位值范围不重要重要的是其空间变化。 % 第二部分加入了一个简单的正弦型像差。 % 构造真实复振幅场 U_true A_obj .* exp(1i * phi_true);实操心得在定义相位phi_true时要特别注意其空间频率不能超过奈奎斯特频率即每个采样点间隔内相位变化不能超过π否则会产生混叠导致仿真失真。对于快速变化的相位需要确保采样足够密集。2. 角谱传播函数用于自由空间衍射模拟function Uout AngularSpectrumPropagation(Uin, lambda, dx, dy, z) % Uin: 输入场 % lambda: 波长 % dx, dy: 输入面采样间隔 % z: 传播距离 % Uout: 输出场 [Ny, Nx] size(Uin); % 生成空间频率坐标 fx linspace(-1/(2*dx), 1/(2*dx), Nx); fy linspace(-1/(2*dy), 1/(2*dy), Ny); [FX, FY] meshgrid(fx, fy); % 计算角谱传递函数 H exp(1i * 2*pi/lambda * z .* sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 对于频率分量满足 (lambda*fx)^2 (lambda*fy)^2 1 的部分sqrt结果为虚数对应倏逝波通常直接置零或做衰减处理。 H(sqrt((lambda*FX).^2 (lambda*FY).^2) 1) 0; % 进行角谱传播 Uin_f fft2(fftshift(Uin)); % 将输入场频谱零频移到角落便于与H相乘 Uout_f Uin_f .* H; Uout ifftshift(ifft2(Uout_f)); % 反变换并移回中心 end注意角谱传播是精确的标量衍射模型但计算时要注意传递函数H中根号内的项为负的情况对应高频倏逝波处理不当会引入数值不稳定。通常的做法是将其滤除置零如上面代码所示。这相当于一个低通滤波是物理合理的。3. GS算法迭代核心循环max_iter 200; threshold 1e-6; error_history zeros(max_iter, 1); % 初始猜测随机相位 phi_guess 2*pi * rand(N, N); U_est A_obj .* exp(1i * phi_guess); % 已知的测量面振幅模拟得到已加噪声 A_measured sqrt(I_measured_noisy); for iter 1:max_iter % 1. 正向传播 (使用角谱传播或FFT需与正向模拟一致) U_prop_est AngularSpectrumPropagation(U_est, lambda, dx, dx, z); % 2. 施加测量面振幅约束 A_prop_est abs(U_prop_est); U_prop_est_constrained A_measured .* exp(1i * angle(U_prop_est)); % 3. 反向传播 U_est_back AngularSpectrumPropagation(U_prop_est_constrained, lambda, dx, dx, -z); % 注意距离取负 % 4. 施加物面振幅约束 U_est_new A_obj .* exp(1i * angle(U_est_back)); % 5. 计算误差测量面振幅误差 error_history(iter) sum(sum((A_prop_est - A_measured).^2)) / sum(sum(A_measured.^2)); % 6. 判断收敛 if error_history(iter) threshold fprintf(在 %d 次迭代后收敛。\n, iter); break; end % 7. 更新估计 U_est U_est_new; end phi_recovered angle(U_est);实操心得反向传播函数理论上应该是正向传播的逆运算。对于角谱传播传递函数是exp(-i*...)即距离取负-z。对于FFT模型4f系统反向传播就是IFFT。务必保持正向模拟和迭代恢复中使用完全相同的传播模型和参数否则算法可能无法收敛到正确解。4. 影响算法性能的关键因素与调优策略GS算法原理简洁但实际应用中其性能受到多种因素影响。理解这些因素并学会调优是完成高质量相位恢复的关键。4.1 初始相位猜测GS算法本质上是一个非线性优化过程初始值会影响收敛速度和最终结果尤其是在解不唯一或存在噪声的情况下。随机相位最常用的方法。2*pi*rand()生成均匀分布的随机相位。优点是简单有助于算法跳出局部极小值。缺点是每次运行结果可能略有不同收敛所需的迭代次数不稳定。常数相位零相位zeros(N, N)。作为最简单的猜测在物体相位变化平缓且信噪比较高时可能收敛很快。但如果真实相位变化剧烈零相位猜测可能引导算法收敛到一个错误的解局部极小值。先验信息如果对物体相位有部分了解例如知道它是一个缓慢变化的像差可以将其作为初始猜测能显著加速收敛并提高准确性。策略在仿真中可以尝试多种初始猜测观察收敛结果的一致性。在实际实验中如果没有先验信息通常采用随机相位并多次独立运行算法选取误差最小的结果或对多次结果进行平均以增强鲁棒性。4.2 振幅约束的施加与支持域Support Constraint在基础GS算法中我们对物面的振幅约束是严格的A_obj。但在很多情况下我们只知道物体的大致范围支持域而不知道其内部精确的振幅分布。例如在晶体学中只知道分子的大致形状。紧支持域定义一个二值掩膜support在物体可能存在的区域为1其他区域为0。物面约束变为U_est_new (A_obj .* exp(1i*angle(U_est_back))) .* support (U_est_back .* (1-support)) * beta。这里beta是一个松弛因子0beta1用于在非支持域缓慢衰减场而不是直接置零这有助于提高收敛性这就是著名的混合输入输出HIO算法的核心思想。松弛振幅约束对于物面振幅如果不是精确已知可以施加一个范围约束而不是固定值。策略如果你的项目涉及非均匀振幅物体或支持域信息强烈建议实现HIO或其变种如OSP、RAAR算法它们比纯GS算法具有更强的抗噪声能力和更广的收敛域。4.3 噪声的影响与处理实验测量的光强I_measured必然包含噪声散粒噪声、读出噪声等。噪声会破坏振幅约束的有效性导致GS算法振荡或不收敛到正确解。现象在噪声较大时误差曲线error_history不会平稳下降到零而是在某个值附近波动。恢复的相位图中会出现明显的随机斑点或条纹。处理方法多次独立运行取平均用不同的随机初始相位多次运行GS算法将最终恢复的相位进行平均可以有效抑制随机噪声带来的波动。正则化在迭代过程中加入正则化项惩罚相位的高频变化即假设相位是平滑的。例如在每次迭代后对恢复的相位进行轻微的高斯滤波。使用更鲁棒的算法如前所述的HIO算法对噪声的容忍度比GS更高。数据预处理对测量的光强图像进行平滑滤波或降噪处理如小波去噪、非局部均值去噪但要注意不能过度平滑而损失细节。4.4 收敛性与停滞问题GS算法可能陷入停滞stagnation即误差不再下降但恢复的结果明显不正确。这通常是因为算法陷入了局部极小值。诊断观察误差曲线如果曲线在迭代几十次后早早地进入平台期且恢复的相位与真实相位相差甚远可能就是陷入了停滞。解决方案改变初始猜测换一个随机种子重新开始。引入扰动在迭代过程中每隔一定次数人为地对当前估计的相位加入一个小的随机扰动帮助算法跳出局部极小。这被称为“扰动GS”或“随机梯度”思想的一种简单实现。切换算法当GS停滞时切换到HIO算法迭代几十次然后再切换回GS往往能打破僵局。这种GS-HIO交替迭代的策略非常有效。多平面相位恢复如果条件允许采集多个不同传播距离或离焦量下的光强图像。对每个平面的光强同时施加约束大大增加了算法的约束条件能显著改善收敛性和唯一性。这就是著名的多平面迭代相位恢复技术。5. 仿真结果分析、可视化与报告撰写完成了算法实现和迭代如何呈现你的工作成果至关重要。一份优秀的课程报告或项目总结需要清晰、直观地展示仿真过程和结论。5.1 结果可视化技巧相位显示直接使用imagesc(phi_recovered)显示相位由于相位通常被包裹在[-π, π]或[0, 2π]区间会呈现出明显的2π跳变条纹。为了更直观地看到相位的连续变化需要进行相位解包裹。% 使用MATLAB的 unwrap函数适用于1D对2D效果有限或第三方解包裹算法 % 一个简单的2D解包裹思路质量图引导 % phi_wrapped angle(U_est); % 包裹相位 % 可以使用基于最小二乘的解包裹算法如 unwrapPhaseMesh (需要自己实现或找工具箱) % 或者直接显示相位梯度如x和y方向的差分梯度图不受2π跳变影响能反映相位的局部变化。 dx_phase diff(phi_recovered, 1, 2); % x方向差分 dy_phase diff(phi_recovered, 1, 1); % y方向差分对于课程作业如果真实相位本身变化范围不超过2π也可以直接比较包裹后的相位。但务必在报告中说明这一点。误差分析图迭代误差曲线semilogy(1:iter, error_history(1:iter))。使用对数坐标可以更清楚地看到误差前期的快速下降和后期的缓慢收敛。相位残差图imagesc(wrapToPi(phi_recovered - phi_true))。显示恢复相位与真实相位的差值包裹到[-π, π]。理想的恢复结果应该是一张均匀的灰度图差值接近0。剖面线对比在相位图上选取一条水平或垂直线绘制真实相位和恢复相位沿该线的变化曲线进行直接对比。line_index N/2; % 中间一行 figure; plot(x, phi_true(line_index, :), b-, LineWidth, 1.5, DisplayName, 真实相位); hold on; plot(x, phi_recovered(line_index, :), r--, LineWidth, 1.5, DisplayName, 恢复相位); xlabel(位置 (m)); ylabel(相位 (rad)); legend; title(相位剖面线对比);定量指标计算% 1. 均方根误差 (RMSE) - 需解包裹或处理跳变 phase_diff phi_recovered - phi_true; % 处理2π跳变将差值映射到[-π, π]区间 phase_diff_wrapped mod(phase_diff pi, 2*pi) - pi; rmse sqrt(mean(phase_diff_wrapped(:).^2)); fprintf(相位RMSE: %.4f rad\n, rmse); % 2. 相关系数 (Correlation Coefficient) C corrcoef(phi_true(:), phi_recovered(:)); correlation C(1,2); fprintf(相位相关系数: %.4f\n, correlation); % 3. 对于振幅恢复可以计算相对误差 amp_recovered abs(U_est); amp_error norm(amp_recovered(:) - A_obj(:)) / norm(A_obj(:)); fprintf(振幅相对误差: %.4f\n, amp_error);5.2 探索性仿真实验设计为了深入理解算法不要只满足于一个成功运行的例子。设计一系列对比实验能让你的报告脱颖而出噪声水平影响实验固定其他参数改变添加到I_measured上的高斯噪声方差例如从0到0.1观察恢复相位的RMSE和迭代收敛曲线如何变化。绘制“噪声水平-恢复误差”关系图。初始猜测敏感性实验对同一个物体和噪声水平分别用随机相位、零相位和一种有偏的初始猜测例如使用一个离焦的球面波相位运行GS算法。比较它们的最终误差、收敛速度和恢复结果的主观质量。支持域约束实验如果物体是一个简单形状如矩形尝试在物面约束中引入支持域。比较使用精确振幅约束、宽松支持域约束HIO算法以及无支持域信息时算法的恢复能力。可以尝试恢复一个在支持域外也有微弱信号的物体。欠采样与混叠实验故意减少采样点数N或者增加相位物体的空间频率使其超过奈奎斯特频率。观察混叠如何导致恢复失败并解释原因。5.3 项目报告与代码组织心得最后将你的工作整理成一份可交付的作业或项目报告。代码组织建议使用脚本.m脚本作为主程序调用多个函数文件.m函数。将核心功能模块化create_test_object.m,forward_propagation.m,gs_algorithm.m,hio_algorithm.m,plot_results.m。使用清晰的变量名添加充分的注释。特别是对于光学参数和坐标定义注释其物理意义和单位。在主脚本开头集中定义所有可调参数方便进行不同的实验。使用savefig或exportgraphics将生成的图表高质量保存便于插入报告。报告撰写要点引言简要说明相位恢复的意义、GS算法的背景和本项目目标。原理清晰阐述GS算法的数学原理和迭代步骤最好配以流程图。方法详细描述你的MATLAB仿真模型包括光学设置、参数、噪声模型以及你实现的GS算法细节是否包含HIO如何处理支持域。结果与分析这是核心部分。展示关键实验结果图原始物体、测量光强、恢复相位、误差曲线、剖面对比等。结合图表分析算法性能讨论前面提到的关键因素初始猜测、噪声、支持域的影响。展示你设计的探索性实验的结果和结论。讨论与结论总结GS算法的优缺点你在实现过程中遇到的主要挑战和解决方案以及对未来改进的设想如尝试更先进的算法Fienup算法、PIE、ePIE等。参考文献引用关键的教材、论文或网络资源。通过这样一个从理论到实践从实现到分析从基础到拓展的完整过程你完成的就不仅仅是一个课程作业而是一个扎实的光学计算研究小课题。这份经验对你理解更复杂的计算成像算法乃至从事相关领域的研究或开发工作都将是一块重要的基石。本文还有配套的精品资源点击获取