模拟哈特曼波前探测器

这是模拟哈特曼波前探测器的另一份代码，这里不需要设置微透镜的参数。

波前 ==> 光点，使用的是 离散一维DFT对应的fft(，)，而不是 二维FFT对应的 fft2()，最大的优点是运算速度极快，其原因是 不需要通过双循环分别对 某个子透镜形成光点处理，利用并行运算直接得到整个微透镜的最终形成的光点阵。

也可以参考OOMAO中的哈特曼波前探测的模拟源码，源码见：

【GitHub - rconan/OOMAO: Object-Oriented, Matlab & Adaptive Optics】

clc;clear all; close all
% function wavePrgted = propagateThrough1(obj,src_in)
%{
for a given wave of size [n, n],   full width at half maximum (FWHM)
nyquistSampling = 1, means 2 pixels per fhwm obtained by setting fftPad = 2,
FHWM 占据两个像素，可通过fftPad = 2 完成设置
fieldStopSize should be set to default 单子透镜分辨率 = n/nLenslet;
        
nyquistSampling = 0.5, means 1 pixels per fhwm obtained by setting fftPad = 1, 
fieldStopSize should be set to default n/nLenslet/nyquistSampling/2

下面用的是 fft(A,) 而不是使用fft(fftshift(A) , ),
fft(),操作使用 有用的信息在四周，而不是中央，故
v((0:nLensletImagePx-1)+centerIndex-halfLength) ...
                    = false; 
半宽高度的范围就是在四周，v设置mask的时候就中心部分的非半宽高度区域内的数据 清除掉 

%}   

nLenslet            = 10;           % 微透镜个数
nLensletImagePx     = 10;           % 单个子透镜分辨率
nLensletsImagePx    = nLenslet * nLensletImagePx;

%% 给出一个湍流波前，必须是弧度单位。 波长为1.215um，==》 rad,乘以系数0.1934

real_coe = [ ...
    -0.0804, 1.5517,-1.7531, 0.7195, 0.3470, 0.1209, 0.0477, 0.4413, 0.0409,-0.0088...
    -0.1640, 0.0846,-0.0254, 0.0734,-0.1656,-0.2289, 0.1348,-0.0130, 0.0624,-0.1272...
    -0.0083, 0.0949, 0.0321, 0.0319, 0.0067,-0.0526,-0.0635,-0.0137,-0.0004, 0.0381...
    -0.0893,-0.0089, 0.0125,-0.0202,-0.0010,-0.0220];

num = numel(real_coe);
Znm = zeros(nLensletsImagePx,nLensletsImagePx,num);       % 初始化像差函数
for i= 1:num
    Znm(:,:,i) = Zernike(i,nLensletsImagePx);% 调用Zernike多项式函数
end
WF = reshape(Znm,[nLensletsImagePx*nLensletsImagePx,num]) * real_coe';
WF = reshape(WF,[nLensletsImagePx,nLensletsImagePx]);     % 系数和对应多项式的线性叠加后的结果

%% 给出一个mask，也就是振幅pupil
x                   = linspace(-1,1,nLensletsImagePx);
[X,Y]               = meshgrid(x,x);
rho                 = sqrt(X.*X + Y.*Y);
rho(rho>1)          = 0;
rho(rho ~=0)        = 1;
mask                = rho * 14.9; % 为啥乘以14.9，我也不清楚哦
val                 = bsxfun( @times , mask , exp(1i*WF) ); % 复振幅

figure(2);
subplot(131);imagesc(WF);colormap(jet);colorbar;title('zernike生成的波前,单位rad')
subplot(132);imshow(val, []);colormap(jet);colorbar;title('振幅.*exp(i * 弧度下的相位)')
subplot(133);imagesc(angle(val));colormap(jet);colorbar;title('angle作用的对象的单位 必须 是弧度')

%% 波前 == > PSF，满足奈奎斯特采样，需要设置像面大小
[nLensletsWavePx,nLensletsWavePxNGuideStar,nWave] ...
                    = size(val);
nLensletsWavePxNGuideStar ... 
                    = nLensletsWavePxNGuideStar*nWave;
val                 = reshape(val,nLensletsWavePx,nLensletsWavePxNGuideStar);
nLensletWavePx      = nLensletsWavePx ./ nLenslet;
nLensletArray       = nLensletsWavePxNGuideStar/nLensletsWavePx;
fftPad              = 2;                        % 采样，nyquistSampling = 1
nOutWavePx          = nLensletWavePx * fftPad;  % 采样前 单个子透镜形成的PSF在像平面坐标区域
evenOdd             = rem(nLensletWavePx,2);
if ~rem(nOutWavePx,2) && evenOdd
    nOutWavePx = nOutWavePx + evenOdd;
end
nOutWavePx          = max(nOutWavePx,nLensletWavePx);
nLensletSquareWavePx ...
                    = nLensletWavePx * nLenslet^2 * nLensletArray;
