简介：面向计算全息（CGH）与光学实验研究的MATLAB开源函数库，整合了菲涅尔衍射传播、迭代傅里叶/菲涅尔算法、透镜与轴锥透镜相位生成、图像边界扩展与相位图保存显示等常用操作，适合从事全息显示、数字光学和图像重建研究的学生与工程师，可显著简化相关编程流程。压缩包共128个文件，其中36个m源码为算法核心，72个png与12个bmp提供实验图像与结果示例，另有license、md说明文档等，整体仅17.05MB，目录结构直观，便于按需调用。截至目前已有1283人学习使用。包内从getLens、getAxilens相位构建，到propDFFT、propTFFT传播计算，再到iterFourier、iterFresnel迭代优化和recFourier、recFresnel重构验证的完整链条，并配有showLCOS、savePhase等实验输出工具，可支撑全息图生成、数值重构和实验验证的完整流程，是一套实用的CGH算法工具箱。

1. 代码补全的定位：一个专治散装脚本的 CGH 函数库

碰过衍射相关课题的人大多经历过这样的场景：代码堆里躺着一堆 MATLAB 脚本，同一个菲涅尔衍射算法写了三份，改一个波长参数要在好几个文件里来回翻。散装脚本最坑的地方不是算法难写，而是没有清晰的边界——传播、量化、重建揉在一个循环里，想抽出来复用需要人肉对照论文公式。HoloLib_MATLAB 就是为了补全这一环而整理的函数库：一套面向计算机生成全息图（CGH）的 MATLAB 实现，系统开源的方式，把菲涅尔衍射、GS 迭代、相位量化和数值重建统一成带清晰接口的函数。输入目标灰度图像，输出对应的相位全息图，附带演示脚本可以直接跑通。适合做衍射光学、全息显示和光计算仿真的开发者，也适合所有被 MATLAB 散装脚本拖到深夜的人。

2. 架构与技术选型：为什么 MATLAB 适合做 CGH 计算链路

2.1 为什么选 MATLAB：二维矩阵运算就是 CGH 的母语

CGH 的计算核心是复振幅在空间频域上的传播，每一步都离不开 FFT、矩阵元素级乘法和相位角度运算。MATLAB 的特点就是把这类二维矩阵运算压到了运算符级别：论文里写 U_out = FFT(U_in · H)，代码里几乎可以照抄这句话。相比用 C 调 FFTW 或者用 Python 手拼 numpy 变换，MATLAB 做算法验证时心智负担最小，这也是很多衍射光学课题默认用它的原因。

HoloLib_MATLAB 的定位不是替代研究平台，而是把一类散装脚本整理成结构统一、函数边界清晰的小型库。小型库最忌讳过度设计：函数少、依赖少，一个目录放得下，全部用函数文件组织而不引入类和全局状态。理由很简单——库很小，函数式接口足够表达所有调用关系；同时现场调试的时候，类封装反而会拦住你看中间结果的路径。另一个实际原因是，散装脚本缺少的从来不是算子，而是调用边界。函数文件天然把这个边界画出来了。

至于为什么选择系统开源，而不是打包成黑盒工具，是因为 CGH 的参数高度依赖实验环境：同一套代码，换一个空间光调制器（SLM），像素尺寸和量化位深都要变。不开源的话，使用者只能对着配置文件猜，出了问题也无法定位。把源码摊开，每个函数头部写清楚输入输出和量纲，二次开发的人能按自己的硬件调整，这才是这类库应该有的态度。

2.2 目录结构：core、tools、scripts 的分层

HoloLib_MATLAB/ ├── core/ # 核心算法函数 │ ├── cghGSA.m # GS 迭代相位恢复 │ ├── cghFresnelRef.m # 传递函数法菲涅尔衍射 │ ├── cghFresnelImp.m # 冲激响应法菲涅尔衍射 │ ├── cghQuantize.m # 相位量化 │ └── cghReconstruct.m # 数值重建 ├── tools/ # 公共工具 │ ├── imgLoadDouble.m # 图像读取并归一化到 double [0,1] │ └── metrics.m # 重建质量评估 (corr2 / RMSE) ├── scripts/ # 可执行示例 │ ├── demo_generate.m # 从目标图生成相位全息图 │ └── demo_reconstruct.m # 重建并评估 └── output/ # 运行结果输出目录

