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.\$

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

$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.');endn  = 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 == 0deltam = 1;elsedeltam = 0;endNorm    = sqrt(2*(n+1)/(1+deltam));% 归一化因子 Rnm_rho = zeros(size(rho));        % 初始化径向多项式for s = 0:(n-abs(m))/2Rnm_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))); % 径向多项式endif m < 0Znm = -Norm.*Rnm_rho.*sin(m.*theta); % m<0时候的zernike多项式elseZnm = Norm.*Rnm_rho.*cos(m.*theta);  % m>0时候的zernike多项式endZnm = 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:numZnm(:,:,i) = Zernike(i,60);
end
figure(5);imagesc(Znm(:,:,5));colormap(jet);colorbar

3.4 结果部分

运行代码后即可出现图像哦。

参考文献

[1] NOLL R J. Zernike polynomials and atmospheric turbulence [J]. JOsA, 1976, 66(3): 207-11.

Zernike多项式的Matlab代码相关推荐

Matlab读Zygo干涉仪面形数据并进行37阶Zernike拟合（附Matlab代码）
本代码具有以下两个主要功能: 一.Matlab读Zygo干涉仪/轮廓仪.dat文件中的面形数据通过在Matlab中开发大量的反编译代码,我们成功地在Matlab中提取到Zygo .dat二进制文件中 ...
【数字信号处理】卷积编程实现 ( Matlab 卷积和多项式乘法 conv 函数 | 使用 matlab 代码求卷积并绘图 )
文章目录一.Matlab 卷积和多项式乘法 conv 函数二.使用 matlab 代码求卷积并绘图一.Matlab 卷积和多项式乘法 conv 函数 Matlab 文档地址 : https:// ...
基于神经网络多项式插值的图像超分辨重构研究-附Matlab代码
⭕⭕ 目录 ⭕⭕ ✳️ 一.引言 ✳️ 二.基于单帧图像的超分辨率重构技术 ✳️ 2.1 最近邻域插值法 ✳️ 2.2 双线性插值法 ✳️ 2.3 双三次插值法(Keys'插值) ✳️ 三.神经网络 ...
matlab数值分析拟合实例,数值分析函数拟合matlab代码.doc
数值分析函数拟合matlab代码.doc 第一题MATLAB代码用SPLINE作图XI0204060810YI098092081064038X10012Y1NEWTON3XI,YI,X源代码见M文件Y ...
5G NR OFDM链路层仿真及Matlab代码实现（1）：LDPC信道编译码之5G Tollbox中相关函数使用介绍
目录 0. 写在前面 1. 信道编码 2. 经过信道(❤) 3. LDPC译码 0. 写在前面这是OFDM链路仿真系列的第一部分. 第一次写博客,望前辈们多多指教!!! 本来想着一口气写完整个过程, ...
matlab函数无限冲激响应滤波器,MATLAB代码有限冲激响应（FIR）滤波器和无限冲激响应（IIR）滤波器...
MATLAB有限冲激响应(FIR)滤波器和无限冲激响应(IIR)滤波器设计附MATLAB代码摘要文章设计了一个数字信号处理仿真实验,产生一个信号,其频率成分为f1和f2,并对其进行理想采样,采样 ...
【路径优化】基于帝企鹅算法求解TSP问题（Matlab代码实现）
目录 1 帝企鹅算法 2 旅行商问题(TSP) 3 运行结果 4 参考文献 5 Matlab代码实现 1 帝企鹅算法帝企鹅优化算法(emperor penguin optimizer,EPO)是Ga ...
【图像增强】基于Step和Polynomial 滤波实现图像增强附matlab代码
1 内容介绍本文提出的Step过滤和多项式过滤方法的实现提供对部分曲线结构的增强. 该方法应用于地球物理图像,以识别线性模式存在于考古遗址 [1] 和故障检测 [2] 中的地下建筑结构.该方法有 ...
【滤波估计】基于双卡尔曼滤波实现soc和soh联合估计附matlab代码
1 内容介绍对电动汽车电池管理系统进行电池状态估计非常重要准确充电,实现电池模型参数的在线更新.在本文中,开路电压的估计转换为开路电压拟合的估计参数,快速时变参数开路电压被转换成几个慢时变参数.提出 ...

Zernike多项式的Matlab代码