2026/9/5 12:59:52
FEATURED · 精选文章 
简介本资源是一套面向光学工程、信息光学及计算成像方向初学者与高校实验教学的数字全息仿真实验MATLAB实现方案聚焦数字全息图生成、零级像抑制、波前再现等核心原理的编程验证。压缩包共2个文件1个BMP格式原始全息图数据、1个holographic.m主程序脚本总大小3.42MB结构精炼便于快速运行与代码剖析其中MATLAB脚本完整涵盖图像读取、傅里叶变换、空域高斯滤波去零级、逆变换及衍射再现全流程可直接用于课堂演示或课后复现。已有2423人学习下载适用于《信息光学》《计算全息》课程实验环节帮助学习者打通“光学原理—数值建模—MATLAB实现—图像分析”的完整链路切实掌握从干涉记录到三维物场重建的关键技术细节与调试逻辑。1. 项目概述从理论到屏幕的数字全息之旅数字全息仿真实验听起来像是光学实验室里高深莫测的玩意儿但实际上它正是一扇连接经典光学理论与现代计算成像的绝佳窗口。简单来说这项目就是用MATLAB这把“数字瑞士军刀”在计算机里完整地模拟一套全息记录与再现系统。你不用真的去搭建昂贵且娇贵的光学平台不用担心激光器的稳定性更不用在暗房里小心翼翼地处理全息干板。所有过程从生成模拟的物光波到模拟参考光干涉形成全息图再到最终的数字再现全部在代码和矩阵运算中完成。这解决了什么问题对于学生和研究者它降低了学习全息原理的门槛让你可以直观地、可重复地观察每一个参数变化对最终成像的影响。对于工程师它成为了一个强大的设计验证工具可以在实际搭建光路前预先仿真不同系统配置比如改变光波长、记录距离、探测器像素尺寸的成像效果节省大量时间和成本。无论你是光学工程专业的学生想深入理解《信息光学》课本里的公式还是从事计算成像、显微成像或三维显示研发的工程师需要快速验证一个新想法这个基于MATLAB的数字全息仿真实验都能提供一个清晰、可控且功能强大的沙盘。2. 仿真实验的核心思路与框架设计数字全息仿真的核心思路是对物理全息过程的严格数学建模和离散化计算。整个过程可以清晰地拆解为三个核心阶段对应着三个主要的MATLAB函数模块。2.1 第一阶段模拟物光波前生成全息记录的是物光波的振幅和相位信息。在仿真中我们首先要“创造”一个虚拟的物体及其发出的光波。这里的关键在于如何用数学描述一个复杂的光场。最直接的方法是采用角谱传播理论。我们假设物体是一个二维的透射或反射率分布图例如一个简单的字母“A”的图片。这个分布图可以看作是一个平面上的复振幅分布其中振幅代表物体的透射或反射强度初始相位可以设为零或一个随机相位板用于模拟粗糙表面。然后我们需要计算这个初始平面光场传播一定距离即记录距离d0后到达全息记录平面如CCD靶面的复振幅分布。这个过程通过角谱传播函数实现其本质是求解标量衍射的积分方程在频域里它表现为一个传递函数的乘积运算。注意为什么不直接用菲涅尔衍射或卷积法角谱理论在数学上是最严格的标量衍射近似只要采样满足奈奎斯特频率它对任何距离的传播计算都是准确的避免了菲涅尔近似在极近场时的误差。这对于构建一个基础扎实的仿真框架至关重要。在MATLAB中这意味着我们需要对物体的二维矩阵进行二维快速傅里叶变换2D-FFT乘以一个对应于传播距离的相位传递函数exp函数构成再进行逆傅里叶变换。这个传递函数H是仿真精度的心脏其表达式为H exp(1i*2*pi*d0/lambda * sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2))其中fx,fy是空间频率坐标lambda是光波长。这里就涉及到第一个关键参数选择如何根据模拟的物理尺寸和像素数正确构建这个频率坐标网格。2.2 第二阶段全息图记录干涉模拟得到物光波O(x,y)后我们需要模拟它与参考光R(x,y)的干涉。参考光通常模拟为平面波或球面波。平面波最简单其复振幅可表示为R Ar * exp(1i * 2*pi/lambda * (sin(theta_x)*x sin(theta_y)*y))其中Ar是振幅常设为1theta_x和theta_y是参考光的倾斜角。这个倾斜角引入了载频对于后续的离轴全息分离衍射级至关重要。两者干涉后记录平面的光强分布即为全息图I_hologram abs(O R).^2。这里得到的I_hologram是一个实数值矩阵模拟了CCD相机记录到的强度信息。它丢失了光波的相位但编码了物光波的振幅和相位信息于干涉条纹中。实操心得参考光与物光的光强比IR/IO是一个需要仔细调节的参数。比值太大参考光过强全息图条纹对比度低再现像信噪比差比值太小物光过强可能导致干涉条纹超过探测器的动态范围产生非线性畸变。通常将这个比值设置在3:1到10:1之间进行仿真尝试观察再现效果。2.3 第三阶段数字全息再现这是从全息图中“解压”出物体信息的过程。数字再现的核心是模拟参考光照射全息图后的衍射过程。最常用的方法是菲涅尔衍射法尽管生成用角谱但再现常用菲涅尔近似因为计算更直观。再现过程在数学上表示为U_recon IFFT2( FFT2(I_hologram .* R_conj) .* H_prop )。这里R_conj是模拟再现照明光通常与参考光共轭即conj(R)用于抵消记录时的倾斜相位使像回到中心。H_prop是菲涅尔衍射的传递函数形式为exp(1i*pi/(lambda*d_recon)*(fx.^2fy.^2))其中d_recon是再现距离通常等于记录距离d0。计算得到的U_recon是一个复矩阵其振幅abs(U_recon)即为再现物体的强度像其相位angle(U_recon)包含了物体的三维形貌信息。对于离轴全息在频谱上会存在三个分离的衍射级零级、正负一级我们需要通过频域滤波提取出包含物体信息的那个一级衍射项再进行逆传播以获得清晰的再现像。3. 关键参数解析与MATLAB实现细节一个仿真能否成功、结果是否物理可信完全取决于一系列关键参数的正确设置和匹配。这些参数构成了连接数字世界与物理世界的桥梁。3.1 空间采样与模拟尺度这是最容易出错的地方。在MATLAB中一切都是以像素为单位的离散数组。我们必须为这些像素赋予物理尺寸。像素尺寸delta这模拟的是CCD相机像元的物理大小例如6.45e-6 m6.45微米。它决定了仿真系统的空间截止频率。网格大小Nx, Ny这是图像矩阵的行列数如1024 x 1024。总模拟的物理尺寸为Lx Nx * delta。波长lambda模拟激光的波长如氦氖激光的632.8e-9 m。记录距离d0物体平面到记录平面的距离。这个距离不能随便设必须满足菲涅尔近似或角谱传播的采样条件以避免混叠。它们之间的约束关系由采样定理决定。对于角谱传播需要满足d0 delta * Lx / lambda以避免频域混叠。在编程时我们首先根据lambda、delta和期望的视场Lx来估算最大允许的d0或者先确定d0再反推所需的delta。MATLAB实现时构建坐标网格的代码至关重要lambda 632.8e-9; % 波长 delta 6.45e-6; % 像素尺寸 N 1024; % 像素数 L N * delta; % 总物理尺寸 % 空间坐标 x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 频率坐标 fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy);注意使用meshgrid生成网格并且频率坐标的构建方式这是后续所有FFT运算的基础。3.2 参考光设计与载频控制对于离轴全息参考光倾斜角的选择直接决定了全息图频谱中各级次的分离程度。参考光波矢在x方向的投影为k_x 2*pi/lambda * sin(theta_x)。在频谱上这会使得物光信息即1级的中心从零频点移动到(f_x0, f_y0) (sin(theta_x)/lambda, sin(theta_y)/lambda)。为了在再现时能完美分离出1级需要满足分离条件f_x0必须大于物光频谱的带宽B约等于物体尺寸除以lambda*d0的1.5倍即f_x0 1.5 * B。否则各级频谱会重叠产生串扰。采样条件f_x0 B必须小于奈奎斯特频率1/(2*delta)否则会发生混叠。在MATLAB中我们通过调整theta_x来满足这些条件。通常先估算物光带宽B然后设置f_x0 2 * B左右再反推theta_x asin(lambda * f_x0)。3.3 相位解包裹与像质评价数字全息再现得到的是包裹相位值域在[-π, π]对于测量物体三维形貌需要进行相位解包裹。MATLAB中有unwrap函数但对于噪声大或不连续的相位图需要更稳健的算法如最小二乘法、质量图导引法。仿真中因为数据干净一维或二维的unwrap通常就够用。评价再现像质量除了主观观察常用客观指标均方误差MSE比较再现像振幅与原始物体图像的差异。峰值信噪比PSNR基于MSE计算值越高越好。结构相似性SSIM从亮度、对比度、结构三方面评价图像相似性更符合人眼感知。在仿真中我们可以通过计算这些指标定量分析不同噪声水平、不同参数误差对成像质量的影响。4. 完整MATLAB仿真流程与代码实现下面我们将上述思路整合成一个可运行的、模块化的MATLAB仿真示例。我们将模拟一个简单的方形孔径作为物体进行离轴菲涅尔全息记录与再现。4.1 步骤一初始化参数与创建物体%% 1. 参数初始化 clear; close all; clc; % 物理参数 lambda 632.8e-9; % 波长单位米 (He-Ne激光) k 2 * pi / lambda; % 波数 delta 6.45e-6; % CCD像素尺寸单位米 N 1024; % 像素数 (假设为正方形) L N * delta; % 总模拟尺寸单位米 d0 0.5; % 记录距离单位米 (需满足采样条件) % 参考光参数 (离轴角) theta_x 0.5 * pi / 180; % x方向倾斜角单位弧度 (0.5度) theta_y 0; % y方向无倾斜 Ar 1.0; % 参考光振幅 % 坐标网格 x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 频率坐标 (用于角谱传播) fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy); %% 2. 创建模拟物体 % 生成一个方形孔径 obj_size 2e-3; % 物体尺寸2mm obj double(abs(X) obj_size/2 abs(Y) obj_size/2); % 可以添加相位信息模拟一个倾斜的相位物体 phase_obj 0.5 * pi * X / (obj_size/2); % 线性相位倾斜 U_obj obj .* exp(1i * phase_obj); % 物体平面复振幅 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, abs(U_obj)); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物体振幅分布‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, angle(U_obj)); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物体相位分布包裹‘);这段代码定义了所有核心物理参数并创建了一个带有线性相位变化的方形物体。坐标网格的构建是后续所有运算的基石。注意我们将单位从米转换到毫米进行显示更符合视觉习惯。4.2 步骤二角谱传播与全息图记录%% 3. 角谱传播计算物体到记录平面的光场 % 角谱传递函数 H_as exp(1i * 2*pi*d0/lambda * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 对物体场进行FFT乘以传递函数再IFFT U_obj_fft fft2(fftshift(U_obj)); % 注意fftshift将零频移到中心与我们的频率坐标匹配 U_rec_fft U_obj_fft .* H_as; U_rec ifftshift(ifft2(U_rec_fft)); % ifftshift将结果移回标准顺序 % 物光波在记录平面的振幅 Ao abs(U_rec); % 为了模拟实际情况可以给物光添加一个衰减使其强度弱于参考光 Ao Ao / max(Ao(:)) * 0.3; % 归一化后调整相对强度 %% 4. 生成参考光波并干涉记录全息图 % 生成平面参考光波带有离轴角 R Ar * exp(1i * k * (sin(theta_x)*X sin(theta_y)*Y)); % 记录平面总光场 U_total U_rec .* (Ao ./ abs(U_rec)) R; % 保持U_rec的相位但使用调整后的振幅Ao % 全息图光强分布 I_hologram abs(U_total).^2; % 显示全息图 subplot(1,3,3); imagesc(x*1e3, y*1e3, I_hologram); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘记录的全息图‘);这里有几个关键点fftshift与ifftshift由于我们构建的频率坐标FX,FY是以零频为中心的所以在对空间域信号做FFT前需要用fftshift将信号零频也移到中心与传递函数对齐。运算完成后再用ifftshift移回来。物光强度调整通过Ao Ao / max(Ao(:)) * 0.3将物光峰值振幅设为参考光振幅的0.3倍大致符合IR/IO ≈ 10:1的强度比以获得高对比度干涉条纹。全息图I_hologram是模拟CCD实际采集到的数据它是一个实值矩阵丢失了相位信息但包含了重建所需的一切。4.3 步骤三数字再现与像分离%% 5. 数字全息再现 % 5.1 频域滤波分离衍射级 I_hologram_fft fft2(I_hologram); I_hologram_fft_shifted fftshift(I_hologram_fft); % 将零频移到中心以便观察 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(log(1 abs(I_hologram_fft_shifted))); axis image; colormap(‘jet‘); title(‘全息图频谱对数显示‘); xlabel(‘空间频率 f_x‘); ylabel(‘空间频率 f_y‘); % 可以观察到三个亮斑中心是零级两侧是正负一级。 % 创建滤波器提取1级 [fxx, fyy] meshgrid(1:N, 1:N); % 估算1级中心位置对应参考光载频 f0_x round(N/2 sin(theta_x) * d0 / (lambda * delta)); % 近似计算 f0_y round(N/2); filter_radius 50; % 滤波器半径需小于载频与零级的距离 % 生成圆形带通滤波器 filter_mask double((fxx - f0_x).^2 (fyy - f0_y).^2 filter_radius^2); % 应用滤波器 I_filtered_fft I_hologram_fft .* filter_mask; subplot(1,3,2); imagesc(filter_mask); axis image; title(‘频域滤波器‘); subplot(1,3,3); imagesc(log(1 abs(fftshift(I_filtered_fft)))); axis image; colormap(‘jet‘); title(‘滤波后的频谱1级‘); % 5.2 菲涅尔衍射法再现 % 构建菲涅尔衍射传递函数再现距离为-d0即共轭再现 d_recon -d0; % 负号表示反向传播 H_fresnel exp(1i * pi/(lambda * d_recon) * (FX.^2 FY.^2) * (delta^2 * N^2)); % 注意这里FX,FY是归一化频率需要转换为实际频率并考虑离散采样效应 % 更标准的写法是使用空间坐标构建传递函数 % H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 使用空间坐标构建传递函数更直观 H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 再现过程滤波后的全息图乘以共轭参考光再进行菲涅尔衍射 R_conj conj(R); % 共轭参考光用于消除载频 U_temp ifft2(I_filtered_fft) .* R_conj; % 回到空域并消除倾斜相位 % 菲涅尔衍射积分通过卷积计算先FFT乘传递函数再IFFT U_recon_fft fft2(U_temp) .* fftshift(H_fresnel); % 注意传递函数需要fftshift对齐 U_recon ifft2(U_recon_fft); % 提取再现像的振幅和相位 amp_recon abs(U_recon); phase_recon angle(U_recon); figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, amp_recon); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像振幅‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, phase_recon); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像相位包裹‘);这一步是整个仿真的核心。频域滤波是关键操作滤波器的大小和位置直接影响再现像的质量和分辨率。滤波器半径filter_radius需要足够大以包含全部物体频谱信息但又不能太大以至于包含零级或其他级的成分这需要根据全息图频谱图手动调整或通过算法自动估计。注意事项菲涅尔衍射传递函数H_fresnel的构建有两种常见方式一种在频率域使用FX, FY一种在空间域使用X, Y。两者在数学上等价但离散化计算时要注意坐标缩放因子。使用空间域形式exp(1i*k/(2*d)*(X.^2Y.^2))通常更直观且不易出错但计算量稍大。在仿真中我们更关注正确性因此推荐空间域形式。4.4 步骤四相位解包裹与结果分析%% 6. 相位解包裹与结果分析 % 相位解包裹 (使用MATLAB内置的unwrap对于仿真简单相位通常有效) phase_unwrapped unwrap(phase_recon, [], 1); % 先按行解包裹 phase_unwrapped unwrap(phase_unwrapped, [], 2); % 再按列解包裹 % 去除倾斜相位因为我们模拟的物体相位本身就是倾斜的这里减去一个平面拟合值作为演示 % 实际上这一步在定量相位测量中用于消除系统误差。 subplot(1,3,3); imagesc(x*1e3, y*1e3, phase_unwrapped); axis image; colormap(‘jet‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像相位解包裹后‘); colorbar; %% 7. 像质评价 % 裁剪出中心区域与原始物体进行比较 crop_ratio 0.3; % 裁剪比例 crop_N round(N * crop_ratio); center_idx N/2 (-crop_N/2 : crop_N/2-1); center_idx round(center_idx); amp_original_crop abs(U_obj(center_idx, center_idx)); amp_recon_crop amp_recon(center_idx, center_idx); % 归一化 amp_original_crop amp_original_crop / max(amp_original_crop(:)); amp_recon_crop amp_recon_crop / max(amp_recon_crop(:)); % 计算均方误差(MSE)和峰值信噪比(PSNR) mse mean((amp_original_crop(:) - amp_recon_crop(:)).^2); max_val 1; % 归一化后最大值为1 psnr 10 * log10(max_val^2 / mse); fprintf(‘图像质量评价\n‘); fprintf(‘ 均方误差 (MSE): %.4e\n‘ mse); fprintf(‘ 峰值信噪比 (PSNR): %.2f dB\n‘ psnr); % 显示对比 figure(‘Position‘ [100 100 800 400]); subplot(1,2,1); imagesc(amp_original_crop); axis image; colormap(‘gray‘); title(‘原始物体裁剪后‘); subplot(1,2,2); imagesc(amp_recon_crop); axis image; colormap(‘gray‘); title(‘再现像裁剪后‘); sgtitle(sprintf(‘PSNR %.2f dB‘ psnr));相位解包裹是获取连续相位分布的必要步骤。MATLAB的unwrap函数对仿真生成的、噪声低的相位图效果很好。但在实际实验数据中由于噪声、阴影和相位跳变可能需要更复杂的算法如phase_unwrap工具箱中的算法。像质评价环节让我们能定量评估仿真系统的性能。PSNR值越高说明再现像与原始物体越接近。在理想仿真中无噪声参数完美匹配PSNR可以非常高60 dB。通过引入噪声或参数误差可以观察PSNR如何下降从而理解系统对各因素的敏感度。5. 仿真中的典型问题、调试技巧与进阶应用即使按照上述流程初学者在仿真中仍会遇到各种问题。下面是一些常见“坑”及其排查思路。5.1 问题一再现像一片模糊或出现鬼影可能原因1频谱滤波不彻底零级或共轭像干扰。排查仔细检查全息图的频谱图log(1abs(fftshift(fft2(I_hologram)))。你是否能看到三个明显分离的亮斑如果零级和1级靠得太近说明参考光载频theta_x太小。解决增大参考光倾斜角theta_x重新计算。确保f_x0 1.5 * B。排查检查你应用的频域滤波器。用imagesc(filter_mask)显示滤波器看其位置是否准确覆盖了1级频谱且没有包含零级中心。解决调整滤波器的中心坐标(f0_x, f0_y)和半径filter_radius。可以尝试先手动在频谱图上选取区域。可能原因2再现距离d_recon设置错误。排查再现距离理论上应等于记录距离d0共轭再现。如果使用菲涅尔衍射法尝试微调d_recon的值观察再现像是否变得清晰。可以写一个循环让d_recon在d0附近微小变化寻找图像最清晰的点聚焦。解决使用自动聚焦算法。常用方法是定义一个清晰度评价函数如图像梯度平方和遍历一系列d_recon取函数值最大的距离作为最佳再现距离。5.2 问题二再现像边缘出现周期性条纹或振铃效应可能原因频谱泄露与吉布斯现象。分析当物体是理想的矩形陡峭边缘时其频谱是无限的sinc函数。我们用有限大小的频域滤波器去截断它相当于在空域与一个sinc函数卷积导致边缘出现振荡。解决对原始物体加窗在生成物体U_obj时对其振幅分布乘以一个缓变的窗函数如高斯窗、汉宁窗使边缘平滑过渡。增大滤波器尺寸适当增加filter_radius包含更多高频分量但要注意不要引入其他级的干扰。使用更优的滤波器将圆形二值滤波器改为高斯衰减滤波器即filter_mask exp(-((fxx-f0_x).^2(fyy-f0_y).^2)/(2*sigma^2))可以平滑截断减少振铃。5.3 问题三计算速度慢特别是对大尺寸图像分析角谱传播和菲涅尔衍射涉及大量FFT运算N1024时很快但当N4096或更大时计算和内存消耗显著增加。优化技巧使用单精度如果精度要求可接受将数据转换为单精度single。U_obj single(U_obj);FFT在单精度下更快内存减半。预计算传递函数如果参数不变可以将H_as或H_fresnel计算一次并保存避免在循环中重复计算。利用GPU如果MATLAB安装了Parallel Computing Toolbox且拥有支持CUDA的NVIDIA GPU可以使用gpuArray将数据转移到GPU上计算。FFT在GPU上对大规模数据有巨大加速。U_obj_gpu gpuArray(U_obj); H_as_gpu gpuArray(H_as); U_rec_gpu ifft2(fft2(U_obj_gpu) .* H_as_gpu); U_rec gather(U_rec_gpu); % 将结果取回CPU减少不必要的可视化在调试完成后关闭中间的图形显示 (close all;)或使用set(0,‘DefaultFigureVisible‘,‘off‘)禁止图形弹出可以节省大量时间。5.4 进阶应用引入噪声与像差仿真一个更贴近现实的仿真需要引入噪声和像差。添加噪声模拟CCD读出噪声、散粒噪声等。SNR_dB 20; % 信噪比 I_hologram_noiseless I_hologram; signal_power mean(I_hologram(:).^2); noise_power signal_power / (10^(SNR_dB/10)); noise sqrt(noise_power) * randn(size(I_hologram)); % 高斯白噪声 I_hologram I_hologram_noiseless noise; I_hologram(I_hologram 0) 0; % 确保强度非负通过改变SNR_dB可以研究噪声对再现像质量PSNR的影响。引入像差模拟光学系统的不完美如球差、彗差、像散等。这可以在角谱传递函数H_as或参考光波前R上乘以一个像差相位板W。% 例如引入初级球差 r2 (X.^2 Y.^2) / (L/2)^2; % 归一化半径 W_spherical 2 * pi / lambda * 0.1e-6 * r2.^2; % 0.1微米的球差 H_as_aberrated H_as .* exp(1i * W_spherical);观察像差如何导致再现像模糊、变形从而理解像差校正如数字相位补偿的重要性。通过这个完整的MATLAB数字全息仿真框架你不仅能够复现教科书中的经典现象更能将其作为一个灵活的工具箱用于探索更复杂的全息成像问题如相移全息、彩色全息、显微全息等为真正的光学实验或工程应用打下坚实的理论和实践基础。本文还有配套的精品资源点击获取
MATLAB数字全息仿真：从角谱传播到离轴全息再现的完整实践
简介本资源是一套面向光学工程、信息光学及计算成像方向初学者与高校实验教学的数字全息仿真实验MATLAB实现方案聚焦数字全息图生成、零级像抑制、波前再现等核心原理的编程验证。压缩包共2个文件1个BMP格式原始全息图数据、1个holographic.m主程序脚本总大小3.42MB结构精炼便于快速运行与代码剖析其中MATLAB脚本完整涵盖图像读取、傅里叶变换、空域高斯滤波去零级、逆变换及衍射再现全流程可直接用于课堂演示或课后复现。已有2423人学习下载适用于《信息光学》《计算全息》课程实验环节帮助学习者打通“光学原理—数值建模—MATLAB实现—图像分析”的完整链路切实掌握从干涉记录到三维物场重建的关键技术细节与调试逻辑。1. 项目概述从理论到屏幕的数字全息之旅数字全息仿真实验听起来像是光学实验室里高深莫测的玩意儿但实际上它正是一扇连接经典光学理论与现代计算成像的绝佳窗口。简单来说这项目就是用MATLAB这把“数字瑞士军刀”在计算机里完整地模拟一套全息记录与再现系统。你不用真的去搭建昂贵且娇贵的光学平台不用担心激光器的稳定性更不用在暗房里小心翼翼地处理全息干板。所有过程从生成模拟的物光波到模拟参考光干涉形成全息图再到最终的数字再现全部在代码和矩阵运算中完成。这解决了什么问题对于学生和研究者它降低了学习全息原理的门槛让你可以直观地、可重复地观察每一个参数变化对最终成像的影响。对于工程师它成为了一个强大的设计验证工具可以在实际搭建光路前预先仿真不同系统配置比如改变光波长、记录距离、探测器像素尺寸的成像效果节省大量时间和成本。无论你是光学工程专业的学生想深入理解《信息光学》课本里的公式还是从事计算成像、显微成像或三维显示研发的工程师需要快速验证一个新想法这个基于MATLAB的数字全息仿真实验都能提供一个清晰、可控且功能强大的沙盘。2. 仿真实验的核心思路与框架设计数字全息仿真的核心思路是对物理全息过程的严格数学建模和离散化计算。整个过程可以清晰地拆解为三个核心阶段对应着三个主要的MATLAB函数模块。2.1 第一阶段模拟物光波前生成全息记录的是物光波的振幅和相位信息。在仿真中我们首先要“创造”一个虚拟的物体及其发出的光波。这里的关键在于如何用数学描述一个复杂的光场。最直接的方法是采用角谱传播理论。我们假设物体是一个二维的透射或反射率分布图例如一个简单的字母“A”的图片。这个分布图可以看作是一个平面上的复振幅分布其中振幅代表物体的透射或反射强度初始相位可以设为零或一个随机相位板用于模拟粗糙表面。然后我们需要计算这个初始平面光场传播一定距离即记录距离d0后到达全息记录平面如CCD靶面的复振幅分布。这个过程通过角谱传播函数实现其本质是求解标量衍射的积分方程在频域里它表现为一个传递函数的乘积运算。注意为什么不直接用菲涅尔衍射或卷积法角谱理论在数学上是最严格的标量衍射近似只要采样满足奈奎斯特频率它对任何距离的传播计算都是准确的避免了菲涅尔近似在极近场时的误差。这对于构建一个基础扎实的仿真框架至关重要。在MATLAB中这意味着我们需要对物体的二维矩阵进行二维快速傅里叶变换2D-FFT乘以一个对应于传播距离的相位传递函数exp函数构成再进行逆傅里叶变换。这个传递函数H是仿真精度的心脏其表达式为H exp(1i*2*pi*d0/lambda * sqrt(1 - (lambda*fx).^2 - (lambda*fy).^2))其中fx,fy是空间频率坐标lambda是光波长。这里就涉及到第一个关键参数选择如何根据模拟的物理尺寸和像素数正确构建这个频率坐标网格。2.2 第二阶段全息图记录干涉模拟得到物光波O(x,y)后我们需要模拟它与参考光R(x,y)的干涉。参考光通常模拟为平面波或球面波。平面波最简单其复振幅可表示为R Ar * exp(1i * 2*pi/lambda * (sin(theta_x)*x sin(theta_y)*y))其中Ar是振幅常设为1theta_x和theta_y是参考光的倾斜角。这个倾斜角引入了载频对于后续的离轴全息分离衍射级至关重要。两者干涉后记录平面的光强分布即为全息图I_hologram abs(O R).^2。这里得到的I_hologram是一个实数值矩阵模拟了CCD相机记录到的强度信息。它丢失了光波的相位但编码了物光波的振幅和相位信息于干涉条纹中。实操心得参考光与物光的光强比IR/IO是一个需要仔细调节的参数。比值太大参考光过强全息图条纹对比度低再现像信噪比差比值太小物光过强可能导致干涉条纹超过探测器的动态范围产生非线性畸变。通常将这个比值设置在3:1到10:1之间进行仿真尝试观察再现效果。2.3 第三阶段数字全息再现这是从全息图中“解压”出物体信息的过程。数字再现的核心是模拟参考光照射全息图后的衍射过程。最常用的方法是菲涅尔衍射法尽管生成用角谱但再现常用菲涅尔近似因为计算更直观。再现过程在数学上表示为U_recon IFFT2( FFT2(I_hologram .* R_conj) .* H_prop )。这里R_conj是模拟再现照明光通常与参考光共轭即conj(R)用于抵消记录时的倾斜相位使像回到中心。H_prop是菲涅尔衍射的传递函数形式为exp(1i*pi/(lambda*d_recon)*(fx.^2fy.^2))其中d_recon是再现距离通常等于记录距离d0。计算得到的U_recon是一个复矩阵其振幅abs(U_recon)即为再现物体的强度像其相位angle(U_recon)包含了物体的三维形貌信息。对于离轴全息在频谱上会存在三个分离的衍射级零级、正负一级我们需要通过频域滤波提取出包含物体信息的那个一级衍射项再进行逆传播以获得清晰的再现像。3. 关键参数解析与MATLAB实现细节一个仿真能否成功、结果是否物理可信完全取决于一系列关键参数的正确设置和匹配。这些参数构成了连接数字世界与物理世界的桥梁。3.1 空间采样与模拟尺度这是最容易出错的地方。在MATLAB中一切都是以像素为单位的离散数组。我们必须为这些像素赋予物理尺寸。像素尺寸delta这模拟的是CCD相机像元的物理大小例如6.45e-6 m6.45微米。它决定了仿真系统的空间截止频率。网格大小Nx, Ny这是图像矩阵的行列数如1024 x 1024。总模拟的物理尺寸为Lx Nx * delta。波长lambda模拟激光的波长如氦氖激光的632.8e-9 m。记录距离d0物体平面到记录平面的距离。这个距离不能随便设必须满足菲涅尔近似或角谱传播的采样条件以避免混叠。它们之间的约束关系由采样定理决定。对于角谱传播需要满足d0 delta * Lx / lambda以避免频域混叠。在编程时我们首先根据lambda、delta和期望的视场Lx来估算最大允许的d0或者先确定d0再反推所需的delta。MATLAB实现时构建坐标网格的代码至关重要lambda 632.8e-9; % 波长 delta 6.45e-6; % 像素尺寸 N 1024; % 像素数 L N * delta; % 总物理尺寸 % 空间坐标 x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 频率坐标 fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy);注意使用meshgrid生成网格并且频率坐标的构建方式这是后续所有FFT运算的基础。3.2 参考光设计与载频控制对于离轴全息参考光倾斜角的选择直接决定了全息图频谱中各级次的分离程度。参考光波矢在x方向的投影为k_x 2*pi/lambda * sin(theta_x)。在频谱上这会使得物光信息即1级的中心从零频点移动到(f_x0, f_y0) (sin(theta_x)/lambda, sin(theta_y)/lambda)。为了在再现时能完美分离出1级需要满足分离条件f_x0必须大于物光频谱的带宽B约等于物体尺寸除以lambda*d0的1.5倍即f_x0 1.5 * B。否则各级频谱会重叠产生串扰。采样条件f_x0 B必须小于奈奎斯特频率1/(2*delta)否则会发生混叠。在MATLAB中我们通过调整theta_x来满足这些条件。通常先估算物光带宽B然后设置f_x0 2 * B左右再反推theta_x asin(lambda * f_x0)。3.3 相位解包裹与像质评价数字全息再现得到的是包裹相位值域在[-π, π]对于测量物体三维形貌需要进行相位解包裹。MATLAB中有unwrap函数但对于噪声大或不连续的相位图需要更稳健的算法如最小二乘法、质量图导引法。仿真中因为数据干净一维或二维的unwrap通常就够用。评价再现像质量除了主观观察常用客观指标均方误差MSE比较再现像振幅与原始物体图像的差异。峰值信噪比PSNR基于MSE计算值越高越好。结构相似性SSIM从亮度、对比度、结构三方面评价图像相似性更符合人眼感知。在仿真中我们可以通过计算这些指标定量分析不同噪声水平、不同参数误差对成像质量的影响。4. 完整MATLAB仿真流程与代码实现下面我们将上述思路整合成一个可运行的、模块化的MATLAB仿真示例。我们将模拟一个简单的方形孔径作为物体进行离轴菲涅尔全息记录与再现。4.1 步骤一初始化参数与创建物体%% 1. 参数初始化 clear; close all; clc; % 物理参数 lambda 632.8e-9; % 波长单位米 (He-Ne激光) k 2 * pi / lambda; % 波数 delta 6.45e-6; % CCD像素尺寸单位米 N 1024; % 像素数 (假设为正方形) L N * delta; % 总模拟尺寸单位米 d0 0.5; % 记录距离单位米 (需满足采样条件) % 参考光参数 (离轴角) theta_x 0.5 * pi / 180; % x方向倾斜角单位弧度 (0.5度) theta_y 0; % y方向无倾斜 Ar 1.0; % 参考光振幅 % 坐标网格 x (-N/2 : N/2-1) * delta; y x; [X, Y] meshgrid(x, y); % 频率坐标 (用于角谱传播) fx (-N/2 : N/2-1) / (N*delta); fy fx; [FX, FY] meshgrid(fx, fy); %% 2. 创建模拟物体 % 生成一个方形孔径 obj_size 2e-3; % 物体尺寸2mm obj double(abs(X) obj_size/2 abs(Y) obj_size/2); % 可以添加相位信息模拟一个倾斜的相位物体 phase_obj 0.5 * pi * X / (obj_size/2); % 线性相位倾斜 U_obj obj .* exp(1i * phase_obj); % 物体平面复振幅 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, abs(U_obj)); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物体振幅分布‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, angle(U_obj)); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘物体相位分布包裹‘);这段代码定义了所有核心物理参数并创建了一个带有线性相位变化的方形物体。坐标网格的构建是后续所有运算的基石。注意我们将单位从米转换到毫米进行显示更符合视觉习惯。4.2 步骤二角谱传播与全息图记录%% 3. 角谱传播计算物体到记录平面的光场 % 角谱传递函数 H_as exp(1i * 2*pi*d0/lambda * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 对物体场进行FFT乘以传递函数再IFFT U_obj_fft fft2(fftshift(U_obj)); % 注意fftshift将零频移到中心与我们的频率坐标匹配 U_rec_fft U_obj_fft .* H_as; U_rec ifftshift(ifft2(U_rec_fft)); % ifftshift将结果移回标准顺序 % 物光波在记录平面的振幅 Ao abs(U_rec); % 为了模拟实际情况可以给物光添加一个衰减使其强度弱于参考光 Ao Ao / max(Ao(:)) * 0.3; % 归一化后调整相对强度 %% 4. 生成参考光波并干涉记录全息图 % 生成平面参考光波带有离轴角 R Ar * exp(1i * k * (sin(theta_x)*X sin(theta_y)*Y)); % 记录平面总光场 U_total U_rec .* (Ao ./ abs(U_rec)) R; % 保持U_rec的相位但使用调整后的振幅Ao % 全息图光强分布 I_hologram abs(U_total).^2; % 显示全息图 subplot(1,3,3); imagesc(x*1e3, y*1e3, I_hologram); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘记录的全息图‘);这里有几个关键点fftshift与ifftshift由于我们构建的频率坐标FX,FY是以零频为中心的所以在对空间域信号做FFT前需要用fftshift将信号零频也移到中心与传递函数对齐。运算完成后再用ifftshift移回来。物光强度调整通过Ao Ao / max(Ao(:)) * 0.3将物光峰值振幅设为参考光振幅的0.3倍大致符合IR/IO ≈ 10:1的强度比以获得高对比度干涉条纹。全息图I_hologram是模拟CCD实际采集到的数据它是一个实值矩阵丢失了相位信息但包含了重建所需的一切。4.3 步骤三数字再现与像分离%% 5. 数字全息再现 % 5.1 频域滤波分离衍射级 I_hologram_fft fft2(I_hologram); I_hologram_fft_shifted fftshift(I_hologram_fft); % 将零频移到中心以便观察 figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(log(1 abs(I_hologram_fft_shifted))); axis image; colormap(‘jet‘); title(‘全息图频谱对数显示‘); xlabel(‘空间频率 f_x‘); ylabel(‘空间频率 f_y‘); % 可以观察到三个亮斑中心是零级两侧是正负一级。 % 创建滤波器提取1级 [fxx, fyy] meshgrid(1:N, 1:N); % 估算1级中心位置对应参考光载频 f0_x round(N/2 sin(theta_x) * d0 / (lambda * delta)); % 近似计算 f0_y round(N/2); filter_radius 50; % 滤波器半径需小于载频与零级的距离 % 生成圆形带通滤波器 filter_mask double((fxx - f0_x).^2 (fyy - f0_y).^2 filter_radius^2); % 应用滤波器 I_filtered_fft I_hologram_fft .* filter_mask; subplot(1,3,2); imagesc(filter_mask); axis image; title(‘频域滤波器‘); subplot(1,3,3); imagesc(log(1 abs(fftshift(I_filtered_fft)))); axis image; colormap(‘jet‘); title(‘滤波后的频谱1级‘); % 5.2 菲涅尔衍射法再现 % 构建菲涅尔衍射传递函数再现距离为-d0即共轭再现 d_recon -d0; % 负号表示反向传播 H_fresnel exp(1i * pi/(lambda * d_recon) * (FX.^2 FY.^2) * (delta^2 * N^2)); % 注意这里FX,FY是归一化频率需要转换为实际频率并考虑离散采样效应 % 更标准的写法是使用空间坐标构建传递函数 % H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 使用空间坐标构建传递函数更直观 H_fresnel exp(1i * k/(2*d_recon) * (X.^2 Y.^2)); % 再现过程滤波后的全息图乘以共轭参考光再进行菲涅尔衍射 R_conj conj(R); % 共轭参考光用于消除载频 U_temp ifft2(I_filtered_fft) .* R_conj; % 回到空域并消除倾斜相位 % 菲涅尔衍射积分通过卷积计算先FFT乘传递函数再IFFT U_recon_fft fft2(U_temp) .* fftshift(H_fresnel); % 注意传递函数需要fftshift对齐 U_recon ifft2(U_recon_fft); % 提取再现像的振幅和相位 amp_recon abs(U_recon); phase_recon angle(U_recon); figure(‘Position‘ [100 100 1200 400]); subplot(1,3,1); imagesc(x*1e3, y*1e3, amp_recon); axis image; colormap(‘gray‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像振幅‘); subplot(1,3,2); imagesc(x*1e3, y*1e3, phase_recon); axis image; colormap(‘hsv‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像相位包裹‘);这一步是整个仿真的核心。频域滤波是关键操作滤波器的大小和位置直接影响再现像的质量和分辨率。滤波器半径filter_radius需要足够大以包含全部物体频谱信息但又不能太大以至于包含零级或其他级的成分这需要根据全息图频谱图手动调整或通过算法自动估计。注意事项菲涅尔衍射传递函数H_fresnel的构建有两种常见方式一种在频率域使用FX, FY一种在空间域使用X, Y。两者在数学上等价但离散化计算时要注意坐标缩放因子。使用空间域形式exp(1i*k/(2*d)*(X.^2Y.^2))通常更直观且不易出错但计算量稍大。在仿真中我们更关注正确性因此推荐空间域形式。4.4 步骤四相位解包裹与结果分析%% 6. 相位解包裹与结果分析 % 相位解包裹 (使用MATLAB内置的unwrap对于仿真简单相位通常有效) phase_unwrapped unwrap(phase_recon, [], 1); % 先按行解包裹 phase_unwrapped unwrap(phase_unwrapped, [], 2); % 再按列解包裹 % 去除倾斜相位因为我们模拟的物体相位本身就是倾斜的这里减去一个平面拟合值作为演示 % 实际上这一步在定量相位测量中用于消除系统误差。 subplot(1,3,3); imagesc(x*1e3, y*1e3, phase_unwrapped); axis image; colormap(‘jet‘); xlabel(‘x (mm)‘); ylabel(‘y (mm)‘); title(‘再现像相位解包裹后‘); colorbar; %% 7. 像质评价 % 裁剪出中心区域与原始物体进行比较 crop_ratio 0.3; % 裁剪比例 crop_N round(N * crop_ratio); center_idx N/2 (-crop_N/2 : crop_N/2-1); center_idx round(center_idx); amp_original_crop abs(U_obj(center_idx, center_idx)); amp_recon_crop amp_recon(center_idx, center_idx); % 归一化 amp_original_crop amp_original_crop / max(amp_original_crop(:)); amp_recon_crop amp_recon_crop / max(amp_recon_crop(:)); % 计算均方误差(MSE)和峰值信噪比(PSNR) mse mean((amp_original_crop(:) - amp_recon_crop(:)).^2); max_val 1; % 归一化后最大值为1 psnr 10 * log10(max_val^2 / mse); fprintf(‘图像质量评价\n‘); fprintf(‘ 均方误差 (MSE): %.4e\n‘ mse); fprintf(‘ 峰值信噪比 (PSNR): %.2f dB\n‘ psnr); % 显示对比 figure(‘Position‘ [100 100 800 400]); subplot(1,2,1); imagesc(amp_original_crop); axis image; colormap(‘gray‘); title(‘原始物体裁剪后‘); subplot(1,2,2); imagesc(amp_recon_crop); axis image; colormap(‘gray‘); title(‘再现像裁剪后‘); sgtitle(sprintf(‘PSNR %.2f dB‘ psnr));相位解包裹是获取连续相位分布的必要步骤。MATLAB的unwrap函数对仿真生成的、噪声低的相位图效果很好。但在实际实验数据中由于噪声、阴影和相位跳变可能需要更复杂的算法如phase_unwrap工具箱中的算法。像质评价环节让我们能定量评估仿真系统的性能。PSNR值越高说明再现像与原始物体越接近。在理想仿真中无噪声参数完美匹配PSNR可以非常高60 dB。通过引入噪声或参数误差可以观察PSNR如何下降从而理解系统对各因素的敏感度。5. 仿真中的典型问题、调试技巧与进阶应用即使按照上述流程初学者在仿真中仍会遇到各种问题。下面是一些常见“坑”及其排查思路。5.1 问题一再现像一片模糊或出现鬼影可能原因1频谱滤波不彻底零级或共轭像干扰。排查仔细检查全息图的频谱图log(1abs(fftshift(fft2(I_hologram)))。你是否能看到三个明显分离的亮斑如果零级和1级靠得太近说明参考光载频theta_x太小。解决增大参考光倾斜角theta_x重新计算。确保f_x0 1.5 * B。排查检查你应用的频域滤波器。用imagesc(filter_mask)显示滤波器看其位置是否准确覆盖了1级频谱且没有包含零级中心。解决调整滤波器的中心坐标(f0_x, f0_y)和半径filter_radius。可以尝试先手动在频谱图上选取区域。可能原因2再现距离d_recon设置错误。排查再现距离理论上应等于记录距离d0共轭再现。如果使用菲涅尔衍射法尝试微调d_recon的值观察再现像是否变得清晰。可以写一个循环让d_recon在d0附近微小变化寻找图像最清晰的点聚焦。解决使用自动聚焦算法。常用方法是定义一个清晰度评价函数如图像梯度平方和遍历一系列d_recon取函数值最大的距离作为最佳再现距离。5.2 问题二再现像边缘出现周期性条纹或振铃效应可能原因频谱泄露与吉布斯现象。分析当物体是理想的矩形陡峭边缘时其频谱是无限的sinc函数。我们用有限大小的频域滤波器去截断它相当于在空域与一个sinc函数卷积导致边缘出现振荡。解决对原始物体加窗在生成物体U_obj时对其振幅分布乘以一个缓变的窗函数如高斯窗、汉宁窗使边缘平滑过渡。增大滤波器尺寸适当增加filter_radius包含更多高频分量但要注意不要引入其他级的干扰。使用更优的滤波器将圆形二值滤波器改为高斯衰减滤波器即filter_mask exp(-((fxx-f0_x).^2(fyy-f0_y).^2)/(2*sigma^2))可以平滑截断减少振铃。5.3 问题三计算速度慢特别是对大尺寸图像分析角谱传播和菲涅尔衍射涉及大量FFT运算N1024时很快但当N4096或更大时计算和内存消耗显著增加。优化技巧使用单精度如果精度要求可接受将数据转换为单精度single。U_obj single(U_obj);FFT在单精度下更快内存减半。预计算传递函数如果参数不变可以将H_as或H_fresnel计算一次并保存避免在循环中重复计算。利用GPU如果MATLAB安装了Parallel Computing Toolbox且拥有支持CUDA的NVIDIA GPU可以使用gpuArray将数据转移到GPU上计算。FFT在GPU上对大规模数据有巨大加速。U_obj_gpu gpuArray(U_obj); H_as_gpu gpuArray(H_as); U_rec_gpu ifft2(fft2(U_obj_gpu) .* H_as_gpu); U_rec gather(U_rec_gpu); % 将结果取回CPU减少不必要的可视化在调试完成后关闭中间的图形显示 (close all;)或使用set(0,‘DefaultFigureVisible‘,‘off‘)禁止图形弹出可以节省大量时间。5.4 进阶应用引入噪声与像差仿真一个更贴近现实的仿真需要引入噪声和像差。添加噪声模拟CCD读出噪声、散粒噪声等。SNR_dB 20; % 信噪比 I_hologram_noiseless I_hologram; signal_power mean(I_hologram(:).^2); noise_power signal_power / (10^(SNR_dB/10)); noise sqrt(noise_power) * randn(size(I_hologram)); % 高斯白噪声 I_hologram I_hologram_noiseless noise; I_hologram(I_hologram 0) 0; % 确保强度非负通过改变SNR_dB可以研究噪声对再现像质量PSNR的影响。引入像差模拟光学系统的不完美如球差、彗差、像散等。这可以在角谱传递函数H_as或参考光波前R上乘以一个像差相位板W。% 例如引入初级球差 r2 (X.^2 Y.^2) / (L/2)^2; % 归一化半径 W_spherical 2 * pi / lambda * 0.1e-6 * r2.^2; % 0.1微米的球差 H_as_aberrated H_as .* exp(1i * W_spherical);观察像差如何导致再现像模糊、变形从而理解像差校正如数字相位补偿的重要性。通过这个完整的MATLAB数字全息仿真框架你不仅能够复现教科书中的经典现象更能将其作为一个灵活的工具箱用于探索更复杂的全息成像问题如相移全息、彩色全息、显微全息等为真正的光学实验或工程应用打下坚实的理论和实践基础。本文还有配套的精品资源点击获取 
