如何用最新版MATLAB解一元二次方程

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('BJUT3.19.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 = 27e-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为什么这个代码别人电脑可以运行出来，我的电脑运行不出来，完全一样的代码，输入数据也一样，是哪里有问题