全息成像中基于菲涅尔法的仿真实现

U(x′,y′)=iλzeikz​∬U0​(x,y)ei2zk​(x′2+y′2)e−i2πλzx′x+y′y​dxdy

import numpy as np
from scipy.fft import fftshift, ifftshift, fft2, ifft2

def fresnel_propagate(hologram, wavelength, distance, pixel_size):
    """
    菲涅尔衍射传播子程序
    :param hologram: 输入全息图 (2D numpy array)
    :param wavelength: 波长 (m)
    :param distance: 传播距离 (m)
    :param pixel_size: 像素尺寸 (m)
    :return: 传播后的复振幅分布
    """
    M, N = hologram.shape
    dx, dy = pixel_size
    k = 2 * np.pi / wavelength
    
    # 计算菲涅尔传播因子
    fx = fftshift(np.fft.fftfreq(N, d=dx))
    fy = fftshift(np.fft.fftfreq(M, d=dy))
    FX, FY = np.meshgrid(fx, fy)
    H = np.exp(1j * k * distance * (FX**2 + FY**2) / (2 * z))
    
    # 执行傅里叶变换
    hologram_fft = fftshift(fft2(ifftshift(hologram)))
    propagated = hologram_fft * H
    field = np.real(ifftshift(ifft2(propagated)))
    
    return field

def generate_hologram(obj_wave, ref_wave):
    """
    生成干涉全息图
    :param obj_wave: 物光波复振幅 (2D array)
    :param ref_wave: 参考光波复振幅 (2D array)
    :return: 全息图强度分布
    """
    intensity = np.abs(obj_wave + ref_wave)**2
    return intensity

lambda_ = 632.8e-9  # 波长 (He-Ne激光)
z_record = 0.1      # 记录距离 (m)
z_recon = 0.3       # 再现距离 (m)
pixel_size = 4.6e-6 # 像素尺寸 (μm)
size = 512          # 图像尺寸

# 生成测试物体（分辨率板）
x = np.linspace(-size/2, size/2, size) * pixel_size
X, Y = np.meshgrid(x, x)
obj_wave = np.exp(-(X**2 + Y**2)/(2*(100e-6)**2))  # 高斯光束

# 参考光（倾斜平面波）
alpha = np.radians(1)  # 倾斜角1°
ref_wave = np.exp(1j * k * z_record) * np.exp(1j * k * (X * np.cos(alpha) + Y * np.sin(alpha)))

hologram = generate_hologram(obj_wave, ref_wave)

# 菲涅尔逆传播
reconstructed = fresnel_propagate(hologram, lambda_, z_recon, pixel_size)

# 分离实像与共轭像
mask = np.fft.fftshift(np.fft.fft2(hologram)) > 0.5 * np.max(np.abs(np.fft.fft2(hologram)))
clean_recon = np.where(mask, reconstructed, 0)

def spectral_filtering(hologram, cutoff_freq):
    """频域滤波抑制噪声"""
    fft_holo = fftshift(fft2(ifftshift(hologram)))
    fft_holo[(np.abs(fft_holo) < cutoff_freq)] = 0
    return np.real(ifftshift(ifft2(fft_holo)))

from skimage.restoration import unwrap_phase

def phase_reconstruction(intensity):
    """从强度分布恢复相位"""
    phase = np.angle(fftshift(fft2(ifftshift(intensity))))
    unwrapped_phase = unwrap_phase(phase, log=True)
    return unwrapped_phase

# 全息成像仿真主程序
import numpy as np
import matplotlib.pyplot as plt

# 参数设置
lambda_ = 632.8e-9
pixel_size = 4.6e-6
size = 512
z_record = 0.1
z_recon = 0.3

# 生成测试物体
x = np.linspace(-size/2, size/2, size)*pixel_size
X, Y = np.meshgrid(x, x)
obj_wave = np.exp(-(X**2 + Y**2)/(2*(100e-6)**2))

# 生成参考光
alpha = np.radians(1)
ref_wave = np.exp(1j*k*z_record) * np.exp(1j*k*(X*np.cos(alpha)+Y*np.sin(alpha)))

# 记录全息图
hologram = generate_hologram(obj_wave, ref_wave)

# 菲涅尔再现
reconstructed = fresnel_propagate(hologram, lambda_, z_recon, pixel_size)

# 显示结果
plt.figure(figsize=(12,6))
plt.subplot(121), plt.imshow(np.abs(reconstructed), cmap='gray'), plt.title('Reconstructed Amplitude')
plt.subplot(122), plt.imshow(np.angle(reconstructed), cmap='hsv'), plt.title('Reconstructed Phase')
plt.show()