简介：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 前先要把这两个平面的振幅约束对应清楚：
| 应用场景 | 物面已知振幅 | 靶面已知振幅 | 传播算子 |
|---|---|---|---|
| 计算全息（CGH） | SLM 入射光振幅，常取均匀 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, iterations=300, alpha=1.0, seed=42): rng = np.random.default_rng(seed) phase = rng.uniform(0.0, 2.0 * np.pi, size=target_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.5～0.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, iterations=200, alpha=1.0, seed=0) 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 算法时最常动的参数和对应的排查思路：
| 参数 | 典型取值范围 | 主要影响 | 出问题时先查什么 |
|---|---|---|---|
| iterations | 100～500 | 收敛进度与平台期位置 | 误差 50 次后是否还有下降趋势 |
| alpha | 0.6～1.0 | 迭代稳定性与高频细节 | 误差震荡时调小到 0.7 |
| seed | 0～1000 | 局部最优解质量 | 跑 5 个 seed 选相关系数最高者 |
| 网格尺寸 | 128～1024 | 靶面采样率与混叠 | 出现周期鬼影时加倍网格 |
| object_amp 复杂度 | 均匀/高斯/support | 可用自由度总量 | 物面像素数不宜小于靶面的 2 倍 |
其中物面自由度是容易忽略的一个坑：物面振幅固定后，真正可调的只剩每像素一个相位自由度，要满足靶面 N 个强度约束，物面像素数至少得是靶面采样点数的 2 倍，否则系统自由度不足，误差平台会明显抬高。现象就是目标图案细节越多、重建越糊，这不是迭代不够，而是物理上就不可能。
4.3 两次停滞时最该检查的环节
第一次检查能量归一化。如果忘了乘以sqrt(N)，每次靶面振幅替换都会把频谱压到目标值附近一个系统性偏小的水平，误差后期会进入高频小幅度震荡而不是平滑下降。最快验证方式是打印循环中np.abs(spec).mean()与target_amp.mean() * norm，两者应接近。用了np.fft.fft2(..., norm="ortho")的话，norm要改成 1，这个细节很容易在复制代码时被忽略。
第二次检查物面约束与目标图案的尺寸匹配。物面只有 64×64、目标图案却要求细腻的 256 级灰阶文字，自由度差太多。常见处理是增大物面分辨率或在目标图案中叠加随机相位扩散（把单一目标点变成一个弥散斑）以增加有效采样面积。若两者都不匹配，任何 GS 变体都救不回来。
5. 计算全息中的加权 GS算法：均匀性权重与最后一公里调试
标准 GS 在计算全息里有一个典型问题：重建图案亮区过亮、暗区过暗，整体均匀性差。这是因为迭代只盯着振幅误差最小化，没有考虑人眼或探测器对暗区误差更敏感。业内常用加权 GS（Weighted 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 章的相关系数评估，跑 3～5 个随机初始化、选相关系数最高的一轮结果，是日常做计算全息相位设计比较稳的工作流。
本文还有配套的精品资源，点击获取
