简介：HoloLib_MATLAB 是一套面向计算全息（CGH）学习与实验的 MATLAB 函数库，适合光学、光电信息及图像处理方向的学生与研究人员使用，用于简化全息图生成、传播与重构的编程流程。压缩包共 128 个文件，约 17.05MB，其中 36 个 .m 脚本构成核心算法，涵盖菲涅尔与傅里叶传播、迭代傅里叶与迭代菲涅尔算法、DOE 元件传播及对应重构函数，并提供轴锥透镜、闪耀光栅、菲涅尔球面透镜等相位生成工具；72 个 png 与 12 个 bmp 图像作为实验素材，另有 license、md 说明及工程配置文件辅助使用。库内还包含路径初始化、边界扩充、矩阵扩展、图像拼接、相位保存与 LCOS 显示等辅助函数，便于快速搭建完整实验链路。目前已有 1285 人学习下载，适合希望跳过底层重复实现、直接验证全息算法效果的读者参考。

1. 从一堆散落的 CGH 脚本说起：HoloLib_MATLAB 到底补的是什么

做计算全息的人大概都有过这种体验：电脑里躺着一个叫cgh_test的文件夹，里面是fresnel.m、angular_spectrum.m、GS.m、plot_phase.m七八个脚本，每个脚本开头都复制粘贴了同一段生成坐标网格的代码，改一个参数就得挨个文件翻。HoloLib_MATLAB 想解决的正是这件事——它把计算机生成全息图（CGH）里反复用到的函数收拢成一个 MATLAB 函数库，让「生成一张全息图」从散装脚本变成可复用的函数调用。这个库适合两类人：一是刚接触 CGH、想跑通菲涅尔和角谱两条基础链路的同学；二是已经在写全息算法、但受够了重复造轮子的开发者。系统开源的好处是你能直接看到每个函数的输入输出约定，而不是对着黑匣子调参。

2. 技术选型与目录结构：为什么是 MATLAB，库该怎么摆

2.1 选 MATLAB 而不是 Python 的现实理由

CGH 的日常工作是矩阵运算加可视化：一次传播就是一次 FFT，一次相位编码就是一次复数取角，中间还要不停imagesc看相位包裹成什么样。MATLAB 在这三件事上的手感是现成的——fft2、fftshift、angle、imagesc一行到位，不用额外装绘图栈，也不用担心数组下标从 0 还是 1 开始。对一个函数库来说，这意味着每个函数的实现可以短到十几行，读者能一眼看完整个数据流。

代价也要说清楚：MATLAB 的循环慢，所以库里的函数必须写成向量化形式，凡是能用矩阵整体算的绝不写for。这反过来约束了函数设计——所有空间域操作都接收完整二维矩阵，而不是逐点回调。如果你后续想迁到 Python，这套「输入矩阵、输出矩阵、中间不落盘」的约定几乎可以一比一翻译成 NumPy，这也是我建议新手先读函数签名再读实现的原因。

2.2 一个函数库该有的目录骨架

原始项目正文没有给出目录清单，但按这类 CGH 函数库的通行组织方式，合格的结构大致是这样（我按自己拆过的同类库补全，命名保持通用）：

HoloLib_MATLAB/ ├── core/ # 物理传播与变换核心 │ ├── propFresnel.m # 菲涅尔传播 │ ├── propAngularSpec.m # 角谱传播 │ └── fftShift2.m # 二维频谱搬移封装 ├── cgh/ # 全息图生成算法 │ ├── cghGS.m # Gerchberg-Saxton 迭代 │ ├── cghDirect.m # 直接相位编码 │ └── cghLohmann.m # 迂回相位编码 ├── utils/ # 通用工具 │ ├── makeGrid.m # 坐标网格生成 │ ├── padTo.m # 补零到指定尺寸 │ └── normComplex.m # 复振幅归一化 ├── viz/ # 可视化 │ ├── showPhase.m # 相位图（含包裹处理） │ └── showIntensity.m # 强度图 ├── examples/ # 可运行示例脚本 │ ├── demo_fresnel.m │ └── demo_gs.m └── README.md

这个划分的关键在于core和cgh分离：传播只负责「光怎么走」，编码只负责「相位怎么落到像素上」，两者通过复振幅矩阵对接。好处是换传播模型不影响编码算法，换编码算法也不动传播代码。utils里最容易被低估的是makeGrid.m——坐标网格一旦在多个函数里各写一遍，采样间距和中心点定义迟早对不上，这是 CGH 里最隐蔽的一类错误来源。

