1.利用Matlab求某个zernike模式的PSF（点扩散函数）图像

clear all;close all;clc
% 1.空域内的参数
psf_sampling      = 0.5e-6;              % focal plane sampling in meters
lambda            = 0.6328e-6;           % 波长，[米]
N                 = 256;                 % 采样点
aperture_diameter = 0.0254;              % 孔径[m]; 1 英寸
focal_length      = 5*aperture_diameter; % 焦距[m]
RMS_SA            = 1;                   % RMS spherical aberration content in waves% 2.定义频域内的参数，为将要进行的傅里叶变换做准备
delta_fx          = 1/(psf_sampling*N);  % 频域内相邻采样点的间距
% 波长*焦距*频率间距 = FTP{瞳孔函数}=h(xi,yi);  fx = xi/λ/f; 故：xi = λffx;
x_pupil           = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length; % 像面上的x轴范围
[X_pupil,Y_pupil] = meshgrid(x_pupil);             % 像面上的积分区域
R_pupil           = sqrt(X_pupil.^2 + Y_pupil.^2); % 直角坐标换成极坐标，极径
Theta             = atan2(X_pupil,Y_pupil);        % 极角
R_norm            = R_pupil/(aperture_diameter/2); % 归一化
assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');% figure(1);imagesc(R_norm);colormap(jet);colorbar;% 3.求出极坐标下某个Zernike模式的符号表达式
j     = 5; % 第j和模式
[n,m] = Noll_to_Zernike(j);
[Zer_pol{1,1},Zer_pol{2,1},...Zer_pol{3,1},Zer_pol{4,1}] = Zernike_string_polynomial(n,-m); % 调用求导后极坐标下的Zernike多项式
W     = num2str(['RMS_SA *', Zer_pol{4,1}]); % 乘一个系数RMS_SAW1    = W;
x                 = linspace(-1,1,90);
[X_pupil,Y_pupil] = meshgrid(x,x);     % 像面上的积分区域
r     = sqrt(X_pupil.^2 + Y_pupil.^2); % 直角坐标换成极坐标，极径
t     = atan2(X_pupil,Y_pupil);        % 极角
W1    = eval(W1) .*(r<=1) ;            % 带入具体的值 得到畸变波前的相位值r     = R_norm;
t     = Theta;
W     = eval(W);        % 带入具体的值 得到畸变波前的相位值
W(R_norm>1) = 0;
W     = W';
Tit   = ['第', num2str(j),'个Zernike模式'];
% figure(2);imagesc(W);colormap(jet);colorbar;title(Tit)E           = exp(sqrt(-1)*2*pi*W);  % Complex amplitude 电磁波的表达形式
E(R_norm>1) = 0;                     % Impose aperture size
x_1 = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure(3);imagesc(x_1*1e6,x_1*1e6,angle(E)/(2*pi));
colormap(jet);colorbar;title('Wavefront Phase(waves)');% 4. 求点扩散函数Create point-spread function
psf     = abs(fftshift(fft2(ifftshift(E)))).^2; % PSF计算公式，
psf     = psf/sum(psf(:));          % 归一化
x_psf   = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure(4);imagesc(x_psf*1e6,x_psf*1e6,psf*1e6);colormap(jet);colorbar;
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));figure(16);subplot(121);imagesc(W1');colormap(jet);colorbar;title(Tit)
subplot(122);imagesc(x_psf*1e6,x_psf*1e6,psf*1e6);colormap(jet);colorbar
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));

任选两个测试结果如下：

2.利用Matlab求固定像差的PSF（点扩散函数）图像

光学中任意像差都可以通过zernike各个模式的线性叠加而成，即各个模式与对应的系数相乘再求和。

