简介：面向全息显示、计算光学成像等领域的开源MATLAB函数库，旨在降低计算机生成全息图（CGH）算法的编程门槛，适合从事数字全息、衍射光学元件设计的研究者与工程师。库内封装了从相位生成到传播重构的完整流程，包括透镜相位、闪耀光栅相位、轴锥透镜相位计算，以及迭代傅里叶/菲涅尔算法、两次/三次傅里叶菲涅尔传播方法等，并附初始化路径、边界扩充、图像拼接、相位保存及LCOS显示等实用工具。压缩包共128个文件，以36个.m源文件为核心，配有72张示例效果图、12个BMP测试图像和说明文档、许可证等，整体大小约17MB，source目录内源码组织清晰，便于按模块学习与二次开发。目前已有1283人学习下载，借助示例图片可快速验证算法输出，是入门计算全息编程和搭建实验原型的实用参考资料。
1. HoloLib_MATLAB 解决什么：把 CGH 从“跑通”变成“能复用”
做全息显示或者激光投影仿真的人，多半经历过同一条路：手上有目标图像，想得到一张能加载到空间光调制器（SLM）上的相位全息图，于是先搜到一段 GS 算法、又找一段 Fresnel 衍射，跑通一版后换了个参数就全乱，最后在 MATLAB 路径里堆出test1.m、final_v2.m。标题里的 HoloLib_MATLAB 看起来是个函数库，本质上是在补这种“每次重写”的缺：把传播计算、相位恢复迭代、量化输出拆成少而稳的 MATLAB 函数，输入输出约定一致，换波长、换像素间距、换重建距离时只改参数不碰算法。适合正在做相位型 SLM 全息显示、CGH 数值重建，或者想从“复制脚本”转成“维护一个库”的工程师。追求的是能复现、能调参、能不返工。
2. CGH 函数库的核心算子：角谱、菲涅尔与 GS 相位恢复的取舍
2.1 为什么函数库统一用角谱衍射而不默认菲涅尔
CGH 库里最容易被调用错的一步是“衍射传播”。常见 CGH 流程分两步：正向传播从全息面到像面，反向传播从约束后的像面回到全息面。两步能用同一个传播算子最省事，而角谱法在这个点上天然占优，它只是同一个传递函数取正负距离，反向传播就是z变成-z。所以我做这类库时默认只放角谱法一个底层传播函数，菲涅尔只作为可选加速项，不参与默认路径。
角谱法的实现可以写成：对复振幅U做一次 FFT，乘上频域传递函数H，再做一次 IFFT。MATLAB 里多数人会选择自己写，因为官方工具箱偏通用，频繁fft2/ifft2时不如在库里集中统一 FFT 布局和归一化。下面这个版本可以直接进函数库，输入是复振幅矩阵、像素间距和传播距离，输出仍是复振幅矩阵：
function U2 = ang_spectrum(U1, dx, dy, lambda, z) % 角谱衍射传播，正向与反向共用同一个函数 % U1 : 源平面复振幅（二维矩阵，double） % dx,dy : 像素间距，单位与 lambda 一致（通常用米） % lambda: 波长 % z : 传播距离，正向>0，反向<0 [Ny, Nx] = size(U1); fx = (-Nx/2 : Nx/2-1) / (Nx * dx); % 空间频率，单位 1/m fy = (-Ny/2 : Ny/2-1) / (Ny * dy); [FX, FY] = meshgrid(fx, fy); k = 2 * pi / lambda; rad = 1 - (lambda * FX).^2 - (lambda * FY).^2; H = zeros(Ny, Nx); H(rad > 0) = exp(1j * k * z * sqrt(rad(rad > 0))); U2 = ifft2(fft2(U1) .* ifftshift(H)); end说明几个容易踩的细节。fx/fy按-N/2到N/2-1排列，这是把频域放在数组中心；H也是中心坐标顺序。MATLAB 的fft2默认把零频放在左上角，所以乘法前需要ifftshift(H)，把中心坐标挪到 FFT 的阵列顺序。偶数尺寸下它和fftshift等价，但既然函数库要接任意偶数网格，统一用ifftshift更稳。rad > 0之外的倏逝波直接置零，长距离传播时不会出现NaN。
这个函数不重采样、不改变像素数，像素间距由外部传入，原因是 SLM 经常不是正方形像素域：1920×1080 的相位型器件，dx和dy可以不同。很多网上版本把两个方向混成一个L/N，用在矩形 SLM 上会直接把目标图横向比例拉歪。至于角谱、卷积型菲涅尔、单次 FFT 菲涅尔的取舍，我一般按这张表记：
| 传播方案 | 适用距离 | 函数库里的角色 |
|---|---|---|
| 角谱 | 毫米到厘米量级，反向传播直接取负距离 | 默认算子 |
| 卷积型菲涅尔 | 中等距离，满足 Fresnel 条件 | 可选算法，接口保持一致 |
| 单次 FFT 菲涅尔 | 较大距离，输出采样自动变粗 | 不推荐作为默认 |
默认角谱不是因为它在所有距离上都更准，而是因为同一个函数实现正反传播最稳。菲涅尔在近场悄悄失效时很难查，角谱出问题反而是明显的边缘混叠，更容易定位。
2.2 GS 算法的最小实现：随机初相位、振幅约束与迭代终止判断
CGH 的本质问题是相位恢复：知道像面强度，求全息面相位。最常见的迭代是 Gerchberg-Saxton，简称 GS。它的思路很直白：全息面只保留相位，像面只保留目标振幅，反复正传、反传。对刚接触全息仿真的人来说，GS 像个小黑匣子，其实核心逻辑就十行。
function phi = gs_kinoform(target, dx, dy, lambda, z, iter, init) % GS 算法生成纯相位全息图 % target : 目标强度图，double，范围 0~1，尺寸与网格一致 % dx,dy : 像素间距；lambda 波长；z 重建距离 % iter : 迭代轮数；init 取 'rand' 或 'zeros' [Ny, Nx] = size(target); if strcmpi(init, 'rand') phi = 2 * pi * rand(Ny, Nx); else phi = zeros(Ny, Nx); end amp = sqrt(max(target, 0)); % 强度转振幅 for k = 1:iter u = exp(1j * phi); % 全息面：纯相位约束 img = ang_spectrum(u, dx, dy, lambda, z); % 正传到像面 img = amp .* exp(1j * angle(img)); % 像面：只换振幅 back = ang_spectrum(img, dx, dy, lambda, -z); % 反传到全息面 phi = angle(back); end phi = mod(phi, 2 * pi); end需要理解两件事。第一，target是强度分布，不是振幅分布，所以代码里用sqrt转成振幅。有人把uint8图像直接当成振幅传进来，重建结果会明显过曝，因为 255 和sqrt(255)差了平方关系。第二，全息面每次迭代都强制exp(1j * phi)，这是纯相位约束，对应相位型 SLM 的物理能力；如果目标器件支持复振幅调制，约束方式才需要改。
迭代轮数不是越多越好。纯亮点或二值图案，一般 10 到 20 轮就能看；连续灰度照片，40 到 60 轮足够；超过 100 轮之后 PSNR 往往不再涨，反而容易出现背景颗粒。所以函数库里最好把iter当独立参数暴露出来，不要在函数内部写死。
初相位用rand会让每次结果不同。实验要可复现时，调用前先rng(42, 'twister')；想要工程稳定性，保留上一轮的相位作为下一轮初值，这也是我后来最常用的“后悔药”做法。GS 收敛到局部最优是常态，不是 bug，对比算法好坏时必须锁初相位。
2.3 传播函数接口设计：dx、dy、波长、距离四个入参就够了
函数库最怕的是接口风格漂移。有的函数收L，有的收dx，有的直接读全局变量，到最后连自己都记不住哪个单位是米、哪个是毫米。我一般会把所有传播类函数统一成ang_spectrum(U, dx, dy, lambda, z)这个签名，四个物理参数全部用米，相位一律用弧度，复数场一律用 double 矩阵。
% 推荐用 struct 集中管理参数，避免调用时把 lambda 和 z 写反 cfg.dx = 3.74e-6; % 像素间距 cfg.dy = 3.74e-6; cfg.lambda = 532e-9; % 波长 cfg.z = 0.05; % 传播距离 U2 = ang_spectrum(U1, cfg.dx, cfg.dy, cfg.lambda, cfg.z); back = ang_spectrum(U2, cfg.dx, cfg.dy, cfg.lambda, -cfg.z);这个接口还有一个好处：反向传播只需要把z取负，不用再写一个ang_spectrum_inverse，函数库少一个函数就少一类调用错误。如果某天要换成 Fresnel 卷积版本，函数名和入参保持不变，只改函数体。量化函数、强度转振幅函数也应该走同一套约定，返回值要么是 double 弧度，要么是 uint8 灰度，不要在同一流程里混用。
3. 本地复现最小流程：从目标灰度图生成 8 位相位全息图
3.1 最小目录结构与可复制的调用脚本
拿到函数库之后，第一步不是急着看算法，而是把目录结构固定下来。我习惯按职责分三层：底层传播函数、上层相位恢复、应用层示例脚本。下面这种结构足够一个单人维护的小库使用：
holo_lib/ ang_spectrum.m gs_kinoform.m quantize_phase.m demos/ make_kinoform.mang_spectrum.m只做衍射传播，gs_kinoform.m只做相位恢复，quantize_phase.m只做相位到灰度映射。示例脚本单独放一个目录，避免每次跑实验都去改库文件。最小可复现脚本可以写成：
cfg.dx = 3.74e-6; % SLM 像素间距，参考常见相位型器件 cfg.dy = 3.74e-6; cfg.lambda = 532e-9; % 绿光波长 cfg.z = 0.05; % 重建距离，先按 5 cm 跑通 I = imread('target.png'); I = im2double(rgb2gray(I)); % 彩色先转灰度 I = imresize(I, [512 512]); % 目标尺寸与网格一致 phi = gs_kinoform(I, cfg.dx, cfg.dy, cfg.lambda, cfg.z, 40, 'rand'); K = quantize_phase(phi, 8); % 8 位相位图，0~255 imwrite(K, 'kinoform_8bit.png');配套的量化函数建议单独写，不要在每个脚本里复制一遍。量化不是简单的四舍五入，它涉及边界和 SLM 灰度响应，见第 4 章。先把最小版本跑通，再用离线重建验证，后面调参才有依据。
function Q = quantize_phase(phi, bits) % 相位量化：把 0~2pi 相位映射到 uint8 灰度 levels = 2^bits; Q = uint8(round(mod(phi, 2*pi) / (2*pi) * (levels - 1))); end3.2 关键参数怎么设：SLM 像素间距、目标尺寸与重建距离
三个参数决定重建结果的基本形态。像素间距dx直接决定全息图能“容纳”多少角度信息：间距越大，可记录的高频角越小，重建画面尺寸也越大。波长lambda影响传递函数的相位曲率，SLM 有设计波长，不要为了仿真好看随便换。重建距离z决定图像大小和角谱采样是否安全，距离越长，H的相位变化越剧烈，混叠风险越高。
| 参数 | 常见取值 | 改大/改小的影响 |
|---|---|---|
| dx / dy | 3.2 ~ 4.5 μm | 变大则重建画面变大，混叠风险也变大 |
| lambda | 405 / 532 / 633 nm | 偏离 SLM 设计波长时，相位深度和亮度都不对 |
| z | 5 ~ 20 cm | 变大则图像变大，角谱采样更容易破 |
| iter | 20 ~ 60 | 太少不收敛，太多可能放大背景噪声 |
| 目标尺寸 | 与网格一致或略小 | 小目标建议补零画布，别硬拉伸 |
这里最容易翻车的是“目标图像分辨率直接等于全息图分辨率”的假设。SLM 是 1920×1080 不代表你一定要生成 1920×1080 的全息图。先用 512×512 跑通算法，再放到 SLM 中央区域，周围相位填 0，能少踩很多暗角、边缘错位的坑。
3.3 灰度图归一化的坑与“补零画布”做法
灰度图归一化看起来简单，实际坑不少。im2double只能把 uint8 变成 0 到 1 的 double，不等于它已经是适合 GS 的目标强度。如果图像整体偏暗，最大值只有 0.6，建议先除以最大值再进 GS，让最大强度为 1；如果整张图就是暗调风格，强行归一化反而会提亮背景，需要按内容决定。
另一个高频问题是长宽比。目标图是 4:3，SLM 是 16:9，直接imresize会把画面压扁。常见做法是补零画布：在黑色画布中央放目标图，目标图周围留至少 64 到 128 像素的空白。这个空白不是装饰，它给衍射加了“边界保护”，能明显减少 FFT 周期延拓造成的重影。
% 补零画布：目标放中央，周围黑色背景留距离 canvas = zeros(1024, 1024); [h, w] = size(I); y0 = floor((1024 - h)/2) + 1; x0 = floor((1024 - w)/2) + 1; canvas(y0:y0+h-1, x0:x0+w-1) = I; target = canvas / max(canvas(:));注意这里的max(canvas(:))是全局最大值，不是局部。黑色背景区域本来就该是 0，不需要抬高。如果你用imresize把 200×200 的文字图强行放大到 1024×1024，边缘会产生磨砂感，那是因为插值把硬边变成了斜边，GS 的迭代会把这种模糊当成目标去拟合，结果更糊。补零画布比拉伸更接近真实 SLM 的利用方式。
4. 避坑向：重建伪影、散斑噪声与三个最让人头疼的参数
4.1 现象：重建有块状伪影，迭代越跑越花
第一次用 GS 生成灰度照片全息图，常见的结果是：目标图案能认出来，但背景像磨砂玻璃，黑色区域全是颗粒，迭代到 50 轮后这些颗粒不但没消失，反而更碎、更密。这不是算法写错，是黑色背景被当成了硬约束。GS 把背景振幅约束成 0，随机相位在像面上互相干涉，最后形成散斑噪声。
解决方式有两种。第一种是补零画布，把目标区域缩小到整个网格的 1/4 到 1/2，背景作为保护带，重建时只看中心 ROI，四周噪声不评分。第二种是给背景一个很小的松弛幅度，而不是 0：
% 像面约束：背景不压到 0，给一点松弛，散斑会均匀很多 amp_target = sqrt(max(target, 0)); mask = target > 0.05; % 目标区域 mask img = ((0.9 * mask + 0.3 * (~mask)) .* amp_target + ... 0.7 * (~mask) .* abs(img)) .* exp(1j * angle(img));这段代码只作为示意见，具体权重要按你的画布比例调。常见调法是背景松弛从 0.2 到 0.4，太低仍然有条纹，太高会让背景发灰。散斑完全消掉不可能，纯相位全息图天然带散斑，目标是把颗粒控制在背景区域而不是目标边缘。
4.2 现象：四周出现环形纹和重复网格
如果重建图左侧出现一条和右侧对称的亮带，或者整个画面像被裁成了两半再拼起来，问题大概率出在角谱采样率。角谱传递函数在频域边缘相位变化很快，z越大，相位曲率越陡，固定像素数的网格就越难满足奈奎斯特条件。我常用的经验参考是：
% 粗略判断角谱是否安全：z_safe 只是一个数量级参考 z_safe = min([Nx * dx^2 / lambda, Ny * dy^2 / lambda]) / 1.5; fprintf('当前 z=%.4f m，经验安全参考 %.4f m\n', cfg.z, z_safe);当z明显超过这个量级，最直接的解法不是换菲涅尔，而是加大网格并补零。同一物理尺寸下，512 网格变 1024 网格，安全距离参考大约提升 4 倍。补零之后全息图的实际有效区域不变，但频域采样变密，混叠自然下降。对于 1080P 级 SLM，生成全息图前先做一次padarray到 2048×2048，对最终重建质量有明显改善。
需要注意，这个z_safe是工程经验值，不是严格界线。它帮你判断“当前设置有没有可能破采样”，真正确认还要看重建图边缘有没有高频折叠。出现网格重影时，先别调迭代轮数，先降z或加网格，往往一次就好。
4.3 现象：8 位量化后亮度漂移
GS 跑出来的phi是 double 弧度，加载到 SLM 前必须转成 8 位灰度。这里的坑在round和floor的差别。用floor会让整个相位系统平均偏掉约半个 LSB，肉眼看着只是亮度略低，但灰度照片的暗部层次会丢。更隐蔽的问题是相位在 2π 边界：接近 2π 的像素四舍五入到 255，接近 0 的像素到 0，两者在物理上是同一个相位，但 SLM 的灰度响应不是理想线性，边界上会压出一条暗线。
我现在的量化习惯是：最高一级留出来做校准余量，避免把相位顶到 2π 极限：
% 保留最高一级作为 SLM 校准余量，减少 wrap 边界伪影 Q = uint8(round(mod(phi, 2*pi) / (2*pi) * (2^bits - 2)));如果你手头有 SLM 厂商标定的相位灰度表，量化的终点应该用标定表的灰阶边界，而不是理论的0~2π。在没有标定数据之前，先不要过度纠结量化精度，8 位量化带来的亮度损失通常在可接受范围，真正影响大的是相位调制深度不够。
4.4 自查顺序：从 FFT 布局一步步排查
遇到重建结果不对，我先查传播函数，再查迭代约束，最后查量化。传播函数最怕 FFT 布局错位。最少验证方式是让z=0，角谱函数必须原样返回输入：
U0 = zeros(4, 4); U0(2,2) = 1; U1 = ang_spectrum(U0, 1e-6, 1e-6, 532e-9, 0); assert(max(abs(U1(:) - U0(:))) < 1e-12, 'z=0 自检未通过');这个断言通过了，说明fft2、ifft2、ifftshift的串联没问题。然后我会在 GS 里临时把iter设成 1，观察重建图是变亮还是变花。如果第一次迭代就出现明显重影，多半是传播方向或坐标问题；如果迭代后完全看不出目标，再看target是否被误当振幅输入。最后才查量化。这个顺序能省掉至少一半的排错时间，比我早期到处打印矩阵靠谱得多。
5. 进阶用法：多平面、动态全息与 SLM 量化的工程处理
5.1 多平面全息：多深度交替投影的迭代写法
单平面 GS 跑稳之后，很多人会想做一个“多个深度同时成像是”的效果。常见做法是交替投影：同一个相位图分别传播到不同距离，在每个距离上都用该平面的目标振幅去约束，转完一圈再进入下一轮迭代。它的本质还是 GS，只是把误差分散到多个平面。下面是一个最小多深度版本：
function phi = gs_multi_plane(targets, ds, dx, dy, lambda, iter) % targets : H x W x P 的强度数组，P 为深度数 % ds : 每个深度对应的传播距离 [Ny, Nx, P] = size(targets); phi = 2 * pi * rand(Ny, Nx); amps = sqrt(max(targets, 0)); for k = 1:iter for p = 1:P img = ang_spectrum(exp(1j*phi), dx, dy, lambda, ds(p)); img = amps(:,:,p) .* exp(1j * angle(img)); back = ang_spectrum(img, dx, dy, lambda, -ds(p)); phi = angle(back); end phi = mod(phi, 2*pi); end end这个写法有两个直觉上的坑。第一个是平面数量别贪多，2 到 3 个深度还能看，5 个以上基本散成噪声，因为一个纯相位面能携带的信息有限。第二个是距离差不能太小，两个深度相差小于 1 cm 时，传播函数对相位的区分度不够，迭代会在两个面之间震荡，最后两个面都模糊。多平面 CGH 的指标应该是“每个平面的重建相关性”，而不是总能量，单独一个平面的信噪比会被其他平面分走。
5.2 动态全息：预生成序列比实时迭代更可控
实时动态 CGH 最容易踩的坑是想在播放时每帧跑 40 轮 GS。MATLAB 里这样做的帧率基本不可用，而且每帧初相位互相独立，画面会出现明显闪烁。工程上有两条路：要么预生成整段相位序列存成视频或 uint8 数组，运行时只切换帧；要么用上一帧相位作为当前帧初值，迭代轮数降到 8 到 10，用时间相关性换稳定性。
% 连续帧生成：上一帧相位作为下一帧初值，能明显抑制闪烁 prev_phi = gs_kinoform(frame1, cfg.dx, cfg.dy, cfg.lambda, cfg.z, 40, 'rand'); for i = 2:num_frames phi_i = gs_kinoform(frame_i, cfg.dx, cfg.dy, cfg.lambda, cfg.z, 10, ... 'rand'); % 实际函数需扩展 init_phase 参数 prev_phi = phi_i; end5.3 量化与 SLM 匹配：8 位相位延迟的补偿习惯
SLM 的相位延迟和灰度之间几乎都不是线性关系。同一个灰度值，在不同波长下的实际相位延迟也不同。函数库里量化函数只能做“理论映射”，真正上光学平台前必须接标定表。常见做法是厂商提供 LUT 文件，加载后把理论灰阶查成实际灰阶：
% LUT 示例：gray 100 对应实际相位 2.1 rad，不是线性值 lut = load('my_slm_lut.mat'); idx = round(mod(phi, 2*pi) / (2*pi) * 255) + 1; Q = uint8(lut.gray_axis(idx)); % 查表后才是 SLM 需要的灰度没有标定表时，至少保留“相位调制深度可能小于 2π”的假设。把 phi 乘一个系数压到器件实测范围，比直接按 2π 量化更安全。这个习惯能避免在光学平台上看到“怎么调都亮不起来”之后才回头改量化代码。
6. 离线验证法：用一屏数值重建判断这次改动值不值得留
6.1 三个指标：相关性、PSNR 与衍射效率的判定顺序
生成相位图后，我做的第一件事不是上光路，而是用同一个ang_spectrum做数值重建。离线重建能过滤掉一大部分算法问题，成本几乎为零，光路调试则要花掉整个下午。判指标时按顺序来：先看相关性，再看 PSNR，最后看衍射效率。
% 数值重建与指标计算（ROI 内评估） replay = abs(ang_spectrum(exp(1j*phi), cfg.dx, cfg.dy, cfg.lambda, cfg.z)).^2; R = replay / max(replay(:)); T = target / max(target(:)); rho = corr(R(:), T(:)); % 相关性，灰度照片可参考 0.6+ mask = target > 0.05; % ROI mask psnr_val = 10 * log10(1 / mean((R(mask) - T(mask)).^2)); eta = sum(R(mask)) / sum(replay(:)); % ROI 能量占比，参考 0.3~0.5 fprintf('corr=%.3f psnr=%.2f dB eta=%.2f\n', rho, psnr_val, eta);相关性高但 PSNR 低，说明轮廓对、暗部崩，问题常在背景噪声；PSNR 高但相关性一般，说明灰度整体偏移，问题常在归一化；衍射效率突然掉下来，先怀疑量化边界和 SLM 调制深度。这三个指标不完美，但足够在改动参数后给出“留还是不留”的结论。我现在的习惯是，任何新版本先进库跑 z=0 自检，再跑一遍这段数值重建，把 corr 和 eta 记下来；忘了记录就把改动直接覆盖掉，这种返工比丢代码还心疼。希望帮到你。
本文还有配套的精品资源，点击获取
