简介：HoloLib_MATLAB 是一套面向计算全息（CGH）学习与实验的 MATLAB 函数库，适合光学、光电信息及相关专业的学生与研究人员使用，用于简化计算全息编程流程、快速搭建相位图生成与传播仿真实验。压缩包共 128 个文件，约 17.05MB，以 36 个 .m 源码函数为主体，辅以 png、bmp 等示例图像与说明文档，源码集中在 source 目录，涵盖路径初始化、图像边界扩充与矩阵扩展、轴锥透镜与闪耀光栅及菲涅尔透镜相位生成、迭代傅里叶与迭代菲涅尔算法、两次与三次傅里叶变换的菲涅尔传播、DOE 元件传播、计算重构、图像拼接、相位保存与 LCOS 显示等模块，函数职责清晰、便于按需调用与二次开发。已有 1285 人学习下载，可作为计算全息入门练习与算法验证的参考工具库。

1. 拿到一份 CGH 函数库，先别急着跑 demo

做光学仿真的人大概都有过这种经历：想验证一个相位恢复算法，翻遍手头资料，重建脚本、传播核、误差评价各写一遍，跑出来的重建图却和文献对不上。问题往往不在算法本身，而在传播距离的单位、采样间隔的定义、相位包裹的处理这些底层细节上。HoloLib_MATLAB 就是冲着这个痛点来的——一套面向计算机生成全息图（CGH）的 MATLAB 函数库，把球面波、平面波、菲涅尔与角谱传播、GS 与 Gerchberg-Saxton 类迭代、相位编码这些常用环节封装成可调用函数。它适合做光学、显示、衍射器件方向的研究生和工程师，也适合需要快速搭一个 CGH 验证环境的人。系统开源意味着你能直接看到每个函数的实现，而不是对着黑匣子猜参数。

2. 函数库的骨架：传播核与迭代引擎怎么分工

拿到一个 CGH 库，第一件事不是跑示例，而是搞清楚它把哪些物理环节拆成了独立函数。HoloLib_MATLAB 的组织逻辑基本遵循「输入场 → 传播 → 目标面约束 → 反传 → 更新」这条闭环，每个环节对应一组函数。理解这个分工，后面调参才知道该动哪一层。

2.1 三类传播函数的选择依据

CGH 里最核心的运算是光场传播，不同传播距离和采样条件对应不同近似。这个库通常提供三类：菲涅尔衍射（近场到中场）、角谱传播（任意距离、采样严格）、球面波卷积（适合大尺寸快速计算）。选错传播核是新手最常见的翻车点——用菲涅尔近似去算毫米级近场，或者用角谱却在频域采样不足导致混叠。

判断依据可以记成一条经验线：传播距离 z 远大于孔径尺寸平方除以波长时，菲涅尔近似够用；z 与孔径尺寸可比时，老老实实上角谱。角谱的代价是必须做补零，否则循环卷积会污染结果。库里的角谱函数一般会暴露padFactor或类似的补零参数，这个参数不是摆设，后面避坑章节会专门说。

2.2 迭代引擎的输入输出约定

GS 类迭代是 CGH 生成的主力。库里的迭代函数通常接受初始相位、目标振幅、传播算子、迭代次数作为输入，返回相位分布和误差曲线。关键在于它约定的是「目标面振幅约束」还是「目标面强度约束」——两者在误差评价上差一个平方根，混用会让收敛曲线看起来很好但重建图发灰。

我一般会先跑一个最小闭环确认约定：

% 最小闭环：确认迭代函数的输入输出约定 lambda = 532e-9; % 波长，单位米 pitch = 8e-6; % 像素间距，单位米 N = 512; % 采样点数 z = 0.3; % 传播距离，单位米 % 构造目标振幅（简单矩形孔径，便于肉眼判断） target = zeros(N); target(N/4:3*N/4, N/4:3*N/4) = 1; % 初始相位：随机相位能打破对称性，避免迭代卡在局部极小 phase0 = 2*pi*rand(N); % 调用库内迭代函数（函数名以实际库为准） [phaseOut, errCurve] = gs_iterate(target, phase0, lambda, pitch, z, 50); % 检查误差曲线是否单调下降，不降说明约定或参数有问题 figure; plot(errCurve); xlabel('迭代次数'); ylabel('误差');

这段代码的重点不在跑通，而在验证三件事：误差曲线是否在前几次迭代明显下降、phaseOut的取值范围是[-pi, pi]还是[0, 2*pi]、传播距离单位是米还是毫米。这三个约定任何一个对不上，后面所有结果都是错的。参数上，pitch和lambda必须同单位，N取 2 的幂是为了 FFT 效率，z的选择要让目标面落在传播核的有效区内。

2.3 相位编码与重建评价

