超表面全息显示

西猴柿

已于 2024-12-23 16:01:15 修改

阅读量1.8k

点赞数 10

分类专栏： Matlab 文章标签： 1024程序员节算法

于 2024-10-24 22:50:43 首次发布

本文链接：https://blog.csdn.net/qq_60088512/article/details/143220125

版权

Matlab 专栏收录该内容

5 篇文章

订阅专栏

本文记录的是一个超表面全息显示的案例。首先是G-S算法生成相位全息图，得到超表面单元的转角分布，再进行数值模拟得到全息显示图像。虽然是一个简单的例子，但是本人一个摸索的过程，对于初学者来说或许会有一定的帮助。如果读者有问题欢迎留言。

1. G-S算法原理

G-S 算法是计算生成全息中常用的优化相位分布的迭代方法，算法流程见图1。首先是初始化， $\left | A_{t} \right |$ 作为目标图像的振幅信息，同时施加一个随机相位RAP来模拟漫反射，这两个信息组成了成像面的复振幅 $A_{0}$ ：

$A_{0}=\left|A_{t}\right| \exp \left(j \varphi_{0}\right)$

在给定初始相位分布下，依次迭代输入和输出的振幅分布，迭代过程有四个步骤：

Step 1：第k次迭代中，对像平面的复振幅 $A_k$ 进行逆傅里叶变换（IFT）计算，得到全息面上的复振幅 $H_k$

$H_{k}=\left|H_{k}\right| \exp \left(j \phi_{k}\right)=\operatorname{IFT}\left(A_{k}\right), k=1,2, \ldots, n$

Step 2：提取 $H_k$ 的相位，并对 $H_k$ 的振幅归一化，得到纯相位全息图（POH） $H'_k$

$H'_k=exp(j\phi '_k)=FT(H'_k)$

Step 3：对POH进行傅里叶计算重构成像面的复振幅 $A'_k$

$A'_k=\left | A'_k \right |exp(j\varphi'_k)=FT(H'_k)$

Step 4：将重构振幅 $\left | A'_k \right |$ 替换为目标振幅 $\left | A_t \right |$ ，生成下一次迭代的输入复振幅 $A_{k+1}$

$A_{k+1}=\left | A_t\right | \cdot \frac{A'_k}{\left | A'_k \right | }$

通常使用均方根误差（MSE）或相关系数（CC）来控制误差。通过计算重建振幅与目标振幅的差来控制迭代过程，当MSE或CC满足一定的预设条件，迭代停止。

图1 相位恢复算法流程图。（图片来源：Light Sci Appl. 2022;11(1):231.）

2. 相位恢复结果

%%%%%%%% G-S相位恢复算法  %%%%%%%%%%
%% 抽样
addpath('dataset')
image0 = imread('Letter_O.png');
image0 = rgb2gray(image0); 
image0 = im2double(image0);
figure
imshow(image0)
title('输入图像')

M=80; N=80;                                            
image0 =imresize(image0, [M,N]); % 指定新的尺寸
figure
imshow(image0)
title('抽样图像')

%% 生成RAP
PI = 3.14159;
image1 = image0;
phase = 2i*PI*rand(M,N);    
image1 = image1.*exp(phase);

%% 计算编码
image2 = ifft2(ifftshift(image1));    % 第一次二维IFT，频域——>时域

% 迭代过程
for t=1:1:1100
    imgangle = angle(image2);                       % 取位相
    image = exp(1i * imgangle);
    image = fftshift(fft2(image));                  % FT 还原，复数图像——>频域
    imgabs = abs(image)/max(max(abs(image)));       % 归一化
    sim = corrcoef(image0,imgabs);                  % 相关系数
    if sim(1,2) >= 0.999   % 满足条件，跳出循环
        break
    else
        imgangle = angle(image2);
        image2 = exp(1i*imgangle);
        image3 = fftshift(fft2(image2));
        imgabs = abs(image3)/max(max(abs(image3)));
        imgangle = angle(image3);
        image3 = exp(1i*imgangle);
        image3 = image3.*(image0+rand(1,1)*(image0-imgabs));
        image2 = ifft2(ifftshift(image3));  
    end
end

%% 输出位相全息图 
imgangle = angle(image2);           % 取位相（弧度制）
image4 = exp(1i*imgangle);
image4 = fftshift(fft2(image4));    % 重建结果
imgabs = abs(image4)/max(max(abs(image4)));

figure
imshow(imgangle);
title('位相全息图')

figure
imshow(imgabs);
title('重建图')

%% 存储位相全息图
 level = 4;                         % 台阶数
 res = (imgangle+PI)/(2*PI/level);
 phi = (round(res)*360/level)/2;    % 转角
 figure
 imshow(phi,[])
 title('相位板');

% name_save = 'phase_only';
% dlmwrite([name_save,'.txt'], phi, 'delimiter', '\t' );