点扩散函数理论

第33卷第4期 Vol_33No．4 红外与激光 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 Infrared and Laser Engineerin 2004年 8月 Aug．2004 点扩散函数的一维数值计算及其MATLAB实现 江月松 ，邱志伟 ，李 铮 (1．北京航空航天大学 电子信息工程学院 光电与信息工程系，北京 100083； 2．北京 东方计量技术研 究所 ，北京 100080) * 摘 要：通过应用复指数函数与 Bessel函数的展开关系和将径向对称的光瞳函数展开成方位 角的Fourier级数的方法，可以将衍射受限光学成像系统的点扩散函数的二维 Fourier变换的计算 转换为一维 Fourier变换和一维 Hankel变换的计算。并借助 MATLAB软件在计算机上实现 了 编程计算。与二维 Fourier变换计算方法相比，一维计算可以方便而清晰地获得像平面上任一确 定方向的点扩散强度分布。与 Fortran、Basic和 C等编程计算语言相比，应用 MATLAB语言编程 计算，程序语言简洁，且大大降低了使用者对数学基础和计算机语言知识的要求，是进行科学与工 程研究的高效 工具。 关键词：点扩散函数； 一维计算； Hankel变换 ； MATLAB软件 中图分类号：TN21；TP31 文献标识码：A 文章编号：1007-2276(2004)04—0405-04 One-dimensional digital computation and realization on M ATLAB of point spread function Jiang Yue—song ，Qiu Zhi—wei ，Li Zheng (1．Department of()ptoelectronics and Information，School of Information Electrical Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100083，China；2．Institute of Beiiing Orient Measuring Technique．Belling 100080，China) Abstract：The computation Of tWO—dimensional Fourier transform for point spread function of diffraction limited pupil can be transformed to the computation of one—dimensional Fourier and Hankel transform by use of the relationship between exponential function and Bessel function and property of an expansion of a optical symmetrical pupil function into azimuthal Fourier series． The detailed numerical procedure based on one—dimensional digital algorithm is completed on MATLAB．In comparison with the classic two—dimensional Fourier transform algorithm ．the one— dimensional method of Fourier and Hankel transform can obtain intensity samples along arbitrary radial direction at the image plane．And compared with Fortran，Basic and C language etc．，the M ATLAB program has advantages of succinct language，great reductions of demand of math— ematic basis and computer language，and it is a efficient tool for scientific and engineering re— search． Key words：Point spread function； One—dimensional computation； Hankel transform ； M ATLAB software 收稿日期：2003—09—02； 修订日期：2003—10—1 5 基金项目：武器装备预研基金资助项目(51401020104HK0106) 作者简介：江月松(1959一)，男，江苏淮安人，教授，博士，主要从事红外成像、综合孔径成像技术与系统、光电制导与对抗等方面的研究。 维普资讯 http://www.cqvip.com 406 红外与激光工程 第33卷 0 引 言 点扩散函数是评价光学系统成像质量的基本工 具，对于衍射受限成像系统，点扩散 函数就是计算 Fraunhofer衍射积分的值 。]，通常在实际工程应用 中，很难求出衍射积分的精确解析式，因此计算点扩 散函数的一般做法是计算含有像差的光瞳函数的二 维快速 Fourier(FFT)变换作数值近似计算。二维 FFT数值计算的好处是可以从总体上给出像面上点 扩散函数强度分布概况，但要精确地观察某一方向的 点扩散函数的强度分布却比较困难。为此，本文研究 了在某一确定的方位角的情况下点扩散函数的一维 数值计算 问题，研究结果表明，在利用复指函数与 Bessel函数之间的展开关系和径向对称光瞳函数可 以展开为 Fourier级数的特性后，可以通过相继的一 维 Fourier变换和 Hankel变换计算来实现点扩散函 数的计算。实际计算表明，一维计算可以清晰地得出 像平面上某一确定方向的点扩散函数的强度分布情 况，同样可以作为工程研究的一个有效方法。此外， 对于实际的编程计算，若用传统的 Basic、Fortran和 c语言编程，既需要熟练地掌握所有语言的语法及编 程技巧，还需要对有关的算法有深刻的了解，这样要 耗去大量的人力和物力，影响工作效率。MATLAB 科学与工程计算软件包的功能强、效率高 ’ ，集应 用程序和图形于一体，表述形式简洁，不仅大大降低 了科学研究和工程设计工作者的语言知识要求，而且 还可以在计算机上直接输出计算结果和精美的图形 显示，已成为当前流行的计算软件。本文的一维计算 方法在 MATLAB上进行了实例编程计算，计算结果 以图形显示， 证明 住所证明下载场所使用证明下载诊断证明下载住所证明下载爱问住所证明下载爱问 了 MAT【 AB软件是进行科学与工 程研究的高效工具。 l 点扩散函数与 Hankel变换 空不变衍射受限成像系统的点扩散函数 h( ， Y )与出射光瞳函数 P(x，．y)的关系可表示为 ： (2i~Yi)一C『』 exp[磕( + )Id：cdy= F{P(Adi．~，Ad )} (1) 式中 ( ，Y )是像平面上的点；C一1／( 。d。d )，d。 是物面到光瞳面的距离，d 是 出瞳平面到像平面的 距离； 是波长； 一寺 ，Y一 ； ，Y是光瞳面上的 “^ “^ 点； 是虚数单位。公式(1)表明，点扩散函数是光瞳 函数的二维 Fourier变换。对于一般情况，很难求得 公式(1)积分的解析计算，通常是通过编程求出二维 离散的快速 Fourier变换作为解析解的近似。而计 算离散的二维 FFT就需在每个指标上连续地计算原 函数的一维 FFT ]，得出的结果是像平面上点扩散 函数强度的总的分布情况，要想清晰地观察某一特定 方向的点扩散函数的强度分布情况，还要做进一步的 工作。但值得注意的是，在实际的光学系统中，出射 光瞳(如透镜、光阑等)的形状通常是圆形的，因此公 式(1)中的光瞳函数通常是径向对称(或圆对称)的， 可以应用这种对称性来研究光瞳函数的二维 Fourier 变换 。 在公式(1)中，将直角坐标用极坐标表示。设 0 和 分别是光瞳面内和像平面中的方位，并令 — rcos0，Y—rsin0， 一pcos~，Y 一psin9，则公式(1)的 积分为： rr h(1D， )一CIIP(r， )K(r，O；p， )rdrdO (2) 式中 K(r， ；lD， )一exp[j2nprcos( 一 )] (3) 应用 Bessel函数将 K(r， ；lD， )展开，有 ： K(r，O；P， )一∑ J (2nrp)e川州 一J。(2nrlD)+ 2∑J”J (2nrp)cosEn(~一 )] (4) 式中 J (2nrp)是第一类 阶 Bessel函数。对于径 向对称的光瞳函数 P(r， )是方位角 0的周期函数， 因此 P(r， )可以展开成如下的 Fourier级数： P(r， )一∑c (r)e 枷 (5) 应用公式(4)和公式(5)，可以将公式(2)积分号内的 乘积项 P(r， )K(r，O；p， )表示为： P(r， )K(r， to)一∑ ∑ c (r)-， (2nr lD)× exp{ [枷 +”(to一 )]} (6) 将公式(6)代入公式(2)可得到： 维普资讯 http://www.cqvip.com 第4期 江月松等：点扩散函数的一维数值计算及其 MATI AB实现 4O7 h(p， )一C∑ ∑ H (J0) ( ) 式中 H (J0)和 ( )分别为： (J0)一j c，，I(r) (2丌rp)rdr 和 ( )一 e，， 下eJ c一 d 一{2丌 ．I 10 其中公式(8)是 阶 Hankel变换。。。。 的正交性，所以公式(7)可简化为： 得光瞳函数展开式(5)中的 Fourier系数 C， (r)。再由 公式(8)求得 H (J0)，然后由公式(12)求得点扩散函 数 h(J0， )，最后求得像面上点扩散函数的强度分布 Ih(J0， )I 。根据离散快速 Fourier算法 j一和 MAT一 ⋯ LAB软件中序列或向量元素下标从 1开始记数的考 虑，公式(5)中 Fourier变换系数 C (r)的离散值可写 为 · ： h(p。 )一21rC∑ e H (J0) (10) 对于实际的数值计算， 取有限值。公式(10)表明： 点扩散函数可以用一系列相继阶的 Hankel变换的 和来近似表示 ，它适用于一般的衍射受限像差的光瞳 函数。光瞳函数在极坐标下表示为 ： P(r， )一 丁(r， )exp[j志W(r， )]一 』exp[ w(r， )]光瞳内 (11) 1 0 光瞳外 式中 志一2丌／ ，志是波数；丁(r， )是振幅透射函数； w(r， )是波前像差函数，包括纵向和横向的离焦量。 造成像差的因素有多种，如透镜本身造成的波面畸 变、镜头内部的不均匀、表面的玷污、传输途中的大气 不稳定等因素。在实际问题中遇到的像差的光瞳函 数总是满足 P(r， )一P(r，一 )_9]，即 P(r， )是 0的 偶函数，由 Fourier 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 得知，公式(5)中的 Fourier 展开系数满足 C (r)一C一 (r)，则公式(5)的 Fourier 级数只由余弦项组成。此外由 Bessel函数的性质 (or)一(一1) J， ( )可知有 H一 一(一1) H ( )，这样点扩散函数公式(10)可表示为： h(J0， )一 2 7rCH。(J0)+4丌c J cos( )H (J0) i= 1 (12) 可见，由此得 到的点扩散 函数可 以由光瞳 函数的 Hankel变换求得。 2 数值计算及其实例 通常难以从解析上求得公式(8)的 Hankel变 换 ，必须采取离散数值计算的方法，为此必须首先求 苣 rP11 P1 2 1 I P21 P22 一 l j LPN1 PN2 P1L P 2L j P Nt (m — I，2，⋯ ，M ) (13) 式中 M 为在进行实际数值计算时所取的前 M 项， P 一P(r ， )，rf、 分别为： fr，一 r0e 卜 (Z一 1，2，⋯，N) 1 一( )q cq 2，⋯ 4 式中 N、L分别为光瞳面径向的离散采样点数和角 向采样点数；a为径 向采样间隔；rn为光瞳面矢径的 初始值。 计算出公式(13)中的 C (r)后便可计算公式(8) 中的 H (J0)，为此，采用 A．E．Siegman在 1977年 提出的准快 速 Hankel变换方 法来计 算。令 r— r。e ，J0一 e ，P0为像平面上矢径方 向采样的初始 值。这样，公式(8)的 Hankel变换积分可化为下列 形式的互相关积分： j三r ( )一I ( )j ( + )d (15) 式中 f ( )一rC， (r)，其对应的逆变换函数为 自 =pH (J0)；j ( + )=2~arpJ， (2 7rrp)。可以看出， 公式(15)是互相关 函数，根据互相关定理 ，它可 以根据一维快速 Fourier变换来计算，在实际的数值 计算中，计算是针对像平面上一系列的离散采样点逐 点进行的，即： 冉， ( ， )一FFr ~FFT[-f (x)~FFTJ ( + ， )]) (16) 式中 与 ， 相对应的径向离散值为： — expE(”一1)口]，(”一 1，2，⋯，N)(17) 在对公式(16)进行了 M 次、N 项的准快速 Hankel 维普资讯 http://www.cqvip.com 4O8 红外与激光工程 第33卷 变换之后，得到矩阵[H ]M ，而 H 一H (1D )。对 于不同方向，设置确定的 ‰，可以得到一系列的 COS (rnq： )值，应用公式(12)最终计算出在 个采样点上 的点扩散函数的强度 l h l 一l h(1D 。)l ，在实际的 工程应用中， 可以根据具体的需要来确定，从而计 算出不同方位的点扩散函数的强度分布情况。计算参 量 _N、r0、 和a可根据参考文献[10]的方法来确定。 将上述点扩散函数的一维计算过程在 MAT— LAB 6．0软件上进行了编程计算，计算了两种像差 的情况，一种是侧向离焦像差，光瞳函数为 P(r，臼)一 exp(ikW11 rcos0)(W11—0．5A)；另一种是彗差，光瞳 函数为 P(r，臼)一exp[ik(w11+w31 r。+wj1 r )reos0] (Wl1—0．5A、W31一一3．12 ，Wj1—2．60A)。计算 中 所用的参量是：N一128，L一16，M一16，rc】一0．08， 一1．0，a一0．0198。图 1是方位角 —0。时两种像差 的点扩散函数的情况，图 2是 —45。的情况，具体计 算时，只需在程序中设置一个 值，对应的计算结果 便会以图形方式显示出来。由图 1和图 2可见，不同 即 啷 Radius ofdifiraction disk on imaging plane／arb．unit 图 1 一0。的点扩散函数 Fig．1 Point spread function at 一0。 Radius of diffraction disk on imaging plane／arb．unit 图 2 =45。的点扩散函数 Fig．2 Point spread function at∞一45。 方位的点扩散函数的强度分布差别相当大，在实际的 工程应用中，可以根据需求计算任一方位的点扩散函 数的强度分布情况，为解决实际问题提供重要的参考 依据。 3 结 论 根据 Bessel函数的性质，将衍射受限光学成像 系统的点扩散函数二维 Fourier变换的计算变换为 一 维 Fourier变换和一维 Hankel变换的数值计算， 并用当今流行的科学与工程计算软件对侧向离焦和 彗差两种像差进行了实际编程计算，与二维 Fourier 变换数值计算方法相比，一维方法可以获得所需要的 任意方位的点扩散函数的强度分布，而二维计算在此 方面却有一定的限制。此外，本文的方法也可应用于 调制传递函数计算、无线电天线阵列分析、光干涉分 析以及光谱分析等领域。计算表明，使用 MATLAB 进行科学和工程计算，程序语言简洁，且大大降低了 使用者对数学基础和计算机语言知识的要求。 参考文献： [1] 苏显渝，李继掏．信息光学[M]．北京：科学出版社，1999．34 84． [2] Barakat R．The calculation of integrals encountered in optical diffraction theory，in The Computer in Optical Research[M]． New York：Springer—Verlag，1980．35—80． [3] 许波，刘征．Mat~b工程数学应用[M]．北京：清华大学出版社，2000． [4] 王末然．MATI AB 6．0与科学计算[M]．北京：电子工业 出版 社 ，2001． [5] Press W H，Teukolsky S A，Vetterling W T，et a1．Numerical recipes in C，the art of scientific computing．2nd ed．[M]．Cam— bridge：Cambridge University Press．1992． [6] 奚定平．贝塞尔函数[M]．北京：高等教育出版社；柏林：施普林 格，1998．62—92． [7] Jack D Gaskil1．Linear systems，Fourier transforms and optics [M]．New York：John Wiley＆ Sons，1978． [8] Born M，Wolf E．Principles of optics．7th edition[M]．London： Cambridge University Press．1999． [9] Robert R Shannon，James C Wyant．Applied optics and optical engineering，XI[M]．New York：Academic Press，Inc，1992． [10] Siegman A E．Quasi fast Hankel transform[J]．Optics letters． 1997，1(1)：l3 15． 维普资讯 http://www.cqvip.com