对于一组Zernike系数为A = [0.0108, -0.2564, -0.2082, -0.2147, 0.0780, -0.0405, 0.1682, 0.0239, -0.1976, 0.0085, 0.0606, 0.0450,0.0642, -0.0527, 0.0230, 0.0003, -0.0797, -0.0044, 0.0610, 0.0504, -0.1023, -0.0480, -0.0335, -0.0223,-0.0082, 0.0038, -0.0008, -0.0328, 0.0128, -0.0128, -0.0147, 0.0201, 0.0473, -0.0083, -0.0072, -0.0207,0.0247, -0.0128, 0.0159, 0.0112, -0.0026, -0.0253, -0.0218, -0.0249, -0.0131, 0.0034, 0.0125, -0.0152,-0.0143, -0.0065, -0.0342, -0.0079, 0.0221, 0.0035, -0.0118, -0.0024, -0.0091, 0.0128, -0.0079, -0.0007,0.0177, 0.0084, -0.0023, 0.0225, 0.0003, -0.0151, -0.0068, -0.0013, 0.0187, 0.0257, 0.0215, -0.0004,-0.0173, 0.0034, -0.0096, 0.0098, 0.0330, -0.0141, -0.0024, 0.0005];其与Zernike多项式各个模式相乘、求和后的像差为：

有了像差后，就可以求其PSF了。

psf_sampling      = 0.5e-6;              % focal plane sampling in meters
lambda            = 0.6328e-6;           % 波长，[米]
N                 = 256;                 % 采样点
aperture_diameter = 0.0254;              % meters; 1 inch
focal_length      = 5*aperture_diameter; % 焦距【meters】
RMS_SA            = 1;                   % RMS spherical aberration content in waves% 2.定义频域内的参数，为将要进行的傅里叶变换做准备
delta_fx          = 1/(psf_sampling*N);  % 频域内相邻采样点的间距
% 波长*焦距*频率间距 = FTP{瞳孔函数}=h(xi,yi);  fx = xi/λ/f; 故：xi = λffx;
x_pupil           = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length; % 像面上的x轴范围
[X_pupil,Y_pupil] = meshgrid(x_pupil);             % 像面上的积分区域
R_pupil           = sqrt(X_pupil.^2 + Y_pupil.^2); % 直角坐标换成极坐标，极径
Theta             = atan2(X_pupil,Y_pupil);        % 极角
R_norm            = R_pupil/(aperture_diameter/2);
assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');% 3.根据Zernike多项式 产生畸变波前的相位值WW  = [];
coe = A; % 把系数复制过来啊for i = 1 : 80     % 80个Zernike模式[n,m] = Noll_to_Zernike(i); % 调用代码见第一篇和第二篇文章哦[~,~,~, Zer_poll{4,1}] = Zernike_string_polynomial(n,-m); % 调用极坐标下的Zernike多项式WW = num2str([WW, '(' ,num2str(coe(i)),') *(' ,Zer_poll{4,1},' )' ,' + ']);clear Zer_polll
end
W                 = num2str(['RMS_SA *', WW , '0']); % 修正最后一个模式后面的 + 号
W1                = W;
x                 = linspace(-1,1,90);
[X_pupil,Y_pupil] = meshgrid(x,x);                 % 像面上的积分区域
r                 = sqrt(X_pupil.^2 + Y_pupil.^2); % 直角坐标换成极坐标，极径
t                 = atan2(X_pupil,Y_pupil);        % 极角
W1                = eval(W1) .*(r<=1) ;            % 带入具体的值 得到畸变波前的相位值
% figure(15);imagesc(W1');colormap(jet);colorbar;
r                 = R_norm;
t                 = Theta;
W                 = eval(W);        % 带入具体的值 得到畸变波前的相位值
W(R_norm>1)       = 0;
W                 = W';
E                 = exp(sqrt(-1)*2*pi*W);  % Complex amplitude 电磁波的表达形式
E(R_norm>1)       = 0;                     % Impose aperture size
x_1               = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure(2);imagesc(x_1*1e6,x_1*1e6,angle(E)/(2*pi));
colormap(jet);colorbar;title('Wavefront Phase(waves)');% 4. 求点扩散函数Create point-spread function
psf               = abs(fftshift(fft2(ifftshift(E)))).^2; %
psf               = psf/sum(psf(:));   % Normalize to unity energyx_psf             = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure(3);imagesc(x_psf*1e6,x_psf*1e6,psf*1e6);colormap(jet);colorbar;
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));figure(16);subplot(121);imagesc(W1');colormap(jet);colorbar; title('像差，单位um')
subplot(122);imagesc(x_psf*1e6,x_psf*1e6,psf*1e6);colormap(jet);colorbar
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));

