函数库的设计与实现)
简介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) - (N1)/2) * dx; y ((1:M) - (M1)/2) * dx; [X, Y] meshgrid(x, y); end(N1)/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); endangle(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/240, N/2-40:N/240) 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像素间距4e-6 ~ 10e-6 m视场变大最高空间频率降低分辨率提高易出现混叠z传播距离0.1 ~ 1 m目标面光斑变大、边缘振铃明显近场近似失效菲涅尔误差增大nIter迭代次数20 ~ 100误差下降但收益递减重建粗糙重影明显M、N尺寸256 ~ 1024计算量按平方增长采样不足高频信息丢失这张表的使用方式是「先固定三个、只动一个」。比如先定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)-(N1)/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以为是毫米结果光斑糊成一片。原因库内所有长度单位统一为米z300意味着传播 300 米远超近场范围。解决在脚本开头用注释标明单位或者封装一个mm2m的小函数强制转换我现在的习惯是所有物理量变量名带单位后缀如z_m、dx_m从命名上杜绝混用。6. 二次开发把 GS 换成带松弛的变体并验证库跑通之后最值得动手改的就是 GS 的收敛质量。标准 GS 的误差曲线很快变平重建有散斑。一个低成本改造是引入松弛因子把反向传播后的相位和上一轮相位做加权混合相当于给迭代加一点惯性能跳出部分局部极小。function phase cghGSRelaxed(targetAmp, lambda, z, dx, nIter, alpha) % 带松弛的 GSalpha 取 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 的重复劳动省下来把时间花在真正需要调的算法上也希望帮到你。本文还有配套的精品资源点击获取