wavePrgted          = zeros(nOutWavePx,nLensletSquareWavePx);
val                 = val ./ nOutWavePx;

%% get fftPhasor
nOutWavePx      = nLensletWavePx * fftPad;    % Pixel length of the output wave
evenOdd         = rem(nLensletWavePx,2);
if ~rem(nOutWavePx,2) && evenOdd
    nOutWavePx  = nOutWavePx + evenOdd;
end
nOutWavePx      = max(nOutWavePx,nLensletWavePx);

% shift the intensity of half a pixel for even sampling
[u,v]           = ndgrid((0 : (nLensletWavePx-1)) .* (~rem(nLensletWavePx,2)-nOutWavePx)./nOutWavePx);
fftPhasor       = repmat( exp(-1i.*pi.*(u+v)) , nLenslet, nLenslet );
% figure(3);imshow(fftPhasor,[]);colormap(jet);colorbar;
  
if isempty(fftPhasor)
    % Shape the wave per columns of lenslet pixels
    % 为lensletd 每列像素 塑造波形
    val             = reshape(val,nLensletWavePx,nLensletSquareWavePx);
    u               = any(val); % Index of non-zeros columns
    wavePrgted(:,u) = fftshift(fft(val(:,u),nOutWavePx),1);
    % Select the field of view
    v               = [];
    if nOutWavePx   > nLensletImagePx
        centerIndex = ceil((nOutWavePx+1)/2);
        halfLength  = floor(nLensletImagePx/2);
        v           = true(nOutWavePx,1);
        v((0:nLensletImagePx-1)-halfLength+centerIndex) ...
                    = false;
    elseif nOutWavePx<nLensletImagePx
        error('lensletArray:propagateThrough:size','The computed image is smaller than the expected image!')
    end
    wavePrgted(v,:) = [];
    % Back to transpose 2D
    val             = reshape( wavePrgted ,...
                                nLensletsImagePx,nLensletWavePx*nLenslet*nLensletArray).';
    figure(2);imshow((val),[]);colormap(jet);colorbar
    % Shape the wave per rows of lenslet pixels
    val             = reshape(val,nLensletWavePx,nLensletsImagePx*nLenslet*nLensletArray);
    u               = any(val); % Index of non-zeros columns
    wavePrgted      = zeros(nOutWavePx,nLensletsImagePx*nLenslet*nLensletArray);
    wavePrgted(:,u) = fftshift(fft(val(:,u),nOutWavePx),1);
    wavePrgted(v,:) = [];

else
    %% 
    val             = val.*repmat(fftPhasor,1,nLensletArray);
    val             = reshape(val,nLensletWavePx,nLensletSquareWavePx);
    u               = any(val); % Index of non-zeros columns
    wavePrgted(:,u) = fft(val(:,u),nOutWavePx);
    
    v               = []; 
   
    
    if nOutWavePx   > nLensletImagePx
        centerIndex = ceil((nOutWavePx+1)/2) + rem(nLensletWavePx,2);
        halfLength  = floor(nLensletImagePx/2);
        v           = true(nOutWavePx,1);
        v((0:nLensletImagePx-1)+centerIndex-halfLength) ...
                    = false; % fft()操作后 能量集中在四周，此时要把中心区域的数据去除掉（对应半宽高度之外的区域）
    elseif nOutWavePx<nLensletImagePx
        error('lensletArray:propagateThrough:size','The computed image is smaller than the expected image!')
    end
        wavePrgted(v,:) = [];
        % Back to transpose 2D
        val             = reshape( wavePrgted ,...
                                nLensletsImagePx,nLensletWavePx*nLenslet*nLensletArray).';
        % Shape the wave per rows of lenslet pixels
        val             = reshape(val,nLensletWavePx,nLensletsImagePx*nLenslet*nLensletArray);
        u               = any(val); % Index of non-zeros columns
        wavePrgted      = zeros(nOutWavePx,nLensletsImagePx*nLenslet*nLensletArray);
        wavePrgted(:,u) = fft(val(:,u),nOutWavePx);
        wavePrgted(v,:) = []; % fft()操作后 能量集中在四周，此时要把中心区域的数据去除掉（对应半宽高度之外的区域）
end
% Back to transpose 2D
wavePrgted  = reshape(wavePrgted,nLensletsImagePx*nLensletArray,nLensletsImagePx).';
[n,m]       = size(wavePrgted);
wavePrgted  = reshape(wavePrgted,[n,m/nWave,nWave]);
psf         = wavePrgted.*conj(wavePrgted);

figure(4);
subplot(121);imagesc(WF);colormap(jet);colorbar;title('平面波')
subplot(122);imagesc(psf);colormap(jet);colorbar;title('对应的光点阵图像')

左侧是像差，右侧是光点阵。