core 目录只放算法，函数内部不读文件、不保存图片，所有数据都通过参数传入传出，这样每个函数都能独立测试。tools 目录负责输入输出和评估，imgLoadDouble 专门把 uint8 图像统一转成 double 且范围固定在 [0,1]，避免用户手滑把 0~255 的整型数据直接丢进迭代；metrics 里放的是相关系数和均方根误差的计算，用于量化评估重建质量。scripts 是组合层，只负责把 core 和 tools 按正确顺序串起来，这层存在的意义是展示参数应该怎么填。

core 文件名统一加 cgh 前缀是刻意的。MATLAB 的 addpath 会把整个目录里的所有函数暴露到全局命名空间，如果库里的函数叫 fresnel 或者 gs，大概率会撞上其他工具箱里的同名函数。MATLAB 的同名函数以路径顺序中最靠前的那个为准，这种冲突是黑匣子级的疑难杂症——明明调用的是自己的函数，却跑出别人的结果。加前缀是低成本防碰撞方式，代价只是多敲几个字母。

core 里同时保留 cghFresnelRef 和 cghFresnelImp 两个衍射实现，不是冗余。传递函数法（TF）在近到中距离传播时稳定、实现简单，冲激响应法（IR）在中到远距离更抗振铃。两个都放，是为了让使用者在不同传播距离下做交叉验证：同一组参数跑出的重建像，如果两种方法给出的焦点位置一致，说明传播模型没有选错。

2.3 一条调用链：目标图、相位、量化、重建

target = imgLoadDouble('my_target.png'); % 归一化到 [0,1] phase = cghGSA(target, 50, 532e-9, 3.74e-6, 0.15); % 50 次迭代 phaseQ = cghQuantize(phase, 8); % 8 级量化 recon = cghReconstruct(phaseQ, 532e-9, 3.74e-6, 0.15);

四个函数串起来就是完整的 CGH 主链路。532e-9 是常用绿色激光波长，3.74e-6 是常见 SLM 像素尺寸，0.15 是桌面尺度的传播距离，三个量纲不匹配，重建像会散焦或者出栅格条纹。这套参数里最容易改错的是波长和像素尺寸的单位——论文里通常写纳米和微米，MATLAB 里必须换算成米，否则相位因子直接差六个数量级。

这张表里的映射关系是从这套库的代码注释里提炼出来的，也是我后来调参时反复对照的基准。值得注意的一点：有些 GS 实现会把目标振幅开平方再丢进迭代，HoloLib 的做法是不在调用层让用户猜这一步，所有归一化和振幅约束都在 cghGSA 内部完成，调用链看起来干净，也避免同一套预处理逻辑在多个 demo 里重复出现。

3. 核心模块拆解：从目标图到相位图的四段旅程

3.1 cghGSA：GS 迭代相位恢复

function phase = cghGSA(targetAmp, iterations, lambda, pixelSize, z) % cghGSA: GS 迭代相位恢复 % targetAmp 目标振幅矩阵, 建议归一化到 [0,1] % iterations 迭代次数, 一般 30~80 % lambda 波长, 单位 m % pixelSize SLM 像素尺寸, 单位 m % z 传播距离, 单位 m % phase 输出相位全息图, 范围 [-pi, pi] targetAmp = max(targetAmp, 0); % 振幅约束不能为负 phase = rand(size(targetAmp)) * 2 * pi; % 随机初始相位 for k = 1:iterations u1 = targetAmp .* exp(1i * phase); % 像面复振幅 uHolo = cghFresnelRef(u1, lambda, pixelSize, z); % 正传到全息面 phaseH = angle(uHolo); % 记录全息面相位 uBack = cghFresnelRef(exp(1i * phaseH), ... % 纯相位约束后逆传播 lambda, pixelSize, -z); phase = angle(uBack); % 更新像面相位 end phase = mod(phase + pi, 2 * pi) - pi; % 卷绕到 [-pi, pi] end