2.3 函数签名约定：先定接口再写实现

库能不能被复用，八成取决于接口是否统一。我一般会强制两条约定：第一，所有传播函数的签名形如U_out = propXxx(U_in, lambda, z, dx)，其中U_in是复振幅矩阵，lambda波长，z传播距离，dx像素间距；第二，所有函数不读全局变量、不弹窗、不自动画图，画图交给viz层。这样你在脚本里串十个函数也不会被中途弹出的 figure 打断。

function U_out = propAngularSpec(U_in, lambda, z, dx) % 角谱传播：频域传递函数法 % U_in : 输入复振幅 (M x N) % lambda: 波长 (m) % z : 传播距离 (m)，可为负 % dx : 像素间距 (m) % U_out : 输出复振幅 (M x N) [M, N] = size(U_in); % 频率坐标，注意 fftshift 后零频在中心 fx = (-N/2:N/2-1) / (N*dx); fy = (-M/2:M/2-1) / (M*dx); [FX, FY] = meshgrid(fx, fy); % 传递函数，倏逝波直接置零 H = exp(1j * 2*pi*z * sqrt(max(0, 1/lambda^2 - FX.^2 - FY.^2))); U_out = ifft2(ifftshift(fftshift(fft2(U_in)) .* H)); end

这段是整个库的地基，逐行看：fx、fy用(-N/2:N/2-1)/(N*dx)生成，保证零频落在数组中心，和后面fftshift的用法严格配套；max(0, ...)是防止根号下出现负数，对应倏逝波分量，物理上它们不传播，置零是标准处理；最后一行先fft2再fftshift把零频搬到中心，乘完传递函数再ifftshift搬回去做ifft2。参数里z允许为负，负距离就是反向传播，做相位恢复时经常用到。dx的单位必须是米，混用毫米是新手最常见的翻车点，后面避坑章节会专门讲。

3. 核心模块拆解：从坐标网格到一张全息图

3.1 坐标网格与采样：所有误差的源头

CGH 里绝大多数「结果不对」最后都能追到采样上。makeGrid.m的职责是把「物理尺寸、像素数、像素间距」三者关系固定下来，避免每个函数各自推导。

function [X, Y] = makeGrid(M, N, dx) % 生成以中心为原点的物理坐标网格 % M, N : 行数、列数 % dx : 像素间距 (m) % X, Y : 物理坐标矩阵 (m) x = ((1:N) - (N+1)/2) * dx; y = ((1:M) - (M+1)/2) * dx; [X, Y] = meshgrid(x, y); end

(N+1)/2这个写法是为了让偶数尺寸时原点落在两个中心像素之间、奇数尺寸时落在中心像素上，两种尺寸都不偏半个像素。如果你自己写(0:N-1)*dx，原点就跑到了左上角，后面所有基于「中心对称」的相位图都会整体偏移，而且偏移量随尺寸变化，极难排查。参数dx同时决定了空间带宽：采样间距越大，能表示的最高空间频率越低，角谱传播里高频被截断，边缘会出现振铃。

3.2 Gerchberg-Saxton 迭代：相位恢复的主循环

GS 算法是 CGH 里最经典的迭代编码，库里的cghGS.m承担「给定目标振幅，求相位分布」的任务。核心是四步循环：正传播、替换振幅、反传播、保留相位。

function phase = cghGS(targetAmp, lambda, z, dx, nIter) % Gerchberg-Saxton 相位恢复 % targetAmp: 目标面振幅 (M x N)，实数非负 % nIter : 迭代次数 phase = 2*pi*rand(size(targetAmp)); % 随机初相位，打破对称 for k = 1:nIter U_slm = targetAmp .* exp(1j*phase); % 全息面复振幅 U_obj = propAngularSpec(U_slm, lambda, z, dx); % 传到目标面 U_obj = targetAmp .* exp(1j*angle(U_obj)); % 强制目标振幅 U_back = propAngularSpec(U_obj, lambda, -z, dx); % 反向传回 phase = angle(U_back); % 只保留相位 end end

第一行随机初相位不是可有可无的：全零初相位会让迭代卡在对称解上，重建出现中心对称的重影，这是 GS 最典型的失败模式。循环里targetAmp被用了两次——一次在全息面当振幅约束（实际 SLM 上振幅恒定），一次在目标面当目标约束，这正是 GS 的物理假设。nIter一般取 20 到 100，再多收益递减，因为 GS 会陷入局部极小，误差曲线早就平了。想改善就得上带松弛因子的变体，或者换加权 GS，这属于二次开发范畴，第 5 章展开。

