波前重构技术

1. 点扩散函数（PSF）总能量与光瞳、波前的严格理论推导

1.1. 问题定义与物理模型

1.1.1. 系统描述

在光学成像系统中，我们考虑一个理想的衍射受限系统（或带有像差的系统）：

• 光瞳掩模 $M(x,y)$：描述光瞳几何形状的函数，在通光孔径内部为 $1$，外部为 $0$。
  • 对于圆形光瞳，$M(x,y) = \begin{cases} 1, & x^2+y^2 \leq R^2 \\ 0, & \text{otherwise} \end{cases}$
• 波前（像差） $W(x,y)$ 或记为 $B$：描述光瞳平面上的相位分布，单位通常为波长 $\lambda$ 或弧度。
• 点扩散函数（PSF） $I(u,v)$：理想点光源经过光学系统后在像平面上形成的强度分布。

1.1.2. 光瞳函数（Pupil Function）

光瞳平面上的复振幅透过率函数（简称光瞳函数）定义为：

$\boxed{P(x,y) = M(x,y) \cdot e^{i k W(x,y)}}$

其中：

• $k = \frac{2\pi}{\lambda}$ 为波数
• $M(x,y)$ 为振幅调制（这里只有 $0$ 或 $1$）
• $e^{i k W(x,y)}$ 为相位调制

注意：波前 $W(x,y)$ 可以是任意函数（离焦、像差、随机相位等），只要它是实函数。

1.2. 连续域严格推导

1.2.1. 从光瞳到像平面

根据标量衍射理论和夫琅禾费近似，像平面（焦平面）上的复振幅分布 $U(u,v)$ 是光瞳函数 $P(x,y)$ 的傅里叶变换（相差一个常数相位因子和坐标缩放）：

$U(u,v) = \mathcal{F}\{P(x,y)\} = \iint_{-\infty}^{+\infty} P(x,y) \, e^{-i 2\pi (xu + yv)} \, dx \, dy$

像平面上的强度分布（即 PSF）为：

$I(u,v) = |U(u,v)|^2 = U(u,v) \cdot U^*(u,v)$

1.2.2. PSF 总能量的定义

PSF 的总能量（或总强度）定义为对整个像平面积分：

$E_{\text{PSF}} = \iint_{-\infty}^{+\infty} I(u,v) \, du \, dv = \iint_{-\infty}^{+\infty} |U(u,v)|^2 \, du \, dv$

1.2.3. 帕塞瓦尔定理（Parseval’s Theorem）

傅里叶变换满足著名的帕塞瓦尔能量守恒定理：

对于任意平方可积函数 $f(x,y)$ 及其傅里叶变换 $F(u,v) = \mathcal{F}\{f(x,y)\}$，有：

$\iint |F(u,v)|^2 \, du \, dv = \iint |f(x,y)|^2 \, dx \, dy$

将 $P(x,y)$ 作为 $f(x,y)$，$U(u,v)$ 作为 $F(u,v)$，直接应用帕塞瓦尔定理：

$E_{\text{PSF}} = \iint |U(u,v)|^2 \, du \, dv = \iint |P(x,y)|^2 \, dx \, dy$

1.2.4. 计算 $|P(x,y)|^2$

由于 $P(x,y) = M(x,y) \cdot e^{i k W(x,y)}$，我们有：

$|P(x,y)|^2 = P(x,y) \cdot P^*(x,y) = M(x,y) \cdot e^{i k W(x,y)} \cdot M(x,y) \cdot e^{-i k W(x,y)}$

由于 $M(x,y)$ 是实函数（$0$ 或 $1$），$M^* = M$：

$|P(x,y)|^2 = M^2(x,y) \cdot \underbrace{e^{i k W} \cdot e^{-i k W}}_{=1} = M^2(x,y) = M(x,y)$

关键观察：相位项 $e^{i k W(x,y)}$ 的模恒为 $1$，因此在 $|P|^2$ 中完全抵消。

1.2.5. 连续域最终结论

代入得：

$\boxed{E_{\text{PSF}} = \iint M(x,y) \, dx \, dy = A_{\text{pupil}}}$

其中 $A_{\text{pupil}}$ 是光瞳孔径的几何面积。

对于圆形光瞳：

$E_{\text{PSF}} = \pi R^2$

其中 $R$ 为光瞳半径。

1.3. 离散域（数值仿真）严格推导

在实际仿真中，我们使用离散傅里叶变换（DFT/FFT）在有限网格上计算 PSF。

1.3.1. 离散化模型

设：

• 光瞳平面采样为 $N \times N$ 网格，采样间隔为 $\Delta x = \Delta y = d$
• 离散光瞳函数：$P_{m,n} = M_{m,n} \cdot e^{i k W_{m,n}}$，其中 $m,n = 0,1,\dots,N-1$
• 离散光瞳掩模：$M_{m,n} \in \{0, 1\}$

1.3.2. 二维 DFT 定义

像平面复振幅（DFT）：

$U_{k,l} = \sum_{m=0}^{N-1} \sum_{n=0}^{N-1} P_{m,n} \cdot e^{-i \frac{2\pi}{N}(mk + nl)}$

PSF 强度：

$I_{k,l} = |U_{k,l}|^2$

1.3.3. 离散帕塞瓦尔定理

二维 DFT 满足离散帕塞瓦尔定理：

$\sum_{k=0}^{N-1} \sum_{l=0}^{N-1} |U_{k,l}|^2 = N^2 \sum_{m=0}^{N-1} \sum_{n=0}^{N-1} |P_{m,n}|^2$

归一化说明：上式对应 numpy.fft.fft2 的默认定义（正变换无 $1/N$ 归一化）。若使用 fftshift、ifft2 等其他归一化约定，比例常数会变化，但比例关系不变。

1.3.4. 计算 $|P_{m,n}|^2$

同理：

$|P_{m,n}|^2 = |M_{m,n}|^2 \cdot |e^{i k W_{m,n}}|^2 = M_{m,n}^2 \cdot 1 = M_{m,n}$

因为 $M_{m,n} \in \{0,1\}$，所以 $M_{m,n}^2 = M_{m,n}$。

1.3.5. 离散域最终结论

设 $N_M$ 为光瞳掩模中值为 $1$ 的像素总数，即：

$N_M = \sum_{m=0}^{N-1} \sum_{n=0}^{N-1} M_{m,n}$

则：

$\sum_{k,l} |U_{k,l}|^2 = N^2 \sum_{m,n} M_{m,n} = N^2 \cdot N_M$

因此 PSF 总能量（像素强度之和）为：

$\boxed{E_{\text{PSF, discrete}} = \sum_{k=0}^{N-1} \sum_{l=0}^{N-1} I_{k,l} = N^2 \cdot N_M}$

在仿真中，如果你固定网格大小 $N$ 和光瞳掩模 $M$，则无论波前 $W_{m,n}$ 如何变化，PSF.sum() 恒等于 $N^2 \cdot N_M$。

若考虑实际的物理单位（像素面积 $\Delta u \cdot \Delta v$），总能量还应乘以像素面积，但这不影响"波前不影响总能量"的结论。

1.4. 物理直观理解

1.4.1. 相位是"搬运工"，不是"生产者"

• 光瞳掩模 $M$ 决定了有多少光能通过系统，它是能量的"源头"。
• 波前 $W$（相位）只是决定了这些光在像平面上如何干涉、如何重新分布。

可以把光瞳想象成一个房间里的多盏灯：

• $M$ 决定开了多少盏灯（总功率）。
• $W$ 决定每盏灯的开关时机和相位关系（干涉图样）。

无论这些灯如何干涉，总功率不变，只是有的地方更亮、有的地方更暗。

1.4.2. 像差的影响

• 理想平面波（$W=0$）：PSF 为艾里斑，能量集中在中心。
• 离焦/像差（$W \neq 0$）：中心能量降低，旁瓣升高，能量向外部扩散。
• 随机相位（强像差）：PSF 变成散斑图样，能量极度分散。

但在以上所有情况下，对全像平面求和，总能量完全相同。

1.5. 特殊情况讨论

1.5.1. 如果光瞳掩模不是严格的 0/1？

如果 $M(x,y)$ 包含灰度值（如部分透过率 $0 < M < 1$），结论仍然成立：

$E_{\text{PSF}} = \iint |M(x,y)|^2 \, dx \, dy$

只是此时 $|P|^2 = |M|^2$，不再是简单的 $M$。

1.5.2. 如果光瞳函数包含振幅衰减？

如果光瞳内存在非均匀的振幅衰减 $A(x,y)$（例如高斯照明），则：

$P(x,y) = A(x,y) \cdot M(x,y) \cdot e^{i k W(x,y)}$

此时：

$|P|^2 = |A|^2 \cdot M$

总能量取决于 $|A|^2$ 在光瞳内的积分，仍然与相位 $W$ 无关。

1.5.3. 数值仿真中能量不守恒的排查

如果你在实际仿真中发现改变波前后 PSF.sum() 变化了，请检查：

可能原因	排查方法
PSF 被截断	增大计算网格 $N$，确保旁瓣没有溢出边界
FFT 频谱泄漏	确保光瞳边缘在网格内清晰采样，避免混叠
归一化不一致	不同参数下使用了不同的常数因子乘在 PSF 上
数值下溢/上溢	检查是否有 float32 溢出，建议用 float64
光瞳定义变化	确认 $M$ 在不同条件下确实完全相同

1.6. 总结

问题	答案
PSF 总能量由什么决定？	仅由光瞳掩模 $M$（的平方积分/求和）决定
波前 $W$ 影响总能量吗？	不影响。无论 $W$ 如何，总能量严格守恒
连续域表达式	$E = \displaystyle\iint M(x,y) \, dx \, dy = A_{\text{pupil}}$
离散域表达式	$E = N^2 \displaystyle\sum_{m,n} M_{m,n}$（无归一化 FFT）
物理本质	相位只重新分配能量空间分布，不改变总能量

1.7. 实测数据验证：曝光时间固定下，不同像差对 PSF 能量的影响

1.7.1. 实验设计

为验证理论推导在实际采集数据中的适用性，我们从 dataset_split_txts/test.txt 中随机抽取了 2000 个样本（先 Shuffle 再随机挑选，确保像差分布均匀）。

• 数据来源：光学系统实测采集的 NPZ 文件（image + zernike）
• 曝光条件：采集期间曝光时间严格固定（WFS=500, PD1/PD2/PD3=50）
• 能量计算：直接读取 npz 原始数据，不经过 dataset.py 中的背景扣除和归一化预处理，分别计算：
  • 在焦能量：$E_{\text{in}} = \sum \text{image}[:,:,0]$
  • 离焦能量：$E_{\text{def}} = \sum \text{image}[:,:,1]$
  • 总能量：$E_{\text{tot}} = E_{\text{in}} + E_{\text{def}}$
• 像差度量：Zernike 多项式系数的 RMS 值

1.7.2. 像差覆盖范围

注：我们对 test.txt 中的 全部 39960 个样本 进行了扫描确认。该数据集的 63 维 Zernike 已经去除了 piston 和 tilt，全部为有效像差模式。

波前 RMS 计算：对于正交归一的 Zernike 多项式，波前 RMS = $\sqrt{\sum_i Z_i^2}$（平方求和再开方）。

统计量	波前 RMS sqrt(∑Z²)	波前峰谷值 P-V	最大单项系数 max|Z_i|
最小值	0.0211	0.0136	0.0073
最大值	0.9255	0.8202	0.6239
均值	0.4295	0.3601	0.2776
P90	0.8898	—	—
P95	0.9093	—	—
P99	0.9199	—	—

1.7.3. 能量统计结果

指标	在焦能量	离焦能量	总能量
均值	4793.56	4813.41	9606.97
标准差	4.80	2.19	5.62
变异系数 (CV)	0.10%	0.05%	0.06%

变异系数（CV = std/mean）是衡量相对波动性的关键指标。 总能量的 CV 仅为 0.06%，说明在 2000 个不同像差样本中，总能量几乎恒定不变。

1.7.4. 能量与像差的相关性

能量类型	与波前 RMS sqrt(∑Z²) 的 Pearson 相关系数
在焦能量	+0.2184（弱正相关）
离焦能量	-0.8639（强负相关）
总能量	-0.1501（弱负相关，接近无关）

解读：

• 离焦能量随像差增大而显著下降（能量从离焦面扩散/转移）。
• 但在焦能量有微弱上升，两者相互补偿。
• 总能量与像差几乎无关（$r = -0.15$），验证了相位只重新分配能量、不改变总量的理论结论。

1.7.5. 分组对比（按像差四分位数）

像差等级	样本数	总能量均值	总能量标准差	CV
低像差（$\leq$ P25）	500	9604.21	7.97	0.08%
中像差（P25 ~ P75）	1000	9610.13	2.94	0.03%
高像差（$\geq$ P75）	500	9603.40	2.29	0.02%
高 vs 低相对差异	—	-0.008%	—	—

在 test.txt 实际像差范围内（波前 RMS 从 0.02λ 到 0.93λ，动态范围约 43 倍），总能量变化不到 0.01%。

1.7.6. 可视化

1.7.6.1. 散点图：能量 vs 波前 RMS

• 左图：在焦能量随像差略有上升（$r=0.22$）。
• 中图：离焦能量随像差显著下降（$r=-0.86$）。
• 右图：总能量基本水平，无系统性趋势（$r=-0.15$）。

1.7.6.2. 总能量 vs 波前 P-V 及最大单项系数

• 左图：总能量 vs 波前峰谷值 P-V（$r = -0.16$）。
• 右图：总能量 vs 最大单项系数（$r = -0.16$）。

两种辅助度量下，总能量均无系统性趋势。

1.7.6.3. 总能量分布直方图

总能量呈极窄的分布，均值 9606.97，半高宽极小。

1.7.6.4. 分组箱线图

低、中、高像差三组的能量中位数几乎重合，箱线高度差异极小。

1.7.7. 实测结论

理论预测	实测结果	一致性
光瞳不变时，PSF 总能量与波前（像差）无关	2000 个样本，总能量 CV = 0.06%，与波前 RMS 相关系数仅 -0.15	✅ 完全符合
相位只改变能量空间分布	在焦能量↑，离焦能量↓，两者补偿	✅ 完全符合
曝光时间固定保证入射光量恒定	总能量在不同像差下差异 < 0.01%	✅ 得到验证

1.7.8. 进一步讨论：实测中观察到的能量"补偿"现象

理论上说，总能量应该与像差严格无关（相关系数为 0）。但实测中我们发现：

• 在焦能量随像差增大有微弱上升（$r = +0.22$）
• 离焦能量随像差增大显著下降（$r = -0.86$）
• 两者相互补偿，使总能量保持恒定

这种"补偿"并非理论上的能量转移，而更可能是相机探测非线性导致的测量效应：

1.7.8.1. 像差小 → 中心过曝（能量截断）

当像差很小时，PSF 能量高度集中在焦平面中心。如果中心像素的光强超过了相机的满阱容量（full well capacity），就会发生饱和/过曝。过曝区域的像素值被截断在最大值（如 4095 for 12-bit），导致：

• 在焦图像测得的能量偏低（真实的中心能量被截断了）
• 随着像差增大，能量扩散，中心不再过曝，测得的在焦能量反而"回升"

这与我们观察到的"在焦能量随像差增大而上升"的趋势一致。

1.7.8.2. 像差大 → 信号淹没于噪声

当像差很大时，PSF 能量被严重稀释到更大的空间范围内，单位面积的光强显著下降。此时：

• 信号可能接近甚至低于相机的读出噪声 + 暗电流噪声基线
• 即使做了背景扣除，弱信号区域的信噪比（SNR）仍然很低
• 离焦图像的能量更容易被噪声"吞没"，导致测得的离焦能量偏低

这与我们观察到的"离焦能量随像差增大而显著下降"（$r = -0.86$）高度吻合。

1.7.8.3. 总结

像差大小	在焦通道	离焦通道	对总能量的影响
小像差	中心过曝 → 测得能量偏低	能量集中 → 测得正常	两者部分抵消
大像差	能量扩散 → 过曝消失，测得能量回升	能量极度扩散 → 被噪声淹没，测得能量偏低	两者再次抵消

这解释了为什么总能量能保持如此惊人的稳定（CV = 0.06%）：不是物理上能量在在焦/离焦之间完美转移，而是两种非线性效应（过曝截断 vs 噪声淹没）恰好相互补偿了。

如果要更严格地验证理论，需要确保：

1. 无过曝：像差小时使用更短曝光或衰减片，避免中心饱和
2. 无噪声淹没：像差大时延长曝光或提高增益，确保弱信号高于噪声基线
3. 或者使用具有更高动态范围（HDR）的探测器

光瞳是"量"，波前是"形"。光瞳决定有多少光，波前决定光长什么样。只要光瞳不变，点扩散的总能量就永恒不变。

2. SHFP（子孔径高频相位）层：频率分区与数学构造

2.1. 核心问题：哈特曼传感器能采到什么频率？

2.1.1. 采样定理的一维推导

哈特曼波前传感器（Hartmann WFS）是 10×10 的子孔径网格。我们先从一维理解：

• 在光瞳直径上，有 N_H = 10 个采样点（子孔径中心）
• 根据奈奎斯特采样定理，能无混叠地重建的最高空间频率为：

$f_{Nyquist}^{Hartmann} = \frac{N_H}{2} = 5 \text{ cycles/pupil}$

这意味着：

• 频率 ≤ 5 c/pupil 的波前变化 → 哈特曼可以分辨
• 频率 > 5 c/pupil 的波前变化 → 哈特曼无法分辨（混叠或丢失）

2.1.2. 一维图示

光瞳（-1 ~ +1，归一化坐标）
|----|----|----|----|----|----|----|----|----|----|
  s1   s2   s3   s4   s5   s6   s7   s8   s9  s10
  ← 10个子孔径中心采样点 →

波前分量：
  3 c/pupil (低频)   : ~~~~  (哈特曼可以采)
  8 c/pupil (中频)   : ~~~~~~~~  (哈特曼混叠)
  15 c/pupil (高频)  : ~~~~~~~~~~~~~~~~  (完全丢失)

2.1.3. 51×51 光瞳网格上的频率映射

在 51×51 的像素网格上做 FFT 时，频率与 FFT 索引的关系：

Centered 格式（零频在中间，可视化常用）：

• 索引 k（相对于中心零频 k=25）的实际频率：

$f(k) = k - 25 \quad \text{(cycles/pupil)}$

• k = 20（距中心 5）→ 频率 = 5 c/pupil（哈特曼奈奎斯特极限）
• k = 15（距中心 10）→ 频率 = 10 c/pupil（SHFP中频目标）
• k = 10（距中心 15）→ 频率 = 15 c/pupil（SHFP高频上限）
• k = 0 或 k = 50 → 频率 = ±25 c/pupil（51网格奈奎斯特极限）

重要：必须理解 ifftshift 的作用

torch.fft.ifft2 要求输入频谱的零频在索引 (0,0)，这是标准 FFT 格式：

标准格式：    [0, 1, 2, ..., 25, -25, ..., -2, -1]  (实际频率)
              ↑零频
              
Centered格式：[ -25, ..., -2, -1, 0, 1, 2, ..., 25]
                        ↑零频在中间

torch.fft.ifftshift 的作用就是将 “零频居中” 的可视化格式转换为 “零频在 (0,0)” 的标准格式。

对于 N=51，中心索引 k=25（距中心 0 像素）对应零频。索引 k=45 距中心 20 像素，但在标准格式中：

• ifftshift 后，k=45 的位置会映射到标准格式的 (45 - 25 = 20)，即 -6 c/pupil（因为 20 > 25.5，实际频率 = 20 - 51 = -6）
• 不加 ifftshift，ifft2 会把 k=45 当成 +45 c/pupil（远超奈奎斯特极限），输出像素噪声

2.1.4. 频率分区总表

区域	频率范围 (c/pupil)	Centered格式索引	物理意义	SHFP作用
低频	0 ~ 5	k = 20~30	哈特曼+Zernike已覆盖	不学
中频	5 ~ 15	k = 10~20, 30~40	哈特曼丢失，Zernike弱覆盖	主要学习目标
高频	15 ~ 25	k = 0~10, 40~50	接近网格奈奎斯特，信噪比差	谨慎学习

2.2. SHFP 层的目标：学习中频（5~15 c/pupil）

2.2.1. 为什么选中频而非高频？

1. 过曝的物理本质：

   • 大像差样本的灰圈/光晕是子孔径尺度的动态效应
   • 子孔径大小约 5 像素（51/10 ≈ 5.1），子孔径内的 1~2 个周期的相位变化对应全光瞳上的 5~10 c/pupil
   • 这正是哈特曼完全丢失、Zernike 基函数难以精细表达的区域

2. 为什么不学 >15 c/pupil？

   • 51 网格的奈奎斯特是 25 c/pupil，但 15~25 区域已经接近像素极限
   • 高频分量能量通常很小，且容易过拟合噪声
   • 实验观察到的光晕展宽主要由中频分量（5~15）贡献

2.2.2. 为什么避开 <5 c/pupil（低频）？

• Zernike 63 项（到第 10~11 径向阶）已经能很好地表示 0~5 c/pupil 的波前
• SHFP 如果学习低频，会与 Zernike 发生参数冲突（competing gradients）
• 之前的方形实现 n_freq=12 覆盖了 0~6 c/pupil，导致与 Zernike 重叠
• 改进后的环形实现只覆盖 5~15 c/pupil，与 Zernike 互补

2.3. SHFP 层的数学构造步骤

