用MATLAB实现Zernike多项式计算PSF与MTF的完整指南

简介：这是一份面向光学设计、图像处理及光学成像系统分析学习者的MATLAB资源包，聚焦点扩散函数（PSF）与调制传递函数（MTF）的计算，以及基于Zernike多项式的波前像差模拟。压缩包共18个文件、约1.09MB，内容以9个MATLAB脚本为主，涵盖Zernike多项式生成、PSF/MTF求解、波前像差分析等核心算法；另含1个PPT理论讲解、1个HTML文档（详细介绍Zernike多项式在人眼波前像差描述中的应用）以及4个GIF示意图，便于对照理解。已有2210人学习下载。资源既提供了可直接运行的示例脚本，也配有图文解释，适合想通过数值实验掌握PSF/MTF与Zernike像差关系的初学者，或需要快速搭建波前模拟程序的研究人员，可帮助深入理解像差对成像质量的影响并完成相关仿真任务。

1. 用Zernike多项式算PSF/MTF，到底在算什么

拿到一组像差系数，最直接的问题是：这个系统到底会把一个点光源糊成什么样？Zernike多项式负责把波前像差拆成一个个模式，点扩散函数PSF给出点光源经过系统后的能量分布，调制传递函数MTF再把这种“糊”翻译成不同空间频率的对比度损失。这套计算点扩散函数资源等于是把“像差→波前→PSF→MTF”的完整链路用MATLAB脚本走了一遍。

资源里的zernike.m、zernike_exam.m、WaveAberrationPSF.m、WaveAberrationMTF.m正好对应链路的四个节点，还附带一个关于人眼波前像差的网页文档和PPT讲稿。适合光学设计、图像复原、计算摄影方向的人作为参考实现，也适合刚接触波动光学仿真的读者照着改参数跑通流程。

2. Zernike多项式的波前拟合：模式编号、径向多项式与MATLAB实现

2.1 为什么把波前像差展开成Zernike多项式

波前像差通常写成光程差函数W(x,y)，直接存一个二维矩阵当然可以，但没法回答问题：“这个系统主要是球差还是彗差？” Zernike多项式的价值在于它是一组定义在单位圆上的正交基函数，低阶项恰好对应光学设计里常见的赛德尔像差类型。于是任意复杂波前都可以写成：

W(rho,theta) = Σ c_i Z_i(rho,theta)

这里的rho是归一化径向坐标，theta是方位角。正交性使得各阶像差在数学上尽量解耦，RMS（均方根波前误差）可以直接由系数平方和估计，这是直接用网点图或二维相位图不好做到的。

这套资源里的网页文档标题是人眼波前像差描述，人眼像差测量仪输出的也是Zernike系数，所以按照(n,m)双索引约定来理解代码最合适。n是径向度数，m是方位角阶数，对于光学系统，常用的前几项里就有离焦、像散、彗差、球差这些“老朋友”。

2.2 zernike.m：径向多项式的求和实现

在 MATLAB 里实现 Zernike 多项式，核心是径向多项式R_n^m(rho)。常见做法是直接按定义式求和：

function Z = zernike(n, m, rho, theta) % 双索引Zernike多项式，计算单个模式在极坐标网格上的值 % 输入： % n 径向度数，整数 % m 方位角阶数，整数，且 n-|m| 必须为偶数 % rho 归一化极径矩阵，范围 0~1 % theta 方位角矩阵，单位rad % 输出： % Z 与 rho/theta 同尺寸，未归一化的Zernike模式值 assert(abs(m) <= n && mod(n - abs(m), 2) == 0, 'invalid (n,m)'); R = zeros(size(rho)); for k = 0:(n - abs(m)) / 2 % 径向多项式的经典求和公式 num = (-1)^k * factorial(n - k); den = factorial(k) * factorial((n + abs(m)) / 2 - k) ... * factorial((n - abs(m)) / 2 - k); R = R + num / den * rho.^(n - 2*k); end if m >= 0 Z = R .* cos(abs(m) * theta); else Z = R .* sin(abs(m) * theta); end end

这里要注意几点。第一，rho.^(n - 2*k)在rho=0且指数为0时得到的就是1，MATLAB 的0^0行为可以放心；但如果网格生成时把圆外坐标也算进去，要先在调用处做口径掩膜，否则圆外数值没有物理意义。第二，factorial.m可以自己写也可以直接调用内置函数，资源里单独放一个factorial.m多半是为了让脚本在旧版MATLAB上也能跑，或者用来计算双阶乘。第三，正负m的约定不同场景可能相反，建议拿到别人的脚本先检查m=1和m=-1的波形方向，再应用到自己的坐标定义中。

2.3 zernike_exam.m 使用时的约定表与常见坑

下面这张表值得贴在显示器边上，它把低阶双索引模式映射到了常见像差名称：

