迭代傅里叶变换算法IFTA：从原理到相位片工程实战

简介：IFTA（迭代傅里叶变换）算法的MATLAB实现资源，面向图像处理与信号处理方向的学生和研究人员，聚焦图像复原、去噪与频谱分析等典型应用。压缩包共含2个文件，一个可直接运行的MATLAB脚本与一幅经典标准测试图像，整体体积仅57KB。代码以快速傅里叶变换为基础，通过迭代方式逐步逼近处理结果，流程覆盖读取图像并转为复数形式、频谱获取、频域迭代处理（如频率选择、阈值调整）、逆变换回空间域以及结果可视化，便于观察算法对频域特性与最终图像效果的影响；代码结构清晰，便于按需修改与调试，适合作为IFTA实现研究的参考。读者可使用配套测试图直接复现算法运行，比较不同迭代阶段的变化，也可在此基础上修改关键策略，将其扩展到更复杂的图像处理任务。已有307人学习参考，是理解迭代傅里叶变换原理并动手实践的轻量级入门资源。

1. 迭代傅里叶不是“傅里叶变换”，而是把目标振幅和相位约束来回投影

拿到一份名为 IFTA.zip 的打包文件，里面通常不是某个新算法，而是迭代傅里叶变换算法（IFTA）的脚本和说明。IFTA 的做法不是改进 FFT，而是在空域和频域之间来回做傅里叶变换，每走一步就强制替换振幅、保留相位，反复逼近目标光强分布。它解决的核心问题是：给定一个期望的衍射图案，反算出相位片或全息图的相位分布。适合做光刻光束整形、激光加工、投影显示和结构光照明的工程师。如果只把 IFTA 当成“更快的傅里叶变换”来读，后面所有参数都会调不明白，因为真正的难点不在变换本身，而在约束的构造方式。

2. IFTA 的迭代骨架：正变换、逆变换与振幅约束的往返

IFTA 家族里最基础、也最容易被重现的是 Gerchberg-Saxton（GS）算法。它的迭代骨架只有四个动作：从谱面做逆傅里叶变换到物面、替换物面振幅、做正傅里叶变换回谱面、替换谱面振幅。所谓“替换”，就是只用目标值替换振幅，相位保持不动。这个流程看起来简单，但它背后的交替投影思想，才是整个迭代傅里叶能收敛的原因。

2.1 正变换与逆变换各自约束什么

在夫琅禾费衍射近似下，透镜后焦面的复振幅与输入面复振幅之间满足二维傅里叶变换关系。这里的物面是目标光强所在的输出平面，谱面是你要加工的相位片平面。IFTA 的每一次迭代，等于在问：如果输出面振幅是目标振幅，输入面应该长什么样？但傅里叶变换成对出现，输入面的振幅也会反过来影响输出面，所以只能在两个域之间反复投影。

正变换（物面到谱面）使用fft2，它把目标振幅和相位拼成的复振幅变换到谱面；谱面约束通常是“纯相位”，也就是把振幅改成 1，只保留相位。逆变换（谱面到物面）使用ifft2，再把谱面相位还原成物面复振幅；物面约束是“目标振幅”，也就是把振幅替换成目标值，保留当前相位。两轮约束合在一起，就是在交替投影到两个集合的交集上。

2.2 最小可运行的 GS 实现

import numpy as np def ifta_gs(target_amp, n_iter=50, seed=0): # target_amp: 目标振幅，二维 float 数组，范围 [0, 1] rng = np.random.default_rng(seed) phase = rng.uniform(-np.pi, np.pi, target_amp.shape) # 谱面初始化为单位振幅加随机相位 spectrum = np.exp(1j * phase) for _ in range(n_iter): # 逆变换到物面 field = np.fft.ifft2(spectrum) # 物面约束：替换振幅，保留相位 field = target_amp * np.exp(1j * np.angle(field)) # 正变换回谱面 spectrum = np.fft.fft2(field) # 谱面约束：纯相位 spectrum = np.exp(1j * np.angle(spectrum)) return np.angle(spectrum)

这段代码是 GS 算法最直接的形态，能跑，但收敛速度慢。注意target_amp的单位：它必须是振幅，而不是 CCD 上读到的强度。如果输入的是灰度图，要先把灰度值开根号再传进来，否则最后衍射出来的是目标图案灰度值的平方。n_iter一般取 20 到 100，超过 200 轮收益很小；seed固定下来，方便对比不同参数时排除随机相位的影响。

