Fienup算法和稀疏约束

1. 原始Fienup算法

原始的Fienup算法是一种用于相位恢复的迭代方法，交替在实域和傅里叶域施加约束。主要步骤如下：

1. 初始化：从相位的初始猜测开始。
2. 傅里叶域约束：用测量到的幅度替换当前傅里叶变换的幅度，同时保留当前的相位。
3. 逆傅里叶变换：执行逆傅里叶变换以获得更新后的实域图像。
4. 实域约束：在实域中应用约束（例如非负性或已知支持）并更新图像。

重复此过程直到收敛。

2. 引入稀疏约束

稀疏约束可以显著提高相位恢复的精度，特别是在已知信号是稀疏的情况下。下面是如何将Fienup算法修改为包括稀疏约束的方法：

1. 初始化：选择相位的初始估计，与原方法类似。
2. 傅里叶域约束：与之前一样，通过保留相位并替换幅度来更新傅里叶域。
3. 逆傅里叶变换：执行逆傅里叶变换以获得实域图像。
4. 稀疏性强制：
   • 投影和阈值处理：将实域图像投影到稀疏基（例如小波基）以获得稀疏系数。
   • 阈值处理：仅保留最大的稀疏系数（对应于最重要的稀疏成分），其余的设为零。
   • 逆投影：使用稀疏系数投影回实域。
5. 迭代：重复上述步骤直到算法收敛。

3. 稀疏Fienup算法的详细步骤

1. 初始化：

• 从初始猜测开始 $z^0[n]=|x[n]|e^{j\varphi [n]}$，其中 $\varphi [n]$ 是随机相位。

2. 傅里叶域更新：

• 计算当前估计 $z^i[n]$ 的傅里叶变换 $Z^i[k]$。

• 更新傅里叶变换幅度：$Z^i[k]\leftarrow |X[k]|e^{j\arg (Z^i[k])}$, 其中 $|X[k]|$​ 是测量到的幅度。

3. 逆傅里叶变换：

   • 计算更新后的 $Z^i[k]$ 的逆傅里叶变换 $z_{\mathrm{inv}}^i[n]$​。

4. 稀疏性强制：

   • 将 $z_{\mathrm{inv}}^i[n]$ 投影到稀疏基 $\Psi$ 以获得稀疏系数 ${\alpha }^i:$ ${\alpha }^i={\Psi }^{-1}{z}_{\mathrm{inv}}^i$。·

   • 对系数 ${\alpha }^i$ 进行阈值处理，仅保留最大的 $k$ 个系数，其余的设为零。

   • 投影回实域：$z_{\mathrm{sparse}}^i=\Psi {\alpha }^i.$​

5. 实域更新：

   • 应用任何附加的实域约束 (例如非负性)。

   • 更新估计：$z^{i+1}[n]=z_{\mathrm{sparse}}^i[n]$。

6. 收敛检查：

   • 检查算法是否收敛(例如，通过测量迭代间的误差变化)。如果未收敛，返回第2

4. 引入稀疏性的优点

• 提高精度：利用信号的稀疏性，算法可以更准确地重建相位并减少解空间。
• 抗噪能力：稀疏表示可以通过聚焦于信号的最重要成分来减轻噪声的影响。
• 效率：稀疏约束减少了非零元素的数量，使问题更易处理，算法可能更快。

5. 算法实现

import numpy as np
import matplotlib.pyplot as plt
from skimage.data import camera, checkerboard
from skimage.transform import resize
import pywt

# 加载skimage库中的示例图像
amplitude_image = camera()
phase_image = checkerboard()

# 如果振幅图像和相位图像的大小不同，调整相位图像的大小
if amplitude_image.shape != phase_image.shape:
    phase_image = resize(phase_image, amplitude_image.shape, anti_aliasing=True)

# 对振幅和相位图像进行填充
padding_width = 50
amplitude_image = np.pad(amplitude_image, pad_width=padding_width, mode='constant', constant_values=0)
phase_image = np.pad(phase_image, pad_width=padding_width, mode='constant', constant_values=0)

# 将图像归一化到[0, 1]范围内
amplitude_image = amplitude_image / np.max(amplitude_image)
phase_image = phase_image / np.max(phase_image)

# 生成复数对象
complex_obj = amplitude_image * np.exp(1j * phase_image)

# 进行傅里叶变换
ft_obj = np.fft.fft2(complex_obj)
abs_ft_obj = np.abs(ft_obj)

# 初始化相位图像
initial_phase = np.angle(np.exp(1j * np.random.rand(*amplitude_image.shape)))