GS 算法是一个从目标像到全息面的往返约束过程。循环里四个操作各对应一步：先搭建带当前相位的像面复振幅，正传播到全息面，得到全息面上的相位分布；SLM 只能调制强度均匀的光，所以用纯相位 exp(1i*phaseH) 代替完整复振幅，这个操作叫全息面约束；再逆传播回像面，取相位作为下一轮迭代的输入。来回循环，相位分布逐步逼近最佳编码，迭代次数越多，重建像越接近目标振幅。

参数 iterations 不是越大越好。我在使用中的经验是，50 次左右就能稳定，超过 100 次相位图上会出现细小而密集的斑点，那是迭代噪声的累积，不是精度提升。随机初始相位这一步最关键也最容易被省略。如果初始相位全零或者用常量，GS 很容易收敛到局部极值，重建像会呈现同心圆状的亮纹，这是典型的迭代锁死现象。常见做法是 rand(size(targetAmp)) * 2 * pi 产生均匀分布的初始相位，相当于给每个像素一个独立的起点，让迭代有足够多的搜索方向。

3.2 cghFresnelRef：传递函数法的菲涅尔衍射

function Uout = cghFresnelRef(Uin, lambda, z, pixelSize) % cghFresnelRef: 传递函数法 (TF) 菲涅尔衍射 % Uin 输入复振幅矩阵 % lambda 波长 (m) % z 传播距离 (m), 负值表示逆传播 % pixelSize 像素尺寸 (m) % Uout 输出复振幅 [M, N] = size(Uin); fx = ((0:N-1) - floor(N/2)) / (N * pixelSize); fy = ((0:M-1) - floor(M/2)) / (M * pixelSize); [FX, FY] = meshgrid(fx, fy); H = exp(-1i * pi * lambda * z * (FX.^2 + FY.^2)); H = ifftshift(H); Uout = ifft2(fft2(Uin) .* H); end

这个函数把菲涅尔衍射近似写成频域传递函数形式，传播距离 z 在指数项里，z 取负值时传递函数变成共轭相位，这就是逆传播的数学基础。关键点在频域采样坐标的构造：fx = ((0:N-1) - floor(N/2)) / (N * pixelSize) 生成从负到正的离散频域坐标，配合 ifftshift(H) 把传递函数中心对齐到 fft2 输出的频谱原点。这个对齐关系如果出错，出来的重建像会错位到四个角上，看起来像图像被切成四块重排。

参数上，pixelSize 是空间采样间隔，N * pixelSize 是采样窗口宽度，频域间隔就是它的倒数。窗口越宽，频域越细，高频成分越丰富；窗口越小，重建像的空间带宽积就越有限。与冲激响应法（cghFresnelImp）相比，TF 法在近场和中距离传播时更稳定，实现也简单，代价是对像素尺寸和距离的配比敏感。传播距离 z 过大，传递函数的高频项会导致欠采样振铃；传播距离过小，频域坐标下的二次相位因子接近常量，看不出衍射效果。第一次跑通这个库的时候，建议先只调 z，确认重建像尺寸和散焦方向，再动其他参数。

3.3 cghQuantize：把连续相位变成 SLM 能用的灰度级

function phaseQ = cghQuantize(phase, levels) % cghQuantize: 连续相位量化到 K 级 % phase 输入相位, [-pi, pi] % levels 量化等级, 典型值 2/4/8/16, 对应 SLM 等效灰度级 % phaseQ 量化后的相位, 仍是 [-pi, pi] 范围内的离散值 phaseW = mod(phase + pi, 2 * pi) - pi; % 卷绕到 [-pi, pi) step = 2 * pi / levels; phaseQ = round(phaseW / step) * step; % 就近量化 phaseQ = min(max(phaseQ, -pi + step), pi - step); end

SLM 或者二元光学元件只能调制有限的相位档位，8 级量化就是把 [-pi, pi] 均匀切成 8 份，每个像素落到最近的档位。量化的代价是衍射效率降低、重建像出现零级光和高阶噪声。levels 取 2 时是二元全息图，衍射效率最低，但制造最方便；levels 取 8 或 16 时重建像质量接近未量化，是实验室里最常见的调节范围。