关于fftshift：理论上，ifft2后坐标原点在数组左上角，目标图案放在中心时，相位片会带一个整体相位倾斜。省掉fftshift不影响收敛，但会影响生成的相位片的周期边界。输出相位片给加工厂时，相位值要能首尾相接，我一般会在最后额外做一次np.fft.fftshift把低频挪到中心，再取angle。

2.3 迭代停止条件：收敛曲线与误差指标

迭代停止不能只看迭代次数，还要看误差。常用指标有三个：归一化均方根误差（NRMSE）、衍射效率、均匀性。

def nrmse(field_amp, target_amp): # 去掉整体能量比例后比较形状差异 a = field_amp / (field_amp.sum() + 1e-12) t = target_amp / (target_amp.sum() + 1e-12) return np.sqrt(np.mean((a - t) ** 2)) / (np.sqrt(np.mean(t ** 2)) + 1e-12)

这里先做能量归一化，是为了避免整体亮度差异主导误差。分母加1e-12防止除零。使用时要在输出平面恢复出的振幅上计算，而不是对傅里叶变换前的复振幅计算。每次迭代后记录这轮的 NRMSE，可以看到曲线先快速下降，再变平；如果曲线锯齿严重，说明约束之间冲突过大，需要调整后面提到的反馈系数。

提示：NRMSE 只能作为参考，不能完全代表实际光学效果。同一个 NRMSE 下，可能是边缘几个像素差异，也可能是低频背景整体偏亮，后者目视更明显。

有了这个最小骨架，IFTA 就能跑通。但跑通和收敛到可用解之间，隔着参数选择的距离。

3. IFTA 的 3 个必调参数：初始相位、反馈系数和衍射距离

GS 类算法是迭代傅里叶家族的核心，但它不是一次参数定终身。实际工程里，我一般只会调三个参数：初始相位的随机种子、物面约束里的反馈系数、以及传播距离对应的衍射模型。调好这三个，90% 的相位片设计问题都能解决。

3.1 初始相位：随机种子决定均匀性

IFTA 要解的相位恢复问题是非凸的，不同的初始相位会收敛到不同的局部极小。GS 算法对初值敏感，不是玄学，而是数学上必然的结果。同一个目标图案，seed=0可能得到均匀性 1.5% 的解，seed=7可能得到 5% 的解。

常见做法是固定跑 3 到 5 个随机种子，选 NRMSE 或均匀性最好的结果，而不是只跑一遍。下面这段代码在迭代结束后返回过程误差，方便做初值对比：

import numpy as np def ifta_gs_with_history(target_amp, n_iter=50, seed=0): rng = np.random.default_rng(seed) phase = rng.uniform(-np.pi, np.pi, target_amp.shape) spectrum = np.exp(1j * phase) history = [] for _ in range(n_iter): field = np.fft.ifft2(spectrum) field_amp = np.abs(field) history.append(nrmse(field_amp, target_amp)) field = target_amp * np.exp(1j * np.angle(field)) spectrum = np.fft.fft2(field) spectrum = np.exp(1j * np.angle(spectrum)) return np.angle(spectrum), history

参数说明里最重要的是seed的选择逻辑：先跑小迭代数，例如 20 轮，观察历史误差，再固定表现最好的种子做完整迭代。这样可以省掉大迭代数下的算力浪费。随机相位的另一个作用是打破对称性，如果你的目标图案是中心对称图形，固定相位会导致迭代陷在对称解里，均匀性很难看。

3.2 反馈系数：从 GS 到加权 IFTA

GS 在物面约束时是硬替换：每一轮都把振幅直接改成目标值。这个做法在迭代初期收敛快，但后期会在两个约束的交界处震荡。常见的改进是引入反馈系数，让新的振幅一部分来自目标，一部分来自上一轮结果：

def ifta_feedback(target_amp, n_iter=100, beta=0.6, seed=1): rng = np.random.default_rng(seed) spectrum = np.exp(1j * rng.uniform(-np.pi, np.pi, target_amp.shape)) for _ in range(n_iter): field = np.fft.ifft2(spectrum) amp = np.abs(field) # 反馈：新的振幅 = (1 - beta) * 目标 + beta * 当前振幅 new_amp = (1 - beta) * target_amp + beta * amp field = new_amp * np.exp(1j * np.angle(field)) spectrum = np.fft.fft2(field) spectrum = np.exp(1j * np.angle(spectrum)) return np.angle(spectrum)

