Zernike多项式的Matlab代码_matlab zernike多项式-CSDN博客

本文链接：https://blog.csdn.net/EyRe1/article/details/120860797

1. 理论介绍：

Zernike多项式【1】（荷兰物理学家弗里茨·泽尔尼克）是定义在单位圆上且满足正交的多项式序列，其在极坐标下可写为：

$\mathop Z\nolimits_j \left( {r,\theta } \right) = \left\{ {\begin{array}{*{20}{c}} {\sqrt {n + 1} \cdot \mathop R\nolimits_n^0 \left( r \right),}&{m = 0}\\ {\sqrt {2\left( {n + 1} \right)} \cdot \mathop R\nolimits_n^m \left( r \right) \cdot \cos \left( {m\theta } \right),}&{m > 0}\\ {\sqrt {2\left( {n + 1} \right)} \cdot \mathop R\nolimits_n^m \left( r \right) \cdot \sin \left( {m\theta } \right),}&{m < 0} \end{array}} \right.\$

其中， $\small \mathop Z\nolimits_j \left( {\rho ,\theta } \right)\$ 是第 j 阶Zernike模式，0≤r≤1，0≤θ≤2π，m、n 分别是多项式的角向级数和径向级数，且满足m≤n；当 n−|m| 是偶数，而径向多项式 $\small {\mathop R\nolimits_n^m \left( r \right)}\$ 定义为：

$\small R^m_n(r) = \sum_{k=0}^{\tfrac{n-m}{2}} \frac{(-1)^k\,(n-k)!}{k!\left (\tfrac{n+m}{2}-k \right )! \left (\tfrac{n-m}{2}-k \right)!} \;r^{n-2\,k}$

2. Zernike多项式的几个性质

2.1 Zernike多项式之间是相关正交的，可以用公式记为：

$\mathop \delta \nolimits_{ij} = \frac{1}{\pi }\int\limits_0^1 {\int\limits_0^{2\pi } {\mathop Z\nolimits_i \left( {r,\theta } \right)\mathop { \cdot Z}\nolimits_j \left( {r,\theta } \right)} } rdrd\theta = \left\{ {\begin{array}{*{20}{c}} {1,{\rm{ }}i = j}\\ {0,{\rm{ }}i \ne j} \end{array}} \right.$

2.2 除平移项(piston模式)外的所有正交多项式的均值为零；

2.3 每个正交多项式(不包括piston模式)的均值为0; 证明如下(利用到的是性质2.1哦)：

$\overline {\mathop Z\nolimits_i } = \frac{1}{A}\int\limits_\sum {\mathop Z\nolimits_i \left( {\rho ,\theta } \right)\rho d\rho d\theta } = \frac{1}{A}\int\limits_\sum {\mathop {\mathop Z\nolimits_1 \left( {\rho ,\theta } \right)Z}\nolimits_i \left( {\rho ,\theta } \right)\rho d\rho d\theta } = 0{\rm{ }}\left( {i \ne 1} \right)$

2.4 波前均值等于平移项(piston模式)的系数

$\begin{array}{l} \overline W \left( {\rho ,\theta } \right) = \frac{1}{A}\int\limits_\sum {W\left( {\rho ,\theta } \right)\rho d\rho d\theta } \\ {\rm{ }} = \frac{1}{A}\int\limits_\sum {\sum\limits_{i = 1}^\infty {\mathop a\nolimits_i } \mathop Z\nolimits_i \left( {\rho ,\theta } \right)} \rho d\rho d\theta \\ {\rm{ }} = \frac{1}{A}\left( {\int\limits_\sum {\mathop {\mathop a\nolimits_1 Z}\nolimits_1 \left( {\rho ,\theta } \right)} \rho d\rho d\theta + \int\limits_\sum {\sum\limits_{i = 2}^\infty {\mathop a\nolimits_i } \mathop Z\nolimits_i \left( {\rho ,\theta } \right)} \rho d\rho d\theta } \right)\\ {\rm{ }} = \mathop a\nolimits_1 + \frac{1}{A}\int\limits_\sum {\sum\limits_{i = 2}^\infty {\mathop a\nolimits_i } \mathop Z\nolimits_i \left( {\rho ,\theta } \right)} \mathop Z\nolimits_1 \left( {\rho ,\theta } \right){\rm{ }}\rho d\rho d\theta \\ {\rm{ }} = \mathop a\nolimits_1 \end{array}$

2.5 每个标准正交多项式都有一个最小方差；

2.6 波前方差为各多项式系数的平方之和，不包括平移项的系数；

$\begin{array}{l} \mathop \sigma \nolimits^2 = \frac{1}{A}{\int\limits_\sum {\left| {W\left( {\rho ,\theta } \right) - \overline W \left( {\rho ,\theta } \right)} \right|} ^2}\rho d\rho d\theta \\ {\rm{ }} = \frac{1}{A}{\int\limits_\sum {\left| {\sum\limits_{i = 1}^J {\mathop a\nolimits_i } \mathop Z\nolimits_i \left( {\rho ,\theta } \right) - \mathop a\nolimits_1 } \right|} ^2}\rho d\rho d\theta \\ {\rm{ }} = \frac{1}{A}\int\limits_\sum {\left( {\mathop a\nolimits_1^2 - 2\mathop a\nolimits_1 \sum\limits_{i = 1}^J {\mathop a\nolimits_i } \mathop Z\nolimits_i \left( {\rho ,\theta } \right) + \sum\limits_{i = 1}^J {\sum\limits_{j = 1}^J {\mathop a\nolimits_i \mathop a\nolimits_j \mathop Z\nolimits_i \left( {\rho ,\theta } \right)\mathop Z\nolimits_j \left( {\rho ,\theta } \right)} } } \right)} \rho d\rho d\theta \end{array}$

$\begin{array}{l} = \mathop a\nolimits_1^2 - \frac{2}{A}\int\limits_\sum {\left( {\mathop a\nolimits_1^2 + \sum\limits_{i = 2}^J {\mathop a\nolimits_1 \mathop {\mathop a\nolimits_i Z}\nolimits_1 \left( {\rho ,\theta } \right)} \mathop Z\nolimits_i \left( {\rho ,\theta } \right)} \right)\rho d\rho d\theta } + \frac{1}{A}\sum\limits_{i = 1}^J {\sum\limits_{j = 1}^J {\mathop a\nolimits_i \mathop a\nolimits_j } } \int\limits_\sum {\mathop Z\nolimits_j \left( {\rho ,\theta } \right)\mathop Z\nolimits_j \left( {\rho ,\theta } \right)} \rho d\rho d\theta \\ {\rm{ }} = \mathop { - a}\nolimits_1^2 + \frac{1}{A}\sum\limits_{i = 1}^J {\sum\limits_{j = 1}^J {\mathop a\nolimits_i \mathop a\nolimits_j } } \mathop \delta \nolimits_{i,j} \\ {\rm{ }} = \sum\limits_{j = 2}^J {\mathop a\nolimits_i^2 } \end{array}\$

由于本人能力有限，只有做到这个程度啊。

那么，第 j 阶模式与n和m是否有相应的联系呢？答案是肯定的咯，一般来说给定任意一个q，可以求出n和m的值。但是n和m的值与Zernike多项式的排序有关系，常见的排序方式分为Noll 序列（应用于大气湍流），OSA（人眼像差）以及Fringe三种，无论怎么排序都不需要Zernike多项式的值哦。

另外，需要注意的是：该多项式是定义在单位圆内，实际操作的时候需要加一个mask函数以屏蔽圆外的数据哦。主要还是圆外的数据比圆内的数据大的多的多啊。

${\rm{mask}} = \sqrt {{x^2} + {y^2}} \le 1\$

3. Matlab仿真部分：

3.1 定义Noll排序方式

Noll中的j，m 和 n的值分别如下所示

对于j = 1:20
j = [1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20]
NOLL下的n和m的值为：
m = [0, 1,-1, 0,-2, 2,-1, 1,-3, 3, 0, 2,-2, 4,-4, 1,-1, 3,-3, 5]
n = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5]

而OSA下的值为：
m = [0,-1, 1,-2, 0, 2,-3,-1, 1, 3,-4,-2, 0, 2, 4,-5,-3,-1, 1, 3]
n = [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5]

用代码实现Noll排序，根据Zernike多项式的序号j，返回n和m的值：

function [n,m] = Noll_to_Zernike(j)
% 序号j从1开始；j = 1 对于的是piston模式
    if(j<1) 
        error('Noll indices start at 1.');
    end
    n  = 0;
    m  = 0;
    j1 = j-1;
    while (j1 > n)
        n  = n + 1;
        j1 = j1 - n;
        m  = (-1)^j * (mod(n,2) + 2*floor((j1+mod(n+1,2))/2));
    end
end

用代码实现OSA排序，根据Zernike多项式的序号j，返回n和m的值

function [n,m] = OSA_to_Zernike(j)
 % 注意 这里的j是从0开始，j = 0对应的是piston模式，
 % 如果要从开始则传入的参数为j-1，此时可以与Noll_to_Zernike匹配
     n = floor(sqrt(2*j+1)+0.5)-1;
     m = 2*j-n*(n+2);
end

3.2 Zernike多项式的主体部分

function Znm= Zernike(j,res)
% 参数j：  Zernike多项式的序号
% 参数res：Zernike多项式的分辨率
    x           = linspace(-1,1,res);
    [x,y]       = meshgrid(x,x);
    [theta,rho] = cart2pol(x,y);% 由(x,y)换算(r,theta)

    [n,m]       = Noll_to_Zernike(j); % 调用函数，返回n和m的值
    % [n,m]     = OSA_to_Zernike(j-1)
    if m == 0
       deltam = 1;
    else
       deltam = 0;
    end
    Norm    = sqrt(2*(n+1)/(1+deltam));% 归一化因子 
    Rnm_rho = zeros(size(rho));        % 初始化径向多项式
    for s = 0:(n-abs(m))/2
        Rnm_rho = Rnm_rho+(-1)^s.*prod(1:(n-s))*rho.^(n-2*s)/(prod(1:s)*...
        prod(1:((n+abs(m))/2-s))*prod(1:((n-abs(m))/2-s))); % 径向多项式
    end
    if m < 0
       Znm = -Norm.*Rnm_rho.*sin(m.*theta); % m<0时候的zernike多项式
    else
       Znm = Norm.*Rnm_rho.*cos(m.*theta);  % m>0时候的zernike多项式
    end
    Znm = Znm.*(rho<=1); % mask = (rho<=1)，只保留单位圆内的数据
 end

3.3 主函数

clc;clear all; close all
% 单独生成某一个模式
res  = 60;                % 分辨率设置为60
Z    = Zernike(3,res);    % 第3个模式（）像散）
figure(1);imagesc(Z);colormap(jet);colorbar

% 生成前60个模式，且分辨率也为60
num = 60;
Znm = zeros(res,res,num);
for i= 1:num
    Znm(:,:,i) = Zernike(i,60);
end
figure(5);imagesc(Znm(:,:,5));colormap(jet);colorbar