# Gerchberg-Saxton (GS) 算法
def gerchberg_saxton(amplitude, initial_phase, num_iterations, epsilon):
    z = amplitude * np.exp(1j * initial_phase)
    errors = []
    for i in range(num_iterations):
        # 正向傅里叶变换
        Z = np.fft.fft2(z)

        # 计算误差并检查终止条件
        error = np.sum((np.abs(Z) - abs_ft_obj) ** 2)
        errors.append(error)
        if error <= epsilon:
            print(f'GS algorithm converged after {i + 1} iterations with error {error:.6f}')
            break

        # 在傅里叶空间中施加振幅约束
        Z = abs_ft_obj * Z / np.abs(Z)
        # 逆傅里叶变换
        z_prime = np.fft.ifft2(Z)
        # 在实空间中施加振幅约束
        z = amplitude * z_prime / np.abs(z_prime + 1e-9)

    final_phase = np.angle(z)
    return final_phase, errors

# Hybrid Input-Output (HIO) 算法
def hio(amplitude, initial_phase, num_iterations, beta, epsilon):
    z = amplitude * np.exp(1j * initial_phase)
    errors = []
    for i in range(num_iterations):
        # 正向傅里叶变换
        Z = np.fft.fft2(z)

        # 计算误差并检查终止条件
        error = np.sum((np.abs(Z) - abs_ft_obj) ** 2)
        errors.append(error)
        if error <= epsilon:
            print(f'HIO algorithm converged after {i + 1} iterations with error {error:.6f}')
            break

        # 在傅里叶空间中施加振幅约束
        Z = abs_ft_obj * Z / np.abs(Z)
        # 逆傅里叶变换
        z_prime = np.fft.ifft2(Z)
        # 在实空间中施加振幅约束，并引入反馈机制
        mask = (amplitude != 0)
        z = np.where(mask,
                     amplitude * z_prime / np.abs(z_prime + 1e-9),
                     z - beta * z_prime)

    final_phase = np.angle(z)
    return final_phase, errors

# 稀疏Fienup算法
def sparse_fienup(amplitude, initial_phase, num_iterations, sparsity_level, epsilon):
    z = amplitude * np.exp(1j * initial_phase)
    errors = []

    for i in range(num_iterations):
        # 正向傅里叶变换
        Z = np.fft.fft2(z)

        # 计算误差并检查终止条件
        error = np.sum((np.abs(Z) - abs_ft_obj) ** 2)
        errors.append(error)
        if error <= epsilon:
            print(f'Sparse Fienup algorithm converged after {i + 1} iterations with error {error:.6f}')
            break

        # 在傅里叶空间中施加振幅约束
        Z = abs_ft_obj * Z / (np.abs(Z) + 1e-9)
        # 逆傅里叶变换
        z_prime = np.fft.ifft2(Z)

        # 在小波基上施加稀疏约束
        coeffs = pywt.wavedec2(z_prime, 'haar', mode='periodization')
        coeffs_flat, coeff_slices = pywt.coeffs_to_array(coeffs)
        threshold = np.sort(np.abs(coeffs_flat))[-sparsity_level]
        coeffs_flat[np.abs(coeffs_flat) < threshold] = 0
        coeffs_thresholded = pywt.array_to_coeffs(coeffs_flat, coeff_slices, output_format='wavedec2')

        # 逆小波变换
        z_sparse = pywt.waverec2(coeffs_thresholded, 'haar', mode='periodization')

        # 应用幅度约束
        z = amplitude * z_sparse / (np.abs(z_sparse) + 1e-9)

    final_phase = np.angle(z)
    return final_phase, errors

# 参数设置
num_iterations = 50
sparsity_level = 250  # 调整此参数以适应你的稀疏性需求
epsilon = 1e-6
beta = 0.9

# 应用GS算法
gs_phase, gs_errors = gerchberg_saxton(amplitude_image, initial_phase, num_iterations, epsilon)

# 应用HIO算法
hio_phase, hio_errors = hio(amplitude_image, initial_phase, num_iterations, beta, epsilon)

# 应用稀疏Fienup算法
sparse_fienup_phase, sparse_fienup_errors = sparse_fienup(amplitude_image, initial_phase, num_iterations,
                                                          sparsity_level, epsilon)

# 设置绘图的字体和字号
font = {'family': 'serif',
        'weight': 'normal',
        'size': 20}
plt.rc('font', **font)

# 显示结果
plt.figure(figsize=(18, 6))

plt.subplot(1, 4, 1)
plt.title('Original Phase')
plt.imshow(phase_image, cmap='gray')
plt.axis('off')

plt.subplot(1, 4, 2)
plt.title('GS Phase')
plt.imshow(gs_phase, cmap='gray')
plt.axis('off')

plt.subplot(1, 4, 3)
plt.title('HIO Phase')
plt.imshow(hio_phase, cmap='gray')
plt.axis('off')

plt.subplot(1, 4, 4)
plt.title('Sparse Fienup Phase')
plt.imshow(sparse_fienup_phase, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()

# 绘制误差下降曲线
plt.figure(figsize=(12, 6))
plt.plot(gs_errors, label='GS Errors')
plt.plot(hio_errors, label='HIO Errors')
plt.plot(sparse_fienup_errors, label='Sparse Fienup Errors')
plt.xlabel('Iteration')
plt.ylabel('Error')
plt.yscale('log')
plt.title('Error Convergence')
plt.legend()
plt.show()

6. 总结

总体感觉没有啥提升