標準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?1{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;endif 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) < tolerancefprintf('在 %d 次迭代后收斂\n', iter);errors = errors(1:iter);break;end% 每10次迭代顯示進度if mod(iter, 10) == 0fprintf('迭代 %d, 誤差 = %.4e\n', iter, errors(iter));endend% 提取恢復的相位recovered_phase = angle(current_field);final_image = current_field;% 可視化結果visualize_results(spatial_intensity, fourier_intensity, ...recovered_phase, final_image, errors);
endfunction 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相位恢復算法 youwenfan.com/contentcsb/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_wavelengthscurrent_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算法在實時成像和大規模數據處理中仍然發揮著重要作用。