量化前必须先卷绕相位，这一步很容易被当成本是多余的。连续相位直接 round 会在 2π 边界处产生跳变，这些跳变在重建里表现为随机散斑。用 mod(phase + pi, 2*pi) - pi 处理之后，跳变点只在相位本身不连续的位置出现。边界那句 min(max(...)) 是为了防止 round 后正好落在 +pi 的边界值，否则 mod 卷绕会让它弹回 -pi 方向，产生一个孤立的离谱像素。

3.4 cghReconstruct：数值重建验证

function recon = cghReconstruct(phaseH, lambda, pixelSize, z) % cghReconstruct: 纯相位全息图的数值重建 % phaseH 量化后的相位全息图 % lambda 波长 (m) % pixelSize 像素尺寸 (m) % z 重建传播距离 (m), 与生成时一致 % recon 重建振幅图 U = exp(1i * phaseH); % 纯相位输入 recon = abs(cghFresnelRef(U, lambda, pixelSize, -z)); % 逆传播 end

生成时用 cghFresnelRef 加 -z 把像面传播回全息面，重建时要把全息图重新传播到目标像面，所以这里再传一次 -z，等于做了一次完整的逆传播。recon 取 abs 是因为我们关心强度分布，相位信息在探测器上不可见。数值重建和真实光路的对应关系是：全息面上的紫外光照射后，传播到 z 处，探测器里看到的强度分布就是这个 abs 结果的平方，一般比较振幅就足够了。

这个函数的设计刻意保持单薄，它只有逆传播加取模两行，不做任何后处理。原因是后处理一旦放进 core，用户就很难分辨重建像的缺陷来自算法本身还是来自滤波美化。lambda、pixelSize、z 三个参数和生成时完全一致的情况下，重建像的灰度分布应接近原始目标图，偏差来源主要是迭代次数和量化等级，用评估指标就能定位问题在哪个环节。

4. 运行与调试：从一个空目录到一张相位全息图

4.1 环境准备与路径加载

% 把整个库加入搜索路径 addpath(genpath('HoloLib_MATLAB'));

addpath 只加根目录是不够的，因为函数分散在 core 和 tools 两个子目录里，必须用 genpath 递归展开子目录才能一次找齐。换机器运行时，注意库路径里不要带中文和空格，MATLAB 在部分版本下对带空格路径的处理并不稳定，直接报“找不到函数”的错，看起来却一脸无辜。

运行环境方面，核心函数只用了 fft2、ifft2、meshgrid、mod、angle 这些基础函数，没有调用任何专业工具箱，所以对 MATLAB 版本不太挑，常规近几年的版本都能直接跑。唯一要注意的是，如果目标图像分辨率超过 2048×2048，内存偏小的机器在 GS 迭代时会有明显等待，建议先把图像缩放到 1024 左右再跑。缩放时顺便把宽高裁成奇数，这个细节在后面的避坑章节会展开。

4.2 生成相位全息图的主流程

% demo_generate.m target = imgLoadDouble('my_target.png'); % 任意灰度图 phase = cghGSA(target, 50, 532e-9, 3.74e-6, 0.15); phaseQ = cghQuantize(phase, 8); mkdir('output'); save(fullfile('output', 'phase_8level.mat'), 'phaseQ'); imshow(mod(phaseQ + pi, 2 * pi) / (2 * pi), []);

目标图建议选 512×512 或 1024×1024 的灰度图，内容对比度最好明显一些，比如白底黑字的图片，便于肉眼看重建像的轮廓是否清晰。保存之前先检查 phaseQ 的 min 和 max，正常情况应该大致落在 -pi 到 pi 的边界附近；如果两个极值都接近零，说明相位分布过度集中，重建像会缺少高频细节。

最后一行 imshow 很多人会写成直接 imshow(phaseQ)，结果看到全黑或者全白，以为是程序出错了。实际原因只是显示范围问题：phaseQ 的取值范围在 [-pi, pi]，而 imshow 默认按 [0,255] 映射，必须先平移缩放到 [0,1]。相位图里那些不规则的灰度条纹就是全息图的结构本身，而不是噪声。

4.3 数值重建与质量评估

% demo_reconstruct.m load(fullfile('output', 'phase_8level.mat'), 'phaseQ'); recon = cghReconstruct(phaseQ, 532e-9, 3.74e-6, 0.15); corrV = corr2(imgLoadDouble('my_target.png'), recon); fprintf('重建相关系数: %.4f\n', corrV); imshowpair(imgLoadDouble('my_target.png'), recon, 'montage');

