针对数字全息相位展开及再现的MATLAB实现，核心可分为两个主要步骤：相位展开（将包裹的相位恢复为连续相位）和全息再现（重建物体波前）。
结合了相位展开和全息再现的完整MATLAB代码框架
%% 1. 加载或生成模拟全息图与包裹相位
% 假设已有全息图 'hologram' 和包裹相位 'wrappedPhase'
% 若无数据，可用以下模拟示例：
[M, N] = [512, 512]; % 图像尺寸
[X, Y] = meshgrid(1:N, 1:M);
objectPhase = 0.5*((X-N/2).^2 + (Y-M/2).^2)/1000; % 模拟物体相位（抛物面）
hologram = exp(1i*objectPhase); % 理想全息图（无载频）
wrappedPhase = angle(hologram); % 包裹相位 [-π, π]%% 2. 相位展开（以质量引导路径积分法为例）
% 计算相位导数方差作为质量图
[dx, dy] = gradient(wrappedPhase);
qualityMap = sqrt(dx.^2 + dy.^2); % 简单的相位导数质量图
unwrappedPhase = unwrapPhaseQualityGuided(wrappedPhase, qualityMap);% 辅助函数：质量引导相位展开
function unwrapped = unwrapPhaseQualityGuided(wrapped, quality)[M, N] = size(wrapped);unwrapped = zeros(M, N);processed = false(M, N);% 找到最高质量点作为起始种子[~, idx] = max(quality(:));[startY, startX] = ind2sub([M, N], idx);queue = [startX, startY]; % 待处理像素队列while ~isempty(queue)% 处理队列中第一个像素（此处简化，实际需优先队列）x = queue(1, 1); y = queue(1, 2);queue(1, :) = [];processed(y, x) = true;% 四邻域neighbors = [x+1, y; x-1, y; x, y+1; x, y-1];for n = 1:4nx = neighbors(n, 1); ny = neighbors(n, 2);if nx>=1 && nx<=N && ny>=1 && ny<=M && ~processed(ny, nx)% 相位差diff = wrapped(ny, nx) - wrapped(y, x);% 调整到 [-π, π] 并补偿 2π 跳变jump = round(diff / (2*pi));unwrapped(ny, nx) = unwrapped(y, x) + diff - jump*2*pi;queue = [queue; nx, ny];processed(ny, nx) = true;endendend
end%% 3. 全息再现（角谱法）
lambda = 632.8e-9; % 波长 (He-Ne激光)
pixelSize = 3.45e-6; % 相机像素尺寸
z = 0.01; % 再现距离% 生成角谱传递函数
kx = 2*pi*(-N/2:N/2-1)/(N*pixelSize);
ky = 2*pi*(-M/2:M/2-1)/(M*pixelSize);
[KX, KY] = meshgrid(kx, ky);
k = 2*pi/lambda;
% 传递函数（角谱法）
H = exp(1i*z*sqrt(k^2 - KX.^2 - KY.^2));
H(isnan(H)) = 0; % 处理 evanescent wave% 再现过程
objectWave = hologram; % 全息图复振幅（此处为模拟，实际需处理）
recField = fftshift(fft2(objectWave)); % 频域
recField = recField .* H; % 角谱传播
recField = ifft2(ifftshift(recField)); % 回到空域%% 4. 结果显示与分析
figure('Position', [100, 100, 1200, 400]);% 原始包裹相位
subplot(1,4,1);
imagesc(wrappedPhase);
colormap jet; colorbar; title('包裹相位');
axis image;% 展开相位
subplot(1,4,2);
imagesc(unwrappedPhase);
colormap jet; colorbar; title('展开相位');
axis image;% 再现强度
subplot(1,4,3);
recIntensity = abs(recField).^2;
imagesc(recIntensity);
colormap gray; colorbar; title('再现强度');
axis image;% 再现相位
subplot(1,4,4);
recPhase = angle(recField);
imagesc(recPhase);
colormap jet; colorbar; title('再现相位');
axis image;
要点详解
1. 相位展开算法选择
相位展开是数字全息中的核心难题，针对不同噪声水平和数据特性，可选择不同算法：
| 算法类型 | 适用场景 | MATLAB实现建议 |
|---|---|---|
| 路径跟踪法 | 噪声较少、连续相位场 | unwrap函数（一维）、枝切法 |
| 质量引导法 | 噪声中等、有孤立区域 | 使用相位导数方差或相关系数作为质量图 |
| 最小范数法 | 噪声较大、需要全局优化 | 离散余弦变换(DCT)求解泊松方程 |
| 深度学习法 | 复杂噪声、实时处理 | 训练UNet等网络（需大量数据） |
实用建议：对于初学者，可从MATLAB自带的unwrap函数（处理每行/列）开始
2. 全息再现方法比较
数字全息再现主要有三种数值方法：
% 菲涅耳衍射法（近场近似）
H_fresnel = exp(1i*k*z) * exp(1i*k/(2*z)*(X.^2 + Y.^2));% 卷积法（适用于中距离）
H_conv = exp(1i*k*sqrt(z^2 + X.^2 + Y.^2));% 角谱法（最精确，无近似）
H_angular = exp(1i*z*sqrt(k^2 - KX.^2 - KY.^2)); % 如上文代码
角谱法最精确但计算量大，菲涅耳法速度快但需满足近似条件。选择时需权衡精度与效率。
参考代码 计算离散傅里叶变换dft www.3dddown.com/cna/84733.html
3. 处理真实数据的注意事项
-
载频去除：实际离轴全息图有载频条纹，需频谱滤波：
F = fftshift(fft2(hologram)); % 设计滤波器选择+1级频谱 filter = fspecial('gaussian', [M,N], 20); F_filtered = F .* filter; objectSpectrum = circshift(F_filtered, [-carrier_y, -carrier_x]); -
相位解缠后处理：
% 去除倾斜（常见于离轴系统） [X,Y] = meshgrid(1:N, 1:M); A = [X(:), Y(:), ones(M*N,1)]; b = unwrappedPhase(:); coeffs = A\b; flattenedPhase = unwrappedPhase - (coeffs(1)*X + coeffs(2)*Y + coeffs(3));