迭代出来的相位是连续分布，真正上器件或上 SLM 前要做相位编码——把连续相位量化到有限级数，或者转成振幅全息的编码图案。库里的编码函数一般支持二值、多阶和误差扩散几种模式。量化级数越低，重建噪声越大，但器件实现越简单，这是典型的工程折中。

重建评价环节容易被忽略。库里通常提供 RMSE、衍射效率、信噪比几个指标，但它们的计算前提不同：RMSE 要在归一化振幅上算，衍射效率要区分「目标区能量」和「总能量」。我见过不少人直接拿原始强度算 RMSE，数值看着漂亮，实际重建图对比度一塌糊涂。

3. 从目标图到相位图：一次完整复现的四个步骤

这一章走一遍完整流程，从读入目标图到输出可用的相位全息图。步骤按实际动手顺序排，每步给出可抄的代码和参数说明。假设你已经把库加入 MATLAB 路径，函数名以库内实际命名为准，逻辑不变。

3.1 目标图预处理与采样对齐

目标图的像素尺寸必须和传播网格对齐，否则传播核的采样间隔和目标图对不上，重建会出现缩放错位。常见做法是把目标图缩放到N×N，并确认物理尺寸N*pitch与你的器件孔径一致。

% 目标图预处理：灰度化、归一化、尺寸对齐 img = imread('target.png'); if size(img,3) == 3 img = rgb2gray(img); % 转灰度，CGH 一般只用振幅信息 end img = im2double(img); img = imresize(img, [N N]); % 对齐到传播网格尺寸 target = img / max(img(:)); % 归一化到 [0,1]，避免能量尺度混乱

这里imresize的顺序有讲究：先缩放到N×N再归一化，而不是反过来，因为插值会改变极值。归一化用max(img(:))而不是固定 255，是为了兼容不同位深的输入图。如果目标图本身带相位信息（比如复振幅目标），这一步要改成分别处理振幅和相位通道。

3.2 传播核参数配置

传播核的参数决定了整个仿真的物理可信度。以角谱为例，关键参数是补零倍数和频域坐标构造方式。补零不足会导致循环卷积混叠，补零过度则浪费算力且可能引入数值噪声。

