使用deformable mirror 拟合zernike多项式

一、目的：使用DM拟合zernike多项式

Deformable mirror（DM）可以用于补偿哈特曼波前探测器（或者其他探测器）测得的像差，假定该像差为zernike多项式中的某些项，若要补偿掉这部分的像差，则DM就应该拟合出对应的面形。因此，本文将介绍如何在仿真中实现它。

使用工具介绍：MATLAB，Object-oriented Matlab adaptive optics（OOMAO）【GitHub - rconan/OOMAO: Object-Oriented, Matlab & Adaptive Optics】

后者是一个adaptive optics （AO）仿真平台，可以完成望远镜系统中的AO校正像差的过程。

二、做法：

1. 测响应函数

       当 DM 表面被驱动器推动时，其实际产生的面形称为其响应函数。即单独对一个驱动器施加恒定电压，使其镜面产生相应的形变后，平行光波经镜面反射再由 波特曼波前探测器（SHWFS）测得畸变波前的斜率值或 Zernike 系数。而整个 DM 的响应函数的测量，则通过对所有驱动器重复相同的操作，最终确定出恒定输入（电压值）和输出（斜率值或 Zernike 系数）之间的静态的线性映射。 可以参照【变形镜的Matlab仿真实现】这一篇进行理解。

2.AO 系统的校正原理

步骤1，t 时刻，相位为 φi(x, y, t) 的畸变波前入射到DM上；

步骤2，光束从DM反射并入射到 SHWFS，SHWFS测的波前斜率；

步骤3，在测得的斜率响应矩阵作用下，t时刻的波前斜率测量值 被转换成  施加到DM上的电压值；

步骤4，DM根据电压拟合出对应的面形，来抵消畸变波前中的一部分相位（PID算法），计算当前波前的RMS值。若RMS值 > 某一阈值，则重复 步骤1-2-3-4；

步骤5，RMS值是一种收敛的曲线，且已经达到了最低效果，此时校正完毕。

      使用AO系统拟合zernike多项式的原理也是如此，之所以要拟合zernike多项式的原因是：评价DM用于拟合高阶zernike多项式的品质，即能否保留下高频信息。

3.代码实现

使用21x21的DM，先测响应矩阵函数

clc;clear all ;close all
%{
    目的：用DM拟合zernike多项式;
    先测DM的的响应函数，根据根据哈特曼测得的斜率==> 驱动电压 
        ==> 施加到DM产生补偿面形。
%}
atm = atmosphere(photometry.V,0.15,30,...
    'altitude',[0,4,10]*1e3,...
    'fractionnalR0',[0.7,0.25,0.05],...
    'windSpeed',[5,10,20],...
    'windDirection',[0,pi/4,pi]);

nLenslet = 20;            
nPx      = nLenslet * 10; 
tel      = telescope(3.6,...
            'fieldOfViewInArcMin',2.5,...
            'resolution',nPx,...
            'samplingTime',1/300);
cam      = imager('diameter',tel.D, 'fieldStopSize',20,'nyquistSampling',8); 
ngs      = source('wavelength',photometry.J,'magnitude',5); 
wfs      = shackHartmann(nLenslet,nPx,0.75);
ngs      = ngs.*tel*wfs*cam;
wfs.INIT    
+wfs;
figure(1);imagesc(cam.frame);colormap(jet);colorbar;title('理想PSF')
bif       = influenceFunction('monotonic',50/100); 
nActuator = nLenslet + 1;       % 某一方向上的 驱动器个数 比 子透镜个数多1
dm        = deformableMirror(nActuator,...
                'modes',bif,...
                'resolution',nPx,...
                'validActuator',wfs.validActuator);

name = ['DMIM_',num2str(nActuator),'nActuators','.mat'];
if exist (name,'file') == 0
    fprintf('DM的响应函数不存在，计算中...\n')
    DM_mode = [];
    Sx      = [];
    figure(2)
    for i = 1:length(dm.coefs)
        dm.coefs(i)    = 0.3038*1e-7;
        ngs            = ngs.*tel *dm*wfs*cam; % ngs.*tel *dm*wfs*cam;
        Sx(:,i)        = wfs.slopes;
        DM_mode(:,:,i) = dm.surface;
        dm.coefs(i) = 0;
        ngs         = ngs.*tel *dm*wfs*cam;
    end
    Slope   = Sx / 0.3038 ;
    save(name, 'Slope','DM_mode')
else
    fprintf('响应函数已经存在，不计算\n')
    load(name)
end

figure(5);
subplot(121);imagesc(Slope);colorbar;colormap(jet)
subplot(122);imagesc(pinv(Slope));colorbar;colormap(jet)

再进行拟合实验

DMvot_total     = []; % 拟合多项式需要的电压
modes           = 35;
for j           = 1:35
    zernn       = zernike(2:modes+1,'resolution',nPx);
    C           = zeros(modes,1);
    C(j)        = 0.1 * 1e-6;       % 拟合第i个Zenrike
    zernn.c     = C;                
    ngs         = ngs.* zernn*wfs*cam ; 
    WF_ini      = ngs.meanRmPhase;  
    loopGain    = 0.030;
    row         = 2;
    column      = 3;
    figure(3)
    subplot(row,column,1);imagesc(WF_ini);colorbar;colormap(jet);title('待拟合出的WF，单位rad')
    subplot(row,column,4);imagesc(cam.frame);colormap(jet);colorbar;title('未用AO前的PSF图像')
    RMS_ini = std(WF_ini);
    i           = 1;
    RMS         = [];
    DM_vot      = [];
    flag        = 1;  
    while (flag == 1)
        ngs          = ngs * dm * wfs * cam; 
        dm.coefs     = 0.4*dm.coefs - loopGain* pinv(Slope) * wfs.slopes *1e-6; % 单位 要到um级别
        DM_vot(:,i)  = dm.coefs;

        subplot(row,column,5);imagesc(cam.frame);colormap(jet);colorbar;title('经过AO后校正后的PSF')
        subplot(row,column,2);imagesc(ngs.meanRmPhase);colorbar;colormap(jet);title('经过AO后校正残余波前')
        subplot(row,column,[3,6]);imagesc(WF_ini - ngs.meanRmPhase);colorbar;colormap(jet);title('拟合出的多项式')
        numm = ['第',num2str(i),'次的RMS = ',num2str( std(ngs.meanRmPhase(:)) )]
        RMS(:,i) = std(ngs.meanRmPhase(:));
        if ( i > 30 && (RMS(:,i-1) - RMS(:,i) < 1e-3) )
            flag = 0;
        end
        drawnow;
        pause(0.1)
        i = i+1;
    end
    DMvot_total(:,j) = sum(DM_vot,2) * 1e6 ;
end

三、结果

1. 拟合倾斜模式的情况

2.  拟合离焦模式情况

3. 其他模式，自己从github上下载OOMAO包，使用上述代码即可运行了哦。