0.632um波长下的测试结果如下：

而1.215um波长下的测试结果如下，影响的是其观测效果：

Zernike多项式和其PSF对应关系如下图所示。

如何根据光学中像差(相位)求出其点扩散函数相关推荐

编写两个函数，分别求10个元素数组的最大和最小值的下标，并在main函数中运行，求出最大值和最小值之差
#include<stdio.h> int Max(int a[10]) {int i, m=0, max = a[0];//定义max的初始值为a[0]for(i=1;i<10;i ...
c语言数组输入4个学生3门课成绩,编程题从键盘输入4个学生和3门课的成绩至数组中，并求出每个学生3门课的平均成绩。...
满意答案 yao713 2016.05.27 采纳率:53% 等级:8 已帮助:1213人 #include void main(){ public static void main(Strin ...
c语言realloc函数中写啥,求大神解惑realloc函数，谢谢！
该楼层疑似违规已被系统折叠隐藏此楼查看此楼 #include #include void f(int * q,int i) { int j; for(j=0;j scanf("%d&quo ...
C语言输入一个整数n，求出其阶乘
输入一个整数n,求出其阶乘 #include <iostream> #include <cstdio> using namespace std; int main() {int ...
80x86汇编语言循环结构找出最小的偶数并在屏幕上显示求出数组的平均值显示在屏幕上
题目1 写一个完整的80X86汇编语言程序:键盘输入15个数据(转换成数值,存储到一维数组中,数值的长度为字),找出最小的偶数并在屏幕上显示,若没有偶数则显示"没有偶数!". .d ...
matlab怎么求两个数的和,matlab怎么求出两个函数的交点
matlab中,两个自变量的函数怎么求最大值(急!1) 需求:利用matlab求解二元函数y=f(x1,x2)=(339-0.01*x1-0.003*x2)*x1+(399-0.004*x1-0.01 ...
python求极限函数_SymPy库常用函数
标签: 简介 SymPy是一个符号计算的Python库.它的目标是成为一个全功能的计算机代数系统,同时保持代码简洁.易于理解和扩展.它完全由Python写成,不依赖于外部库.SymPy支持符号计算. ...
ACMNO.39 分解质因数求出区间[a,b]中所有整数的质因数分解。蓝桥杯训练！
题目描述求出区间[a,b]中所有整数的质因数分解. 输入输入两个整数a,b. 输出每行输出一个数的分解,形如k=a1*a2*a3...(a1< =a2< =a3...,k也是从小到大 ...
从10W个数中随机抽走2个数，求出那两个数是多少
这道题目是从51js论坛上看到的,链接在这里>> 题目大意是: 从1到10w(共10w个数)中随机抽走2个数,然后打乱剩下的数的顺序,问如果从这剩下的数中快速的找出抽走的是哪2个数? 我想 ...
pku 1486 求出二分匹配图中的必须边
开始楞是没看懂意思,E文让我很纠结... 要判断一条边是否为二分图中必须边,方法如下: 1.先求出原图的任意最大匹配 2.对二分图某一边的所有点,删去其当前的匹配边.删的过程不是简单的将原图设为不连通 ...

如何根据光学中像差(相位)求出其点扩散函数

1.利用Matlab求某个zernike模式的PSF（点扩散函数）图像

2.利用Matlab求固定像差的PSF（点扩散函数）图像

如何根据光学中像差(相位)求出其点扩散函数相关推荐

最新文章

热门文章