% 角谱传播核配置 padFactor = 2; % 补零倍数，2 倍是常用起点 Npad = N * padFactor; dx = pitch; % 空域采样间隔 fx = (-Npad/2:Npad/2-1) / (Npad*dx); % 频域坐标，单位 1/m [FX, FY] = meshgrid(fx, fx); % 角谱传递函数：注意 evanescent 波要截断，否则出现虚数 arg = 1 - (lambda*FX).^2 - (lambda*FY).^2; H = exp(1i*2*pi*z/lambda*sqrt(max(arg,0))); % 传播函数：补零 → FFT → 乘传递函数 → IFFT → 裁剪 propagate = @(U) crop(ifft2(fft2(padarray(U,[N*(padFactor-1)/2 ... N*(padFactor-1)/2]) .* H), N);

max(arg,0)这一步是血泪经验：角谱传递函数在arg<0时对应倏逝波，物理上快速衰减，数值上会变成虚数导致结果发散。截断为 0 是标准处理。padFactor=2对多数目标够用，如果重建边缘出现周期性伪影，加到 3 或 4 再试。crop是库内或自写的裁剪函数，把补零后的结果裁回N×N。

3.3 迭代生成与收敛监控

迭代环节要监控的不只是误差曲线，还有相位分布的直方图。相位过度集中意味着迭代陷入局部极小，重建图会出现明显伪影。

% 迭代生成相位，带收敛监控 maxIter = 100; errCurve = zeros(maxIter,1); phase = 2*pi*rand(N); % 随机初相位 for k = 1:maxIter U0 = target .* exp(1i*phase); % 目标面复振幅 Uz = propagate(U0); % 正向传播 Uz = target .* exp(1i*angle(Uz)); % 目标面振幅约束 Uback = propagate(conj(Uz)); % 反向传播（共轭技巧） phase = angle(Uback); % 更新相位 errCurve(k) = norm(abs(Uz(:)) - target(:)) / norm(target(:)); end

反向传播用conj加正向传播实现，是 GS 的经典技巧，避免单独写反传核。误差用归一化 L2 范数，数值在 0.05 以下通常认为收敛可接受。如果 100 次迭代后误差还在 0.2 以上，先别加迭代次数，回头检查传播核的补零和单位。

3.4 相位量化与结果导出

连续相位要量化到器件支持的级数。量化方式直接影响重建质量，误差扩散比直接取整好，但计算量略大。

% 相位量化到 8 级，并导出 levels = 8; phaseQ = round((phase + pi) / (2*pi) * (levels-1)) / (levels-1) * 2*pi - pi; % 重建验证 Ufinal = propagate(target .* exp(1i*phaseQ)); recon = abs(Ufinal); recon = recon / max(recon(:)); % 导出为 PNG，注意归一化到 0-255 imwrite(uint8(255 * (phaseQ + pi) / (2*pi)), 'phase_hologram.png');

量化级数levels按你的 SLM 或加工能力定，8 级是常见起点。导出前把相位映射到[0,255]灰度，是因为多数器件驱动按灰度读相位。重建验证这一步不能省，量化后的重建图要和量化前对比，差异过大说明级数不够或目标图对比度太高。

4. 避坑与排查：五个让重建图对不上的原因

这一章按「现象 → 原因 → 解决」列五条实际踩过的坑。每条都对应前面章节的某个参数或约定，出问题时按顺序排查能省不少时间。

4.1 重建图整体缩放错位

现象：重建出的图案形状对，但尺寸比目标大一圈或小一圈。原因：传播距离z的单位和目标图物理尺寸不匹配，或者pitch用了像素间距却按毫米填。解决：统一所有长度单位为米，用N*pitch反算物理孔径，和器件参数核对一遍。菲涅尔传播对z敏感，角谱对pitch敏感，先确认用的是哪个核。

4.2 边缘出现周期性条纹

现象：重建图四周有规律性亮纹或暗纹。原因：角谱传播补零不足，循环卷积把对侧能量卷了回来。解决：把padFactor从 2 提到 3 或 4，观察条纹是否消失。如果目标图本身在边界有强跳变，先做边缘渐变再传播，能进一步压住高频泄漏。

4.3 误差曲线不降反升

现象：迭代前几次误差下降，之后震荡或上升。原因：目标面约束用错了量纲——把强度当振幅约束，或者相位更新时符号反了。解决：确认约束是target .* exp(1i*angle(Uz))而不是sqrt(target)，反传用conj而不是直接ifft。另外检查初相位是否全零，全零初相位在对称目标上容易卡死，换随机相位。

4.4 量化后重建噪声骤增

现象：连续相位重建正常，量化到 4 级或 8 级后噪声明显。原因：目标图动态范围太大，低阶量化无法表达细节。解决：先对目标图做伽马校正压缩动态范围，或者改用误差扩散量化。级数实在上不去，就在迭代时直接把量化步骤嵌进去，让算法自己适应量化误差。

4.5 不同机器结果不一致

现象：同一份代码换台机器跑，重建图有细微差异。原因：随机初相位没固定种子，或者 FFT 实现差异导致浮点误差累积。解决：迭代前rng(0)固定种子，误差评价用相对范数而不是绝对范数。跨平台对比时，把中间相位分布存成.mat逐层比对，定位是哪一步开始分叉。

5. 进阶技巧：把迭代嵌进量化，少走一轮返工

前面流程是「先迭代出连续相位，再量化」，这个顺序在级数低的时候会浪费迭代成果。更稳的做法是把量化算子直接嵌进迭代循环，让算法在每一步都「知道」最终要量化，收敛到的相位天然适配器件级数。代价是每次迭代多一次量化运算，但省掉了量化后重建不达标再返工的时间。

具体改法是在目标面约束之后、反传之前插入量化：

% 量化感知迭代：在循环内嵌量化 levels = 8; for k = 1:maxIter U0 = target .* exp(1i*phase); Uz = propagate(U0); Uz = target .* exp(1i*angle(Uz)); % 关键改动：反传前把相位量化，模拟器件行为 phaseQ = round((angle(Uz)+pi)/(2*pi)*(levels-1))/(levels-1)*2*pi - pi; Uz = abs(Uz) .* exp(1i*phaseQ); Uback = propagate(conj(Uz)); phase = angle(Uback); errCurve(k) = norm(abs(Uz(:)) - target(:)) / norm(target(:)); end

这个改动的逻辑是：量化误差在迭代中被反复「看到」并补偿，最终相位分布对量化更鲁棒。实测在 4 级量化下，量化感知迭代的 RMSE 通常比「先迭代后量化」低一截，重建图对比度也更均匀。参数上，levels要和实际器件一致，别用 8 级迭代却上 4 级器件，那等于白做。

验证方法也简单：同一目标图跑两条路径，一条先迭代后量化，一条量化感知迭代，把两条重建图的 RMSE 和衍射效率列成表对比。如果量化感知那条没有优势，多半是levels填错或者量化公式的相位范围没对齐。

还有一个容易忽略的点：量化感知迭代的误差曲线不会像连续迭代那样单调平滑，会有小幅台阶，这是量化算子引入的离散跳变，属于正常现象，别看到曲线抖就以为发散。判断收敛看整体趋势和最终重建图，不看单点数值。

从那以后我每次搭 CGH 流程，都强制先跑一遍最小闭环确认单位、约定和补零，再上目标图。这套习惯帮我省掉的返工时间，比任何加速技巧都值。希望帮到你。

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