corr2 计算原始目标图和重建振幅图之间的相关系数，0.8 以上通常可以视为视觉上可辨认，8 级量化时的典型值在 0.85~0.95 之间，具体取决于迭代次数。相关系数对整体亮度的线性缩放不敏感，对局部噪声敏感，所以如果重建像整体偏暗但形状对，corr2 仍会比较高；反过来，如果出现棋盘格或条纹，corr2 会明显掉到 0.6 以下，这时候该找的是衍射传播的参数问题，而不是亮度问题。

imshowpair 的 montage 模式把两张图并排显示，适合第一眼判断。更严格的做法是叠加显示并把差值图也打出来，差值图里的亮斑位置就是重建误差集中的区域。质量评估本身没有绝对标准，我的习惯是把 corrV 记录在文件名里，比如 phase_8level_corr0.91.mat，这样后面调参时可以快速对比哪个参数组合更优。

4.4 复现失败时的常见信号

第一次跑这套库时最常遇到的现象有三个。重建像出现棋盘格条纹，通常是 pixelSize 和 z 的取值导致频域混叠，先把 z 缩小一个数量级看是否改善；重建像出现同心圆，是 GS 迭代次数不够或初始相位不是随机分布；重建像只有边缘一圈亮，中间发黑，多半是量化等级太低，0 级衍射光占了主导。判断结果好不好，不要只盯着一幅图，把 target、phase、phaseQ、recon 四个变量的尺寸、量纲、数据范围打印出来对照一遍，80% 的玄学问题都出在这些基础值上。

5. 避坑记录：跑 CGH 函数库最容易翻车的五个现场

这套库我跑了很多轮，踩过的坑比想象中多。下面五条不是从文档里抄来的，是每个都真实消耗过时间才定位出来的问题，按现象、原因、解决的顺序写，方便直接对照排查。

5.1 重建像错位成四块，图像被切碎重排

现象：正确生成的相位全息图做逆传播后，重建像的四个角互相调换，图像像被刀切成了四块又重贴到一起。

原因：fft2 输出的频谱原点在左上角，而构造频域坐标 fx、fy 时把零频放在了中心。没有对传递函数做 ifftshift 的话，频域相乘时高、低频分量位置错位，等于把频谱的四个象限旋转了一圈再反变换回来。

解决：在构造传递函数后执行 H = ifftshift(H)，保证传递函数与 fft2 的频谱排列方式一致。初版代码里我漏过这一步，重建像一直是四块拼接图，检查了一天才发现是 ifftshift 的位置问题。这个函数名的语义容易让人忽略，它处理的是数组排列顺序，和 fftshift 正好是相反操作，用错一个符号就会出这种匪夷所思的画面。

5.2 偶数尺寸矩阵导致中心采样偏差

现象：重建像出现细微错位或中心亮斑偏移，同一套参数改用奇数尺寸后一切恢复正常，换了偶数尺寸又复发。

原因：很多教程里写频域坐标用 fx = (-N/2:N/2-1)，这个写法只对偶数 N 成立。当 N 是奇数时，-N/2 不是整数，序列长度凑不齐 N 个点，零频位置也偏了半格。cghFresnelRef 里用 ((0:N-1) - floor(N/2)) 就是为了同时兼容奇偶尺寸，但只要输入图是偶数边长，采样中心依旧落在两个像素之间。

解决：目标图裁剪成奇数尺寸，比如 513×513 或 1025×1025，零填充也是一种办法，但要记得在重建后把边缘裁掉。实验上必须用偶数尺寸时，可以在传播前先对复振幅做奇数尺寸的零填充，传播完再裁回原尺寸，代价是多一次 fft，精度却稳定很多。

5.3 重建像过暗，亮度差一个数量级

现象：生成相位全息图时一切正常，重建像轮廓清晰但亮度很低，和目标图对比像蒙了一层灰。

原因：振幅约束没做归一化，图像数据是 uint8 范围 0~255，直接参与迭代。迭代过程中误差项被大数值放大，GS 很容易收敛到局部低质量解，重建像的对比度被压扁。

