```
clc; clear; close all;
% 参数设置
M = 1024; N = 1024;
file = '1.mat'; % 输入你的全息图
data = load(file);
fn = fieldnames(data);
holo= double(data.(fn{1}));
holo_amp = sqrt(holo / max(holo(:))); % 归一化幅值
amp_gt = im2double(imread('BJUT.png'));% 原始复振幅的“幅值图”
amp_gt = imresize(amp_gt, [1024, 1024]);
amp_gt = rgb2gray(amp_gt);
amp_gt = amp_gt / max(amp_gt(:));%归一化幅值
lambda = 532e-9; % 波长
dx = 2.9e-6; % 像素尺寸
z = 17e-3; % 传播距离
N_iter = 100; % 迭代次数
% 相位复原
[U_obj, recon_amp, recon_phase] = phase_recovery(holo_amp, lambda, dx, z, N_iter, amp_gt);
% 显示结果
figure;
subplot(1,2,1); imshow(recon_amp, []); title('复原幅值');
subplot(1,2,2); imshow(recon_phase, []); title('复原相位');
% 评估指标
recon_amp_norm = recon_amp / max(recon_amp(:));
mse = mean((amp_gt(:) - recon_amp_norm(:)).^2);
rmse_val = sqrt(mse);
psnr_val = 10 * log10(1 / mse);
fprintf('[评估结果] PSNR = %.2f dB,RMSE = %.4f\n', psnr_val, rmse_val);
%% === 角谱传播函数 ===
function U_out = angular_spectrum(U_in, z, wavelength, dx)
[M, N] = size(U_in);
k = 2 * pi / wavelength;
fx = (-N/2:N/2-1) / (N * dx);
fy = (-M/2:M/2-1) / (M * dx);
[FX, FY] = meshgrid(fx, fy);
H = exp(1j * k * z * sqrt(1 - (wavelength * FX).^2 - (wavelength * FY).^2));
H((wavelength * FX).^2 + (wavelength * FY).^2 > 1) = 0;
U_out = ifftshift(ifft2(fft2(fftshift(U_in)) .* fftshift(H)));
end
%% === 相位复原主函数 ===
function [U_obj, recon_amp, recon_phase] = phase_recovery(holo_amp, wavelength, dx, z, num_iter, amp_gt)
[M, N] = size(holo_amp);
H = holo_amp; % 初始记录面幅值
phase_est = zeros(M, N); % 初始相位设为 0
U_rec = H .* exp(1j * phase_est); % 初始复振幅
for iter = 1:num_iter
U_obj = angular_spectrum(U_rec, -z, wavelength, dx); % 传播到物面
A = abs(U_obj);
phi = angle(U_obj);
A = min(A, 1); % 正吸收约束
U_obj = A .* exp(1j * phi);
% TV 去噪(Split Bregman)
U_obj = TV_denoise(U_obj, 0.02, 5); % 可调参数
% 传播回记录面并施加幅值约束
U_rec = angular_spectrum(U_obj, z, wavelength, dx);
U_rec = H .* exp(1j * angle(U_rec));
% 可视化中间迭代
if mod(iter, 20) == 0
figure(2); imshow(abs(U_obj), []); title(['Iter ' num2str(iter)]); drawnow;
end
end
% 最后再传播一次获得最终复原图
U_obj = angular_spectrum(U_rec, -z, wavelength, dx);
recon_amp = abs(U_obj);
recon_phase = angle(U_obj);
end
%% === Split Bregman TV 去噪
% ===
function U_out = TV_denoise(U_in, lambda, iter_tv)%lambda:TV正则项系数
[m, n] = size(U_in);%初始化
U_out = U_in;
dx = zeros(m, n);
dy = zeros(m, n);
bx = zeros(m, n);
by = zeros(m, n);
mu = 1;
for k = 1:iter_tv
[Ux, Uy] = gradient(U_out);%计算梯度项
dx = shrink(Ux + bx, lambda / mu);%更新 d(软阈值)
dy = shrink(Uy + by, lambda / mu);
div_d_b = divergence(dx - bx, dy - by);%更新 U(解一个波动方程)
U_out = solve_poisson(U_in + mu * div_d_b, mu);
[Ux, Uy] = gradient(U_out);%更新 Bregman 参数
bx = bx + (Ux - dx);
by = by + (Uy - dy);
end
end
function out = shrink(x, t)%软阈值函数,适用于复数
out = max(abs(x) - t, 0) ./ (abs(x) + 1e-8) .* x;
end
function div = divergence(px, py)%计算散度
[m, n] = size(px);
div = zeros(m, n);
div(:,1:end-1) = px(:,1:end-1) - px(:,2:end);
div(:,end) = px(:,end);
div(1:end-1,:) = div(1:end-1,:) + py(1:end-1,:) - py(2:end,:);
div(end,:) = div(end,:) + py(end,:);
end
function U = solve_poisson(rhs, mu)%使用频域方法解泊松方程
[m, n] = size(rhs);
[fx, fy] = meshgrid(0:n-1, 0:m-1);
fx = fx / n;
fy = fy / m;
denom = 1 + mu * ((2 - 2 * cos(2 * pi * fx)) + (2 - 2 * cos(2 * pi * fy)));
U = real(ifft2(fft2(rhs) ./ denom));
end```请帮忙改善一下这个代码