在zernike_exam.m这类演示脚本里，最容易踩的坑是模式编号不一致。有的脚本用单索引Z(1)到Z(37)，有的用(n,m)双索引，Noll 编号和 Fringe 编号又不相同。拿到资源后不要直接改系数，先运行一遍zernike_exam.m，确认输出的波前图里哪个位置对应离焦、哪个位置对应彗差，再开始做自己的拟合。

另外一个容易被忽略的点是归一化半径。rho必须在口径边界处等于1，如果实际仿真区域是正方形网格而光学口径是内切圆，需要先构造pupil = rho <= 1，再把波前矩阵乘上这个掩膜。否则圆外的大数值会把FFT结果污染得很严重，PSF看起来像被加了一层窗函数而不是真实像差。

3. 波前到PSF：夫琅禾费FFT与参数设定

3.1 PSF计算的物理模型与离散化

一个理想点光源经过光学系统后，在像面产生的不再是几何点，而是衍射斑加上像差带来的扩散。空间不变的前提下，非相干成像系统的PSF可以写成瞳孔函数的傅里叶变换模平方：

PSF(x,y) = |F{ P(rho,theta) * exp(i * W_phase) }|^2

其中P是光瞳函数，圆口径内为1，圆外为0。W_phase是波前像差对应的相位，单位要统一。MATLAB里用fft2实现的正是离散傅里叶变换，等价于一种离散化下的夫琅禾费衍射计算。

FFT输出的是一个周期延拓的频谱，所以习惯上fft2之后紧跟fftshift，把零频移到数组中心。反过来，输入光瞳矩阵时要先用ifftshift把坐标原点挪回FFT算法的(1,1)位置。这一步写错，输出的PSF会整体平移半个周期，而且肉眼不容易看出来。

3.2 WaveAberrationPSF.m 的完整实现链路

把上一章的zernike.m接进来，PSF计算脚本可以收敛成下面这个结构：

function psf = computePSF(coeffs, modes, lambda, f, R, N) % 由Zernike系数计算非相干点扩散函数 % coeffs: 系数向量，量纲要与lambda一致 % modes: Nx2矩阵，每行是 (n,m)，与coeffs一一对应 % lambda: 工作波长 % f: 像方焦距 % R: 光瞳半径 % N: 网格尺寸，建议至少256 x = linspace(-R, R, N); [X, Y] = meshgrid(x); rho = hypot(X, Y) / R; theta = atan2(Y, X); pupil = rho <= 1; % 合成波前 W = zeros(size(X)); for k = 1:numel(coeffs) W = W + coeffs(k) * zernike(modes(k,1), modes(k,2), rho, theta); end % 瞳孔复振幅 E = pupil .* exp(1i * (2*pi/lambda) * W); % 夫琅禾费衍射：傅里叶变换取模平方 amp = fftshift(fft2(ifftshift(E))); psf = abs(amp).^2; % 能量归一化，便于后续Strehl比和MTF计算 psf = psf / sum(psf(:)); % 像面像素物理间距：lambda * f / (2*R) dx_psf = lambda * f / (2*R); end

关于参数，这里最值得确认的是系数单位。如果coeffs直接给的是波长数，那么相位计算应该是exp(1i * 2*pi * W)；如果系数单位是微米或纳米，就必须像代码里这样除以lambda。资源里的WaveAberrationPSF.m无论采用哪种写法，只要最终波前值W和波长lambda同单位即可，混用单位会导致相位整体缩小放大几十倍，PSF看起来像完全离焦。

另一个工程细节是网格生成方式。linspace(-R,R,N)简单直观，但首尾都取到了边界点，FFT实际等效的采样长度是N*dx = N*2R/(N-1)，略大于2R。做严格仿真时我一般用dx = 2*R/N; x = -R:dx:R-dx;，避免边界多算一个点。两者差异在低阶像差仿真里通常不影响结论，但当你计算接近衍射极限的MTF时，这种细节会把截止频率偏差百分之几，不值得踩。

3.3 参数表：改哪些量会让PSF发生明显变化

仿真过程中可以用斯特列尔比做快速质量判断：

strehl = max(psf(:)) / max(psfDL(:)); if strehl >= 0.8 disp('接近衍射极限'); else disp('像差明显，需要校正或后处理'); end

这里的psfDL是同样参数下把所有Zernike系数置零得到的理想PSF。斯特列尔比的物理含义是实际峰值强度与衍射极限峰值强度之比，0.8附近对应光学系统常见的“衍射受限”判据。

4. PSF到MTF：频域变换、脚本实现与MTF曲线判读

4.1 OTF与MTF的关系

MTF不是直接对波前做傅里叶变换，而是对PSF做傅里叶变换。PSF经过能量归一化后，其傅里叶变换称为光学传递函数OTF，它的模被称为调制传递函数MTF：

OTF(ξ,η) = F{ PSF(x,y) } MTF(ξ,η) = |OTF(ξ,η)|

MTF在空间频率(0,0)处的值是1，代表零频率对比度不损失；频率越高，MTF越低，对应系统对细密条纹的调制能力变差。之所以不直接看PSF二维图，是因为MTF一维截面能更直观地给出“这个系统能分辨多少线对每毫米”的结论，这也是镜头评测、图像复原中经常用MTF作为核心指标的原因。

