标准GS相位恢复算法详解与MATLAB实现

Gerchberg-Saxton (GS) 算法是一种经典的相位恢复方法，广泛应用于光学成像、衍射成像和全息技术等领域。该算法通过迭代过程从未知相位的强度测量中恢复相位信息。

算法原理

GS算法的核心思想是利用傅里叶变换关系在空间域和频率域之间交替施加约束，逐步逼近真实解。其基本步骤如下：

1. 初始化：生成随机相位作为初始估计
2. 空间域约束：保留估计的相位，替换振幅为已知测量值
3. 傅里叶变换：将更新后的场变换到频域
4. 频率域约束：保留频域相位，替换振幅为已知测量值
5. 逆傅里叶变换：将更新后的场变换回空间域
6. 迭代：重复步骤2-5直到收敛

数学表示为：

gₖ = |f| exp(jφₖ)       // 空间域约束
Gₖ = F{gₖ}              // 傅里叶变换
Gₖ' = |F| exp(j∠Gₖ)     // 频率域约束
gₖ₊₁ = F⁻¹{Gₖ'}         // 逆傅里叶变换

MATLAB实现

function [recovered_phase, final_image] = gerchberg_saxton(...
        spatial_intensity, fourier_intensity, max_iter, tolerance)
    % GERCHBERG_SAXTON 标准GS相位恢复算法
    %
    % 输入参数:
    %   spatial_intensity - 空间域强度测量 (MxN矩阵)
    %   fourier_intensity - 傅里叶域强度测量 (MxN矩阵)
    %   max_iter - 最大迭代次数 (默认: 100)
    %   tolerance - 收敛容差 (默认: 1e-6)
    %
    % 输出参数:
    %   recovered_phase - 恢复的相位信息 (弧度)
    %   final_image - 恢复的复值图像
    
    % 参数检查与默认值设置
    if nargin < 3 || isempty(max_iter)
        max_iter = 100;
    end
    if nargin < 4 || isempty(tolerance)
        tolerance = 1e-6;
    end
    
    % 确保输入强度矩阵大小一致
    if ~isequal(size(spatial_intensity), size(fourier_intensity))
        error('输入强度矩阵必须具有相同尺寸');
    end
    
    % 获取矩阵尺寸
    [M, N] = size(spatial_intensity);
    
    % 步骤1: 初始化随机相位
    rand_phase = 2*pi * rand(M, N);
    current_field = sqrt(spatial_intensity) .* exp(1i * rand_phase);
    
    % 预分配误差数组
    errors = zeros(max_iter, 1);
    
    % 迭代过程
    for iter = 1:max_iter
        % 保存前一次迭代结果
        prev_field = current_field;
        
        % 步骤2: 应用空间域约束
        current_field = sqrt(spatial_intensity) .* exp(1i * angle(current_field));
        
        % 步骤3: 傅里叶变换到频域
        fourier_field = fftshift(fft2(ifftshift(current_field)));
        
        % 步骤4: 应用频率域约束
        fourier_field = sqrt(fourier_intensity) .* exp(1i * angle(fourier_field));
        
        % 步骤5: 逆傅里叶变换回空间域
        current_field = fftshift(ifft2(ifftshift(fourier_field)));
        
        % 计算收敛误差 (均方根误差)
        diff = abs(current_field) - abs(prev_field);
        errors(iter) = sqrt(sum(diff(:).^2) / (M*N));
        
        % 检查收敛条件
        if errors(iter) < tolerance
            fprintf('在 %d 次迭代后收敛\n', iter);
            errors = errors(1:iter);
            break;
        end
        
        % 每10次迭代显示进度
        if mod(iter, 10) == 0
            fprintf('迭代 %d, 误差 = %.4e\n', iter, errors(iter));
        end
    end
    
    % 提取恢复的相位
    recovered_phase = angle(current_field);
    final_image = current_field;
    
    % 可视化结果
    visualize_results(spatial_intensity, fourier_intensity, ...
        recovered_phase, final_image, errors);
end

function visualize_results(spatial_intensity, fourier_intensity, ...
        recovered_phase, final_image, errors)
    % 可视化结果函数
    
    % 创建新图窗
    figure('Position', [100, 100, 1200, 800], 'Name', 'GS相位恢复结果');
    
    % 原始空间域强度
    subplot(2, 3, 1);
    imagesc(spatial_intensity);
    colormap gray; axis image; colorbar;
    title('原始空间域强度');
    
    % 原始傅里叶域强度
    subplot(2, 3, 2);
    imagesc(log(1 + abs(fourier_intensity)));
    colormap gray; axis image; colorbar;
    title('傅里叶域强度 (对数尺度)');
    
    % 恢复的相位
    subplot(2, 3, 3);
    imagesc(recovered_phase);
    colormap hsv; axis image; colorbar;
    title('恢复的相位');
    
    % 恢复的强度
    subplot(2, 3, 4);
    recovered_intensity = abs(final_image).^2;
    imagesc(recovered_intensity);
    colormap gray; axis image; colorbar;
    title('恢复的空间域强度');
    
    % 相位与原始强度叠加
    subplot(2, 3, 5);
    rgb_image = ind2rgb(gray2ind(mat2gray(spatial_intensity), gray(256));
    hsv_image = rgb2hsv(rgb_image);
    hsv_image(:, :, 1) = mat2gray(recovered_phase);  % 相位映射到色相
    hsv_image(:, :, 2) = 1;  % 最大饱和度
    hsv_image(:, :, 3) = mat2gray(spatial_intensity);  % 亮度为原始强度
    rgb_phase = hsv2rgb(hsv_image);
    imshow(rgb_phase);
    title('相位(颜色)与强度(亮度)');
    
    % 收敛曲线
    subplot(2, 3, 6);
    semilogy(errors, 'LineWidth', 2);
    grid on;
    xlabel('迭代次数');
    ylabel('均方根误差 (log scale)');
    title('收敛曲线');
    
    % 恢复图像与原始图像比较
    figure('Name', '恢复效果比较');
    subplot(1, 2, 1);
    imshow(mat2gray(spatial_intensity));
    title('原始图像');
    subplot(1, 2, 2);
    imshow(mat2gray(abs(final_image).^2));
    title('恢复图像');
end

使用示例

% 生成测试图像
[M, N] = deal(256); % 图像尺寸
spatial_intensity = double(imresize(rgb2gray(imread('cameraman.tif')), [M, N]));
fourier_intensity = abs(fftshift(fft2(ifftshift(spatial_intensity)))).^2;

% 添加噪声 (可选)
noise_level = 0.05;
spatial_intensity = spatial_intensity + noise_level * max(spatial_intensity(:)) * randn(size(spatial_intensity));
fourier_intensity = fourier_intensity + noise_level * max(fourier_intensity(:)) * randn(size(fourier_intensity));

% 设置算法参数
max_iterations = 200; % 最大迭代次数
tolerance = 1e-8;    % 收敛容差

% 运行GS算法
tic;
[recovered_phase, final_image] = gerchberg_saxton(...
    spatial_intensity, fourier_intensity, max_iterations, tolerance);
fprintf('计算时间: %.2f 秒\n', toc);

% 计算重建质量指标
original_image = sqrt(spatial_intensity) .* exp(1i * angle(final_image));
mse = mean(abs(final_image(:) - original_image(:)).^2);
psnr = 10 * log10(max(abs(original_image(:)).^2) / mse);
ssim_val = ssim(abs(final_image).^2, spatial_intensity);

fprintf('重建质量指标:\n');
fprintf('  MSE: %.4e\n', mse);
fprintf('  PSNR: %.2f dB\n', psnr);
fprintf('  SSIM: %.4f\n', ssim_val);

参考源码 标准GS相位恢复算法 www.youwenfan/contentcnn/80971.html

算法特点与改进方向

标准GS算法的特点：

1. 简单易实现：算法结构清晰，易于编程实现
2. 收敛性：通常能在50-200次迭代内收敛
3. 局限性：
   • 可能陷入局部极小值
   • 对初始相位敏感
   • 对噪声较为敏感

改进策略：

% 1. 输入输出算法 (IO) - 加速收敛
function [current_field] = input_output_update(current_field, prev_field, spatial_intensity, beta)
    % beta: 松弛参数 (0.5-1.5)
    constraint = sqrt(spatial_intensity) .* exp(1i * angle(current_field));
    current_field = constraint - beta * (prev_field - constraint);
end

% 2. 混合输入输出算法 (HIO) - 避免局部极小
function [current_field] = hybrid_input_output(current_field, prev_field, spatial_intensity, beta)
    mask = abs(current_field) > 0;
    constraint = sqrt(spatial_intensity) .* exp(1i * angle(current_field));
    current_field = constraint .* mask + (prev_field - beta * constraint) .* ~mask;
end

% 3. 共轭梯度法加速 - 提高收敛速度
% 在迭代过程中添加共轭梯度方向

应用场景扩展：

1. 多平面相位恢复：

% 对多个焦平面进行测量
planes = 3; % 焦平面数量
z_positions = [0, 0.1, 0.2]; % 离焦距离

for p = 1:planes
    % 计算传播到第p个平面
    propagated_field = angular_spectrum(current_field, wavelength, z_positions(p));
    % 应用强度约束
    propagated_field = sqrt(intensity_measurements(:,:,p)) .* exp(1i * angle(propagated_field));
    % 传播回原平面
    current_field = angular_spectrum(propagated_field, wavelength, -z_positions(p));
end

2. 部分相干照明：

% 考虑光源的部分相干性
for w = 1:num_wavelengths
    current_field_w = propagate_to_object_plane(source_fields(:,:,w));
    % 应用GS更新
    % ...
    % 相干叠加
    total_field = total_field + current_field_w;
end

实际应用注意事项

1. 数据预处理：

   • 强度数据归一化：intensity = intensity / max(intensity(:))
   • 边缘加窗减少边界效应：hann_window = hanning(M) * hanning(N)'

2. 相位解包裹（恢复后处理）：

% 使用Goldstein算法解包裹相位
unwrapped_phase = GoldsteinUnwrap2D(recovered_phase);

3. 硬件考虑：

   • 考虑CCD采样与奈奎斯特频率关系
   • 校准傅里叶域坐标系统
   • 补偿光学系统的像差

4. 计算优化：

   • 使用GPU加速（gpuArray）
   • 并行计算多波长情况
   • 优化FFT计算（预计算FFT计划）

Gerchberg-Saxton算法作为相位恢复的经典方法，尽管已有50多年历史，但其核心思想仍广泛应用于现代成像技术中。通过结合现代优化方法和计算加速技术，GS算法在实时成像和大规模数据处理中仍然发挥着重要作用。