matlab GS相位恢复算法代码附加注释

以下是Matlab中的GS相位恢复算法代码，附加注释以便理解：

%% GS相位恢复算法
clc; clear; close all;

%% 生成信号
N = 128; % 信号长度
f1 = 10; % 信号频率
f2 = 30;
t = (0:N-1) / N;
s1 = sin(2*pi*f1*t);
s2 = sin(2*pi*f2*t);
s = s1 + s2;

%% 信号FFT
S = fft(s);

%% 相位恢复
phi = angle(S); % 原始相位
phi_est = zeros(size(phi)); % 估计相位
max_iter = 100; % 最大迭代次数
tol = 1e-5; % 收敛精度
for k = 1:max_iter
    for n = 1:N
        phi_est(n) = phi_est(n) + angle(S(n)) - angle(exp(1i*phi_est(n)));
    end
    if norm(phi_est - phi) < tol % 判断是否收敛
        break;
    end
end

%% 信号重构
S_est = abs(S) .* exp(1i*phi_est);
s_est = ifft(S_est);

%% 绘图
subplot(3,1,1);
plot(t, s);
title('原始信号');
xlabel('时间（s）');
ylabel('幅值');

subplot(3,1,2);
plot(t, s_est);
title('相位恢复后的信号');
xlabel('时间（s）');
ylabel('幅值');

subplot(3,1,3);
plot(t, phi, '-r', t, phi_est, '-b');
title('相位恢复');
xlabel('时间（s）');
ylabel('相位');
legend('原始相位', '估计相位');

注释已经标注在代码中，其中最重要的部分是相位恢复的循环部分，即：

for k = 1:max_iter
    for n = 1:N
        phi_est(n) = phi_est(n) + angle(S(n)) - angle(exp(1i*phi_est(n)));
    end
    if norm(phi_est - phi) < tol % 判断是否收敛
        break;
    end
end

这里使用了GS算法迭代计算相位，其中phi_est(n)表示第n个采样点的相位估计值，angle(S(n))表示原始FFT结果中第n个采样点的相位，angle(exp(1i*phi_est(n)))表示估计相位的指数项，exp(1i*phi_est(n))表示将估计相位转换成指数形式。每次迭代都更新估计相位值，直到满足收敛条件为止。最后，通过重构信号的FFT值和估计相位值计算出重构后的信号s_est。