3.3 相位可视化：包裹相位怎么看才不误判

showPhase.m看着简单，但它是调试时最常被打开的函数。相位是-pi到pi的包裹值，直接imagesc会看到锯齿状跳变，新手容易误以为算法崩了。

function showPhase(phase, titleStr) % 显示包裹相位，使用 hsv 色图让 2pi 跳变连续 figure; imagesc(angle(exp(1j*phase))); % 强制回到 [-pi, pi] colormap(hsv); colorbar; axis image off; title(titleStr); end

angle(exp(1j*phase))这一步是幂等的：输入已经在[-pi, pi]就原样返回，输入超出范围就自动折回，省得调用方自己判断。hsv色图的首尾颜色接近，包裹跳变处看起来是连续的，比jet更容易判断相位分布是否合理。axis image保证像素是方的，不然相位图的几何比例会被拉伸，判断光斑形状时会看走眼。

3.4 一个最小可运行示例

把上面几个函数串起来，examples/demo_gs.m的骨架大致如下，这也是你拿到库后第一个该跑的脚本：

% demo_gs.m：生成一个方形光斑的全息图并重建 lambda = 532e-9; % 绿光，单位米 dx = 8e-6; % SLM 像素间距 8 微米 M = 512; N = 512; z = 0.3; % 传播 30 厘米 target = zeros(M, N); target(M/2-40:M/2+40, N/2-40:N/2+40) = 1; % 方形目标 phase = cghGS(target, lambda, z, dx, 50); showPhase(phase, 'GS 全息图相位'); % 验证：用求得的相位重建 U = propAngularSpec(exp(1j*phase), lambda, z, dx); showIntensity(abs(U).^2, '重建强度');

跑通这个脚本，你就验证了「网格—传播—迭代—可视化」整条链路。如果重建出来是一个模糊的方斑而不是清晰方斑，先别怀疑算法，去查dx和z的单位是否统一，以及target的尺寸是否和传播函数里的M、N一致——这两条能解决八成问题。

4. 跑起来与改起来：环境、调用与参数调优

4.1 环境准备与路径挂载

这个库对 MATLAB 版本要求不高，R2016b 之后基本都能跑，因为用到的都是fft2、meshgrid、imagesc这类老函数。不需要额外工具箱，这一点很关键——很多 CGH 库依赖 Image Processing Toolbox 的padarray或fftshift变体，一旦缺工具箱就跑不动。把库放进工作目录后，第一件事是把各子目录加进搜索路径：

addpath(genpath('HoloLib_MATLAB')); % 递归加入所有子目录 savepath; % 可选：持久化，下次不用再加

genpath会把examples也加进去，方便直接调示例。如果你在多个项目间切换，建议不要savepath，而是在每个项目脚本开头调一次addpath，避免不同版本的库互相覆盖——这是我在同时维护两个分支时踩过的坑。

4.2 参数怎么调：一张对照表

调参是 CGH 里最耗时间的环节，把关键参数和影响列清楚能省很多试错：

这张表的使用方式是「先固定三个、只动一个」。比如先定dx和z，只调nIter看误差收敛曲线，找到拐点再停；然后固定nIter去调z。同时动多个参数，你永远不知道是哪个起了作用。

4.3 换传播模型：菲涅尔与角谱的取舍

库同时提供propFresnel.m和propAngularSpec.m，不是冗余。角谱是严格解，近场远场都适用，但要求采样满足带宽条件；菲涅尔是傍轴近似，远场更快，但在大角度或极近距离下会失真。判断标准是菲涅尔数：N_F = a^2/(lambda*z)，a是孔径半宽，N_F远大于 1 时菲涅尔够用，接近 1 或更小时老老实实上角谱。我一般默认用角谱，只有在做大规模远场仿真、嫌 FFT 太慢时才切菲涅尔，切换时记得同步检查采样条件，否则结果对不上还以为是代码 bug。

4.4 从示例到自己的场景

示例跑通后，改造顺序建议是：先换target，把方形换成你要的图案或字母；再换lambda和dx匹配你手上的 SLM 参数；最后才动算法。换target时注意它必须是实数非负矩阵，且尺寸和传播函数一致，否则.*会因维度不匹配直接报错。如果你的目标图案是灰度图，记得先归一化到[0,1]，GS 对振幅绝对尺度不敏感，但对负值会直接产生错误相位。

5. 避坑与排查：五条血泪经验

5.1 重建结果整体偏移半个视场

