简介GS算法Gerchberg-Saxton算法是由Gerchberg与Saxton于1972年提出的经典相位恢复迭代方法主要面向光学成像、数字全息、X射线衍射等需要从强度信息重建波前相位的场景适合正在学习计算光学或傅里叶光学、需要MATLAB可运行范例的研究生与工程师。压缩包内共2个文件包含.m实现源码与.png相位恢复效果示意图整体仅53KB结构轻盈源码中完整呈现物域与频域交替迭代的核心流程。目前已有2135人学习下载适合快速入门与二次开发。对照代码可逐步理解初始化相位分布、施加强度约束、傅里叶正逆变换更新等关键步骤同时借助效果图直观检验恢复精度与收敛表现免去从零推导的繁琐。在此基础上还可扩展硬约束、多平面迭代或与遗传算法等改进策略结合为课程实验、课题预研或实际系统搭建提供便捷起点。1. GS算法只靠振幅找回相位的傅里叶迭代法相机传感器只能记录光强相位在曝光那一刻就被丢弃可计算全息、相干衍射成像、激光整形这些应用偏偏需要把相位“造”出来。GS算法Gerchberg-Saxton是解决这类相位问题最经典的迭代方案把复振幅在物面和靶面之间来回做傅里叶变换每经过一个平面就用已知的振幅约束替换一次振幅、只保留相位循环几百次后得到一个能同时满足两个平面约束的相位分布。整个过程不需要先验样本、不需要训练数据几十行 NumPy 就能从随机相位初始化跑出一个可用结果适合刚接触相位恢复的工程师也适合要评估“标准迭代还能不能压榨出精度”的熟手。2. GS算法为什么能工作傅里叶变换对、交替投影与收敛边界2.1 相位问题如何变成两个集合上的投影光场在任意横截面上都能写成u(x) A(x) exp(i φ(x))的复振幅形式其中 A 是振幅、φ 是相位。在夫琅禾费衍射条件下从物面传播到靶面就是一个二维傅里叶变换U(f) B(f) exp(i ψ(f))。一般场景里已知物面振幅 A激光器输出是高斯分布、空间光调制器入射近似均匀平面波和靶面振幅 B目标光斑或衍射图案灰度开方φ 和 ψ 都未知这就是经典“相位问题”。GS算法把“振幅已知”看成约束集合物面约束要求复振幅满足振幅等于 A靶面约束要求频谱振幅等于 B。一次完整迭代就是从物面集合投影到靶面集合再投影回物面集合。投影操作本身不复杂对复振幅做 FFT把谱的振幅替换成 B、相位原样保留再做 IFFT然后把物面振幅替换成 A、相位继续保留。两个约束被交替满足靶面误差替换后振幅与目标振幅的范数差单调不增。要说明收敛性的边界交替投影在凸集合上保证收敛到交集但“固定振幅”的集合是非凸的GS算法只保证误差不增加不保证找到全局最优解。这就是为什么初始相位用随机值、换随机种子会得到不同结果——不是代码写错而是问题本身让迭代停在某个局部解上。我一般保留目标图案的整体平移、缩放自由度不参与评价只看靶面强度分布是否达标。2.2 物面约束和靶面约束的实际对应关系不同光学任务里“物面”和“靶面”的角色不完全一样用 GS 前先要把这两个平面的振幅约束对应清楚应用场景物面已知振幅靶面已知振幅传播算子计算全息CGHSLM 入射光振幅常取均匀 1目标衍射图案灰度开方FFT夫琅禾费激光光束整形高斯光束振幅分布平顶、环形等目标光强开方FFT透镜焦平面相干衍射成像CDI样品支撑域 support 内的透过率振幅探测器记录的衍射强度开方FFT 过采样自由空间传播入射面复振幅已知部分传播距离 z 后的强度开方角谱法 / 菲涅耳传播核这些场景共用同一个 GS 迭代框架差别只在传播算子。计算全息和焦平面整形直接走 FFT 就够了自由空间传播则需要把 FFT 换成角谱传递函数H(f) exp(i 2π z sqrt(1/λ² - f_x² - f_y²))。换算子时振幅替换逻辑完全不用动牢牢记住“正变换、替换振幅、逆变换、替换振幅”四步即可。3. 用 Python 复现 GS 算法FFT 归一化、迭代循环与首组参数3.1 最小可运行代码与关键参数说明GS 算法用 NumPy 实现只有二十行多点核心是 FFT 的能量归一化。np.fft.fft2默认不归一化正变换后频域能量是空域能量的 N 倍N 为像素总数所以靶面目标振幅要乘以sqrt(N)才能和实际频谱振幅对齐这个因子漏掉会让误差曲线一直高位震荡。import numpy as np def gs(target_amp, object_amp, iterations300, alpha1.0, seed42): rng np.random.default_rng(seed) phase rng.uniform(0.0, 2.0 * np.pi, sizetarget_amp.shape) field object_amp * np.exp(1j * phase) # 物面复振幅初始化 norm np.sqrt(field.size) # FFT 能量放大因子 errors [] for _ in range(iterations): # 物面 - 靶面 spec np.fft.fftshift(np.fft.fft2(field)) target_freq target_amp * norm # 按 FFT 缩放对齐目标振幅 amp np.abs(spec) new_amp alpha * target_freq (1.0 - alpha) * amp spec new_amp * np.exp(1j * np.angle(spec)) err np.linalg.norm(new_amp - target_freq) / np.linalg.norm(target_freq) errors.append(err) # 靶面 - 物面 field np.fft.ifft2(np.fft.ifftshift(spec)) field object_amp * np.exp(1j * np.angle(field)) return np.angle(field), errors这段代码的参数选择直接决定收敛行为iterations默认 300简单圆斑图案 100 次就够复杂灰阶图案给到 500 以上alpha是振幅替换混合系数1.0 是标准 GS 纯替换0.50.9 可以压制早期震荡代价是高频细节变软seed控制随机初始相位固定它可以复现结果换 seed 相当于换一个局部最优解object_amp是物面振幅约束数组尺寸必须与target_amp完全一致通常是全 1、高斯分布或 support 掩模。fftshift的作用是把零频移到数组中心让目标图案放在网格中央对应光轴附近可视化直观。如果不用 shift目标图案设计在数组左上角也可以工作但实际靶面图案相对光轴会整体平移物理含义不够直观。用 shift 时记得逆变换前要ifftshift还原排列顺序否则相位分布会附加一个线性相位倾斜。3.2 跑通第一个圆形光斑用一个 128×128 的网格验证最小闭环目标图案是中心半径 32 像素的圆物面给均匀振幅模拟平面波入射到纯相位元件上。y, x np.mgrid[0:128, 0:128] target (((x - 64) ** 2 (y - 64) ** 2) 32 ** 2).astype(float) object_amp np.ones((128, 128)) phase, errors gs(target, object_amp, iterations200, alpha1.0, seed0) print(final error:, errors[-1])输出误差一般在 0.05 以下把np.exp(1j * phase)再做一次 FFT 得到的靶面强度就是目标圆斑。这里注意我得到的phase是物面调制相位实际光学系统中把它加载到空间光调制器上入射均匀振幅光场后透镜焦平面即可再现目标图案。若目标图案包含明亮背景和精细文字建议把文字渲染成灰度图后先开方再送入target_amp因为 GS 约束的是振幅目标灰度直接作为振幅会导致暗区被过度放大。4. GS算法收敛诊断相关系数、能量归一化与停滞处理4.1 用相关系数判断重建质量不能只盯着误差GS 的误差曲线单调下降但降到平台后就不再动平台高低与目标图案的高频成分占比直接相关。判断结果是否真正可用我一般不看 MSE而是看靶面强度与目标图案的结构相关系数。原因很简单一个所有像素亮度都偏暗 30% 的重建结果MSE 很大但形状完全正确而相关系数能扣掉整体缩放更准确反映“像不像”。def amp_correlation(a, b): a a - a.mean() b b - b.mean() denom np.linalg.norm(a) * np.linalg.norm(b) if denom 0: return 0.0 return float(np.sum(a * b) / denom)使用时把重放强度np.abs(np.fft.fft2(np.exp(1j * phase))) ** 2与目标强度都归一化到 [0,1] 再比较或者统一乘以相同缩放因子。相关系数超过 0.95 对大多数计算全息任务就够用这时再增加迭代次数收益很小应该去调初始相位或改用加权策略。4.2 参数影响与调试方向速查下面这张表总结了我调 GS 算法时最常动的参数和对应的排查思路参数典型取值范围主要影响出问题时先查什么iterations100500收敛进度与平台期位置误差 50 次后是否还有下降趋势alpha0.61.0迭代稳定性与高频细节误差震荡时调小到 0.7seed01000局部最优解质量跑 5 个 seed 选相关系数最高者网格尺寸1281024靶面采样率与混叠出现周期鬼影时加倍网格object_amp 复杂度均匀/高斯/support可用自由度总量物面像素数不宜小于靶面的 2 倍其中物面自由度是容易忽略的一个坑物面振幅固定后真正可调的只剩每像素一个相位自由度要满足靶面 N 个强度约束物面像素数至少得是靶面采样点数的 2 倍否则系统自由度不足误差平台会明显抬高。现象就是目标图案细节越多、重建越糊这不是迭代不够而是物理上就不可能。4.3 两次停滞时最该检查的环节第一次检查能量归一化。如果忘了乘以sqrt(N)每次靶面振幅替换都会把频谱压到目标值附近一个系统性偏小的水平误差后期会进入高频小幅度震荡而不是平滑下降。最快验证方式是打印循环中np.abs(spec).mean()与target_amp.mean() * norm两者应接近。用了np.fft.fft2(..., normortho)的话norm要改成 1这个细节很容易在复制代码时被忽略。第二次检查物面约束与目标图案的尺寸匹配。物面只有 64×64、目标图案却要求细腻的 256 级灰阶文字自由度差太多。常见处理是增大物面分辨率或在目标图案中叠加随机相位扩散把单一目标点变成一个弥散斑以增加有效采样面积。若两者都不匹配任何 GS 变体都救不回来。5. 计算全息中的加权 GS算法均匀性权重与最后一公里调试标准 GS 在计算全息里有一个典型问题重建图案亮区过亮、暗区过暗整体均匀性差。这是因为迭代只盯着振幅误差最小化没有考虑人眼或探测器对暗区误差更敏感。业内常用加权 GSWeighted GS做补偿做法是在每次靶面振幅替换前根据上一轮实际达到的强度分布更新一个逐像素权重偏暗的像素加大权重、偏亮的像素减小权重让算法把更多能量赶去暗区。target_intensity target_amp ** 2 # 目标强度分布 weight np.ones_like(target_intensity) # 均匀性权重初始为 1 gamma 0.8 # 反馈系数典型 0.5~1.0 for _ in range(iterations): spec np.fft.fftshift(np.fft.fft2(field)) achieved np.abs(spec) ** 2 # 当前重建强度 weight * (target_intensity / achieved.clip(1e-12)) ** gamma weight np.clip(weight, 1e-3, 1e3) # 防止个别像素权重爆炸 target_amp_weighted np.sqrt(weight * target_intensity) * norm spec target_amp_weighted * np.exp(1j * np.angle(spec)) field np.fft.ifft2(np.fft.ifftshift(spec)) field object_amp * np.exp(1j * np.angle(field))关键在权重的更新方向target_intensity / achieved大于 1 表示该像素偏暗权重变大下一轮目标振幅被抬高从而驱动更多能量流向这里。gamma控制反馈力度——太小均匀化速度慢太大亮区暗区交替过冲误差曲线出现周期性波动我通常从 0.8 起步。限制权重范围是必要的否则个别初始相位极暗的像素会在一百轮内把权重推到 10⁶ 量级整个相位分布被它带走。这个加权版本只在标准 GS 循环里插入了五行权重更新却能把靶面强度不均匀度从 ±30% 压到 ±5% 左右。配合第 4 章的相关系数评估跑 35 个随机初始化、选相关系数最高的一轮结果是日常做计算全息相位设计比较稳的工作流。本文还有配套的精品资源点击获取 GS相位恢复算法：傅里叶迭代与计算全息Python实战
简介GS算法Gerchberg-Saxton算法是由Gerchberg与Saxton于1972年提出的经典相位恢复迭代方法主要面向光学成像、数字全息、X射线衍射等需要从强度信息重建波前相位的场景适合正在学习计算光学或傅里叶光学、需要MATLAB可运行范例的研究生与工程师。压缩包内共2个文件包含.m实现源码与.png相位恢复效果示意图整体仅53KB结构轻盈源码中完整呈现物域与频域交替迭代的核心流程。目前已有2135人学习下载适合快速入门与二次开发。对照代码可逐步理解初始化相位分布、施加强度约束、傅里叶正逆变换更新等关键步骤同时借助效果图直观检验恢复精度与收敛表现免去从零推导的繁琐。在此基础上还可扩展硬约束、多平面迭代或与遗传算法等改进策略结合为课程实验、课题预研或实际系统搭建提供便捷起点。1. GS算法只靠振幅找回相位的傅里叶迭代法相机传感器只能记录光强相位在曝光那一刻就被丢弃可计算全息、相干衍射成像、激光整形这些应用偏偏需要把相位“造”出来。GS算法Gerchberg-Saxton是解决这类相位问题最经典的迭代方案把复振幅在物面和靶面之间来回做傅里叶变换每经过一个平面就用已知的振幅约束替换一次振幅、只保留相位循环几百次后得到一个能同时满足两个平面约束的相位分布。整个过程不需要先验样本、不需要训练数据几十行 NumPy 就能从随机相位初始化跑出一个可用结果适合刚接触相位恢复的工程师也适合要评估“标准迭代还能不能压榨出精度”的熟手。2. GS算法为什么能工作傅里叶变换对、交替投影与收敛边界2.1 相位问题如何变成两个集合上的投影光场在任意横截面上都能写成u(x) A(x) exp(i φ(x))的复振幅形式其中 A 是振幅、φ 是相位。在夫琅禾费衍射条件下从物面传播到靶面就是一个二维傅里叶变换U(f) B(f) exp(i ψ(f))。一般场景里已知物面振幅 A激光器输出是高斯分布、空间光调制器入射近似均匀平面波和靶面振幅 B目标光斑或衍射图案灰度开方φ 和 ψ 都未知这就是经典“相位问题”。GS算法把“振幅已知”看成约束集合物面约束要求复振幅满足振幅等于 A靶面约束要求频谱振幅等于 B。一次完整迭代就是从物面集合投影到靶面集合再投影回物面集合。投影操作本身不复杂对复振幅做 FFT把谱的振幅替换成 B、相位原样保留再做 IFFT然后把物面振幅替换成 A、相位继续保留。两个约束被交替满足靶面误差替换后振幅与目标振幅的范数差单调不增。要说明收敛性的边界交替投影在凸集合上保证收敛到交集但“固定振幅”的集合是非凸的GS算法只保证误差不增加不保证找到全局最优解。这就是为什么初始相位用随机值、换随机种子会得到不同结果——不是代码写错而是问题本身让迭代停在某个局部解上。我一般保留目标图案的整体平移、缩放自由度不参与评价只看靶面强度分布是否达标。2.2 物面约束和靶面约束的实际对应关系不同光学任务里“物面”和“靶面”的角色不完全一样用 GS 前先要把这两个平面的振幅约束对应清楚应用场景物面已知振幅靶面已知振幅传播算子计算全息CGHSLM 入射光振幅常取均匀 1目标衍射图案灰度开方FFT夫琅禾费激光光束整形高斯光束振幅分布平顶、环形等目标光强开方FFT透镜焦平面相干衍射成像CDI样品支撑域 support 内的透过率振幅探测器记录的衍射强度开方FFT 过采样自由空间传播入射面复振幅已知部分传播距离 z 后的强度开方角谱法 / 菲涅耳传播核这些场景共用同一个 GS 迭代框架差别只在传播算子。计算全息和焦平面整形直接走 FFT 就够了自由空间传播则需要把 FFT 换成角谱传递函数H(f) exp(i 2π z sqrt(1/λ² - f_x² - f_y²))。换算子时振幅替换逻辑完全不用动牢牢记住“正变换、替换振幅、逆变换、替换振幅”四步即可。3. 用 Python 复现 GS 算法FFT 归一化、迭代循环与首组参数3.1 最小可运行代码与关键参数说明GS 算法用 NumPy 实现只有二十行多点核心是 FFT 的能量归一化。np.fft.fft2默认不归一化正变换后频域能量是空域能量的 N 倍N 为像素总数所以靶面目标振幅要乘以sqrt(N)才能和实际频谱振幅对齐这个因子漏掉会让误差曲线一直高位震荡。import numpy as np def gs(target_amp, object_amp, iterations300, alpha1.0, seed42): rng np.random.default_rng(seed) phase rng.uniform(0.0, 2.0 * np.pi, sizetarget_amp.shape) field object_amp * np.exp(1j * phase) # 物面复振幅初始化 norm np.sqrt(field.size) # FFT 能量放大因子 errors [] for _ in range(iterations): # 物面 - 靶面 spec np.fft.fftshift(np.fft.fft2(field)) target_freq target_amp * norm # 按 FFT 缩放对齐目标振幅 amp np.abs(spec) new_amp alpha * target_freq (1.0 - alpha) * amp spec new_amp * np.exp(1j * np.angle(spec)) err np.linalg.norm(new_amp - target_freq) / np.linalg.norm(target_freq) errors.append(err) # 靶面 - 物面 field np.fft.ifft2(np.fft.ifftshift(spec)) field object_amp * np.exp(1j * np.angle(field)) return np.angle(field), errors这段代码的参数选择直接决定收敛行为iterations默认 300简单圆斑图案 100 次就够复杂灰阶图案给到 500 以上alpha是振幅替换混合系数1.0 是标准 GS 纯替换0.50.9 可以压制早期震荡代价是高频细节变软seed控制随机初始相位固定它可以复现结果换 seed 相当于换一个局部最优解object_amp是物面振幅约束数组尺寸必须与target_amp完全一致通常是全 1、高斯分布或 support 掩模。fftshift的作用是把零频移到数组中心让目标图案放在网格中央对应光轴附近可视化直观。如果不用 shift目标图案设计在数组左上角也可以工作但实际靶面图案相对光轴会整体平移物理含义不够直观。用 shift 时记得逆变换前要ifftshift还原排列顺序否则相位分布会附加一个线性相位倾斜。3.2 跑通第一个圆形光斑用一个 128×128 的网格验证最小闭环目标图案是中心半径 32 像素的圆物面给均匀振幅模拟平面波入射到纯相位元件上。y, x np.mgrid[0:128, 0:128] target (((x - 64) ** 2 (y - 64) ** 2) 32 ** 2).astype(float) object_amp np.ones((128, 128)) phase, errors gs(target, object_amp, iterations200, alpha1.0, seed0) print(final error:, errors[-1])输出误差一般在 0.05 以下把np.exp(1j * phase)再做一次 FFT 得到的靶面强度就是目标圆斑。这里注意我得到的phase是物面调制相位实际光学系统中把它加载到空间光调制器上入射均匀振幅光场后透镜焦平面即可再现目标图案。若目标图案包含明亮背景和精细文字建议把文字渲染成灰度图后先开方再送入target_amp因为 GS 约束的是振幅目标灰度直接作为振幅会导致暗区被过度放大。4. GS算法收敛诊断相关系数、能量归一化与停滞处理4.1 用相关系数判断重建质量不能只盯着误差GS 的误差曲线单调下降但降到平台后就不再动平台高低与目标图案的高频成分占比直接相关。判断结果是否真正可用我一般不看 MSE而是看靶面强度与目标图案的结构相关系数。原因很简单一个所有像素亮度都偏暗 30% 的重建结果MSE 很大但形状完全正确而相关系数能扣掉整体缩放更准确反映“像不像”。def amp_correlation(a, b): a a - a.mean() b b - b.mean() denom np.linalg.norm(a) * np.linalg.norm(b) if denom 0: return 0.0 return float(np.sum(a * b) / denom)使用时把重放强度np.abs(np.fft.fft2(np.exp(1j * phase))) ** 2与目标强度都归一化到 [0,1] 再比较或者统一乘以相同缩放因子。相关系数超过 0.95 对大多数计算全息任务就够用这时再增加迭代次数收益很小应该去调初始相位或改用加权策略。4.2 参数影响与调试方向速查下面这张表总结了我调 GS 算法时最常动的参数和对应的排查思路参数典型取值范围主要影响出问题时先查什么iterations100500收敛进度与平台期位置误差 50 次后是否还有下降趋势alpha0.61.0迭代稳定性与高频细节误差震荡时调小到 0.7seed01000局部最优解质量跑 5 个 seed 选相关系数最高者网格尺寸1281024靶面采样率与混叠出现周期鬼影时加倍网格object_amp 复杂度均匀/高斯/support可用自由度总量物面像素数不宜小于靶面的 2 倍其中物面自由度是容易忽略的一个坑物面振幅固定后真正可调的只剩每像素一个相位自由度要满足靶面 N 个强度约束物面像素数至少得是靶面采样点数的 2 倍否则系统自由度不足误差平台会明显抬高。现象就是目标图案细节越多、重建越糊这不是迭代不够而是物理上就不可能。4.3 两次停滞时最该检查的环节第一次检查能量归一化。如果忘了乘以sqrt(N)每次靶面振幅替换都会把频谱压到目标值附近一个系统性偏小的水平误差后期会进入高频小幅度震荡而不是平滑下降。最快验证方式是打印循环中np.abs(spec).mean()与target_amp.mean() * norm两者应接近。用了np.fft.fft2(..., normortho)的话norm要改成 1这个细节很容易在复制代码时被忽略。第二次检查物面约束与目标图案的尺寸匹配。物面只有 64×64、目标图案却要求细腻的 256 级灰阶文字自由度差太多。常见处理是增大物面分辨率或在目标图案中叠加随机相位扩散把单一目标点变成一个弥散斑以增加有效采样面积。若两者都不匹配任何 GS 变体都救不回来。5. 计算全息中的加权 GS算法均匀性权重与最后一公里调试标准 GS 在计算全息里有一个典型问题重建图案亮区过亮、暗区过暗整体均匀性差。这是因为迭代只盯着振幅误差最小化没有考虑人眼或探测器对暗区误差更敏感。业内常用加权 GSWeighted GS做补偿做法是在每次靶面振幅替换前根据上一轮实际达到的强度分布更新一个逐像素权重偏暗的像素加大权重、偏亮的像素减小权重让算法把更多能量赶去暗区。target_intensity target_amp ** 2 # 目标强度分布 weight np.ones_like(target_intensity) # 均匀性权重初始为 1 gamma 0.8 # 反馈系数典型 0.5~1.0 for _ in range(iterations): spec np.fft.fftshift(np.fft.fft2(field)) achieved np.abs(spec) ** 2 # 当前重建强度 weight * (target_intensity / achieved.clip(1e-12)) ** gamma weight np.clip(weight, 1e-3, 1e3) # 防止个别像素权重爆炸 target_amp_weighted np.sqrt(weight * target_intensity) * norm spec target_amp_weighted * np.exp(1j * np.angle(spec)) field np.fft.ifft2(np.fft.ifftshift(spec)) field object_amp * np.exp(1j * np.angle(field))关键在权重的更新方向target_intensity / achieved大于 1 表示该像素偏暗权重变大下一轮目标振幅被抬高从而驱动更多能量流向这里。gamma控制反馈力度——太小均匀化速度慢太大亮区暗区交替过冲误差曲线出现周期性波动我通常从 0.8 起步。限制权重范围是必要的否则个别初始相位极暗的像素会在一百轮内把权重推到 10⁶ 量级整个相位分布被它带走。这个加权版本只在标准 GS 循环里插入了五行权重更新却能把靶面强度不均匀度从 ±30% 压到 ±5% 左右。配合第 4 章的相关系数评估跑 35 个随机初始化、选相关系数最高的一轮结果是日常做计算全息相位设计比较稳的工作流。本文还有配套的精品资源点击获取