这里beta的取值范围是 0 到 1。beta=0是原始 GS，beta=0.6到0.8能明显抑制后期震荡。反馈系数越大，收敛越慢，但跳离局部极小的能力越强。如果你发现误差曲线后期呈等幅振荡，优先调小beta；如果曲线平得太早，说明陷入局部极小，优先调大beta。

加权 IFTA 是反馈思路的延伸：在目标区域和非目标区域使用不同的权重，让零阶或者杂散光区域的能量被惩罚。实际代码里，只需要把target_amp替换成加权的混合目标：

roi = target_amp > 0.05 weight = np.ones_like(target_amp) weight[roi] = 0.8 weight[~roi] = 1.2 # 迭代时把目标振幅改为加权版本 weighted_target = target_amp * weight

权重小于 1 的区域允许更大的振幅偏差，权重大的区域则被收紧。这个技巧在很多商用的衍射光学元件设计软件里叫“ROI weighting”。

3.3 衍射距离：选错模型会收敛到错误解

IFTA 默认基于夫琅禾费衍射，也就是fft2/ifft2直接对应透镜后焦面。如果你要设计的是近场图案，或者相位片到目标面之间存在一段自由空间传播，就必须换成角谱传播。角谱传播的代码并不复杂：

def angular_spectrum(field, d, wavelength, dx): # field: 输入复振幅，二维复数数组 # d: 传播距离，单位与 wavelength 一致 # dx: 采样间距，单位与 wavelength 一致 k = 2 * np.pi / wavelength fy = np.fft.fftfreq(field.shape[0], dx) fx = np.fft.fftfreq(field.shape[1], dx) FX, FY = np.meshgrid(fx, fy) # 角谱传递函数 H = np.exp(1j * k * d * np.sqrt(1 - (wavelength * FX) ** 2 - (wavelength * FY) ** 2)) return np.fft.ifft2(np.fft.fft2(field) * H)

这段代码里，d、wavelength、dx必须统一单位。H是角谱传递函数，它把频谱的每一个频率分量乘以不同的相位延迟。如果sqrt内出现负数，对应倏逝波成分，此时的正确做法是把该频率置零，而不是保留指数衰减项，否则数值上会溢出。

选错模型的表现很典型：误差曲线收敛得很好，但实际打出来的光斑边界模糊或位置偏移。这不是优化问题，而是你的“目标振幅”被放在了错误的参考面上。

参数讲完之后，真正让 IFTA 工程化的，是把迭代中出现的各种反直觉现象逐个排掉。

4. IFTA 跑不出好结果的 4 个坑：频谱折叠、平方根、模型边界和局部极小

很多从 GitHub 或网盘里拿到 IFTA.zip 的人，会觉得代码没错，但结果就是不对。这里的问题往往不在迭代算法本身，而在输入和模型边界。以下四个坑我几乎每次帮人调 IFTA 都会遇到。

4.1 目标频谱碰到边界：采样率不足

迭代傅里叶里的 FFT 默认是周期的，目标图案的高频成分如果超过了奈奎斯特频率，会被折返到低频，形成摩尔纹或中心光斑异常的条纹。判断方法很简单：把目标图案做一次 FFT，检查高频能量是否集中在边界。

def check_spectrum_border(field): spec = np.fft.fftshift(np.fft.fft2(field)) spec_db = 20 * np.log10(np.abs(spec) + 1e-12) # 观察数组四边的能量 h, w = spec_db.shape border = np.concatenate([spec_db[0, :], spec_db[-1, :], spec_db[:, 0], spec_db[:, -1]]) return border.max()

如果这个值只比中心峰值低 20 dB 以内，说明高频已经顶到边界。解决办法不是加密 FFT 点数，而是把目标图案缩小，让图案外至少留 1/4 到 1/2 的空白区域。零填充不是简单补零，而是让目标图案在输出面上“呼吸”的空间。

4.2 把强度图当成振幅图

IFTA 的目标振幅和 CCD 探测器读到的强度不是一回事。光强正比于振幅的平方，所以大多数测试图案拿到手时是强度图（灰度值），直接送进ifta_gs，收敛结果会是目标灰度的平方根效果，表现为图案暗部抬起、亮部过曝。

正确做法是先归一化再开根号：