解决：统一走 imgLoadDouble 把数据转成 double 并归一化到 [0,1]。验证方法很简单，打印重建像的 max 值，如果不到 0.1，优先检查输入范围而不是算法逻辑。这个坑在接外部图片时特别容易踩，png 读进来通常已经是 double，jpg 却往往是 uint8，同一个脚本换了个图片格式就出问题，排查方向完全跑偏。

5.4 更换 SLM 后重建像散焦

现象：同一套脚本，把 pixelSize 从 3.74μm 换成 8μm 后，重建像变得模糊，像相机对焦不准。

原因：pixelSize 直接影响频域采样间隔，也就是全息图的视角范围。8μm 像素对应的高频截止频率比 3.74μm 小，全息图能编码的视角更窄，而 z 没跟着调整，等效于重建平面离开了焦点位置。

解决：按采样定理先估算可分辨的最大空间频率，再设定传播距离。常见做法是保持 z 和 pixelSize 的乘积不变，换硬件时按比例换算。也可以用 cghFresnelImp 的冲激响应法交叉验证，看两种传播模型给出的焦点位置是否一致，如果不一致，说明参数换算还有问题。

5.5 GS 迭代锁死，相位图出现规则同心圆

现象：迭代超过 200 次仍不收敛，生成的相位图出现一圈圈同心圆，重建像像从望远镜里看月亮。

原因：初始相位是常量，比如全零，GS 收敛到一个局部最优解。全零相位相当于所有像素同相起振，迭代过程被低阶轴对称项主导，最终锁死在同心圆条纹上。

解决：用随机相位初始化，并固定随机种子方便调试对比。运行前先检查 phase 初始矩阵的 std，如果 std 为 0，基本就是锁死状态。随机种子固定后，同一组参数每次跑出的全息图完全一致，调参时才能区分是参数变化还是随机性带来的差异。这个习惯我在后面搭多平面的 wrapper 时也保留了下来。

6. 二次开发：把函数库接进另一个光学系统

这套库最值得改造的点，是让传播函数可替换。cghGSA 内部写死了调用 cghFresnelRef，如果我的光学系统用到了 lensless 成像或者需要更高精度的传播模型，就得改 core 里的代码。更稳的做法是给 cghGSA 加一个可选参数，传入函数句柄来决定用哪种传播模型：

function phase = cghGSA(targetAmp, iterations, lambda, pixelSize, z, propagator) % propagator 为可选参数, 传入函数句柄可切换传播模型 if nargin < 6, propagator = @cghFresnelRef; end ... uHolo = propagator(u1, lambda, pixelSize, z); uBack = propagator(exp(1i * phaseH), lambda, pixelSize, -z);

这样原有的五个参数调用完全不受影响，旧脚本不用改一行。想换冲激响应法就传 @cghFresnelImp，想换镜头像差模型就自己写一个同签名的 m 文件传进来。这种扩展方式在 MATLAB 里最轻量，不需要重构任何现有函数。

另一个常用改造方向是多平面全息图。单平面的 cghGSA 只能优化一个距离处的重建像，实际系统里经常需要让不同目标出现在不同深度，比如近处显示图标、远处显示文字。做法是包一层循环，对每个深度分别做一次 GS，再把得到的相位图合并。合并方式可以直接叠加相位角，也可以按权重做复振幅平均，前者简单但会引入串扰，后者效果好一些但要小心振幅归一化。我一般先按相位角叠加做通，视觉可接受后再换成复振幅平均。验证方法还是那套：重建，算 corr2，对比各深度目标图的相关性，确认目标间的串扰在可接受范围。

这类函数库最大的改造价值不在算法本身，而在调用边界。我后来给这套库加自定义损耗函数时，完全没碰 core 里的任何文件，只在外部包了一层 wrapper。那年拆另一套代码时的教训是：不要一上来就改算法，先跑通对方的 demo，再花十分钟按调用参数逆推每个函数的量纲和约束，全部记到函数头部注释里，然后再动手写自己的 wrapper。从那以后，我每次拿到新的函数库，都强制自己走一遍这个流程——黑盒跑 demo、白盒读核心、红笔标量纲，这套流程帮我在 CGH 函数调试上少走了很多弯路。希望帮到你。

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