2.3.1. 步骤 1：定义目标频谱区域（环形中频带）

Centered 频谱（51×51，中心区域示意，距离按欧氏距离 √u²+v² 计算）：
    v
    ↑
    | ...RRRRRRRRRRRRRRRRRRR...
    | ..RRRRRRRRRRRRRRRRRRRRR..
    | .RRRRRRRRRRRRRRRRRRRRRRR.
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRBBBBBRRRRRRRRRR   <- dist=5（哈特曼奈奎斯特边界）
    | RRRRRRRRRBBBBBBBRRRRRRRRR
    | RRRRRRRRBBBBBBBBBRRRRRRRR
    | RRRRRRRRBBBBBBBBBRRRRRRRR
  0-| RRRRRRRRBBBBBBBBBRRRRRRRR ←—— 零频 (u=0, v=0)
    | RRRRRRRRBBBBBBBBBRRRRRRRR
    | RRRRRRRRBBBBBBBBBRRRRRRRR
    | RRRRRRRRRBBBBBBBRRRRRRRRR
    | RRRRRRRRRRBBBBBRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | RRRRRRRRRRRRRRRRRRRRRRRRR
    | .RRRRRRRRRRRRRRRRRRRRRRR.
    | ..RRRRRRRRRRRRRRRRRRRRR..
    | ...RRRRRRRRRRRRRRRRRRR...   <- dist=15（SHFP高频边界）
    ←———————————————— u ———————————————→

    图例: B = 低频盲区（圆形，dist<5，Zernike已覆盖，SHFP不学）
          R = SHFP可学习区域（环形带，5≤dist≤15 c/pupil）
          . = 高频区（dist>15，接近网格奈奎斯特，不学习）

2.3.2. 步骤 2：构造环形频率掩膜

在 PyTorch 中：

N = 51
center = N // 2  # 25
y, x = torch.meshgrid(torch.arange(N), torch.arange(N), indexing='ij')
dist = torch.sqrt((x - center)**2 + (y - center)**2)

# 低频盲区：r < 5（Zernike已覆盖）
low_freq_mask = (dist >= 5).float()

# 高频截断：r > 15（避免过拟合噪声）  
high_freq_mask = (dist <= 15).float()

# SHFP有效区域：环形带 5 <= r <= 15
shfp_mask = low_freq_mask * high_freq_mask  # [51, 51]

2.3.3. 步骤 3：可学习频谱参数

在有效环形区域内放置可学习的复数傅里叶系数：