4.2 psf2mtf函数：把PSF矩阵转成MTF

资源里的WaveAberrationMTF.m本质上就是先调用PSF计算，再做一次FFT并取模。自己写的时候可以封装成下面这样：

function [mtf, fx, fy] = psf2mtf(psf, dx) % 将归一化PSF转为MTF % psf: 上一章得到的点扩散函数矩阵 % dx: PSF像素的物理尺寸，单位与波长/焦距一致 % 输出： % mtf 与psf同尺寸的调制传递函数 % fx/fy 空间频率坐标 psf = psf / sum(psf(:)); % 反中心化 + FFT + 再中心化 otf = fftshift(fft2(ifftshift(psf))); % 取模并除以零频，保证MTF(0,0)=1 mtf = abs(otf); mtf = mtf / mtf(floor(size(mtf,1)/2)+1, floor(size(mtf,2)/2)+1); % 空间频率坐标，单位：周期/长度 N = size(psf, 1); fx = (-N/2 : N/2-1) / (N * dx); fy = fx; end

这段代码里最关键的是ifftshift和fftshift成对使用。PSF能量集中在矩阵中心附近，但MATLAB的FFT认为数组左上角是起始点，所以变换前要把中心变量挪到左上角；变换后频谱中心才在当前(1,1)位置，再用fftshift挪回中心。少一个ifftshift，MTF的相位会多出一个线性斜坡，实部虚部都偏离，但取模后的MTF可能“看起来还正常”，这类bug非常隐蔽。

频率坐标fx = (-N/2 : N/2-1) / (N*dx)的单位取决于dx的单位。如果dx是微米，那么fx单位是周期/微米；要换算成光学设计里常用的lp/mm，再乘1000即可。实际读取曲线时，一般取过中心的水平或垂直切片：

center = floor(size(mtf,1)/2) + 1; mtf_x = mtf(center, center:end); freq_x = fx(1, center:end) * 1000; % 转换为 cyc/mm

注意这里的mtf_x是从零频往单侧取，避免把对称曲线重复画两遍。

4.3 有像差系统的MTF特征与判读

对于口径均匀的圆孔衍射受限系统，非相干MTF截止频率近似为：

fc = 1 / (lambda * F#)

其中F#是像方F数。以lambda=0.55um、F#=10为例，截止频率约182 cyc/mm。仿真时如果MTF没有落在这个范围，先检查坐标换算而不是怀疑代码，因为不少人直接用像素索引当空间频率，得到的曲线数值完全不可读。

不同Zernike模式对MTF的压制方式不一样，常见表现如下：

因此，跑WaveAberrationMTF.m时不要只盯某一条频率线，建议同时画出有像差和衍射极限两条MTF曲线，观察两者相差最大的频率区间。相差集中在中频，说明是离焦或低阶球差主导；相差集中在高频，往往是高阶模式或采样不足造成的伪影。

5. 从Zernike系数反向优化：最小二乘拟合与验证三板斧

5.1 用最小二乘从波前斜率反解Zernike系数

如果手里只有波前相位采样矩阵，而目标是得到一组Zernike系数，最直接的方法是把每个采样点上的模式值拼成设计矩阵A，再对波前向量做线性最小二乘：

mask = pupil(:); M = size(modes, 1); A = zeros(size(W(mask), 1), M); for k = 1:M z = zernike(modes(k,1), modes(k,2), rho, theta); A(:, k) = z(mask); end c = A \ W(mask);

这里有几个实际经验：一是模式数量不要超过采样点数的三分之一，否则矩阵条件数变差，高频模式会和噪声互相竞争；二是拟合前先把波前的活塞项去掉，也就是把均值置零，否则(0,0)项会吸收大部分能量；三是拟合后一定要看重建残差：

Wfit = A * c; rms_resid = sqrt(mean((W(mask) - Wfit).^2));

残差的量级应该远小于波前RMS本身，如果残差偏大，多半是模式不够，或口径偏移导致Zernike正交性被破坏。

5.2 验证三板斧：残差、Strehl比和MTF对比

每次修改系数后，建议依次做三件事。第一，检查波前残差的RMS；第二，用上一章的computePSF计算Strehl比；第三，把MTF与衍射极限曲线画在一起。这三步可以分别暴露拟合问题、能量集中度问题和实际分辨率问题，比单看一张PSF彩图可靠得多。

调试的时候还习惯用单一变量法：在zernike_exam.m里只把某一个系数从0改成0.3，其他保持不变，观察PSF是否出现对应的不对称（彗差）或旋转对称扩散（离焦）。这样能快速确认当前坐标约定和模式编号没有搞错。整套脚本跑通后，再回到真实波前数据做最小二乘拟合，得到的系数才能放心用于像差补偿或图像去卷积。

本文还有配套的精品资源，点击获取