gray = plt.imread("target.png").astype(np.float32) gray = gray / gray.max() target_amp = np.sqrt(gray) # 把强度转成振幅

这个转换的误差非常隐蔽，因为 NRMSE 也会同步下降，很多人误以为迭代成功。验证方法是在仿真里做一次正向衍射，把输出光强和原始灰度图直接对比，如果形状一致但对比度不对，就是振幅和强度的问题。

4.3 近场与远场的标量衍射边界

使用角谱法时，传递函数里的空间频率必须满足物理上限f < 1 / wavelength。当采样间距dx过大，或者图案高频部分过多时，(wavelength * FX) ** 2 + (wavelength * FY) ** 2会大于 1，对应倏逝波。

在数值实现中，需要对传递函数加掩膜：

mask = (1 - (wavelength * FX) ** 2 - (wavelength * FY) ** 2) > 0 H = np.zeros_like(FX) H[mask] = np.exp(1j * k * d * np.sqrt(1 - (wavelength * FX[mask]) ** 2 - (wavelength * FY[mask]) ** 2))

如果不加掩膜，sqrt得到复数值，传播结果会失控，误差曲线可能突然飙升。另一个边界问题是距离太小、采样点不足时，角谱法需要的零填充量急剧增长，一般要求传播距离至少大于采样间隔的数值，否则直接改用近场菲涅尔模型更稳定。

4.4 迭代停滞：用加权惩罚打破局部极小

GS 算法是交替投影，非凸问题决定了它一定会陷入局部极小。实际表现为误差曲线在某个值附近不再下降，加大迭代次数也没有用。除了上一章的反馈系数，更有效的办法是在目标区域里按“误差热点”加权重。误差热点通常出现在目标图案的锐利边缘，因为这些位置的高频分量受限于采样率，无法被相位片完全还原。

def apply_error_weight(target_amp, field_amp, base_weight=1.0, penalty=2.0): # 找出当前振幅偏差最大的 10% 像素 err = np.abs(field_amp - target_amp) threshold = np.quantile(err, 0.9) hot = err > threshold weight = np.ones_like(target_amp) * base_weight weight[hot] = penalty return weight

把权重乘到目标振幅上再迭代，会牺牲部分整体均匀性，但能压住视觉上最明显的边缘过曝。需要说明的是，这个操作不是免费午餐，它只是把误差从热点推到非热点区域，必须在衍射效率和均匀性之间重新做权衡。

排掉这些坑之后，IFTA 的输出已经能看了。下一步是把相位片从数组变成可加工的文件。

5. 从目标图到相位片：IFTA 实战里的 3 个验证技巧

相位片设计和普通图像滤波不一样，最后交付的是相位值，而不是灰度图。这里给三个我每次都会做的步骤，能帮你少返工。

5.1 目标准备：先做能量归一化再开根号

我习惯把目标区域限定在一个圆形或矩形 ROI 内，避免无关背景参与迭代。目标图案先用gray / gray.max()归一化，再开根号，最后把目标区域外的振幅置为 0。这样做的好处是让衍射效率有一个明确的定义基准。

5.2 用收敛曲线判断是否陷入局部极小

把 5 个 seed 的误差曲线画在同一张图上，如果某一轮曲线在相同水平附近停滞，且多 seed 结果差异明显，基本可以断定问题出在目标的高频成分太多。此时优先降频，而不是盲目调参。

5.3 相位片的 2π 周期归一化与 16 位输出

最后得到的相位值要 wrap 到 [0, 2π) 再量化输出。直接用np.angle得到 [-π, π] 会导致相位分布看起来有突变，但 2π 跳变对光学上无影响，可加工设备不一定喜欢这种表示。

phase = ifta_feedback(target_amp, n_iter=100, beta=0.7, seed=3) phase_wrapped = np.mod(phase, 2 * np.pi) # 映射到 16 位灰度 phase_16bit = (phase_wrapped / (2 * np.pi) * 65535).astype(np.uint16)

这里phase_wrapped的值域是 [0, 2π)，再做线性映射到 0 到 65535。保存为 PNG 时要注意像素类型，uint16能保留相位精度，uint8只有 256 级，量化噪声会直接变成衍射光斑里的背景杂散光。输出前也可以做一次轻微的高斯羽化，把相位片的边缘过渡抹平，能有效抑制高阶衍射，但羽化半径不要超过一个采样间隔，否则会糊掉细节。

本文还有配套的精品资源，点击获取