class HighFrequencyPhaseLayer(nn.Module):
    def __init__(self, grid_size=51, min_freq=5, max_freq=15):
        super().__init__()
        self.grid_size = grid_size
        self.min_freq = min_freq
        self.max_freq = max_freq
        
        # 在环形区域内创建可学习参数
        y, x = torch.meshgrid(
            torch.arange(grid_size), torch.arange(grid_size), indexing='ij'
        )
        dist = torch.sqrt((x - grid_size//2)**2 + (y - grid_size//2)**2)
        
        # 注册环形掩膜（不可学习）
        self.register_buffer('ring_mask', 
            ((dist >= min_freq) & (dist <= max_freq)).float())
        
        # 可学习复数系数（只在mask=1的位置有效）
        self.amp_real = nn.Parameter(
            torch.zeros(grid_size, grid_size) * 0.01)
        self.amp_imag = nn.Parameter(
            torch.zeros(grid_size, grid_size) * 0.01)
        
        # 自适应幅度缩放
        self.intensity_scale = nn.Parameter(torch.tensor(0.05))
    
    def forward(self, batch_size, device, complexity):
        # 1. 构造完整频谱（centered格式，零频在中间）
        spectrum = torch.complex(self.amp_real, self.amp_imag)
        spectrum = spectrum * self.ring_mask  # 只保留环形带
        
        # 2. ifftshift：将零频从中心移到 (0,0)
        full_spectrum = torch.fft.ifftshift(spectrum)
        
        # 3. ifft2：从频域合成空间相位屏
        phase_hf = torch.fft.ifft2(full_spectrum).real
        
        # 4. 归一化（保持单位标准差）
        phase_hf = phase_hf / (phase_hf.std() + 1e-6)
        
        # 5. 与波前复杂度自适应耦合
        scale = torch.sigmoid(self.intensity_scale) * complexity.view(-1, 1, 1) * 0.1
        
        return phase_hf.unsqueeze(0).expand(batch_size, -1, -1) * scale

2.3.4. 步骤 4：叠加到 Zernike 波前

# 在 Zernike2PSF_layer.forward() 中
wavefront_rad = zernike_wavefront / wavelength_m * 2 * torch.pi  # [B, 51, 51]

# 计算波前复杂度（用于SHFP自适应幅度）
complexity = torch.norm(zernike_um, dim=1, keepdim=True)  # [B, 1]

# 生成高频相位扰动
hf_phase = self.hf_phase_layer(batch_size, device, complexity)  # [B, 51, 51]

# 叠加：Zernike低频 + SHFP中频
wavefront_rad = wavefront_rad + hf_phase

2.3.5. 步骤 5：物理衍射（fft2）

pupil_complex = effective_amp * torch.exp(1j * wavefront_rad)

# 注意：这里 fft2 是物理衍射，与 SHFP 的 ifft2（傅里叶级数合成）不同
# fft2 输入也要求零频在 (0,0)， pupil_complex 是空间域，自然满足
psf_amplitude = torch.fft.fft2(pupil_complex)
psf_intensity = torch.fft.fftshift(torch.abs(psf_amplitude)**2)

2.4. SHFP 与 zernike_residual 的互补性

模型中已有 zernike_residual（63-D 可学习向量），它与 SHFP 层是否重叠？不重叠，是互补关系。

2.4.1. zernike_residual 的本质

# 63-D 可学习向量，所有样本共享同一残差（全局校准误差）
self.zernike_residual = nn.Parameter(torch.zeros(63, dtype=torch.float32))

使用方式：zernike_um = zernike_um + self.zernike_residual

• 物理意义：校正输入 Zernike 系数的全局系统误差（如哈特曼标定偏差）
• 表示空间：Zernike 域（63 个系数）→ 通过 zernikePol 映射到 51×51 空间网格
• 与输入关系：固定偏移，所有样本加同一个 63 维向量

2.4.2. 一个常见的误解

“63 项 Zernike 在 51×51 网格上的频率上限只有 0~5 c/pupil，所以 zernike_residual 只覆盖低频。”

这个说法不对。 Zernike 多项式是连续正交多项式，高阶模式（如 n=10, m=10）在 51 点网格上采样后，空间振荡完全可以达到 10~20 c/pupil。这些高频不是采样定理意义上的"可分辨频率"，而是混叠后在离散网格上呈现的高频内容。

所以 zernike_residual 在空间域上的频率范围不限于 0~5 c/pupil。

2.4.3. 那为什么 zernike_residual 仍然替代不了 SHFP？

核心区别在于自由度约束和输入相关性：

	zernike_residual	SHFP 层
自由度	63 个系数	640 个频谱像素
基函数形状	Zernike 多项式（全局、有固定解析形式）	傅里叶正弦/余弦（任意形状）
空间局部性	❌ 全局分布	✅ 可局部可全局
输入相关性	❌ 固定（所有样本相同）	✅ 自适应（大像差→大扰动）
物理约束	受 Zernike 正交性约束	自由（仅频率带限）

关键问题：Zernike 是"结构化"的基。

63 项 Zernike 只能生成 63 维子空间 中的波前。这个子空间包含径向多项式 × 角向三角函数的特定组合，不能生成任意的子孔径尺度局部扰动。

具体例子：假设过曝只在左上角第 (2,3) 个子孔径产生了局部相位畸变。Zernike 基函数无法单独"激活"一个子孔径——所有 Zernike 模式都是覆盖整个光瞳的。要用 63 个全局基去拟合一个局部特征，需要大量基函数的组合，且会引入不必要的全局副作用。

而 SHFP 是像素级自由的。它可以直接在 (2,3) 子孔径位置生成局部扰动，不影响其他区域。

2.4.4. zernike_residual 是否已经吸收了部分过曝效应？

有可能。 训练过程中，zernike_residual 会学习到"让平均 PSF 更接近真实"的方向。如果过曝导致的灰圈有一个固定的高频分量（比如所有大像差样本都有类似的子孔径边缘模糊），zernike_residual 的高阶模式可能会部分吸收它。

但从实验现象看：

多跑 epoch 后模型自发全局加宽（blur_kernel/spatial_scale），大像差改善但小像差过拟合

这说明：

1. zernike_residual 不够用——如果它足够，模型不需要靠全局 blur 来补偿
2. 过曝效应主要是动态的（大像差时有、小像差时无），不是固定的
3. 全局参数无法区分"大像差需要宽 PSF、小像差不需要"

这正是引入 SHFP 的根本原因：提供一个输入相关的、局部自由的自由度。

2.4.5. 总结

问题	答案
zernike_residual 有高频内容吗？	✅ 有（高阶 Zernike 在 51 网格上可产生 >5 c/pupil）
它的自由度足够吗？	❌ 不够（只有 63 维，且受 Zernike 形状约束）
它能表示局部子孔径扰动吗？	❌ 不能（Zernike 是全局基）
它是输入相关的吗？	❌ 不是（固定偏移，所有样本相同）
它已经吸收了过曝效应吗？	可能部分吸收了固定的低频/全局分量
SHFP 还有必要吗？	✅ 有必要（处理动态的、局部的、Zernike 表示不了的分量）

3. 哈特曼波前测量中的子孔径间振幅不均匀性及其对聚焦点扩散函数的影响

哈特曼-夏克（Hartmann-Shack）波前传感器通过子孔径焦斑阵列重构入射波前。传统上，子孔径光斑的质心位移被用来解算相位斜率，从而重构波前相位 $\phi(\mathbf{r})$；然而，不同子孔径光斑之间的亮度差异直接反映了瞳孔面上波前振幅 $A(\mathbf{r})$ 的宏观空间分布。本文建立从哈特曼子孔径间光斑不均匀性到波前振幅重构、再到后续光学系统聚焦点扩散函数（PSF）的完整链路，推导了振幅调制对聚焦 PSF 的定量影响，并通过数值模拟验证了不同子孔径振幅分布下的 PSF 畸变规律。

3.1. 问题背景与物理链路

3.1.1. 哈特曼波前传感器的复振幅重构

哈特曼传感器将入射波前分割为 $N \times M$ 个子孔径。对于第 $(m,n)$ 个子孔径，探测器上记录的光强分布为：

$h_{m,n}(\boldsymbol{\rho}) = \left| \mathcal{F}\{ P_{m,n}(\mathbf{r}) \, A(\mathbf{r}) \, e^{i\phi(\mathbf{r})} \}(\boldsymbol{\rho}) \right|^2$

其中 $\mathbf{r}=(x,y)$ 为瞳孔面（波前平面）的横向笛卡尔坐标，$\boldsymbol{\rho}=(u,v)$ 为焦平面（透镜后焦面）的横向笛卡尔坐标；$P_{m,n}(\mathbf{r})$ 为子孔径掩模，$A(\mathbf{r})$ 和 $\phi(\mathbf{r})$ 分别为瞳孔面的振幅与相位。（注：本文中粗体字母 $\mathbf{r}, \boldsymbol{\rho}$ 均表示二维笛卡尔坐标矢量，而非极坐标；标量 $r=|\mathbf{r}|$ 仅在表示径向距离时出现。）

传统波前重构算法仅提取 $h_{m,n}$ 的质心位移，反演局部相位斜率，从而得到相位图 $\phi(\mathbf{r})$。但 $h_{m,n}$ 的总亮度（光斑积分强度）还包含振幅信息：

• 子孔径间光斑亮度差异 → 瞳孔面不同区域的平均振幅 $A(\mathbf{r})$ 不同。

若将哈特曼测得的所有子孔径振幅值拼接，即可得到 pupil 面上的宏观振幅分布图 $A(\mathbf{r})$。本文假设每个子孔径内部振幅均匀，仅考虑子孔径之间的振幅差异。重构出的完整波前应表示为复振幅场：

$\boxed{ U(\mathbf{r}) = A(\mathbf{r}) \, e^{i\phi(\mathbf{r})} } \tag{1}$

其中 $A(\mathbf{r})$ 是一个阶梯函数（piecewise constant），在每个子孔径内为常数。

3.1.2. 从重构波前到聚焦点扩散函数

在自适应光学链路中，哈特曼测得的波前会被送入后续光学系统（如聚焦透镜）。根据傅里叶光学，透镜后焦面的光场正比于瞳孔面复光场的傅里叶变换：

$U_f(\boldsymbol{\rho}) = \mathcal{F}\{ U(\mathbf{r}) \}(\boldsymbol{\rho}) = \mathcal{F}\{ A(\mathbf{r}) e^{i\phi(\mathbf{r})} \}(\boldsymbol{\rho})$

聚焦**点扩散函数（PSF）**定义为焦面强度分布：

$\boxed{ h(\boldsymbol{\rho}) = \left| \mathcal{F}\{ A(\mathbf{r}) e^{i\phi(\mathbf{r})} \}(\boldsymbol{\rho}) \right|^2 } \tag{2}$

式 (2) 建立了完整的物理链路：哈特曼子孔径间光斑亮度不均匀 → 表征为波前振幅 $A(\mathbf{r})$ 的宏观不均匀 → 与相位 $e^{i\phi}$ 共同决定聚焦 PSF。

3.2. 理论推导：子孔径间振幅不均匀如何影响聚焦 PSF

3.2.1. 平面波相位下的纯振幅效应

为剥离相位因素的干扰，考虑一种常见场景：相位已被哈特曼测得并完全校正（或入射光本身为平面波，$\phi(\mathbf{r})=0$）。此时式 (2) 简化为：

$\boxed{ h(\boldsymbol{\rho}) = \left| \mathcal{F}\{ A(\mathbf{r}) \}(\boldsymbol{\rho}) \right|^2 } \tag{3}$

式 (3) 是本文的核心结果：当相位恒定时，聚焦 PSF 完全由瞳孔振幅分布 $A(\mathbf{r})$ 的傅里叶变换模平方决定。 这意味着，即使哈特曼测得"理想平面波"相位，只要子孔径间存在亮度不均匀（即 $A(\mathbf{r}) \neq \text{const}$），聚焦光斑就必然偏离理想衍射极限。

3.2.2. 子孔径间振幅调制的数学描述

将 pupil 面划分为 $M \times M$ 个子孔径，每个子孔径内的振幅为常数 $A_{m,n}$。则 pupil 面振幅可写为：

$A(\mathbf{r}) = \sum_{m,n} A_{m,n} \, P_{m,n}(\mathbf{r})$

其中 $P_{m,n}(\mathbf{r})$ 为第 $(m,n)$ 个子孔径的矩形窗函数。代入式 (3)：

$h(\boldsymbol{\rho}) = \left| \sum_{m,n} A_{m,n} \, \mathcal{F}\{P_{m,n}\}(\boldsymbol{\rho}) \right|^2$

对于尺寸为 $d \times d$ 的矩形子孔径，$\mathcal{F}\{P_{m,n}\} \propto \text{sinc}(u d) \, \text{sinc}(v d)$。因此 PSF 是多个 sinc 函数的相干叠加，加权系数为各子孔径的振幅 $A_{m,n}$。

物理意义：

• 若所有 $A_{m,n}=A_0$（均匀），各 sinc 函数相干叠加后中心主瓣相互增强，旁瓣相消，得到理想窄 PSF；
• 若 $A_{m,n}$ 空间变化，不同位置的 sinc 函数加权不同，破坏相干相消条件，导致旁瓣能量泄漏；
• 若 $A_{m,n}$ 呈周期性变化（如棋盘格），旁瓣会出现明显的周期性结构。

3.2.3. 能量守恒关系

由帕塞瓦尔定理：

$\int |A(\mathbf{r})|^2 d\mathbf{r} = \int h(\boldsymbol{\rho}) d\boldsymbol{\rho} = E_{\text{total}}$

总能量守恒。但振幅不均匀会在空间上重新分配能量：将中心主瓣的能量转移到旁瓣。

定义旁瓣能量占比：

$\eta = \frac{\int_{\text{sidelobe}} h(\boldsymbol{\rho}) d\boldsymbol{\rho}}{E_{\text{total}}}$

对于均匀振幅，理想情况下 $\eta$ 仅受衍射极限限制。对于非均匀振幅，$\eta$ 随子孔径间振幅起伏的剧烈程度单调增加。

3.2.4. 子孔径间振幅不均匀的典型模式与 PSF 对应关系

子孔径振幅分布 $A_{m,n}$	哈特曼观测现象	对聚焦 PSF 的影响
全部相等（Uniform）	所有子孔径光斑亮度一致	理想衍射极限 PSF，旁瓣最低
随机起伏（Random）	子孔径光斑亮度无规则差异	不规则旁瓣，峰值强度下降
中心高、边缘低（Gaussian）	中心子孔径亮，边缘子孔径暗	主瓣略宽，旁瓣被压制
棋盘格交替（Checkerboard）	相邻子孔径亮度交替变化	强烈的周期性旁瓣，能量泄漏严重
部分子孔径丢失（Dropout）	部分子孔径完全无光	PSF 出现十字/星形旁瓣，类似稀疏孔径

3.3. 数值模拟与可视化

采用离散 FFT 模拟式 (3) 的物理过程。设 pupil 面采样为 $N \times N$，划分为 $8 \times 8$ 个子孔径，并统一限制在圆形孔径 $r \le 0.9$ 内（孔径外强制为 0）。相位恒为 0，每个子孔径内振幅均匀，仅改变子孔径间的振幅值 $A_{m,n}$，计算聚焦 PSF $h = |\text{FFT2}\{A\}|^2$。

3.3.1. 模拟参数与子孔径振幅模型

对比五种典型子孔径振幅分布：

1. Uniform：所有子孔径 $A_{m,n}=1$
2. Random subap：每个子孔径随机振幅 $0.2 \sim 1.0$
3. Gaussian subap：子孔径振幅按高斯分布（中心亮、边缘暗）
4. Checkerboard：相邻子孔径振幅 $0.2$ 和 $1.0$ 交替
5. Partial dropout：随机 20% 子孔径振幅置零（模拟遮挡/坏点）

3.3.2. 结果：子孔径振幅分布与 PSF 的完整链路

下图从左至右依次展示了子孔径振幅分布 $A(\mathbf{r})$（带网格线）、PSF 线性放大图、PSF 对数全局图，完整对应式 (3) 的推导链条。

观察结果：

• Uniform：所有子孔径等亮，焦面振幅为理想 sinc 包络，PSF 主瓣最尖锐，旁瓣最低。
• Random subap：子孔径亮度随机起伏破坏了相干相消，对数图中可见明显的无规则旁瓣。
• Gaussian subap：边缘子孔径被抑制，等效于对有限孔径做软边切趾（apodization），旁瓣能量被压到约 3%，为所有情况中最低。
• Checkerboard：相邻子孔径振幅反相交替，引入强高频分量，PSF 出现显著的周期性旁瓣，能量大量泄漏。
• Partial dropout：部分子孔径缺失，等效于稀疏孔径，PSF 出现十字形旁瓣结构，与干涉合成孔径的点扩散函数形态一致。

4. 离焦系数与物理离焦距离的等价性证明及数值验证

推导为什么"在仿真中加 Zernike 离焦项"与"真实光路中移动相机"是完全等价的。

4.1. 光学结构

先画一个简单的示意图：

光瞳平面 (Pupil Plane)          焦平面 (Focal Plane)      离焦平面 (Defocus Plane)
    ↓                                  ↓                          ↓
   [○]  ←——— 距离 f ———→  [·]  ←——— 距离 Δz ———→  [○]
   圆形孔径                  理想聚焦点                离焦光斑

• 光瞳平面：Mask.mat 描述的圆孔，直径 $D = 60\,\text{mm}$
• 焦平面：距离光瞳 $f = 1248.158\,\text{mm}$，这里是光束完美会聚的位置
• 离焦平面：你把相机往前或往后移动了 $\Delta z$，光束在这里还没完全会聚（或已经发散），所以光斑比焦平面大

光瞳平面上的场 $U_p$ 经过透镜后，传播到像面。在仿真中，我们用 FFT（快速傅里叶变换）来计算这个传播过程。

4.2. 方法 A：Zernike 离焦项是怎么工作的

4.2.1. 加离焦项后的相位

在仿真中，你的代码做了这样一件事：在原来的 Zernike 系数基础上，给离焦项额外加一个系数 $C_d$（单位是 $\mu\text{m}$）。

于是光瞳上的波前（以长度为单位）变为：

$W_A(\rho) = C_d \cdot Z_{\text{defocus}}(\rho) = C_d \, (c_2 \rho^2 + c_0)$

对应的相位为：

$\phi_A(r) = \frac{2\pi}{\lambda} \, W_A(\rho)$

把 $\rho = 2r/D$ 代入：

$\begin{aligned} \phi_A(r) &= \frac{2\pi}{\lambda} \, C_d \left(c_2 \cdot \frac{4r^2}{D^2} + c_0\right) \\[8pt] &= \underbrace{\frac{8\pi c_2}{\lambda D^2} \, C_d \, r^2}_{\text{随 } r \text{ 变化的部分}} + \underbrace{\frac{2\pi c_0}{\lambda} \, C_d}_{\text{常数部分}} \end{aligned}$

4.2.2. 3.4 为什么常数相位不影响 PSF？

上式中第二项 $\frac{2\pi c_0}{\lambda} C_d$ 是一个与 $r$ 无关的常数。这意味着整个光瞳上的每一点都额外转动了相同的角度。

由于 PSF 的计算是取模平方：

$\text{PSF} = \bigl|\mathcal{F}\{A \cdot e^{i(\phi + \phi_0)}\}\bigr|^2 = \bigl|e^{i\phi_0} \cdot \mathcal{F}\{A \cdot e^{i\phi}\}\bigr|^2 = \bigl|\mathcal{F}\{A \cdot e^{i\phi}\}\bigr|^2$

常数相位 $e^{i\phi_0}$ 的模是 1，所以会被消掉。我们在比较两种方法时，只需要关心 $r^2$ 的系数是否一致。

4.3. 方法 B：移动相机是怎么改变光场的

4.3.1. 直觉理解

想象你用放大镜聚焦阳光到一张纸上的一个小点。如果你把纸往上移一点（远离放大镜），光点就会变大变模糊。这就是离焦。

在数学上，移动相机 $\Delta z$ 意味着：光波从光瞳出发后，需要传播距离 $f + \Delta z$ 才能到达像面，而不是原来的 $f$。

4.3.2. 严格推导：从惠更斯原理到菲涅耳衍射

4.3.2.1. 第一步：惠更斯原理——每一点都是新的波源

核心问题：光从光瞳传播到像面，中间发生了什么？

荷兰物理学家惠更斯（Christiaan Huygens）在 1678 年提出了一个惊人的想法：

波前上的每一点，都可以看作是一个新的球面波的波源。这些次级球面波向四面八方传播，它们互相叠加，就形成了下一时刻的波前。

光从光瞳传播到像面的过程，本质上就是：光瞳上的每一点发出球面波，这些球面波在像面上叠加干涉。

你可能会问：如果每一点都向四面八方发球面波，平面波不就散开了吗？

这是一个非常经典的问题！答案是：球面波确实向四面八方扩散了，但它们在绝大多数方向上会互相抵消，只在特定方向上相干增强。

用排队扔石头来理解：

想象 100 个人站成一条直线（模拟平面波的波前），每个人同时扔一块石头到水里。每个人产生的波纹都是圆形，向四面八方扩散。

现在问你：在这些人正前方远处的水面上，波纹是什么样的？

答案是：还是一条直线形波纹（平面波）！

为什么？因为在正前方远处的某一点：

• 从左边第 1 个人传来的波纹走了最远的路，相位滞后最多
• 从中间第 50 个人传来的波纹走了最短的路，相位滞后最少
• 从右边第 100 个人传来的波纹走了最远的路，相位滞后最多

但这些波纹到达时，波峰恰好对齐——左边的波峰虽然出发早，但走得远；右边的波峰虽然出发晚，但走得近。最终所有波峰同时到达，叠加成一个更强的平面波。

而在侧向（比如垂直于这条直线的方向）：

• 左边的人和右边的人到侧向某点的距离不同
• 他们发出的球面波到达时，有的波峰、有的波谷
• 一正一负，互相抵消了

数学上的结果是： 无穷多个球面波的相干叠加，恰好重构出了原来的平面波。这正是惠更斯原理的精妙之处——它看起来让波"散开"了，但通过干涉，波实际上保持了原来的传播方向。

回到我们的系统： 光瞳不是一个完整的平面波前，而是一个被圆孔截断的平面波。圆孔内每一点都发出球面波：

• 在正前方（光轴方向），球面波相干增强，形成主光斑
• 在其他方向，由于圆孔边缘的限制，球面波不能完全抵消，形成了衍射环（Airy 环）

这就是光从光瞳传播到像面时产生 PSF 的本质。

4.3.2.2. 第二步：菲涅耳衍射公式——把惠更斯原理写成数学

现在我们有了物理图像：光瞳上的每一点发出球面波，在像面上叠加。接下来我们需要把这个图像翻译成数学公式。

4.3.2.2.1. 一个点发出的球面波是什么样？

在三维空间中，一个点波源发出的球面波，在距离它 $R$ 远的地方，场的形式是：

$\text{球面波} \propto \frac{e^{ikR}}{R}$

其中：

• $k = 2\pi/\lambda$ 是波数
• $e^{ikR}$ 表示相位随距离变化（每走一个波长 $\lambda$，相位转 $2\pi$）
• $1/R$ 表示振幅随距离衰减（能量守恒，球面面积 $4\pi R^2$）

4.3.2.2.2. 从光瞳到像面：距离是多少？

光瞳上的一点 $(x_p, y_p)$ 到像面上的一点 $(x, y)$ 的距离 $R$ 为：

$R = \sqrt{(x-x_p)^2 + (y-y_p)^2 + z^2}$

其中 $z$ 是光瞳到像面的轴向距离（对理想聚焦 $z=f$，对离焦 $z=f+\Delta z$）。

4.3.2.2.3. 傍轴近似（Paraxial Approximation）

在我们的系统中，光瞳直径 $D = 60\,\text{mm}$，焦距 $f = 1248\,\text{mm}$。光线偏离光轴的角度很小（约 $D/(2f) \approx 0.024$ 弧度 $\approx 1.4°$）。

这意味着我们可以做一个近似：$(x-x_p)$ 和 $(y-y_p)$ 远小于 $z$。把 $R$ 用泰勒展开：

$\begin{aligned} R &= z\sqrt{1 + \frac{(x-x_p)^2 + (y-y_p)^2}{z^2}} \\[6pt] &\approx z \left(1 + \frac{(x-x_p)^2 + (y-y_p)^2}{2z^2}\right) \\[6pt] &= z + \frac{(x-x_p)^2 + (y-y_p)^2}{2z} \end{aligned}$

这个近似叫做傍轴近似，它把复杂的球面波传播简化为一个二次函数。

4.3.2.2.4. 把距离代入球面波公式

把 $R \approx z + \frac{(x-x_p)^2 + (y-y_p)^2}{2z}$ 代入球面波：

$\frac{e^{ikR}}{R} \approx \frac{e^{ikz}}{z} \cdot \exp\!\left[i \frac{k}{2z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right]$

因为 $k = 2\pi/\lambda$，所以 $k/(2z) = \pi/(\lambda z)$。上式变为：

$\frac{e^{ikz}}{z} \cdot \exp\!\left[i \frac{\pi}{\lambda z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right]$

这就是从光瞳上一点到像面一点的球面波传播因子。

4.3.2.2.5. 把所有点的贡献叠加起来

根据惠更斯原理，像面上 $(x, y)$ 点的总场，等于光瞳上所有点发出的球面波到达那里的叠加。

用积分表示：

$U(x, y, z) = \frac{e^{ikz}}{i\lambda z} \iint U_p(x_p, y_p) \, \exp\!\left[i \frac{\pi}{\lambda z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right] dx_p dy_p$

这就是菲涅耳衍射公式。逐项解释：

项	含义
$\displaystyle \frac{e^{ikz}}{i\lambda z}$	整体比例因子和相位偏移。$e^{ikz}$ 是传播距离 $z$ 的相位积累；$i$ 表示 $90°$ 相位旋转；$1/(\lambda z)$ 是归一化因子
$\displaystyle U_p(x_p, y_p)$	光瞳平面上的场（振幅 + 相位），也就是每个"次级波源"的强度
$\displaystyle \exp\!\left[i \frac{\pi}{\lambda z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right]$	从 $(x_p, y_p)$ 到 $(x, y)$ 的球面波传播相位。平方项 $(x-x_p)^2$ 来自傍轴近似下的距离公式
$\displaystyle \iint \cdots dx_p dy_p$	对光瞳上所有点积分——把所有球面波叠加起来

所以，公式里的平方项 $\exp[i\frac{\pi}{\lambda z}(\cdots)^2]$ 不是凭空出现的，它正是"球面波传播"在数学上的体现！ 每一点到每一点的距离不同，相位也就不同，所有不同的相位叠加在一起，就产生了衍射图案。

4.3.2.3. 第三步：展开指数项

现在我们有了菲涅耳衍射公式，接下来要做一件关键的数学操作：把指数里的平方项展开。这会让我们看到三个清晰的物理含义。

把指数项展开：

$\begin{aligned} &\exp\!\left[i \frac{\pi}{\lambda z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right] \\[6pt] =& \exp\!\left[i \frac{\pi}{\lambda z}(x^2 + y^2)\right] \cdot \exp\!\left[-i \frac{2\pi}{\lambda z}(xx_p + yy_p)\right] \cdot \exp\!\left[i \frac{\pi}{\lambda z}(x_p^2 + y_p^2)\right] \end{aligned}$

这三项分别代表：

1. $\exp\!\left[i \frac{\pi}{\lambda z}(x^2 + y^2)\right]$：像面上的二次相位（整体因子，可以提到积分外面）
2. $\exp\!\left[-i \frac{2\pi}{\lambda z}(xx_p + yy_p)\right]$：线性相位项——这一项使得积分变成傅里叶变换的形式
3. $\exp\!\left[i \frac{\pi}{\lambda z}(x_p^2 + y_p^2)\right]$：光瞳上的二次相位——这就是菲涅耳衍射与夫琅禾费衍射的区别所在

4.3.2.4. 第四步：加入薄透镜

在我们的系统中，光瞳平面紧贴透镜。透镜的作用是给光波增加一个二次相位：

$\phi_{\text{lens}}(r) = -\frac{\pi}{\lambda f} r^2$

为什么是负号？因为会聚透镜让光波向中心"收拢"，相当于给边缘更多的相位延迟（或者说让中心的相位"超前"）。

所以光瞳平面上的总场为：

$U_p(x_p, y_p) = A(r) \cdot \exp\!\left[i\phi_{\text{aberration}}(r)\right] \cdot \exp\!\left(-i\frac{\pi}{\lambda f}r^2\right)$

4.3.2.5. 第五步：合并两个二次相位——最关键的"抵消"

现在我们做整个推导中最关键的一步：把透镜引入的二次相位和菲涅耳衍射带来的二次相位合并在一起。

4.3.2.5.1. 合并前的两个二次相位

在光瞳平面上，有两个二次相位因子：

1. 菲涅耳衍射带来的（来自传播距离 $z$）：

   $\exp\!\left[i\frac{\pi}{\lambda z}(x_p^2 + y_p^2)\right]$

2. 薄透镜带来的（焦距 $f$）：

   $\exp\!\left(-i\frac{\pi}{\lambda f}(x_p^2 + y_p^2)\right]$

把它们相乘（合并）：

$\exp\!\left(-i\frac{\pi}{\lambda f}(x_p^2 + y_p^2)\right) \cdot \exp\!\left(i\frac{\pi}{\lambda z}(x_p^2 + y_p^2)\right) = \exp\!\left[i\frac{\pi}{\lambda}\left(\frac{1}{z} - \frac{1}{f}\right)(x_p^2 + y_p^2)\right]$

4.3.2.5.2. 焦平面上发生了什么？（$z = f$）

当像面恰好在焦平面时，$z = f$：

$\frac{1}{z} - \frac{1}{f} = \frac{1}{f} - \frac{1}{f} = 0$

代入合并后的相位：

$\exp\!\left[i\frac{\pi}{\lambda} \cdot 0 \cdot (x_p^2 + y_p^2)\right] = \exp(0) = 1$

二次相位完全抵消了！

4.3.2.5.3. 抵消后，积分退化成傅里叶变换

让我们把 $z = f$ 代入完整的菲涅耳衍射公式，看看发生了什么。

原始公式（回忆第二步）：

$U(x, y, z) = \frac{e^{ikz}}{i\lambda z} \iint U_p(x_p, y_p) \, \exp\!\left[i \frac{\pi}{\lambda z}\bigl((x-x_p)^2 + (y-y_p)^2\bigr)\right] dx_p dy_p$

我们已经把指数展开成了三项（第三步）：

$\exp\!\left[i \frac{\pi}{\lambda z}(x^2 + y^2)\right] \cdot \exp\!\left[-i \frac{2\pi}{\lambda z}(xx_p + yy_p)\right] \cdot \exp\!\left[i \frac{\pi}{\lambda z}(x_p^2 + y_p^2)\right]$

现在 $z = f$，并且光瞳场 $U_p$ 已经包含了透镜相位。所以被积函数中的总二次相位（光瞳上的）是：

$\underbrace{\exp\!\left[i\frac{\pi}{\lambda f}(x_p^2 + y_p^2)\right]}_{\text{菲涅耳传播}} \cdot \underbrace{\exp\!\left(-i\frac{\pi}{\lambda f}(x_p^2 + y_p^2)\right]}_{\text{透镜调制}} = 1$

它们精确抵消了！

剩下的部分：

$U(x, y, f) = \frac{e^{ikf}}{i\lambda f} \cdot \exp\!\left[i\frac{\pi}{\lambda f}(x^2 + y^2)\right] \iint U_p^{\text{原始}}(x_p, y_p) \, \exp\!\left[-i\frac{2\pi}{\lambda f}(xx_p + yy_p)\right] dx_p dy_p$

其中 $U_p^{\text{原始}}(x_p, y_p) = A(r) \cdot \exp[i\phi_{\text{aberration}}(r)]$ 是不含透镜相位的光瞳场。

现在看积分号里面的东西：

$\iint U_p^{\text{原始}}(x_p, y_p) \, \exp\!\left[-i\frac{2\pi}{\lambda f}(xx_p + yy_p)\right] dx_p dy_p$

这正是二维傅里叶变换的定义！

回忆傅里叶变换的标准形式：

$\mathcal{F}\{g(x_p, y_p)\} = \iint g(x_p, y_p) \, e^{-i2\pi(f_x x_p + f_y y_p)} dx_p dy_p$

对比可得：

• 空间频率 $f_x = \dfrac{x}{\lambda f}$
• 空间频率 $f_y = \dfrac{y}{\lambda f}$

物理意义：

• 像面上位置 $x$ 越大的点，对应的空间频率 $f_x$ 越高
• 这是因为从光瞳边缘发出的光（高频成分）会聚焦到像面的边缘

4.3.2.5.4. 为什么叫"夫琅禾费衍射"？

历史上，夫琅禾费（Joseph von Fraunhofer）发现：当观察屏距离衍射孔足够远时（远场），衍射图案就是孔径函数的傅里叶变换。

在我们的系统中，透镜起到了"压缩距离"的作用：

• 没有透镜时，你需要让光传播到无穷远才能看到夫琅禾费衍射图案
• 有了透镜，焦平面就是等效的无穷远——因为透镜把所有平行光（不同角度）聚焦到焦平面的不同位置

所以：焦平面上的场 = 光瞳场的傅里叶变换。

这正是你的代码中用 fft2 计算 PSF 的数学基础。

4.3.2.5.5. 直观理解"抵消"

为什么透镜相位和菲涅耳相位会抵消？

• 菲涅耳传播相位 $\left(+\dfrac{\pi}{\lambda z}r^2\right)$：光从光瞳传播到像面时，边缘比中心走更远的路，所以边缘相位更"滞后"。这个相位让波前变得更弯曲。

• 透镜相位 $\left(-\dfrac{\pi}{\lambda f}r^2\right)$：透镜给边缘更大的相位延迟（负号表示延迟），让平面波变得弯曲，准备会聚到焦点。

当 $z = f$ 时：

• 透镜让波前弯曲了 $-\dfrac{\pi}{\lambda f}r^2$
• 传播恰好需要 $+\dfrac{\pi}{\lambda f}r^2$ 来"展平"这个弯曲
• 两者一加一减，完美抵消，波前在焦平面上变成一个完美的点（理想情况下）

4.3.2.6. 第六步：像面移动了 $\Delta z$

现在你把相机移动了 $\Delta z$，像面位置变为 $z = f + \Delta z$。

计算：

$\frac{1}{z} - \frac{1}{f} = \frac{1}{f+\Delta z} - \frac{1}{f} = \frac{f - (f+\Delta z)}{f(f+\Delta z)} = -\frac{\Delta z}{f(f+\Delta z)}$

在傍轴近似下，$\Delta z \ll f$，所以 $f + \Delta z \approx f$，于是：

$\frac{1}{z} - \frac{1}{f} \approx -\frac{\Delta z}{f^2}$

代入相位表达式：

$\exp\!\left[i\frac{\pi}{\lambda}\left(-\frac{\Delta z}{f^2}\right)r^2\right] = \exp\!\left(-i\frac{\pi\Delta z}{\lambda f^2}r^2\right)$

4.3.2.7. 第七步：得到方法 B 的相位

取指数中的相位部分，得到方法 B 的核心表达式：

$\boxed{ \phi_B(r) = -\frac{\pi}{\lambda f^2} \, \Delta z \, r^2 }$

物理意义再解释一遍：

• 当 $\Delta z > 0$（相机远离透镜，像面在焦点后方）时，$\phi_B(r)$ 为负
• 这意味着光瞳中心相对于边缘有额外的相位超前
• 直观理解：如果像面在焦点后方，光束到达那里时还没完全会聚，所以光瞳上的波前应该"更平一些"（相比理想会聚球面波）
• 透镜原本让波前非常弯曲（要会聚到焦点），现在像面远了，所以需要在光瞳上加一个"反弯曲"的相位来补偿

4.4. 核心证明：两种方法等价性的直接比较

4.4.1. 证明思路

我们可以用一种直接的方式来证明等价性：

两种方法都是在光瞳平面上引入了一个二次相位因子，然后对同一个光瞳做完全相同的夫琅禾费衍射。如果这两个二次相位因子相同，那么像面场 $U(x,y,z)$ 必然相同，PSF 也就必然相同。

让我分别写出两种方法的光瞳场，然后直接比较。

4.4.2. 方法 A 的光瞳场

方法 A 直接在光瞳的波前上叠加 Zernike 离焦项：

$U_p^{(A)}(x_p, y_p) = A(r) \cdot \exp\!\left[i\phi_{\text{base}}(r)\right] \cdot \exp\!\left[i\phi_A(r)\right]$

其中：

• $A(r)$ 是光瞳振幅（Mask.mat 的二值圆孔）
• $\phi_{\text{base}}(r)$ 是其他像差（如随机像差）引入的基础相位
• $\phi_A(r) = \dfrac{2\pi}{\lambda} C_d \, Z_{\text{defocus}}(\rho) = \dfrac{8\pi c_2}{\lambda D^2} C_d \, r^2 + \text{const}$ 是 Zernike 离焦项引入的相位

4.4.3. 方法 B 的光瞳场

方法 B 中，光瞳平面上的原始场（不含额外离焦）为：

$U_p^{(B,\text{原始})}(x_p, y_p) = A(r) \cdot \exp\!\left[i\phi_{\text{base}}(r)\right]$

光瞳紧贴透镜，透镜引入相位 $-\dfrac{\pi}{\lambda f}r^2$。光瞳场传播到距离 $z = f + \Delta z$ 的像面时，菲涅耳衍射公式中的传播因子为 $+\dfrac{\pi}{\lambda z}r^2$（见第四步展开后的第三项）。

这两个二次相位合并后（见第五步推导）：

$\exp\!\left(-i\frac{\pi}{\lambda f}r^2\right) \cdot \exp\!\left(i\frac{\pi}{\lambda z}r^2\right) = \exp\!\left[i\frac{\pi}{\lambda}\left(\frac{1}{z} - \frac{1}{f}\right)r^2\right]$

当 $z = f + \Delta z$ 且 $\Delta z \ll f$ 时：

$\frac{1}{z} - \frac{1}{f} \approx -\frac{\Delta z}{f^2}$

所以光瞳场上的总二次相位为：

$\phi_B(r) = -\frac{\pi}{\lambda f^2} \, \Delta z \, r^2$

因此方法 B 的等效光瞳场为：

$U_p^{(B)}(x_p, y_p) = A(r) \cdot \exp\!\left[i\phi_{\text{base}}(r)\right] \cdot \exp\!\left[i\phi_B(r)\right]$

4.4.4. 直接比较两种光瞳场

对比 $U_p^{(A)}$ 和 $U_p^{(B)}$：

$\begin{aligned} U_p^{(A)} &= A(r) \cdot e^{i\phi_{\text{base}}} \cdot \exp\!\left(i\frac{8\pi c_2}{\lambda D^2} C_d \, r^2 + i\,\text{const}\right) \\[8pt] U_p^{(B)} &= A(r) \cdot e^{i\phi_{\text{base}}} \cdot \exp\!\left(-i\frac{\pi}{\lambda f^2} \Delta z \, r^2\right) \end{aligned}$

两者具有完全相同的结构：

• 相同的光瞳振幅 $A(r)$
• 相同的基础相位 $\phi_{\text{base}}$
• 都乘以一个关于 $r^2$ 的二次相位因子

唯一的区别是二次项的系数。如果我们选择合适的 $\Delta z$，使得两个二次项系数相等（允许相差常数）：

$\frac{8\pi c_2}{\lambda D^2} \, C_d = -\frac{\pi}{\lambda f^2} \, \Delta z$

两边约去 $\pi/\lambda$：

$\frac{8 c_2}{D^2} \, C_d = -\frac{1}{f^2} \, \Delta z$

解得：

$\boxed{ \Delta z = -\frac{8 f^2 c_2}{D^2} \, C_d }$

4.4.5. 结论

当 $\Delta z$ 和 $C_d$ 满足上述换算关系时：

$U_p^{(A)}(x_p, y_p) = U_p^{(B)}(x_p, y_p) \cdot e^{i\,\text{const}}$

两者相差一个全局常数相位（来自 Zernike 离焦项的 piston $c_0$）。如 3.4 节所证，全局常数相位不影响 PSF。

由于两种方法的光瞳场完全相同（允许全局相位差），它们经过完全相同的夫琅禾费衍射过程后，必然得到完全相同的像面场 $U(x,y,z)$，从而生成完全相同的 PSF。

这就严格证明了两种方法的等价性。

4.5. 数值仿真验证

4.5.1. 仿真代码做了什么？

我们写了一个 Python 程序 compare_defocus_methods_v2.py，严格按照用户代码中的参数进行仿真：

"""
对比两种离焦仿真方式（修正版）：
  方式A：加Zernike离焦项（用户现有代码逻辑）
  方式B：移动相机 —— 在光瞳上乘以菲涅耳二次相位因子（傍轴近似下严格等价于移动相机），
         然后做夫琅禾费衍射得到离焦PSF。

核心问题：真实光路中通过前后移动相机得到离焦，是否等价于在仿真中直接加Zernike离焦项？
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import os

# ---------------------------------------------------------------------------
# 1. 加载光瞳和Zernike数据（与用户代码一致）
# ---------------------------------------------------------------------------
zernikePol = loadmat("../ReconMatrix_C.mat")["ReconMatrix_C"][:, 2:].astype(np.float64)
Wx = loadmat("../Mask.mat")["Mask"].astype(np.float64)
Wx[np.isnan(Wx)] = 0.0

H, W = Wx.shape
valid_mask = Wx > 0
valid_idx = valid_mask.flatten().nonzero()[0]
N_valid = len(valid_idx)
print(f"光瞳尺寸: {H}x{W}, 有效点: {N_valid}, Zernike阶数: {zernikePol.shape[1]}")

# ---------------------------------------------------------------------------
# 2. 光学参数（与用户代码一致）
# ---------------------------------------------------------------------------
NM_TO_M = 1e-9
UM_TO_M = 1e-6
MM_TO_M = 1e-3

wavelength_m = 530 * NM_TO_M
total_folcal_len_m = 6240.79 * MM_TO_M
ratio_before_pupil = 5
pupil_ccd_focal_len_m = total_folcal_len_m / ratio_before_pupil

WFS_pitch_m = 300 * UM_TO_M
WFS_N = 10
pupil_to_WFS_ratio = 20
pupil_diameter_m = WFS_pitch_m * WFS_N * pupil_to_WFS_ratio

CCD_pixel_pitch_m = 4.5 * UM_TO_M
pad_size = 256

print(f"波长 λ = {wavelength_m*1e9:.1f} nm")
print(f"成像焦距 f = {pupil_ccd_focal_len_m*1e3:.3f} mm")
print(f"光瞳直径 D = {pupil_diameter_m*1e3:.2f} mm")
print(f"F数 = {pupil_ccd_focal_len_m / pupil_diameter_m:.2f}")

# ---------------------------------------------------------------------------
# 3. 确定Zernike离焦项形式并建立 Cd ↔ Δz 换算关系
# ---------------------------------------------------------------------------
defocus_vals = zernikePol[:, 1]
center = ((W - 1) / 2.0, (H - 1) / 2.0)
radius_px = min(center)
y_coords, x_coords = np.indices((H, W))
rho_flat = np.sqrt((x_coords - center[0])**2 + (y_coords - center[1])**2).flatten()[valid_idx] / radius_px
A = np.vstack([rho_flat**2, np.ones_like(rho_flat)]).T
a, b = np.linalg.lstsq(A, defocus_vals, rcond=None)[0]
c2 = a  # 二次项系数
print(f"\nZernike离焦项拟合: Z_defocus = {c2:.6f}*rho^2 + {b:.6f}")
print(f"理论 Noll Z4 = sqrt(3)*(2rho^2-1) = {2*np.sqrt(3):.6f}*rho^2 + {-np.sqrt(3):.6f}")

# 换算公式推导：
# 方式A相位: phi_A = (2pi/lambda) * Cd * c2 * rho^2
# 方式B相位（菲涅耳近似）: phi_B = -(pi/lambda) * (dz/f^2) * r^2
#   其中 r = rho * D/2, 所以 r^2 = rho^2 * D^2/4
#   phi_B = -(pi/lambda) * (dz/f^2) * (D^2/4) * rho^2 = -(pi*D^2 / (4*lambda*f^2)) * dz * rho^2
# 令 phi_A = phi_B:
#   (2pi/lambda) * Cd * c2 = -(pi*D^2 / (4*lambda*f^2)) * dz
#   2 * Cd * c2 = -(D^2 / (4*f^2)) * dz
#   dz = -(8 * f^2 * c2 / D^2) * Cd

def zernike_defocus_to_dz(Cd):
    return - (8.0 * pupil_ccd_focal_len_m**2 / pupil_diameter_m**2) * c2 * Cd

# 测试几个典型值
test_Cds = [(-0.5 * 0.530) * UM_TO_M, -0.1937 * UM_TO_M, -0.530 * UM_TO_M, -1.060 * UM_TO_M]
print("\n--- Zernike离焦系数 ↔ 物理离焦距离换算 ---")
for Cd in test_Cds:
    dz = zernike_defocus_to_dz(Cd)
    print(f"  Cd = {Cd/UM_TO_M:+.4f} um  <=>  dz = {dz*1e3:+.4f} mm")

# ---------------------------------------------------------------------------
# 4. 核心仿真函数
# ---------------------------------------------------------------------------
def simulate_defocus_method_A(zernike_coeffs, Cd_defocus, pad=None):
    """
    方式A：Zernike离焦项法（与用户代码一致）
    """
    zc = zernike_coeffs.copy()  # um单位
    zc[1] += Cd_defocus / UM_TO_M  # Cd_defocus [m] -> um
    zc = zc * UM_TO_M  # 转为m单位

    wavefront_m_flat = zc @ zernikePol.T
    wavefront_m = np.zeros(H * W, dtype=np.float64)
    wavefront_m[valid_idx] = wavefront_m_flat
    wavefront_m = wavefront_m.reshape(H, W)

    phase = wavefront_m / wavelength_m * 2 * np.pi
    pupil_complex = Wx * np.exp(1j * phase)

    if pad is not None and pad > max(H, W):
        M = pad
        pad_h = (M - H) // 2
        pad_w = (M - W) // 2
        padded = np.zeros((M, M), dtype=np.complex128)
        padded[pad_h:pad_h+H, pad_w:pad_w+W] = pupil_complex
        pupil_complex = padded
    else:
        M = H

    U_f = np.fft.fft2(pupil_complex, norm="ortho")
    U_f = np.fft.fftshift(U_f, axes=(-2, -1))
    intensity = np.abs(U_f)**2

    n_eff = Wx.sum()
    max_intensity_theoretical = (n_eff / M)**2
    intensity = intensity / max_intensity_theoretical

    sim_pixel_pitch_m = wavelength_m * pupil_ccd_focal_len_m * H / (M * pupil_diameter_m)
    return intensity, sim_pixel_pitch_m, M


def simulate_defocus_method_B(zernike_coeffs, Cd_defocus, pad=None):
    """
    方式B：菲涅耳近似下的移动相机法。
    在光瞳上乘以二次相位因子 exp(-i * pi/lambda * dz/f^2 * r^2)，
    这严格等价于傍轴近似下将像面沿光轴移动距离 dz。
    """
    # 基础波前（不含额外离焦）
    zc = zernike_coeffs.copy() * UM_TO_M

    wavefront_m_flat = zc @ zernikePol.T
    wavefront_m = np.zeros(H * W, dtype=np.float64)
    wavefront_m[valid_idx] = wavefront_m_flat
    wavefront_m = wavefront_m.reshape(H, W)

    # 基础相位
    base_phase = wavefront_m / wavelength_m * 2 * np.pi
    pupil_complex = Wx * np.exp(1j * base_phase)

    # 计算对应的物理离焦距离
    dz = zernike_defocus_to_dz(Cd_defocus)

    # 构建光瞳平面物理坐标网格（单位：m）
    # 光瞳直径 D，采样 HxH，像素尺寸 = D/H
    dp = pupil_diameter_m / H  # 光瞳像素尺寸
    y_idx, x_idx = np.indices((H, W))
    # 以光瞳中心为原点的物理坐标
    x_p = (x_idx - (W - 1) / 2.0) * dp
    y_p = (y_idx - (H - 1) / 2.0) * dp
    r_sq = x_p**2 + y_p**2

    # 菲涅耳二次相位因子（移动相机 dz 的等价相位）
    # 注意：对于会聚光束，像面在焦点后方 (dz>0) 时，光瞳中心相位滞后于边缘
    # 公式：phi = -pi/(lambda) * dz/f^2 * r^2
    phi_defocus = -np.pi / wavelength_m * (dz / pupil_ccd_focal_len_m**2) * r_sq

    # 乘以二次相位因子
    pupil_complex = pupil_complex * np.exp(1j * phi_defocus)

    # 零填充
    if pad is not None and pad > max(H, W):
        M = pad
        pad_h = (M - H) // 2
        pad_w = (M - W) // 2
        padded = np.zeros((M, M), dtype=np.complex128)
        padded[pad_h:pad_h+H, pad_w:pad_w+W] = pupil_complex
        pupil_complex = padded
    else:
        M = H

    # 夫琅禾费衍射 → 离焦像面
    U_f = np.fft.fft2(pupil_complex, norm="ortho")
    U_f = np.fft.fftshift(U_f, axes=(-2, -1))
    intensity = np.abs(U_f)**2

    n_eff = Wx.sum()
    max_intensity_theoretical = (n_eff / M)**2
    intensity = intensity / max_intensity_theoretical

    sim_pixel_pitch_m = wavelength_m * pupil_ccd_focal_len_m * H / (M * pupil_diameter_m)
    return intensity, sim_pixel_pitch_m, M, dz


# ---------------------------------------------------------------------------
# 5. 对比实验
# ---------------------------------------------------------------------------
zernike_zero = np.zeros(zernikePol.shape[1], dtype=np.float64)
Cd_test = -0.1937 * UM_TO_M

print("\n" + "="*60)
print(f"对比实验: Cd = {Cd_test/UM_TO_M:.4f} um")
print("="*60)

psf_A, pitch_A, M_A = simulate_defocus_method_A(zernike_zero, Cd_test, pad=pad_size)
psf_B, pitch_B, M_B, dz_actual = simulate_defocus_method_B(zernike_zero, Cd_test, pad=pad_size)
print(f"方式A (Zernike离焦): PSF shape={psf_A.shape}, pixel_pitch={pitch_A*1e6:.3f} um")
print(f"方式B (菲涅耳移动):  PSF shape={psf_B.shape}, pixel_pitch={pitch_B*1e6:.3f} um, dz={dz_actual*1e3:.4f} mm")

assert np.isclose(pitch_A, pitch_B)
assert M_A == M_B

# 定量差异
psf_A_norm = psf_A / psf_A.max()
psf_B_norm = psf_B / psf_B.max()
abs_diff = np.abs(psf_A_norm - psf_B_norm)

threshold = 0.01
mask_compare = (psf_A_norm > threshold) | (psf_B_norm > threshold)
mae = np.mean(abs_diff[mask_compare])
rmse = np.sqrt(np.mean(abs_diff[mask_compare]**2))
max_err = abs_diff[mask_compare].max()

print(f"\n差异统计 (mask > {threshold} of max):")
print(f"  MAE  = {mae:.8f}")
print(f"  RMSE = {rmse:.8f}")
print(f"  Max  = {max_err:.8f}")

mae_full = np.mean(abs_diff)
rmse_full = np.sqrt(np.mean(abs_diff**2))
print(f"\n差异统计 (全图):")
print(f"  MAE  = {mae_full:.8f}")
print(f"  RMSE = {rmse_full:.8f}")

# ---------------------------------------------------------------------------
# 6. 可视化
# ---------------------------------------------------------------------------
def add_crosshair(ax, img, color='lime'):
    cy, cx = np.array(img.shape) // 2
    ax.axhline(cy, color=color, lw=0.8, alpha=0.7)
    ax.axvline(cx, color=color, lw=0.8, alpha=0.7)

fig, axes = plt.subplots(2, 4, figsize=(16, 8))

vmax = max(psf_A_norm.max(), psf_B_norm.max())

axes[0, 0].imshow(psf_A_norm, cmap='hot', vmin=0, vmax=vmax)
axes[0, 0].set_title(f"Method A: Zernike defocus\nCd={Cd_test/UM_TO_M:.4f} um")
axes[0, 0].axis('off')
add_crosshair(axes[0, 0], psf_A_norm)

axes[0, 1].imshow(psf_B_norm, cmap='hot', vmin=0, vmax=vmax)
axes[0, 1].set_title(f"Method B: Fresnel propagation\ndz={dz_actual*1e3:.3f} mm")
axes[0, 1].axis('off')
add_crosshair(axes[0, 1], psf_B_norm)

im2 = axes[0, 2].imshow(abs_diff, cmap='jet', vmin=0, vmax=max_err)
axes[0, 2].set_title(f"Absolute diff |A-B|\nMAE={mae:.6f}, Max={max_err:.6f}")
axes[0, 2].axis('off')
add_crosshair(axes[0, 2], abs_diff)

rel_diff = abs_diff / (psf_B_norm + 1e-12)
rel_diff_vis = np.log10(rel_diff + 1e-6)
axes[0, 3].imshow(rel_diff_vis, cmap='jet')
axes[0, 3].set_title("Relative diff log10(|A-B|/B+1e-6)")
axes[0, 3].axis('off')
add_crosshair(axes[0, 3], rel_diff)

# 中心剖面
cy = psf_A_norm.shape[0] // 2
axes[1, 0].plot(psf_A_norm[cy, :], 'r-', lw=1.5, label='Method A (Zernike)')
axes[1, 0].plot(psf_B_norm[cy, :], 'b--', lw=1.5, label='Method B (Fresnel)')
axes[1, 0].set_title(f"Horizontal profile (y={cy})")
axes[1, 0].set_xlabel("pixel")
axes[1, 0].set_ylabel("normalized intensity")
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(abs_diff[cy, :], 'g-', lw=1.5)
axes[1, 1].set_title("Horizontal profile diff")
axes[1, 1].set_xlabel("pixel")
axes[1, 1].set_ylabel("|A-B|")
axes[1, 1].grid(True, alpha=0.3)

cx = psf_A_norm.shape[1] // 2
axes[1, 2].plot(psf_A_norm[:, cx], 'r-', lw=1.5, label='Method A (Zernike)')
axes[1, 2].plot(psf_B_norm[:, cx], 'b--', lw=1.5, label='Method B (Fresnel)')
axes[1, 2].set_title(f"Vertical profile (x={cx})")
axes[1, 2].set_xlabel("pixel")
axes[1, 2].set_ylabel("normalized intensity")
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

# 径向平均
def radial_profile(img):
    y, x = np.indices(img.shape)
    center = np.array(img.shape) // 2
    r = np.sqrt((x - center[1])**2 + (y - center[0])**2).astype(np.int64)
    r_max = min(center[0], center[1])
    radial_mean = np.zeros(r_max)
    for i in range(r_max):
        mask_r = (r == i)
        if mask_r.sum() > 0:
            radial_mean[i] = img[mask_r].mean()
    return radial_mean

rad_A = radial_profile(psf_A_norm)
rad_B = radial_profile(psf_B_norm)
axes[1, 3].plot(rad_A, 'r-', lw=1.5, label='Method A (Zernike)')
axes[1, 3].plot(rad_B, 'b--', lw=1.5, label='Method B (Fresnel)')
axes[1, 3].set_title("Radial average")
axes[1, 3].set_xlabel("radius (pixel)")
axes[1, 3].set_ylabel("normalized intensity")
axes[1, 3].legend()
axes[1, 3].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('comparison_defocus_methods_v2.png', dpi=200, bbox_inches='tight')
print("\nSaved: comparison_defocus_methods_v2.png")

# ---------------------------------------------------------------------------
# 7. 不同离焦量扫掠
# ---------------------------------------------------------------------------
print("\n" + "="*60)
print("Error vs defocus amount:")
print("="*60)

defocus_ums = np.linspace(-1.0, 1.0, 21) * 0.530
results = []
for Cd_um in defocus_ums:
    Cd = Cd_um * UM_TO_M
    psf_A, _, _ = simulate_defocus_method_A(zernike_zero, Cd, pad=pad_size)
    psf_B, _, _, _ = simulate_defocus_method_B(zernike_zero, Cd, pad=pad_size)
    psf_A_n = psf_A / psf_A.max()
    psf_B_n = psf_B / psf_B.max()
    diff = np.abs(psf_A_n - psf_B_n)
    mae_val = diff.mean()
    rmse_val = np.sqrt(np.mean(diff**2))
    results.append((Cd_um, mae_val, rmse_val))
    print(f"  Cd={Cd_um:+.3f} um  MAE={mae_val:.8f}  RMSE={rmse_val:.8f}")

results = np.array(results)

fig2, ax2 = plt.subplots(1, 1, figsize=(8, 5))
ax2.plot(results[:, 0], results[:, 1], 'bo-', label='MAE', markersize=4)
ax2.plot(results[:, 0], results[:, 2], 'rs-', label='RMSE', markersize=4)
ax2.axvline(x=-0.1937, color='lime', linestyle='--', alpha=0.7, label='Trained Cd=-0.1937um')
ax2.set_xlabel("Zernike defocus Cd (um)")
ax2.set_ylabel("PSF difference")
ax2.set_title("Method A vs Method B error")
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_yscale('log')
plt.tight_layout()
plt.savefig('defocus_sweep_error_v2.png', dpi=200, bbox_inches='tight')
print("Saved: defocus_sweep_error_v2.png")

# ---------------------------------------------------------------------------
# 8. 含随机像差的对比
# ---------------------------------------------------------------------------
np.random.seed(42)
zernike_random = np.random.randn(zernikePol.shape[1]) * 0.05  # um
zernike_random[1] = 0.0

Cd_rand = -0.1937 * UM_TO_M
psf_A_rand, _, _ = simulate_defocus_method_A(zernike_random, Cd_rand, pad=pad_size)
psf_B_rand, _, _, _ = simulate_defocus_method_B(zernike_random, Cd_rand, pad=pad_size)

psf_A_rand_n = psf_A_rand / psf_A_rand.max()
psf_B_rand_n = psf_B_rand / psf_B_rand.max()
diff_rand = np.abs(psf_A_rand_n - psf_B_rand_n)
mae_rand = diff_rand.mean()
rmse_rand = np.sqrt(np.mean(diff_rand**2))

print("\n" + "="*60)
print(f"With random aberrations (Cd=-0.1937um): MAE={mae_rand:.8f}, RMSE={rmse_rand:.8f}")
print("="*60)

fig3, axes3 = plt.subplots(1, 3, figsize=(14, 4))
vmax_r = max(psf_A_rand_n.max(), psf_B_rand_n.max())
axes3[0].imshow(psf_A_rand_n, cmap='hot', vmin=0, vmax=vmax_r)
axes3[0].set_title("Method A: Zernike (with aberrations)")
axes3[0].axis('off')
axes3[1].imshow(psf_B_rand_n, cmap='hot', vmin=0, vmax=vmax_r)
axes3[1].set_title("Method B: Fresnel (with aberrations)")
axes3[1].axis('off')
im3 = axes3[2].imshow(diff_rand, cmap='jet')
axes3[2].set_title(f"Absolute diff |A-B|\nMAE={mae_rand:.6f}")
axes3[2].axis('off')
plt.colorbar(im3, ax=axes3[2], fraction=0.046)
plt.tight_layout()
plt.savefig('comparison_with_aberrations_v2.png', dpi=200, bbox_inches='tight')
print("Saved: comparison_with_aberrations_v2.png")

print("\nAll done.")

• 光瞳：Mask.mat 的 $51\times51$ 二值圆孔
• 零填充到 $256\times256$ 以提高采样率
• 波长 $\lambda = 530\,\text{nm}$
• 焦距 $f = 1248.158\,\text{mm}$
• 光瞳直径 $D = 60\,\text{mm}$
• FFT 使用 ortho 归一化（能量守恒）

4.5.2. 纯离焦对比（$C_d = -0.1937\,\mu\text{m}$）

下图对比了两种方法生成的离焦 PSF：

看图说话：

• 左图（方法 A） 和 中图（方法 B） 看起来几乎一模一样，都是一个中心亮斑加上外围的同心圆环
• 右图（差异图） 用蓝色表示差异小、红色表示差异大。可以看到大部分区域是深蓝色，说明两者非常接近

定量统计：

指标	数值	含义
全图 MAE	$0.0019$	平均每像素差异不到 0.2%
有效区域 MAE（>1% max）	$0.0259$	在亮斑区域内平均差异约 2.6%
全图 RMSE	$0.0156$	均方根误差约 1.6%
最大差异	$0.1931$	极个别像素差异约 19%

为什么还有微小差异？

1. 离散采样误差：$51\times51$ 的光瞳网格导致拟合的 $c_2 = 3.3296$ 与理论值 $3.4641$ 有约 $3.9\%$ 偏差。这意味着我们换算的 $\Delta z$ 有微小误差。
2. FFT 数值精度：256 点 FFT 的有限精度会引入约 $10^{-15}$ 量级的舍入误差。
3. 常数相位的数值处理：虽然理论上常数相位不影响 PSF，但在数值计算中，$e^{i \times 3974}$ 这样的大相位旋转可能引入微小误差。

4.5.3. 不同离焦量扫掠

我们把 $C_d$ 从 $-1\lambda$ 变到 $+1\lambda$（步长 $0.05\lambda$），看看误差如何变化：

看图说话：

• 横轴是 Zernike 离焦系数 $C_d$（单位 $\mu\text{m}$）
• 纵轴是两种方法生成 PSF 的差异（对数坐标）
• 在 $C_d = 0$ 处，误差严格为 0（两个都没有离焦，当然一样）
• 随着离焦量增大，误差缓慢上升
• 在整个 $\pm 1\lambda$ 范围内，误差都很小（MAE $< 0.05$）
• 绿虚线标记了你训练后的实际工作点 $C_d = -0.1937\,\mu\text{m}$，此处误差处于较低水平

4.5.4. 加入随机像差后的对比

为了验证等价性在更复杂情况下是否依然成立，我们给光瞳加入了一些随机的小像差（模拟真实系统中除了离焦之外的其他不完美）：

看图说话：

• 左图和方法 B 的 PSF 都呈现出复杂的斑点结构（因为有随机像差）
• 但两者的整体形态高度一致
• 差异图 mostly 是深蓝色，说明差异很小

定量统计：

指标	数值
MAE	$0.0030$
RMSE	$0.0069$

有趣的是，加入随机像差后，两种方法的差异反而更小了（从 MAE $0.026$ 降到 $0.003$）。这是因为随机像差作为"共同背景"，在一定程度上掩盖了两种离焦模型之间的微小数值差异。

5. 光瞳仿真中的傅里叶变换：从周期延拓的数学宿命到波动光学的衍射极限

5.1. 从夫琅禾费衍射到 FFT：物理定律的数值翻译

在物理光学中，计算一个光学系统焦平面上的光场分布，通常始于夫琅禾费衍射积分（Fraunhofer Diffraction Integral）。对于位于透镜前焦面（或光瞳面）的复振幅分布 $P(x,y)$，其后焦面上的复振幅 $U_f(x', y')$ 正比于 $P(x,y)$ 的二维连续傅里叶变换（Continuous Fourier Transform, CFT）：

$U_f(f_x, f_y) \propto \iint_{-\infty}^{+\infty} P(x,y) \, e^{-i 2\pi (f_x x + f_y y)} \, dx \, dy$

其中：

• $P(x,y) = A(x,y) \cdot e^{i W(x,y)}$ 是光瞳函数（Pupil Function）
• $A(x,y)$ 是孔径透过率（光瞳内通常为 1，光瞳外严格为 0）
• $W(x,y)$ 是波前像差（可能包含倾斜 $ax+by$、离焦、像散等）
• $f_x = x'/(\lambda f)$, $f_y = y'/(\lambda f)$ 是空间频率坐标，$\lambda$ 为波长，$f$ 为透镜焦距

这个积分的物理意义非常深刻：焦平面上的每一点，都是整个光瞳面上所有子波源发出的球面波相干叠加的结果。透镜的作用，本质上是把远场的夫琅禾费衍射图样聚焦到了它的后焦平面上，从而把原本需要传播很远距离才能看到的衍射斑，"压缩"到了焦平面附近。

5.1.1. 为什么计算机需要 DFT？

连续傅里叶变换是一个无穷积分，涉及连续函数和无限域。计算机无法直接处理无穷和连续，因此必须进行三重离散化：

物理现实	数值计算	引入的代价
连续空间 $(x,y) \in \mathbb{R}^2$	离散采样网格 $(x_n, y_m)$	采样定理约束：空间细节必须被足够密的网格捕捉
积分限无穷大	有限求和（$N \times N$ 网格）	截断误差：用有限区域近似无穷积分
直接积分 $O(N^4)$ 复杂度	快速傅里叶变换（FFT）$O(N^2 \log N)$	周期延拓假设：DFT 的数学根基

当我们写下：

U_f = torch.fft.fft2(pupil_complex)

我们实际上是在计算二维离散傅里叶变换（2D-DFT）：

$U_f[k,l] = \sum_{n=0}^{N-1} \sum_{m=0}^{N-1} P[n,m] \, e^{-i 2\pi \left(\frac{kn}{N} + \frac{lm}{N}\right)}$

这个公式在数学上是对 CFT 的黎曼和近似。当采样足够密（满足奈奎斯特采样定理）、网格足够大时，DFT 的结果可以收敛到 CFT 的结果。

5.2. DFT 的数学宿命：周期延拓不是物理假设，而是算法根基

5.2.1. 紧支集的物理现实

在真实物理中，光瞳函数 $P(x,y)$ 具有紧支集（Compact Support）：光瞳内（透镜口径范围内）有光场分布，光瞳外严格为 0。

这不是数值近似，而是物理现实。 透镜的边框、镜筒、光阑（Aperture Stop）把光挡住了。光瞳外没有光能参与成像。因此，真实的夫琅禾费衍射积分虽然在形式上写成了 $\int_{-\infty}^{+\infty}$，但由于 $P(x,y)$ 在光瞳外为 0，积分实际上只在有限区域内非零：

$U_f(f_x, f_y) \propto \iint_{\text{光瞳}} P(x,y) \, e^{-i 2\pi (f_x x + f_y y)} \, dx \, dy$

这是一个紧支集函数的连续傅里叶变换。从数学分析可知，紧支集光滑函数的傅里叶变换是整函数（Entire Function），在频域上是无限延伸且无限可微的，没有任何周期性。

5.2.2. DFT 天生的周期性

离散傅里叶变换的公式定义在周期序列上。DFT 从不"认识"有限序列，它只认识周期为 $N$ 的无限序列的一个主值周期：

$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-i 2\pi kn/N}$

这个公式的隐含前提是：输入序列 $x[n]$ 以 $N$ 为周期无限重复。换句话说，DFT 计算的不是"一个有限图像的频谱"，而是这个图像作为无限周期瓷砖阵列的一个周期单元的离散频谱。

这意味着：

• 空间域周期化：你的 $N \times N$ 图像被数学强制理解为在上下左右四个方向无限平铺。
• 频率域周期化：输出的频谱也是周期为 $N$ 的离散序列。

周期延拓不是你可以通过某种参数"关闭"的选项，它是你用 FFT 这个工具时必须接受的数学现实。 FFT 算法之所以能达到 $O(N^2 \log N)$ 的复杂度（而非直接积分的 $O(N^4)$），正是因为它利用了这种周期性结构（通过单位根的周期性简化计算）。

5.2.3. 为什么不能直接"假设光瞳外为 0"？

DFT 的基函数是复指数函数 $e^{-i 2\pi kn/N}$，这些基函数本身是**全局支撑（global support）**的——它们在整个周期上非零。DFT 本质上是把输入信号投影到这组正交周期基上。你无法用一组周期基去"表达"一个非周期边界条件（紧支集），除非通过某种方式把边界推到足够远，使其影响可忽略。

换句话说，当你使用 FFT 时，你已经在数学上"签约"接受了周期世界。紧支集的物理现实，必须在这个周期框架内被重建出来，而不是被直接施加。

5.3. 零填充：在周期框架内重建紧支集物理

既然周期延拓不可消除，那么如何让它对真实物理的干扰最小化？答案就是零填充（Zero-Padding）。

5.3.1. 零填充的数学本质

零填充的操作极其简单：把原本大小为 $D \times D$ 的光瞳函数，嵌入到一个更大的 $N \times N$ 数组（$N > D$）的中央，光瞳外的区域全部填充 0。

# 假设光瞳是 D x D，嵌入到 N x N 的网格中
padded = torch.zeros(N, N, dtype=torch.complex64, device=device)
start = (N - D) // 2
padded[start:start+D, start:start+D] = pupil_complex

# 执行 FFT
U_f = torch.fft.fftshift(torch.fft.fft2(padded, norm='ortho'))

在数学上，这等价于对原始光瞳函数 $P(x,y)$ 进行上采样（Upsampling）——即在更细的频率网格上计算其傅里叶变换。从卷积定理的角度看，零填充后的 DFT 输出，是原始光瞳 DFT 与 sinc 函数的卷积，但由于零填充增加了采样密度，这个 sinc 函数在频域上被压缩，从而更精确地采样了真实的连续频谱。

5.3.2. 为什么零填充能消除周期延拓的伪影？

考虑不做零填充的情况：光瞳填满整个 $N \times N$ 网格。此时 DFT 的周期延拓会在边界处产生硬截断跳变：

周期单元（光瞳填满网格）：
[ ████ ] [ ████ ] [ ████ ]
[ ████ ] [ ████ ] [ ████ ]
       ↑
    边界跳变：光瞳右边缘的值直接拼接左边缘的值

如果光瞳在边界处不是周期连续的（例如圆形光瞳被方形网格截断），这种跳变在数学上是一个不连续函数，其傅里叶变换会产生严重的吉布斯现象（Gibbs Phenomenon）和频谱泄漏（Spectral Leakage）。此时 FFT 计算的不是"单个孔径的衍射"，而是无限大周期孔径阵列（Periodic Aperture Array）的衍射——这在物理上完全错误。

而做了零填充之后：

周期单元（光瞳只占网格中心一小部分）：
[ 000000000000000 ]
[ 000  光瞳  0000 ]
[ 000000000000000 ]
       ↑
    边界跳变：0 与 0 相接，没有不连续

此时，周期延拓的相邻"瓷砖"之间被大量 0 隔开，拼接处是 0 与 0 相接，没有任何跳变不连续。在这个巨大的周期单元内部，光瞳被 0 包围，数学上完美复现了"紧支集"的物理图景。

周期延拓引入的"假邻居"被推到很远的地方，其频谱泄漏主要落在高频区域，而你的有效频带内几乎不受影响。

5.3.3. 零填充不会增加物理信息，但提高数值精度

这是一个常见的误解：零填充能提高分辨率吗？

不能。 零填充不会增加任何物理上不存在的信息，也不会突破衍射极限。光学系统的物理分辨率（瑞利判据）由光瞳口径 $D$ 和波长 $\lambda$ 决定：

$\Delta x \approx 1.22 \frac{\lambda f}{D}$

零填充改变不了这个极限。它改善的是数值采样精度：

• 不做零填充（即使不做0填充，实际上圆形光瞳外还是有0填充的，也会有一个被污染的艾里斑）：FFT 输出可能只采样了 sinc/Airy 函数的 3~4 个点，峰值位置可能被网格量化误差偏移，旁瓣结构扭曲。
• 做零填充（在数学上逼近了紧支集假设）：同一个衍射斑被采样了几十甚至上百个点，峰值位置更准确，旁瓣的对称性和深度更干净，能量守恒更容易满足。

5.3.4. 实验验证：零填充对 PSF 重建的影响

为了直观展示零填充对仿真精度的影响，下面是同一组波前（原始尺寸 $51\times51$）在两种处理方式下的对比结果。上排为实验采集的真实点扩散函数（Orig），下排为对应波前生成的仿真点扩散函数（Sim）。

将波前零填充至 $256\times256$ 后仿真：

不做零填充（原始 $51\times51$）直接仿真：

5.3.5. 填充多少合适？

填充倍数	光瞳占网格比例	适用场景
1×（不填充）	100%	仅用于快速测试，会引入严重混叠
2×	25%	快速估算，混叠基本可控，旁瓣可见
4×	6.25%	标准精度，旁瓣干净，能量守恒良好
8×	1.56%	极高精度，或光瞳边缘相位变化剧烈（如高阶像差）

经验法则：如果你看到焦平面的 PSF 在边缘处有奇怪的振荡，或者总能量在正逆变换后不守恒，通常就是零填充不够。

5.3.6. 归一化的选择

torch.fft.fft2 的 norm 参数决定了正逆变换的归一化方式：

• norm='backward'（默认）：正变换不归一化，逆变换除以 $N^2$。这是传统信号处理用法。
• norm='ortho'（推荐）：正逆变换都除以 $N$（即 $1/\sqrt{N^2}$）。这保证了**帕塞瓦尔定理（Parseval’s Theorem）**成立——空间域和频率域的能量相等。

在光学仿真中，能量守恒通常是重要的物理约束（光瞳内的总功率应等于焦平面上的总功率），因此推荐使用 norm='ortho'。

这里你 pad 了之后并不会影响最终仿真图像的大小，因为 pad 只是让频域的采样点更密（频率分辨率更高），它并不改变 FFT 能分析到的最大频率。真正决定最大频率的是光瞳面的采样率 $\Delta x$——$\Delta x$ 越小，频域覆盖越宽；$\Delta x$ 不变，最大频率就不变。而焦平面的总尺寸正比于这个最大频率，所以无论你 pad 多少，仿真输出的像素数虽然变多了、每个像素变小了，但总的物理视场始终不变。再经过 scale 缩放到物理靶面坐标后，最终图像的大小自然也是一样的。换句话说，pad 改变的是“采样密度”，不是“采样范围”；图像的物理尺寸由采样范围决定，因此与 pad 无关。

量	不 padding（$N_1$）	Padding 4×（$N_2=4N_1$）	是否改变
光瞳采样间隔 $\Delta x$	$\Delta x$	$\Delta x$	❌ 不变
频域总宽度 $1/\Delta x$	$1/\Delta x$	$1/\Delta x$	❌ 不变
焦平面总视场 $\lambda f/\Delta x$	FOV	FOV	❌ 不变
频率步长 $\Delta f = 1/(N\Delta x)$	$1/L_1$	$1/L_2 = 1/(4L_1)$	✅ 变密 4 倍
焦平面像素大小 $\lambda f \cdot \Delta f$	$p_1$	$p_2 = p_1/4$	✅ 变小 4 倍
输出数组像素数	$N_1$	$N_2 = 4N_1$	✅ 变多 4 倍
总物理尺寸 $N \cdot p$	$N_1 p_1 = \text{FOV}$	$N_2 p_2 = 4N_1 \cdot (p_1/4) = \text{FOV}$	❌ 不变

5.4. 无限平面波 vs 有限孔径：几何光学与波动光学的分水岭

5.4.1. 几何光学的直觉

在几何光学（射线光学，Ray Optics）中，一个倾斜的平行平面波被理想透镜汇聚后，会在焦平面上形成一个几何点，其位置仅由倾斜角决定：

$x_f \propto a, \quad y_f \propto b$

如果截取平面波的中间一部分（通过有限口径），几何光学的结论是：焦点位置不变，只是参与聚焦的光线数量减少，因此像点光强减弱。

这个直觉在工程上非常有用，它是波动光学在波长 $\lambda \to 0$ 时的极限。当孔径远大于波长时，衍射斑非常小，肉眼或常规仪器看来就像一个"点"。

5.4.2. 波动光学的回应：截断改变了一切

但在**波动光学（Wave Optics）**中，有限孔径的衍射效应不可忽略。被透镜口径截断的倾斜平面波，其数学形式是：

$P(x,y) = \text{rect}\left(\frac{r}{D}\right) \cdot e^{i(ax+by)}$

根据卷积定理，它的傅里叶变换是孔径函数的变换与平面波变换的卷积：

$\mathcal{F}\{P\} = \mathcal{F}\{\text{rect}\} * \mathcal{F}\{e^{i(ax+by)}\} = \left[D^2 \cdot \text{sinc}(D f_x) \text{sinc}(D f_y)\right] * \delta\left(f_x - \frac{a}{2\pi}, f_y - \frac{b}{2\pi}\right)$

这意味着：理想点被孔径的 sinc 函数"涂抹"开了。焦平面上的光场不是一个点，而是一个中心偏移的 sinc 函数（矩形光瞳）或 Airy 斑（圆形光瞳）。

5.4.3. 为什么截断会改变"聚焦行为"？

从惠更斯-菲涅尔原理（Huygens-Fresnel Principle）看，焦平面上任意一点的光场，是整个孔径上所有子波源发出的球面子波相干叠加的结果。

当你把孔径截断时，你不是简单地"扔掉一些光线"（几何光学视角），而是删除了一大批参与干涉的相干子波源。这些被移除的子波源原本负责在焦平面周围形成精确的相位抵消——它们定义了"理想点"的边界。它们消失后，能量就从中心峰值泄漏到了周围区域，形成旁瓣。

这不是"补充 0 带来的数值干扰"，而是物理衍射——即使你用解析积分而不是 FFT，结果也是 sinc 斑。零填充只是让 FFT 更准确地计算这个 sinc 斑，而不是引入它。

5.4.4. 理想点需要什么条件？

要在物理上真正得到理想的几何点，需要同时满足两个不可实现的条件：

1. 无限大的平面波：波前完全平坦，无边界截断，在全空间 $\mathbb{R}^2$ 上都是 $e^{i(ax+by)}$。
2. 无限大的透镜/光学系统：没有任何口径限制，能接收无限大波前的全部能量。

只要其中任何一个变成有限，结果就立即退化为有限孔径衍射。

所以，"无限平面波聚焦为理想点"是一个数学极限，物理上永远无法实现。 任何真实光学系统（望远镜、镜头、人眼、显微镜）的焦平面上，你看到的永远是衍射斑，只是大小不同：

• 大望远镜（口径数米）：Airy 斑可能只有几微米，但绝不是 0 尺寸。
• 小孔径：衍射非常明显（如针孔相机、手机镜头夜景的星芒）。

5.5. 结语

当我们写下 torch.fft.fft2(pupil_complex) 时，我们实际上同时操作着三个层面：

代码层：一个高效的矩阵运算，利用 FFT 算法在毫秒级完成数百万点的变换。

数学层：一个周期序列的离散傅里叶级数系数计算。DFT 的基函数——复指数 $e^{-i 2\pi kn/N}$——是全局支撑的周期函数，它们天生不理解"边界"和"遮挡"。周期延拓不是 bug，而是这个数学空间的结构常数。

物理层：光波穿过有限孔径后，在远场形成的带有衍射极限的干涉图样。透镜的有限口径就是物理光阑，光瞳外严格为 0。紧支集、有限能量、边界衍射，这些是物理现实强加给我们的约束。

数值光学工程的核心智慧，不在于消除数学与物理之间的张力——这种张力是内禀的、不可消除的——而在于深刻理解张力双方的结构，并用精巧的技巧（零填充、过采样、正交归一化）在数学框架内重建物理真实。

几何光学给了我们简洁而强大的直觉：光线汇聚于一点。波动光学则提醒我们：只要孔径有限、波长非零，世界就永远是模糊的。衍射不是计算的缺陷，不是 FFT 的伪影，而是光的内禀属性——它是光作为波的最诚实告白。

“光在传播时，每一寸边界都在诉说它的存在。你无法截断一束波而不留下痕迹，正如你无法在时空中划定边界而不产生涟漪。”

6. 从连续到离散：傅里叶变换家族全解析

6.1. 连续世界：无限与有限的根本分歧

6.1.1. 无限区间：经典的连续傅里叶变换（CFT）

当信号 $f(x)$ 定义在整个实数域上时，我们使用最标准的傅里叶变换对：

正变换：

$\hat{f}(\xi) = \int_{-\infty}^{\infty} f(x) \, e^{-2\pi i x \xi} \, dx$

逆变换：

$f(x) = \int_{-\infty}^{\infty} \hat{f}(\xi) \, e^{2\pi i x \xi} \, d\xi$

这是理论最完美的情况：时域无限，频域连续。但现实中，我们永远无法处理无限长的信号。

6.1.2. 有限区间：两条截然不同的路径

当信号仅在 $[a, b]$ 上有定义时，数学上出现了分叉：

路径 A：紧支集假设（Truncation） 直接将积分限改为有限区间，假设信号在区间外严格为零：

$\hat{f}(\xi) = \int_{a}^{b} f(x) \, e^{-2\pi i x \xi} \, dx$

• 频谱特征：连续的、解析的（无限可微），呈现 sinc 函数型振荡
• 逆变换结果：若使用完整频谱 $\hat{f}(\xi)$（$\xi \in (-\infty, \infty)$），重构信号在 $[a,b]$ 外精确为 0
• 物理意义：适用于瞬态信号（如脉冲、衰减振动）

关于为什么叫紧支集的解释：函数的支集（函数值不为零的所有点的闭包）本身就是一个紧集 （在 $\mathbb{R}^n$ 中即有界闭集），意味着非零区域被限制在一个有限且封闭的范围内，外部恒为零。

路径 B：周期延拓（Periodic Extension） 将有限区间视为周期信号的一个周期，使用傅里叶级数：

$f(x) = \sum_{n=-\infty}^{\infty} c_n \, e^{i \frac{2\pi}{T}nx}, \quad c_n = \frac{1}{T}\int_{-T/2}^{T/2} f(x) \, e^{-i \frac{2\pi}{T}nx} \, dx$

• 频谱特征：离散的频率点 $\xi_n = n/T$，谐波结构
• 逆变换结果：周期函数，在区间外循环重复原波形
• 物理意义：适用于 inherently 周期的现象（如声波、交流电）

这里两种方式都是有限区间内的积分，只是一个除以了区间长度，一个没有除，为什么就有定义域外默认为0,和定义域外默认为周期延拓的区别？

6.1.2.1. 基函数的内禀属性

两种方法确实都在有限区间上积分，区别似乎只是归一化因子 $1/T$。但这只是单位换算的技术细节，真正决定延拓行为的，是两种方法选用的**基函数（Basis Functions）**具有根本不同的数学属性。

路径 A（傅里叶变换）的基函数：

$e^{-2\pi i x \xi}, \quad \xi \in \mathbb{R}$

这些基函数是非周期的平面波，在 $|x| \to \infty$ 时既不衰减也不归零。当你进行逆变换时：

$f(x) = \int_{-\infty}^{\infty} \hat{f}(\xi) e^{2\pi i x\xi} \, d\xi$

本质上是用无限多个无限延伸的平面波来合成原信号。要使这个合成结果在 $[a,b]$ 之外精确为零（紧支集），需要所有基波在区间外通过相位干涉完美相消——这要求频谱 $\hat{f}(\xi)$ 必须是特定的解析函数（整函数）。换言之，紧支集特性是频谱分析的结果，而非随意假设。

路径 B（傅里叶级数）的基函数：

$e^{i \frac{2\pi}{T} n x}, \quad n \in \mathbb{Z}$

这些基函数天生具有周期 $T$：

$e^{i \frac{2\pi}{T} n (x+T)} = e^{i \frac{2\pi}{T} n x}$

无论系数 $c_n$ 如何选择，它们的线性组合必然满足：

$f(x+T) = f(x)$

因此，当你写出：

$f(x) = \sum_{n=-\infty}^{\infty} c_n e^{i \frac{2\pi}{T} n x}$

周期性已经内禀地蕴含在基函数之中，不是可选项，而是数学必然。

6.1.2.2. 更深层的差异：函数生活在哪里？

问题的关键在于：当你写下有限区间的积分时，你正在对什么对象做积分？

紧支集方法：函数 $f$ 定义在 $\mathbb{R}$（整条实数轴） 上。你写下 $\int_a^b$ 时，实际上是 $\int_{-\infty}^\infty$ 的偷懒写法，因为 $f(x)$ 在 $[a,b]$ 外已经被定义为零。数学表述：$f \in L^2(\mathbb{R})$ 且 $\text{supp}(f) \subseteq [a,b]$。"区间外为零"是前提假设，积分只是验证这个假设的副产品。

周期延拓方法：函数 $f$ 定义在 $\mathbb{T}_T$（周长为 $T$ 的圆环/周期域） 上，或者等价地，是 $\mathbb{R}$ 上满足 $f(x+T)=f(x)$ 的周期函数。你写下 $\int_{-T/2}^{T/2}$ 时，是在圆环上的一个 Fundamental Domain（基本域）上积分，不是在"实数轴的一个区间"上积分。数学表述：$f \in L^2(\mathbb{T}_T)$，这里 $x+T$ 就是 $x$ 本身，没有"区间外"的概念。"周期重复"是定义域的内在结构，不是延拓出来的。

特征	紧支集方法	周期延拓方法
函数生活在哪里	实数轴 $\mathbb{R}$（开集，无限延伸）	圆环 $\mathbb{T}_T$（紧集，首尾粘连）
积分含义	截断整个实数轴（外面已经没东西了）	遍历整个圆环（选一个代表元区间）
$f(0)$ 和 $f(T)$ 的关系	没关系，可能 $f(0)\neq 0$ 而 $f(T)=0$	$f(0) = f(T)$ 是强制性等同，因为它们是同一个点
逆变换结果	支集外为 0	支集外周期性复制

6.1.2.3. 频谱后果的对比

当首尾值不相等或者首尾值不等于0时，周期延拓或者紧支都有可能产生不连续性， 从而在频域产生高频分量，但形态不同。

6.1.2.4. 傅里叶不确定性原理

当你只能观测或采集有限长的一段信号（有限时间窗口 $T$，或有限孔径 $D$）时，数学上，有限窗口等价于给信号乘了一个矩形窗，而矩形窗的傅里叶变换是一个 sinc 函数（$\sin x / x$ 形状）。根据卷积定理，真实频谱会被这个 sinc 函数“涂抹”。

6.1.2.4.1. 连续傅里叶变换（如光学透镜）

透镜的孔径是有限的，相当于空间域的矩形窗。焦平面上的频谱仍然是连续分布的，你理论上可以读取任意频率点的值，不存在离散的“频率间隔”。

但是，真实频谱已经被孔径的 sinc 谱涂抹了：

• 一个理想的点光源（空间 delta 函数）成像后，变成一个艾里斑或 sinc 斑；
• 两个靠得很近的空间频率成分，会因主瓣重叠而糊成一片，无法区分。

表现：没有离散的频率格子，但存在最小可分辨间隔（瑞利判据），其量级约为 $1/D$（$D$ 为孔径）。频谱是连续的，但“模糊”了。

6.1.2.4.2. 离散傅里叶变换（DFT/FFT）

在数值计算中，你不仅时间有限（$T$），还只采了 $N$ 个点。有限窗口同样导致 sinc 涂抹，但 DFT 额外做了一步：它在频域只取 $N$ 个采样点。

sinc 涂抹的主瓣宽度约为 $1/T$，而 DFT 的频点间隔恰好也是：

$\Delta f = \frac{1}{T} = \frac{f_s}{N}$

表现：频域被强制画成了等间距的离散网格。你看到的不是连续的模糊带，而是只能落在这 $N$ 根谱线上的采样值。两根谱线之间的细节被“格子”锁死了；即使你用补零（Zero Padding）让曲线看起来更平滑，也只是插值，真实分辨率仍然由 $1/T$ 决定。

6.1.2.4.3. 为什么频率间隔为$1/T$

有限窗口 $[0, T]$ 内的傅里叶分析，本质上是在问：两个不同频率的复指数信号，在这个区间内是否"独立"？

数学上，"独立"的标准是正交——它们的内积必须为零（才能将其区分开）。

6.1.2.4.3.1. 第一步：定义内积

在区间 $[0, T]$ 上，两个复指数信号的内积定义为：

$\langle e^{i2\pi f_1 t},\, e^{i2\pi f_2 t} \rangle = \int_0^T e^{i2\pi f_1 t} \cdot e^{-i2\pi f_2 t} \, dt = \int_0^T e^{i2\pi (f_1 - f_2)t} \, dt$

令 $\Delta f = f_1 - f_2$，计算积分：

$\langle \mathbf{e}_{f_1}, \mathbf{e}_{f_2} \rangle = \frac{e^{i2\pi \Delta f \, T} - 1}{i2\pi \Delta f}$

利用欧拉公式 $e^{i\theta} - 1 = e^{i\theta/2} \cdot 2i\sin(\theta/2)$，化简为：

$\langle \mathbf{e}_{f_1}, \mathbf{e}_{f_2} \rangle = T \cdot \frac{\sin(\pi \Delta f \, T)}{\pi \Delta f \, T} \cdot e^{i\pi \Delta f \, T}$

6.1.2.4.3.2. 第二步：正交条件

两个基向量正交，当且仅当内积为零：

$\sin(\pi \Delta f \, T) = 0$

这要求：

$\pi \Delta f \, T = k\pi, \quad k = \pm 1, \pm 2, \dots$

即：

$\Delta f = \frac{k}{T}$

6.1.2.4.3.3. 第三步：最小正交间隔

满足正交的最小非零频率间隔对应 $k = 1$：

$\boxed{\Delta f_{\min} = \frac{1}{T}}$

这意味着：在有限窗口 $[0, T]$ 内，只有当两个频率相差 $1/T$ 的整数倍时，它们对应的复指数基函数才是严格正交的。

6.1.2.4.3.4. 第四步：正交性为什么决定分辨率？

傅里叶分析是线性投影：把信号往各个频率基向量上投影，投影系数就是该频率的幅度。

• 若基向量正交（$\Delta f \geq 1/T$），不同频率的能量互不泄漏，可以被独立、唯一地提取。
• 若基向量不正交（$\Delta f < 1/T$），内积 $\frac{\sin(\pi \Delta f T)}{\pi \Delta f T}$ 接近 1，两个基向量近似平行。此时：
  • 10 Hz 的投影里混入了 10.05 Hz 的能量；
  • 10.05 Hz 的投影里也混入了 10 Hz 的能量；
  • 逆问题病态，存在无穷多组系数都能同样好地解释同一段数据。

6.1.2.4.3.5. 第五步：与 DFT 的对应

DFT 选取的频率点正是：

$f_k = \frac{k}{T}, \quad k = 0, 1, \dots, N-1$

这不是人为设定，而是这套基函数在 $[0, T]$ 上唯一能保持正交的频率网格。相邻谱线的间隔恰好是 $1/T$。

低于这个间隔的频率对，因为基向量不正交，DFT 无法将它们独立分辨——能量会泄漏、重叠、混为一谈。

6.1.2.4.4. 小结

有限采样空间触发的 sinc 涂抹，是连续和离散情况下的共同根源。连续 FT 让你看到涂抹后的连续模糊；DFT 则把这个模糊采样成了间距为 $1/T$ 的离散格子——两者都逃不开 $\Delta t \cdot \Delta f \sim 1$ 的约束。

	离散傅里叶变换 (DFT)	连续傅里叶变换 (如光学透镜)
有限窗口的数学作用	对无限信号加矩形窗，取 $N$ 点	对波前加有限孔径 $P(x,y)$
频域效应	与 sinc 函数（Dirichlet核）卷积后，只采样离散频点	与 sinc 函数（或艾里斑）卷积，频域仍是连续的
表现形式	频率间隔（Frequency Spacing） $\Delta f = 1/T$，谱线稀疏	频谱展宽 / 衍射极限 $\delta$ 函数变成 sinc/艾里斑，旁瓣重叠
分辨极限	频率分辨率受限，谱泄漏	瑞利判据：两点刚好可分辨的最小角间距

6.1.3. 从连续到离散：采样与周期化的对偶

当我们将连续信号数字化时，面临两个操作：时域采样 和 频域采样。这两个操作对应着对偶的周期化效应。

6.1.3.1. 时域采样导致频域周期化（DTFT）

对连续信号 $f(t)$ 以间隔 $T_s$ 采样，得到离散序列 $x[n] = f(nT_s)$。其频谱变为：

$X(\omega) = \sum_{n=-\infty}^{\infty} x[n] \, e^{-i\omega n}$

这是离散时间傅里叶变换（DTFT）。注意到 $\omega$ 是连续的，但具有 $2\pi$ 周期性——时域的离散化强制频域变成周期函数（看公式，$\omega + 2\pi$带入，其实和不加$2\pi$是一样的）。

6.1.3.2. 2. 频域采样导致时域周期化（DFT 的核心）

工程中最常用的 DFT，实际上是对 DTFT 的频域采样：

$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-i\frac{2\pi}{N}kn}, \quad k = 0, 1, \ldots, N-1$

DFT 隐含的数学假设：

• 时域周期：$x[n] = x[n + N]$（强制首尾相连，周期为 $N$）
• 频域周期：$X[k] = X[k + N]$

这与"截断"有本质区别！ 如果你有一段信号 $[1, 2, 3, 4]$，DFT 实际处理的是无限序列 $\ldots, 1,2,3,4,1,2,3,4,1,2,3,4, \ldots$，而不是 $\ldots,0,0,1,2,3,4,0,0,\ldots$。

6.1.4. 全家族对比总结

变换类型	时域特征	频域特征	数学假设	适用场景	边界/周期特性
连续傅里叶变换 (CFT)	连续，无限区间 $(-\infty, \infty)$	连续，无限区间	绝对可积	理论分析	无周期假设
有限区间 FT（截断）	连续，有限 $[a,b]$，区间外为 0	连续，无限	紧支集	瞬态信号分析	区间外严格零值
傅里叶级数 (FS)	连续，周期延拓	离散谱线 $\xi_n = n/T$	周期函数	周期振动、谐波分析	区间外周期重复
离散时间 FT (DTFT)	离散采样，无限长	连续，周期 $2\pi$	时域离散	数字滤波器设计	频域自然周期化
离散傅里叶变换 (DFT)	离散 + 周期延拓（$N$ 点）	离散 + 周期（$N$ 点）	双重周期	FFT 实现、数字信号处理	时域首尾相连，循环卷积

7. 目标解耦度量

7.1. M4度量

7.1.1. 定义

设 $I_0(x,y)$ 和 $I_d(x,y)$ 为两帧图像，其离散傅里叶变换（DFT）为：

$S_0(u,v) = \mathcal{F}\{I_0\}, \quad S_d(u,v) = \mathcal{F}\{I_d\}$

M4 定义为功率谱比值：

$M_4(u,v) = \frac{|S_0(u,v)|^2 - |S_d(u,v)|^2}{|S_0(u,v)|^2 + |S_d(u,v)|^2}$

7.1.2. 均匀背景不敏感

在做基于 M4 metric 的无信标波前重构时，我尝试在 dataset.py 里把图像背景（灰度值约 130）减掉，结果测试 RMS 从 0.16 涨到了 0.22。深入排查后发现：不是背景不该减，而是我的数据增强 pipeline 与“绝对灰度值”强耦合，减去背景后，RandomShiftWithPad0、RandomErasing 等操作几乎失效，导致训练正则化不足、泛化变差。

预处理方式	最佳 Test RMS
不减背景 (image / 4096.0)	~0.16
减去背景 130 ((image - 130) / 4096.0)	~0.22

从 log 里可以清晰看到，减去背景的模型在第 6 个 epoch 后就开始过拟合，test_rms 停止下降甚至反弹；而不减背景的模型收敛得更平稳，最终 RMS 明显更低。

这让我非常困惑：背景不是噪声吗？为什么去掉噪声反而变差？

7.1.2.1. 原因剖析

7.1.2.1.1. M4 Metric 本身对均匀背景不敏感

先回到物理本质。M4 metric 的定义是：

$M_4 = \frac{|S_o|^2 - |S_d|^2}{|S_o|^2 + |S_d|^2}$

其中 $S_o$、$S_d$ 分别是两幅图像的 FFT。
一个均匀常数背景 $B$ 的傅里叶变换只集中在零频 $(0,0)$，对非零频率的贡献为 0。因此，如果没有后续的数据增强，减不减背景，M4 输出在绝大多数频点上应该是完全一样的。

我直接用训练数据验证过：只做 ToTensor 和归一化时，两种情况算出的 M4 均值差异仅有 $5\times10^{-6}$，几乎可以忽略。

所以问题不在 M4 层，而在 M4 层之后的数据增强。

7.1.2.1.2. 数据增强与“绝对灰度值”强耦合

我的训练 transform 如下：

transform_train = transforms.Compose([
    transforms.ToTensor(),
    RandomShiftWithPad0(max_shift_h=5, max_shift_w=5, p=0.5),   # 零填充平移
    transforms.ColorJitter(brightness=0, contrast=1, saturation=0, hue=0),
    transforms.RandomErasing(p=0.3, scale=(0.02, 0.33), ratio=(0.3, 3.3), value=0),
])

这三个操作都严重依赖图像的绝对灰度值，而减去背景会彻底改变它们的行为：

① RandomShiftWithPad0 —— 平移后补 0

• 不减背景：图像背景约为 $130/4096 \approx 0.032$。平移后空缺处填充 0，会在图像边缘形成明显的“暗边”，引入强烈的频谱扰动。
• 减去背景：背景被拉到 0 附近。填充 0 等同于填充背景，暗边消失，频谱扰动几乎被抹平。

② RandomErasing(value=0) —— 用 0 擦除一块区域

• 不减背景：在亮背景上打一块“黑斑”，产生显著的局部频域变化。
• 减去背景：擦除区域和周围背景灰度几乎一致，擦了个寂寞，对网络来说相当于没做增强。

③ ColorJitter(contrast=1) —— 以均值为基准拉伸

• 不减背景：均值约为 0.032，对比度拉伸会把“信号+背景”一起放大，保持动态范围变化。
• 减去背景：均值接近 0，拉伸的幅度和效果都大打折扣。

一句话总结：减去背景后，数据增强被“阉割”了，训练样本多样性严重不足，模型学不到足够鲁棒的特征，测试性能自然崩塌。

7.1.2.1.3. 减去的 130 本身也不准确

更雪上加霜的是，我统计了训练集里 200 个样本的灰度分布：

统计量	数值
像素最小值	118.43 ~ 119.01（平均 118.79）
像素均值	122.58

这说明真实的暗电流/偏置（background）应该在 ~118.8 左右，而不是 130（130是背景的最大值）。强行减 130 会在图像上留下约 -11.2 的残余负偏置，导致：

• RandomErasing(value=0) 的“黑斑”变成了相对背景而言的亮斑；
• 图像整体出现负值，可能让 BatchNorm 或 ReLU 的分布发生微妙偏移。

7.1.2.1.4. 论文里说的 “subtract” 不是空域减背景

回头再看参考论文 Phase-diversity wave-front sensor for imaging systems，其中提到的 subtract 出现在噪声分析章节：

“the camera noise can be closely approximated by a decoupled spectrum that will cancel in the numerator. This spectrum can be experimentally determined and also subtracted as a separate term in the denominator.”

注意，这里说的是在 M4 metric 的分母里减去噪声功率谱（denominator subtraction），是在频域对功率谱做的操作，绝不是在空域直接减去一个固定灰度值 130。

7.1.2.2. 实验验证

为了量化增强失效的程度，我写了一个对比脚本：对同一张训练样本施加完全相同的平移（5,5）、对比度抖动（factor=1.5）、随机擦除，然后计算 M4。

# 模拟相同增强后计算 M4
for sub_bg in [False, True]:
    x = apply_transforms(img, sub_bg=sub_bg).unsqueeze(0)
    ...
    m4 = (power_So - power_Sd) / (power_So + power_Sd + 1e-8)
    print(f"sub={sub_bg}: M4 mean={m4.mean():.4f}")

结果：

预处理	M4 均值
不减背景	0.2059
减去背景 130	0.5478

两者 M4 的差异均值达到 -0.34，最大单点差异接近 2.0。 这证实了减去背景操作对数据增广的影响。

这并不能定量度量减去背景对数据增广的影响？为什么？

7.1.2.3. 结论

不是背景不能减，而是你的数据增强策略与“非零背景”共生。
在这个特定的 M4 + ResNet 框架下，保留约 130 的背景偏置，恰好让 RandomShiftWithPad0、RandomErasing 和 ColorJitter 发挥了最强的正则化作用。

深度学习里的很多“玄学”问题，往往根子都在预处理与增强的耦合上。这次踩坑让我意识到：

“去掉背景”这个操作在物理上很干净，但在工程实现里，必须考虑整条数据管道的上下文。

7.1.3. 循环平移不变特性

定理：对于离焦图像对 $I_0$（聚焦）和 $I_d$（离焦），经循环平移后，M4 指标保持不变。

7.1.3.1. DFT 的周期延拓本质

对于长度为 $N$ 的一维离散序列 $x[n]$，其中 $n = 0, 1, \dots, N-1$，其 DFT 定义为：

$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-i 2\pi k n / N}, \quad k = 0, 1, \dots, N-1$

关键洞察：这个公式本身就假设了 $x[n]$ 是周期为 $N$ 的周期信号，即：

$x[n + N] = x[n], \quad \forall n \in \mathbb{Z}$

7.1.3.2. 循环平移 = 周期信号的自然平移

在周期信号框架下，平移 $\Delta$ 定义为：

$x'[n] = x[(n - \Delta) \bmod N]$

这正是 torch.roll 所做的操作。在这个定义下：

平移定理（DFT 版本）： 如果 $x[n] \xrightarrow{\text{DFT}} X[k]$，那么：

$x[(n-\Delta) \bmod N] \xrightarrow{\text{DFT}} X[k] \cdot e^{-i 2\pi k \Delta / N}$

证明：

$\begin{aligned} \sum_{n=0}^{N-1} x[(n-\Delta) \bmod N] \cdot e^{-i 2\pi k n / N} &= \sum_{m=-\Delta}^{N-1-\Delta} x[m \bmod N] \cdot e^{-i 2\pi k (m+\Delta) / N} \quad (\text{令 } m = n-\Delta) \\ &= e^{-i 2\pi k \Delta / N} \sum_{m=0}^{N-1} x[m] \cdot e^{-i 2\pi k m / N} \quad (\text{周期求和不变}) \\ &= X[k] \cdot e^{-i 2\pi k \Delta / N} \end{aligned}$

7.1.3.3. 功率谱不变性

取模平方：

$|X'[k]|^2 = |X[k] \cdot e^{-i 2\pi k \Delta / N}|^2 = |X[k]|^2 \cdot |e^{-i 2\pi k \Delta / N}|^2 = |X[k]|^2$

关键：因为 $e^{-i\theta}$ 的模为 1，功率谱（能量）完全守恒。

7.1.3.4. 对比：零填充破坏了周期假设

零填充平移在时域的数学表达：

$x_{zp}[n] = \begin{cases} x[(n-\Delta) \bmod N] & \text{if } 0 \leq n < N \\ 0 & \text{otherwise} \end{cases}$

但这不是 DFT 框架下的操作！零填充实际上等价于：

7.1.3.4.1. 时域：与矩形窗的乘积

$x_{zp}[n] = x_{shifted}[n] \cdot \Pi[n]$

其中 $\Pi[n]$ 是矩形窗（$n \in [0,N-1]$ 时为 1，否则为 0）。

7.1.3.4.2. 频域：与 Sinc 函数的卷积

根据卷积定理：

$\text{DFT}\{x_{zp}\} = \text{DFT}\{x_{shifted}\} \circledast \text{DFT}\{\Pi\} = (X[k] \cdot e^{-i\phi_k}) \circledast \text{sinc}(k)$

这里 $\text{sinc}(k) = \frac{\sin(\pi k)}{\pi k}$（Dirichlet 核的周期形式）。

结果：功率谱被 smear（涂抹）：

$|X_{zp}[k]|^2 \neq |X[k]|^2$

8. 仿真数据优化

8.1. 让波长"可学习"

在光学模型里，我们经常希望波长能自动调整——比如让神经网络自己优化出最佳波长。 但波长有个硬约束：必须是正数，不能为 0，更不能为负。

下面这段代码展示了一个非常巧妙的做法：

raw_wl = float(init_wavelength_m)
if raw_wl > 100e-9:
    raw_wl = np.log(np.exp(raw_wl) - 1.0)
self.raw_wavelength_m = nn.Parameter(torch.tensor(raw_wl))

@property
def wavelength_m(self):
    return F.softplus(self.raw_wavelength_m) + 1e-12

8.1.1. 编码：把波长变成"暗号"存起来

假设你的初始波长是 1.0，先做了一道"反向计算"：

raw_wl = np.log(np.exp(raw_wl) - 1.0)

这步就是在算暗号：应该存什么数字，才能让后面解码出来恰好等于原始波长？

打个比方：保险箱规定"取出来的时候自动加 1"。你想最终拿到 500，就不能存 500，得存 499。这步就是在算"499"。

8.1.2. 解码：把暗号还原成合法的波长

@property
def wavelength_m(self):
    return F.softplus(self.raw_wavelength_m) + 1e-12

softplus 的作用很简单：无论你输入什么数字，输出永远是正数。

内部暗号	softplus 输出
10	10.0
0	0.69
-5	0.007
-100	约 0

因为内部暗号可以任意变化（训练时梯度下降会不断更新它），但经过 softplus 一过滤，最终波长永远被锁在正数区间。

最后的 + 1e-12 只是再加个保险，确保波长不会变成 0，避免后面做除法时崩溃。

8.1.3. 为什么要绕这个弯？

8.1.3.1. 物理结果要合法

波长不能为负，不能为 0。softplus 在数学上就是一个"单向阀门"：进去的是任意实数，出来的永远是正数。

8.1.3.2. 训练过程要自由

nn.Parameter 意味着这个值会被自动优化。优化过程中参数可能变正、变负、变小、变大——内部暗号随便变，没关系。因为解码器（softplus）会把它"压"成正数。

编码是为了让初始值不被扭曲；解码是为了让输出永远合法。
内部存一个可以随便变的暗号，对外暴露一个永远为正的真实波长——这就是"可学习参数 + 物理约束"的经典解法。

8.2. 平移不变的点扩散函数（PSF）监督损失：从 FFT 幅度谱出发

背景：在真实光学系统中采集的在焦/离焦点扩散函数（PSF）往往因为光轴对准误差、探测器安装公差等原因，并不严格位于图像中心。然而，物理仿真得到的 PSF 默认是以光轴为中心生成的。如果直接在像素空间做 L2/L1 监督，一个微小的平移就会让损失函数爆炸，优化器被迫去"凑"位置而不是真正学习像差和离焦量。

8.2.1. 核心原理：傅里叶变换的平移不变性

8.2.1.1. 平移定理的推导

我们从二维离散傅里叶变换（2D-DFT）的定义出发，一步一步推导。

对一幅尺寸为 $H \times W$ 的离散图像 $I(x,y)$，其 DFT 定义为：

$\mathcal{F}\{I\}(u,v) = \sum_{x=0}^{W-1} \sum_{y=0}^{H-1} I(x,y) \; e^{-j 2\pi \left(\frac{ux}{W} + \frac{vy}{H}\right)}$

其中：

• $(x,y)$ 是空间域像素坐标；
• $(u,v)$ 是频率域坐标；
• $j$ 是虚数单位，满足 $j^2 = -1$。

现在，图像发生空间平移 $(x_0, y_0)$，得到新图像：

$I'(x,y) = I(x-x_0,\; y-y_0)$

将 $I'$ 代入 DFT 定义：

$\mathcal{F}\{I'\}(u,v) = \sum_{x=0}^{W-1} \sum_{y=0}^{H-1} I(x-x_0,\; y-y_0) \; e^{-j 2\pi \left(\frac{ux}{W} + \frac{vy}{H}\right)}$

做变量代换：令 $m = x - x_0$，$n = y - y_0$，则 $x = m + x_0$，$y = n + y_0$。当 $x$ 遍历 $0 \dots W-1$ 时，$m$ 遍历 $-x_0 \dots W-1-x_0$。由于 DFT 隐含周期性边界条件 $I(x+W, y) = I(x, y)$，求和区间的平移不会改变求和结果。因此我们可以将求和变量改回标准区间：

$\mathcal{F}\{I'\}(u,v) = \sum_{m=0}^{W-1} \sum_{n=0}^{H-1} I(m,n) \; e^{-j 2\pi \left(\frac{u(m+x_0)}{W} + \frac{v(n+y_0)}{H}\right)}$

将指数项拆开。利用指数性质 $e^{a+b} = e^a \cdot e^b$：

$\mathcal{F}\{I'\}(u,v) = \sum_{m=0}^{W-1} \sum_{n=0}^{H-1} I(m,n) \; e^{-j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)} \cdot e^{-j 2\pi \left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}$

提取公因子。观察第二项指数 $e^{-j 2\pi \left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}$，它只依赖于频率 $(u,v)$ 和平移量 $(x_0,y_0)$，与求和变量 $(m,n)$ 完全无关。因此它可以被提到双重求和号外面：

$\mathcal{F}\{I'\}(u,v) = \underbrace{\left[\sum_{m=0}^{W-1} \sum_{n=0}^{H-1} I(m,n) \; e^{-j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)}\right]}_{\text{这正是 } \mathcal{F}\{I\}(u,v)} \cdot e^{-j 2\pi \left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}$

于是得到平移定理（Shift Theorem）：

$\boxed{\mathcal{F}\{I'\}(u,v) = \mathcal{F}\{I\}(u,v) \;\cdot\; e^{-j 2\pi \left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}}$

物理意义：空间域的平移 $(x_0, y_0)$，在频率域体现为每个频率分量 $(u,v)$ 都乘以了一个单位复数 $e^{-j 2\pi (\cdots)}$。这个复数的模为 1，因此不改变各频率分量的"强度"；但它带有相位信息，因此会改变各频率分量的"对齐方式"。

8.2.1.2. 为什么幅度不变、相位改变？

我们把原图像在频率 $(u,v)$ 处的傅里叶系数写成极坐标形式（幅度-相位形式）：

$\mathcal{F}\{I\}(u,v) = A(u,v) \cdot e^{j\phi(u,v)}$

其中：

• $A(u,v) = \big|\mathcal{F}\{I\}(u,v)\big| \ge 0$ 是幅度（Amplitude / Magnitude）；
• $\phi(u,v) = \arg\big(\mathcal{F}\{I\}(u,v)\big)$ 是相位（Phase）。

将平移定理代入：

$\mathcal{F}\{I'\}(u,v) = \left[A(u,v) \cdot e^{j\phi(u,v)}\right] \cdot e^{-j 2\pi \left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}$

合并指数项：

$\mathcal{F}\{I'\}(u,v) = A(u,v) \cdot e^{j\left[\phi(u,v) \;-\; 2\pi\left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)\right]}$

8.2.1.2.1. 幅度分析

对上式两边取复数模（Magnitude）。利用复数模的乘法性质 $|z_1 \cdot z_2| = |z_1| \cdot |z_2|$：

$\big|\mathcal{F}\{I'\}(u,v)\big| = \Big|A(u,v) \cdot e^{j\left[\phi(u,v) - 2\pi(\cdots)\right]}\Big|$

$= |A(u,v)| \cdot \underbrace{\left|e^{j\left[\phi(u,v) - 2\pi(\cdots)\right]}\right|}_{\text{（关键一步）}}$

这里用到复分析中的基本结论：对任意实数 $\theta$，都有

$|e^{j\theta}| = \sqrt{\cos^2\theta + \sin^2\theta} = 1$

因此，无论括号里的相位是什么，单位复数的模恒为 1：

$\left|e^{j\left[\phi(u,v) - 2\pi(\cdots)\right]}\right| = 1$

代回上式：

$\big|\mathcal{F}\{I'\}(u,v)\big| = |A(u,v)| \cdot 1 = A(u,v)$

而 $A(u,v) = \big|\mathcal{F}\{I\}(u,v)\big|$，所以：

$\boxed{\big|\mathcal{F}\{I'\}(u,v)\big| = \big|\mathcal{F}\{I\}(u,v)\big|}$

结论：图像在空间域发生任意平移 $(x_0, y_0)$ 后，其傅里叶变换的幅度谱（Magnitude Spectrum）在每个频率上都严格保持不变。

8.2.1.2.2. 相位分析

合并指数后的表达式直接读取相位：

$\arg\Big(\mathcal{F}\{I'\}(u,v)\Big) = \phi(u,v) \;-\; 2\pi\left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)$

即：

$\boxed{\arg\Big(\mathcal{F}\{I'\}\Big) = \arg\Big(\mathcal{F}\{I\}\Big) \;-\; 2\pi\left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)}$

结论：

• 空间平移确实改变了相位谱；
• 在每个频率 $(u,v)$ 上，相位改变量 $\Delta\phi = -2\pi\left(\frac{ux_0}{W} + \frac{vy_0}{H}\right)$ 与频率坐标 $(u,v)$ 成线性关系；
• 相位改变量与图像内容 $I(x,y)$ 本身无关，仅由平移量 $(x_0, y_0)$ 决定。

$\text{平移后} \quad \Rightarrow \quad \begin{cases} \text{幅度谱：} & \big|\mathcal{F}\{I'\}\big| = \big|\mathcal{F}\{I\}\big| & \text{（完全不变）}\\[8pt] \text{相位谱：} & \arg(\mathcal{F}\{I'\}) = \arg(\mathcal{F}\{I\}) - 2\pi\left(\frac{ux_0}{W} + \frac{vy_0}{H}\right) & \text{（线性偏移）} \end{cases}$

这就是平移不变性的严格数学来源。也正因如此，当我们用 FFT 幅度谱作为监督目标时，就完全不用担心真实采集的 PSF 与仿真 PSF 之间有几像素的中心偏移——因为偏移只影响相位，而我们根本不看相位。

• 幅度谱：描述图像中各频率分量的"强度"——即图像中有没有某种粗细的纹理、边缘、光斑形状。它与这些纹理出现在图像的哪个位置无关。
• 相位谱：描述各频率分量的"位置对齐信息"。一旦图像平移，所有频率的相位都会按线性规律发生偏转。

因此，如果我们只拿幅度谱做对比，就天然忽略了"图像整体偏移了多少"这个问题，只关心"图像里的形状对不对"。

8.2.2. 损失函数设计：FFT Log-Magnitude Loss

知道了"幅度谱平移不变"这个性质，接下来的问题是：如何把幅度谱转化为一个鲁棒的、可优化的损失函数？

直接对 $|FFT|$ 做 L1 会遇到一个工程问题：PSF 的中心亮斑幅度极高，外围像差细节的幅度极低。如果不做处理，损失函数会被中心亮斑完全主导，外围的高频像差信息被淹没。

以下是我们在实际代码中采用的四步法：

8.2.2.1. Step 1：背景扣除（Background Subtraction）

真实探测器采集的图像通常带有暗电流或环境杂散光基底。这个恒定背景会在 FFT 的零频（DC 分量）处产生一个巨大的尖峰，掩盖其他频率的差异。

real = real - real.amin(dim=(-2, -1), keepdim=True)
sim  = sim  - sim.amin(dim=(-2, -1), keepdim=True)

8.2.2.2. Step 2：峰值归一化（Peak Normalization）

仿真 PSF 的绝对强度是物理量（W/m²），而真实图像是相机灰度级（ADU）。两者尺度完全不同，必须归一化到同一量级：

real = real / (real.amax(dim=(-2, -1), keepdim=True) + eps)
sim  = sim  / (sim.amax(dim=(-2, -1), keepdim=True) + eps)

这样做保留了 PSF 的相对能量分布（在焦与离焦的相对亮度、光斑能量集中程度），而去除了探测器增益差异的影响。

8.2.2.3. Step 3：FFT 与 Log 压缩

mag_real = torch.log1p(torch.abs(torch.fft.fft2(real, dim=(-2, -1))))
mag_sim  = torch.log1p(torch.abs(torch.fft.fft2(sim,  dim=(-2, -1))))

• torch.fft.fft2：计算二维傅里叶变换。
• torch.abs：取幅度谱，丢弃受平移影响的相位。
• torch.log1p(x)：即 $\ln(1+x)$，对高动态范围进行压缩。中心亮斑的 $10^4$ 倍差异，在 log 域只变成几倍差异，使得外围中低频像差也能获得合理的梯度权重。

8.2.2.4. Step 4：Smooth L1 回归

loss = F.smooth_l1_loss(mag_sim, mag_real)

相比于 MSE，Smooth L1（Huber Loss）对离群值更鲁棒，在频域某些异常频率点上不会导致梯度爆炸。

8.2.3. 注意事项与进阶方向

8.2.3.1. 当前方案的局限

• 旋转与缩放不变性：FFT 幅度谱对平移不变，但对旋转和缩放不是不变的。如果你的光路还存在旋转对准误差或像素尺度误差，需要额外处理（例如先把仿真图插值到估计的像素尺度，或在极坐标频域做监督）。
• 直流分量：虽然做了背景扣除，但如果图像存在大面积非均匀光照（不是恒定背景），低频区域仍可能受影响。

8.2.3.2. 进阶改进

1. 频域加权： 像差引起的光学传递函数（OTF）退化主要集中在中低频。可以构造一个二维频率权重矩阵 $W(u,v)$，给低频更高权重：

   loss = (W * (mag_sim - mag_real)).abs().mean()

2. 互功率谱对齐（Phase Correlation）： 如果平移量极大（接近图像边界），仅靠幅度谱可能不够。可以先计算互功率谱的逆 FFT，得到平移的置信度图，再据此做软对齐（Soft-Argmax）或作为辅助监督。

3. 加窗处理： FFT 假设信号周期性。如果 PSF 图像边缘截断明显（例如能量没有衰减到接近零），可以加 Hann 窗或 Hamming 窗后再做 FFT，减少频谱泄漏。

8.3. 背景噪声导致的"频域Loss 下降但视觉变差"问题分析

核心结论：真实采集图像存在非零背景噪声底板（min-max 后均值约 0.019），而仿真 PSF 背景严格为 0。当前 fft_psf_loss 对频域轮廓敏感，但对背景噪声底板几乎无感知，导致模型可以通过高斯模糊"骗过"频域 loss，却在空间域视觉上严重失真。

8.3.1. 问题现象

观察到以下矛盾现象：

实验	训练 epoch	Test Loss	视觉质量（人眼观察）
no_截断（仅高斯模糊）	7	0.136 ⬇️	更差 ❌
no_高斯（仅截断）	1	0.161	更好 ✅
完整 detector（高斯+截断）	1	0.218	一般

矛盾点：no_截断 的 loss 最低，但视觉上 PSF 与真实数据差异反而更大。

8.3.2. 核心诊断：背景噪声底板缺失

8.3.2.1. 真实采集数据的背景特征

原始数据文件（.npz）中的 image 在送入网络前，经过 dataset.py 的预处理：

image = image - image.amin(dim=(-2, -1), keepdim=True)   # 逐样本扣最小值
image = image / (image.amax(dim=(-2, -1), keepdim=True)) # 逐样本除最大值

预处理后，真实采集图像（focus_real / defocus_real）的像素分布如下：

统计指标	focus_real	defocus_real
严格为 0 的像素数	0~2 个（占比 ~0%）	0~2 个（占比 ~0%）
噪声带 (0.001~0.05)	97.1%	97.1%
中间亮度 (0.05~0.5)	2.8%	2.8%
高亮区 (≥0.5)	0.1%	0.1%
P1 分位数	0.0093	0.0106
P10 分位数	0.0149	0.0152
底部 20% 均值	0.0150	0.0152

关键发现：

• 真实图像预处理后，几乎所有像素（97%）都分布在一个非零的噪声带内，形成一层均匀的"灰色噪声底板"
• 即使是最暗的 1% 像素，亮度也在 0.009~0.011 之间
• 图像中不存在大片严格为 0 的纯黑区域

8.3.2.2. 仿真数据的背景特征

三个实验的仿真 PSF（focus_train / defocus_train）经过相同的 min-max 归一化后：

实验	focus_train 严格为 0	defocus_train 严格为 0	噪声带 (0.001~0.05)
完整 detector	91.7%	92.1%	5.4%
no_截断	94.2%	94.2%	3.8%
no_高斯	91.7%	91.3%	6.8%

关键发现：

• 仿真图像 90%+ 像素严格为 0，背景是纯黑的
• 只有约 4~7% 的像素分布在噪声带，这些通常是高斯模糊或截断产生的过渡像素

8.3.2.3. 背景区（角落）精细对比

选取图像四个角落各 50×50 区域（远离光斑峰值，视为纯背景）：

指标	真实 focus_real 角落	完整 detector 仿真角落
均值	0.0193	0.0000
标准差	0.0039	0.0000
最小值	0.0000	0.0000
最大值	0.0392	0.0000

关键发现：

• 真实背景区有明显的随机波动（std=0.004），这是读出噪声和热噪声的特征
• 仿真背景区完全平坦（std=0），没有任何噪声

8.3.2.4. 误差按区域分解

以 focus 为例，分析 train - real 的差异在不同区域的分布：

区域定义	完整 detector MAE	no_截断 MAE	no_高斯 MAE
背景区 (real < 0.05)	0.0206	0.0205	0.0201
中间区 (0.05 ≤ real < 0.5)	0.1150	0.1131	0.1170
峰值区 (real ≥ 0.5)	0.5475	0.5493	0.5602
整体 MAE	0.0237	0.0235	0.0233

关键发现：

• 三个实验在背景区的误差方向一致且显著为负（仿真比真实暗约 0.017~0.018），说明所有实验都存在"背景缺失"问题
• no_高斯 在背景区的 MAE 最小（0.0201），_std 也最小（0.0148），说明其背景最"干净"但也最偏离真实噪声
• 空间域误差主要由峰值区结构差异主导（MAE _{0.55），但背景区的系统性偏差（}0.02）在视觉上非常突出

8.3.3. 根本原因分析

8.3.3.1. 为什么 fft_psf_loss 会"被骗"？

当前损失函数定义：

def fft_psf_loss(sim_psf, real_psf, eps=1e-8):
    real_fft = fft.fft2(real_psf, dim=(-2, -1))
    sim_fft = fft.fft2(sim_psf, dim=(-2, -1))
    real_mag = torch.log1p(torch.abs(real_fft))
    sim_mag = torch.log1p(torch.abs(sim_fft))
    return F.smooth_l1_loss(sim_mag, sim_mag)

8.3.3.1.1. 实验验证：噪声对 FFT loss 的真实影响

真实数据类型	vs Sim_Clean (清晰)	vs Sim_Blur (模糊)
Clean 艾里斑	0.000	0.006
Clean + 均匀常数背景 (0.015)	0.000	0.006
Clean + 随机背景噪声 (mean=0.015, std=0.004)	0.377	0.383

关键发现：

1. 均匀常数背景经过 min-max 归一化后，对 FFT loss 零影响
2. 随机波动的背景噪声（std=0.004）让同一个艾里斑的 FFT loss 暴涨到 0.377
3. 噪声并没有让"模糊"变得更有优势——无论目标是否带噪声，清晰 PSF 始终比模糊 PSF 的 loss 低约 0.006

8.3.3.1.2. 真正的问题机制：噪声→不规则频谱→信号扭曲

第一步：真实频谱被噪声"污染"

真实数据的频谱 = 信号频谱（艾里斑） + 噪声频谱（宽谱随机波纹）

随机噪声虽然空间域 std 只有 0.004，但它分布在 400×400 = 16 万个像素上。FFT 后，这些扰动不会消失，而是给每一个频率分量都叠加了一个随机偏移。

第二步：log1p 没有"消除"噪声，而是把它变成了一层"不规则的频谱波纹"

噪声在每个频率上的能量不高（被 16 万像素分摊），log1p 把它压缩成一个小的随机偏移量。但这个偏移量分布在成千上万个频率分量上，累积起来形成了一个不规则的频谱基线。

第三步：SmoothL1Loss 强迫模型扭曲信号去匹配噪声形状

当真实频谱的某个频率因为噪声而随机偏高时，模型为了降低总 loss（数千个频率分量的累积），只能通过扭曲信号本身来在该频率上"制造"类似的能量：

• 调整波长（如 no_截断 从 530 nm 推到 306 nm，偏离 42%）
• 调整 defocus（从 -0.265 推到 -0.142）
• 加高斯模糊（改变频谱衰减斜率）

第四步：空间域灾难

频域上"凑出相似衰减曲线"的代价是：空间域的 PSF 变成了一个物理上不合理、结构上严重失真的光斑。艾里环被抹平，光斑尺寸被强行改变，悬浮在纯黑背景上——这就是"loss 降、视觉崩"的根源。

随机背景噪声给真实频谱添加了"无规律的波纹"。FFT loss 强迫模型扭曲信号（波长、离焦、模糊）去贴合这些波纹，从而降低了 loss 数值，但空间域的 PSF 结构被严重破坏。

8.3.3.2. 为什么 no_截断 视觉上更差？

频域匹配 ≠ 空间域匹配

维度	no_截断（7 epoch）	no_高斯（1 epoch）	真实数据
FWHM（半高宽）	9 px（过窄）	24 px（接近）	19 px
第一暗环深度	0.0496	0.0492	0.0196
第一暗环位置	25 px	14 px	32 px
背景	纯黑（0）	纯黑（0）	灰色底板（~0.02）
艾里环清晰度	被高斯模糊严重抹平	较清晰	清晰但叠加噪声

关键机制：

• no_截断 训练 7 epoch 后 loss 从 0.223 降到 0.136，主要不是靠高斯模糊（实验证明模糊反而增加 loss），而是靠物理参数被严重扭曲：
  • 波长：530 nm → 306 nm（偏离 42%，非物理）
  • defocus：-0.265 μm → -0.142 μm
• 模型通过扭曲物理参数来"拟合"被噪声污染的不规则频谱，在空间域产生了一个过窄、过平滑、无艾里环结构的光斑
• 真实数据虽然有噪声，但艾里环结构仍然可辨；no_截断 的仿真却把结构完全抹掉了
• no_高斯 虽然只有 1 epoch，但没有高斯模糊，且物理参数偏离小（波长仅到 518 nm），保留了艾里环结构，虽然背景仍是 0，但至少"轮廓是对的"

8.3.3.3. min-max 归一化如何放大了问题？

当前数据预处理流程：

原始图像（有背景噪声）
    ↓ 逐样本减去最小值
最小值被拉到 0，但其他背景像素仍 > 0
    ↓ 逐样本除以最大值
噪声底板被归一化到约 0.01~0.02 的范围

问题：

• 真实相机的背景是系统性的（暗电流、读出噪声、固定环境光），应该全局扣除
• 逐样本 amin 导致每个样本的背景基准不同，模型难以学到统一的背景
• 仿真 PSF 背景严格为 0，min-max 后仍然是 0，与真实数据的非零底板形成强烈视觉反差

8.3.4. 解决方案

8.3.4.1. 方案一：修改数据预处理（dataset.py）

核心思想：用全局/鲁棒背景估计替代逐样本 amin，让真实数据的背景更接近 0。

import torch
import torch.nn.functional as F
import numpy as np
from torch.utils.data import Dataset

class CustomDataset(Dataset):
    def __init__(self, sample_txt, transform=None, global_bg=None):
        self.samples = np.loadtxt(sample_txt, dtype=np.str_)
        self.transform = transform
        self.global_bg = global_bg  # 若提供全局背景值（原始尺度）

    def __len__(self):
        return len(self.samples)

    def __getitem__(self, idx):
        while True:
            try:
                data = np.load(self.samples[idx])
                image = data["image"]  # shape 假设为 [H, W, C] 或 [C, H, W]
                zernike = data["zernike"]

                if self.transform:
                    image = self.transform(image)

                # ========== 背景扣除策略（三选一）==========

                # 策略 A：全局背景扣除（推荐）
                # 预先计算所有样本角落区域的中位数作为全局背景
                if self.global_bg is not None:
                    image = image - self.global_bg
                    image = F.relu(image)  # 防止负值

                # 策略 B：逐样本角落中位数扣除（无需预计算）
                else:
                    h, w = image.shape[-2], image.shape[-1]
                    corner_size = min(h, w) // 8  # 取 1/8 边长作为角落
                    corners = torch.cat([
                        image[..., :corner_size, :corner_size].flatten(),
                        image[..., :corner_size, -corner_size:].flatten(),
                        image[..., -corner_size:, :corner_size].flatten(),
                        image[..., -corner_size:, -corner_size:].flatten()
                    ])
                    bg = corners.median()  # 中位数比最小值更鲁棒
                    image = image - bg
                    image = F.relu(image)

                # 归一化：只除以最大值（背景已扣除，min 接近 0）
                image = image / (image.amax(dim=(-2, -1), keepdim=True) + 1e-8)

                return image, zernike
            except BaseException as e:
                print(e)
                import time
                time.sleep(1)

为什么用角落中位数？

• 光斑通常位于图像中心，四个角落远离光斑，可视为纯背景
• 中位数比最小值更鲁棒，不受单个坏点/热像素影响
• 全局背景扣除保证了样本间背景基准一致

8.3.4.2. 方案二：给 Detector 加背景噪声模拟

在 generative_M4_model.py 中修改 DetectorResponseLayer，添加可学习的背景偏移和噪声。

8.3.4.3. 方案三：组合损失函数（空间域 + 频域）

当前仅用 fft_psf_loss，对背景不敏感。建议增加空间域损失项：

from torch import fft
import torch.nn.functional as F

def combined_psf_loss(sim_psf, real_psf, alpha=0.5, eps=1e-8):
    """
    组合损失：频域 loss（平移不变）+ 空间域 loss（背景敏感）

    Args:
        sim_psf:  [B, 2, H, W] 仿真 PSF
        real_psf: [B, 2, H, W] 真实 PSF
        alpha:    频域 loss 权重（0~1），建议从 0.7 开始调试
    """
    # ===== 频域 loss（保留平移不变性）=====
    real_fft = fft.fft2(real_psf, dim=(-2, -1))
    sim_fft = fft.fft2(sim_psf, dim=(-2, -1))
    real_mag = torch.log1p(torch.abs(real_fft))
    sim_mag = torch.log1p(torch.abs(sim_fft))
    fft_loss = F.smooth_l1_loss(sim_mag, sim_mag)

    # ===== 空间域 loss（对背景噪声敏感）=====
    # Smooth L1 比 MSE 更鲁棒，对异常值（如过曝像素）不敏感
    spatial_loss = F.smooth_l1_loss(sim_psf, real_psf)

    # ===== 可选：背景区域加权 loss =====
    # 让模型更关注背景区的匹配（防止背景被忽略）
    # 定义背景掩码：真实图像中亮度 < 0.05 的区域
    bg_mask = (real_psf < 0.05).float()
    bg_loss = F.smooth_l1_loss(sim_psf * bg_mask, real_psf * bg_mask)

    # 组合
    total_loss = alpha * fft_loss + (1 - alpha) * spatial_loss + 0.1 * bg_loss

    return total_loss, {
        'fft_loss': fft_loss.item(),
        'spatial_loss': spatial_loss.item(),
        'bg_loss': bg_loss.item()
    }

参数调试建议：

阶段	alpha	说明
初期	0.7	以频域为主，保证整体轮廓正确
中期	0.5	空间域与频域同等重要
后期	0.3	以空间域为主，精细化背景匹配

8.3.5. 验证实验

8.3.5.1. 方案1验证试验，中位数背景扣除效果（没有重新训练，仅测试）

以下为实际测试数据：将修改后的 dataset.py（角落中位数背景扣除）分别应用到三个实验配置中，不重新训练模型，直接运行 test_optical_param.py 对比新旧 dataset 的 raw data。

8.3.5.1.1. 实验设置

实验	训练模型	Dataset 预处理	说明
旧-完整 detector	旧 checkpoint（1 epoch）	逐样本 amin	基线
旧-无截断	旧 checkpoint（7 epoch）	逐样本 amin	基线
旧-无高斯	旧 checkpoint（1 epoch）	逐样本 amin	基线
新-完整 detector	旧 checkpoint（1 epoch）	角落中位数扣除	验证
新-无截断	旧 checkpoint（7 epoch）	角落中位数扣除	验证
新-无高斯	旧 checkpoint（1 epoch）	角落中位数扣除	验证

8.3.5.1.2. 真实数据背景分布变化

旧 dataset（逐样本 amin）：

• 底部 10% 像素 = 0.0149
• 底部 20% 均值 = 0.0150
• 严格为 0 的像素 = ~0%
• 角落背景均值 = 0.0093

新 dataset（角落中位数）：

• P1/P5/P10 = 0.000000
• 严格为 0 的像素 = 49.3%
• 微弱残余带 (0.001~0.01) = 44.5%
• 角落背景均值 = 0.000884（旧值的 1/10）

8.3.5.1.3. 按区域分解的改善来源

区域	旧 MAE	新 MAE	变化	解读
背景区 (real < 0.05)	~0.010	~0.002	↓ 80% ✅	中位数扣除直接消除了背景系统偏差
中间区 (0.05~0.5)	~0.120	~0.125	基本持平	中位数扣除不影响信号结构
峰值区 (real ≥ 0.5)	~0.500	~0.500	基本持平	峰值结构不受预处理影响

8.3.5.1.4. 关键发现

1. 改善 100% 来自背景区，峰值区不变

   • 中位数扣除只影响背景，不影响艾里斑结构
   • 峰值区 MAE 不变是因为模型未重训（旧 checkpoint 的物理参数没变）

2. 无高斯在新 dataset 下仍然是最佳配置

   • 新-无高斯 Focus MAE = 0.00328（三实验中最低）
   • 再次验证：高斯模糊是不必要的

3. 仿真数据背景仍为 0，真实数据仍有微弱残余

   • 新 dataset 下仿真背景 = 0.000000
   • 真实背景 = 0.000884，仍有约 44.5% 像素在 0.001~0.01
   • 这是因为角落中位数可能略低于某些暗区像素，扣除后仍有少量残余

4. 如果重新训练模型，峰值区可能也会改善

   • 旧模型是在"高背景"数据上训练的，参数是为了匹配旧数据优化的
   • 新数据背景更干净，如果重训，模型可能学到更准确的物理参数

8.3.5.2. 方案1验证试验，中位数背景扣除效果（重新训练1 epoch）

实验设置：三个配置（完整 detector / 无截断 / 无高斯）均使用新的 dataset.py（角落中位数背景扣除）重新训练 1 个 epoch，然后运行测试。

8.3.5.2.1. 训练日志对比

Epoch 0 Loss

实验	旧 Dataset Train	旧 Dataset Test	新 Dataset Train	新 Dataset Test	Train 降幅	Test 降幅
完整 detector	0.2458	0.2182	0.2032	0.1785	-17.3%	-18.2%
无截断	0.2512	0.2234	0.2083	0.1834	-17.1%	-17.9%
无高斯	0.1712	0.1611	0.1415	0.1340	-17.3%	-16.8%

关键发现：

• 新 dataset 下，所有配置的初始 loss 都降低了约 17%，说明中位数背景扣除显著降低了数据与模型的初始失配
• 无高斯仍然是 loss 最低的（Test=0.1340），且与其他两个配置的差距保持（完整/无截断约 0.18）

8.3.5.2.2. 物理参数收敛对比

训练后模型参数

参数	旧-完整	新-完整	旧-无截断	新-无截断	旧-无高斯	新-无高斯
波长 λ	464 nm	464 nm	403 nm ❌	466 nm ✅	518 nm	535 nm ✅
defocus	-0.184 μm	-0.185 μm	-0.154 μm	-0.187 μm	-0.290 μm	-0.301 μm
focus_defocus	-0.009 μm	-0.011 μm	-0.038 μm	-0.016 μm	-0.062 μm	-0.066 μm
detector σ	1.13	1.13	0.60	1.13	1.21	1.21
detector saturation	0.92	0.92	1.00	1.00	1.04	1.03

关键发现：

1. 无高斯的波长最接近真实值：534.8 nm，仅偏离初始值 530 nm 约 0.9%。

2. 无截断的 σ 回归正常：旧 dataset 下 σ 被迫降到 0.60 来减弱高斯模糊，新 dataset 下回归至 1.13（接近初始 1.2）。

3. defocus 系数趋于一致：完整/无截断在新 dataset 下 defocus 都收敛到约 -0.186 μm（接近初始 -0.265 μm 的合理范围），而旧无截断是 -0.154 μm。

8.3.5.2.3. 高斯模糊与波长的耦合理论分析

耦合机制:

高斯模糊（卷积）和波长变化在频域上是耦合的，两者可以相互补偿。

频域效应对比：

操作	空间域	频域效应	对高频的影响
波长 λ↓	艾里斑变小	频谱展宽	增加
高斯模糊 σ↑	光斑展宽	频谱收缩	减少

关键洞察：

• 波长减小 → 频谱高频能量增加
• 高斯模糊增强 → 频谱高频能量减少
• 两者产生相反的频域效应

实验验证：不同 (λ, σ) 组合是否能产生相似频谱?

以参考 (λ=530nm, σ=0) 为基准，计算各组合的 fft_psf_loss：

(λ, σ)	vs 参考的 fft_loss	解读
(530, 0)	0.0000	参考
(464, 0)	0.0108	仅改变波长
(464, 0.6)	0.0057 ⭐	λ↓ + σ↑ 组合，loss 比单纯 λ↓ 更低
(464, 1.2)	0.0441	λ↓ + σ↑↑，过度补偿
(530, 1.2)	0.0436	仅加高斯模糊

关键发现：

• λ=464nm + σ=0.6 的 loss（0.0057）比 λ=464nm + σ=0（0.0108）更低
• 这说明：高斯模糊确实可以补偿波长变化带来的频谱差异
• 模型可以通过 (λ↓, σ↑) 的组合来逼近目标频谱，而无需学到真正的物理参数

对 detector 设计的启示

如果高斯模糊和波长是耦合的，同时保留两者会导致参数解空间退化，模型学不到真正的物理参数。

建议：

1. 去掉高斯模糊（已验证无高斯配置在所有指标上最优）
2. 固定波长（如果激光波长已知为 530nm，直接锁死 raw_wavelength_m 的梯度）
3. 如果必须保留高斯模糊，应同时限制 λ 的学习范围（如 500~560nm），防止模型用 λ 来补偿 σ

8.3.5.2.4. PSF 空间域匹配度分析

⚠️ 重要说明：以下分析基于逐像素 MAE。由于采集的 PSF 和仿真的 PSF 之间存在未配准的平移（峰值偏移平均约 1.3~6.5 像素，最大达 13~33 像素），峰值区的逐像素 MAE 会被平移误差污染。

但以下分析仍然有效：

• 整体 MAE：97% 像素为背景区，背景均匀，平移对整体 MAE 影响极小
• 背景区 MAE：平移几乎不影响均匀背景的统计
• 峰值区 MAE：受平移污染，仅供趋势参考；已补充平移不变度量（NCC）验证

整体 MAE

实验	旧 Focus MAE	新 Focus MAE（未重训）	新 Focus MAE（重训 1 epoch）	vs 旧改善
完整 detector	0.01169	0.00377	0.00377	-67.7%
无截断	0.01171	0.00316	0.00378	-67.7%
无高斯	0.01131	0.00325	0.00324	-71.3%

实验	旧 Defocus MAE	新 Defocus MAE（重训）	vs 旧改善
完整 detector	0.01434	0.00453	-68.4%
无截断	0.01434	0.00455	-68.3%
无高斯	0.01405	0.00431	-69.3%

关键发现：

• 整体 MAE 改善约 70%，几乎全部来自 dataset 预处理，重训 1 epoch 对空间域 MAE 几乎没有额外贡献
• 这是因为 1 epoch 训练主要优化频域 loss，而空间域结构由模型初始化决定

按区域分解（Focus）

区域	旧-完整	新-完整（重训）	旧-无截断	新-无截断（重训）	旧-无高斯	新-无高斯（重训）
背景 (<0.05)	0.01020	0.00238	0.01020	0.00238	0.00987	0.00190
中间 (0.05~0.5)	0.1209	0.1267	0.1149	0.1256	0.1142	0.1202
峰值 (≥0.5)	0.4357	0.4347	0.5488	0.4909	0.4434	0.4367

⚠️ 峰值区 MAE 受平移污染，以下补充平移不变度量验证。

平移不变验证（NCC 归一化互相关）：

实验	NCC (init)	NCC (train)	训练后改善
完整 detector	0.4973	0.4985	+0.0012
无截断	—	—	—
无高斯	—	—	—

平移不变验证（对齐后 MAE）：

实验	对齐前 MAE	对齐后 MAE	平移导致误差占比
完整 detector	0.003830	0.003502	~8.6%

关键发现：

1. 背景区 MAE 降低 80%：这是 dataset 改动的直接效果，与是否重训无关，且不受平移影响

2. 峰值区存在平移，但训练后仍有改善趋势（无截断 Peak MAE 0.549→0.491）

   • 无截断改善最大：旧 checkpoint 训练 7 epoch 参数扭曲，新 dataset 下重训 1 epoch 后参数回归
   • 平移不变 NCC 验证：训练后仅改善 0.0012，说明 1 epoch 对结构改善确实有限

3. 中间区 MAE 反而略有上升（旧→新：0.121→0.127）

   • 真实数据背景被扣得更干净后，中间亮度区域边界更清晰，仿真模型尚未充分调整

8.3.5.2.5. FWHM 对比

实验	真实 FWHM	仿真 FWHM	偏差
完整 detector	19 px	37 px	+95% ❌
无截断	19 px	39 px	+105% ❌
无高斯	19 px	22 px	+16% ✅

关键发现：

• 无高斯的 FWHM 最接近真实（22 vs 19 px），完整/无截断因高斯模糊严重过宽
• FWHM 基于峰值位置计算，不受平移影响

8.3.5.2.6. 能量集中度

实验	真实中心 41×41 能量占比	仿真中心 41×41 能量占比
完整 detector	48.1%	64.8%
无截断	48.1%	64.7%
无高斯	48.1%	64.4%

关键发现：

• 所有仿真 PSF 的能量都过于集中在中心（比真实高约 16 个百分点）
• 仿真背景为 0，真实数据扣除后仍有微弱残余背景（corner mean=0.0009），会分散能量
• 无高斯的集中度最接近真实（64.4% vs 64.8%）

8.3.5.2.7. M4 精度分析

实验	旧 M4 init	旧 M4 train	新 M4 init	新 M4 train	Train 改善
完整 detector	0.7048	0.7138	0.6797	0.6905	-3.3%
无截断	0.6155	0.6518	0.6154	0.6290	-3.5%
无高斯	0.7277	0.7097	0.7043	0.6849	-3.5%

关键发现：

• 新 dataset 下 M4 MAE 普遍比旧 dataset 低 3~4%
• 无高斯训练后 M4 MAE = 0.6849，是三实验中最低的
• 无高斯是唯一一个训练后 M4 比 init 改善的配置（0.7043 → 0.6849）

8.3.5.2.8. 综合结论

中位数背景扣除的效果（重训验证）

维度	效果
训练收敛速度	Epoch 0 loss 降低 17%，初始匹配更好
物理参数合理性	波长不再被扭曲到非物理值（无截断 403→466nm，无高斯 518→535nm）
PSF 空间匹配	整体 MAE 降低 70%，背景区降低 80%
峰值结构	无截断 Peak MAE 降低 12%（参数回归），无高斯降低 1.6%
M4 精度	提升 3~4%，无高斯最优

三个配置的最终排名（新 dataset + 1 epoch）

排名	配置	优势	劣势
⭐ 1	无高斯	Test loss 最低（0.134），MAE 最低（0.00324），波长最接近 530nm，FWHM 最接近真实，M4 最优	背景仍为 0，与真实有微弱差距
2	完整 detector	参数稳定	FWHM 过宽（37px），loss 比无高斯高 33%
3	无截断	Peak 区改善最大（重训后）	FWHM 最宽（39px），高斯模糊 unnecessary

关键认知更新

1. 高斯模糊完全不必要：无高斯在所有指标上都是最优或接近最优，且没有高斯模糊的副作用（FWHM 失真）。

2. 1 epoch 重训对空间域改善有限：整体 MAE 的 70% 改善来自 dataset 预处理本身，重训 1 epoch 主要改善了频域 loss 和峰值区结构（尤其是无截断）。

3. 下一步应继续训练：当前仅 1 epoch，峰值区和中间区仍有优化空间。建议用新 dataset + 无高斯配置训练 10~20 epoch。

你的目标是"仿真 PSF 与真实 PSF 在空间域上尽可能一致"，而不仅仅是"频域 loss 最小"。当前 fft_psf_loss 是一个好的正则化项（保证平移不变性），但不应是唯一的监督信号。空间域的 SmoothL1Loss 直接惩罚像素级差异，包括背景噪声底板，是达到"视觉一致"的必经之路。

8.4. Interpolate 模式选择分析

在 Zernike2PSF_layer 中，仿真生成的 PSF（像素间距约 2.2μm）需要通过 F.interpolate 下采样到真实 CCD 的像素尺度（4.5μm）。下采样方式直接影响仿真 PSF 的空间形态，进而影响物理参数（波长 λ、离焦量 defocus）的学习精度。

模式	数学本质	物理意义
bilinear	双线性插值取点	从高密度网格中插值采样
area	区域像素平均	像素积分（真实相机行为）

8.4.1. 为什么 area 更物理？

真实相机的每个像素是一个 4.5μm × 4.5μm 的物理面积，接收的是该区域内总光强：

像素值 = ∫∫ PSF(x,y) dx dy  （在像素物理面积内）

这等价于对理想 PSF 做一个 box filter（低通滤波）。mode='area' 在缩小图像时，输出像素 = 输入区域内像素的平均值，正好模拟这一物理过程。

而 bilinear 只是从高密度采样点中插值取点，没有积分效应，旁瓣能量被压缩到亚像素级别，不可见。

8.4.2. 实验对比

8.4.2.1. 测试条件

• 数据集：相同的训练/测试划分
• 训练目标：fft_psf_loss + 0.1 × m4_loss
• 训练时长：2 epoch
• 唯一变量：interpolate 模式（bilinear vs area）

8.4.2.2. 结果对比

指标	bilinear (temp_test_noise)	area (temp_test_arae)	变化
λ	535 nm	522 nm	更接近真实 530nm
defocus	-0.301 μm	-0.291 μm	更接近真实 -0.265μm
训练 loss	0.143	0.134	相当
M4 MAE	0.676	0.646	降低 4.4%
Init Avg Loss	0.140	0.117	降低 16%

8.4.2.3. 遇到的问题：旁瓣仍然不明显

尽管 area 模式改善了能量分布，但视觉上仿真 PSF 的"余光"/旁瓣仍然远不如真实数据明显。

原因分析：

1. 下采样比例太大：scale = 0.49（约 1/2），旁瓣结构被压缩到 1-2 个像素宽
2. 理想衍射旁瓣本身就极弱：艾里斑第一亮环强度 ≈ 中心峰值的 1.75%
3. 真实相机的"余光"来源复杂：像素积分 + 光学像差 + 散射光 + 读出噪声底板