现象：重建光斑不在中心，而是偏到一角，且偏移量随图像尺寸变化。原因：坐标网格用了(0:N-1)*dx，原点在左上角，而传播函数里的频率坐标是按中心为零频生成的，两者不一致。解决：统一用makeGrid.m的((1:N)-(N+1)/2)*dx写法，所有函数只调它，不自己推导坐标。

5.2 相位图出现同心圆环状伪影

现象：showPhase出来的相位图有一圈圈同心环，重建强度中心有亮斑。原因：targetAmp里有零值区域，GS 迭代时这些点的相位是随机的，反向传播后在全息面形成环状干涉。解决：给目标振幅加一个极小的底数（如1e-6），或者在零值区做相位平滑，避免随机相位参与反传。

5.3 换尺寸后 FFT 结果和预期对不上

现象：M、N从 512 改成 500 后，传播结果出现明显不对称。原因：偶数尺寸下fftshift的零频位置和奇数尺寸不同，(-N/2:N/2-1)的写法只对偶数严格成立。解决：要么统一用偶数尺寸，要么在频率坐标生成时用floor处理，并在文档里注明只支持偶数——库的约定比通用性更重要。

5.4 迭代误差不降反升

现象：nIter加到 200，重建质量反而比 50 次差。原因：GS 是单调下降的，但如果你在循环里手动改了相位（比如加了随机扰动做相位多样性），扰动幅度过大会破坏收敛。解决：检查循环体内是否有额外操作，扰动幅度控制在0.1*pi以内，并记录每次迭代的误差，画曲线确认单调性。

5.5 单位混用导致传播距离完全错误

现象：z填了300（以为是毫米），结果光斑糊成一片。原因：库内所有长度单位统一为米，z=300意味着传播 300 米，远超近场范围。解决：在脚本开头用注释标明单位，或者封装一个mm2m的小函数强制转换；我现在的习惯是所有物理量变量名带单位后缀，如z_m、dx_m，从命名上杜绝混用。

6. 二次开发：把 GS 换成带松弛的变体并验证

库跑通之后，最值得动手改的就是 GS 的收敛质量。标准 GS 的误差曲线很快变平，重建有散斑。一个低成本改造是引入松弛因子，把反向传播后的相位和上一轮相位做加权混合，相当于给迭代加一点惯性，能跳出部分局部极小。

function phase = cghGSRelaxed(targetAmp, lambda, z, dx, nIter, alpha) % 带松弛的 GS，alpha 取 0.5~0.9 phase = 2*pi*rand(size(targetAmp)); for k = 1:nIter U_slm = targetAmp .* exp(1j*phase); U_obj = propAngularSpec(U_slm, lambda, z, dx); U_obj = targetAmp .* exp(1j*angle(U_obj)); U_back = propAngularSpec(U_obj, lambda, -z, dx); phaseNew = angle(U_back); % 相位差折回 [-pi,pi] 再加权，避免 2pi 跳变污染混合 dphi = angle(exp(1j*(phaseNew - phase))); phase = phase + alpha * dphi; end end

关键在dphi那一行：直接alpha*phaseNew + (1-alpha)*phase在相位跨越±pi时会把跳变当成真实差异，混合出错误结果；先做angle(exp(1j*...))把差值折回主值区间，再加权，才是正确的相位插值。alpha从 0.5 试起，太大（接近 1）退化成标准 GS，太小收敛慢，0.7 左右通常是甜点。

验证方法不能只看图。我一般同时记录两个指标：目标面振幅误差err = norm(abs(U_obj) - targetAmp, 'fro') / norm(targetAmp, 'fro')，以及全息面相位的光滑度（相邻像素相位差的均值）。前者看重建质量，后者看相位是否可被 SLM 稳定显示——相位跳变太剧烈，实际器件上会出现串扰。把alpha从 0.3 扫到 0.9，各跑 50 次迭代，画两条曲线，选误差最低且相位光滑度可接受的那个点。

再往深走可以换初始化：随机相位换成二次相位（透镜相位），能让 GS 更快收敛到类高斯分布的目标；或者引入加权 GS，在目标面按区域给不同权重，优先保证中心区域质量。这些改动都只动cghGS一个文件，传播和可视化层完全不用碰，这正是第 2 章把core和cgh分开的价值。

从那以后我每次改迭代类算法，都强制先跑一遍误差曲线再动可视化，图好看不代表误差降了，反过来误差降了图一定不会太差。希望这个库能帮你把 CGH 的重复劳动省下来，把时间花在真正需要调的算法上，也希望帮到你。

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