繆正祥，李重光，張中恒

（昆明理工大學(xué)理學(xué)院，昆明650500）

基于角譜理論的頻移插值法衍射計(jì)算

繆正祥，李重光*，張中恒

（昆明理工大學(xué)理學(xué)院，昆明650500）

在衍射分布的離散計(jì)算中，當(dāng)衍射距離增大時(shí)，滿足抽樣定理的抽樣點(diǎn)數(shù)會(huì)隨之增加，為了解決補(bǔ)零算法帶來的計(jì)算量大、計(jì)算機(jī)內(nèi)存溢出等問題，提出了一種頻移插值法，在不改變抽樣點(diǎn)數(shù)的情況下，獲得更為豐富的頻譜成分，計(jì)算每次插值后的光場(chǎng)分布，然后進(jìn)行線性疊加，從而得到完整的衍射分布。結(jié)果表明，該算法可以避開衍射距離增大必須增加抽樣點(diǎn)數(shù)這一離散計(jì)算限制條件，和傳統(tǒng)補(bǔ)零算法相比，計(jì)算量極大地減小，并降低了對(duì)計(jì)算機(jī)內(nèi)存的要求。

衍射；光柵；角譜；頻移；插值

引 言

衍射問題是光學(xué)中遇到的最困難的問題之一，在衍射理論中很少存在某種意義上可以認(rèn)為是嚴(yán)格的解［1-2］。根據(jù)角譜理論，光在自由空間中由衍射屏到觀察屏的傳播過程，在頻域中等效為一個(gè)理想低通濾波器［3］，在有限的觀察面上，只能接收部分頻率的衍射光波，離散計(jì)算時(shí)，需要選擇合適的離散點(diǎn)數(shù)對(duì)空間頻率離散化，當(dāng)衍射距離增大時(shí)，為獲得更為豐富的頻譜成分，必須增加抽樣點(diǎn)數(shù)以減小抽樣頻率間隔，但這樣計(jì)算量會(huì)成倍增加［4］。

本文中提出一種頻移插值方法，對(duì)物光波的頻譜進(jìn)行多次頻移計(jì)算，在不顯著增加計(jì)算量的前提下獲取更小頻率間隔下的頻譜。頻移過程是一個(gè)對(duì)物光波頻譜插值的過程，經(jīng)頻移插值后，即縮小了抽樣頻率間隔又不需要增加抽樣點(diǎn)數(shù)，確保衍射計(jì)算中不丟失太多的頻譜成分，實(shí)現(xiàn)衍射分布的正確計(jì)算。

1 空間頻率分析

以平面波為例，討論1維坐標(biāo)下，觀察面上各點(diǎn)的空間頻率隨空間位置的變化關(guān)系。如圖1所示，衍射面和觀察面相互平行，衍射面中心和觀察面中心在同一直線上，S為衍射面寬度，D為有效觀察面寬度，觀察面到衍射面的距離為z0，θ為平面光波與x軸的夾角。

Fig.1 Schematic of wave propagation

觀察平面上的光場(chǎng)分布可以看作許多不同方向傳播的單色平面波的線性疊加，每一平面波的振幅和相位都取決于相應(yīng)的角譜［5-6］，因此，觀察面上任意點(diǎn)沿x方向的空間頻率可表示為：

Fig.2 Spatial frequency changing with spatial position

圖2 中，波長(zhǎng)λ＝632.8nm，衍射面寬度S＝2mm，觀察面尺寸D＝10mm，“*”部分為z＝0.5mm時(shí)的空間頻率隨空間位置的變化曲線，“－”部分為z＝20mm時(shí)的空間頻率隨空間位置的變化曲線。衍射距離較小時(shí)，在觀察面尺寸范圍內(nèi)，能保留較多的高頻分量。隨著衍射距離的增大，高頻分量急劇減少，要完整地描述衍射分布，必須減小抽樣頻率間隔，這就意味著抽樣點(diǎn)數(shù)必須增加。通過頻移插值法不僅可以減小抽樣頻率間隔，而且能在不增加抽樣點(diǎn)數(shù)的情況下保留下更多頻譜成分［7］。

2 頻移插值算法

離散數(shù)值計(jì)算時(shí)，須對(duì)變量x和y離散化處理，取離散間隔分別為Δx，Δy。衍射面上沿x，y方向的抽樣點(diǎn)數(shù)分別為M和N，光波在z＝0處的復(fù)振幅分布u（r，t；0）＝u（rΔx，tΔy；0），根據(jù)離散傅里葉變換的定義，物光波的傅里葉頻譜為：

式中，Δfx，Δfy為x，y方向的抽樣頻率間隔；m，r＝根據(jù)衍射的角譜計(jì)算公式，衍射的角譜傳遞函數(shù)G（z0）可表示為：

衍射光波的頻譜為：

隨著衍射距離的增大，在已定尺寸的觀察面上，用同樣離散條件所能描述別的高頻分量減小，只有通過補(bǔ)零增加抽樣點(diǎn)數(shù)才能更好的滿足耐奎斯特抽樣定理。補(bǔ)零后將M擴(kuò)大P倍，變?yōu)镻×M，N擴(kuò)大Q倍，變?yōu)镼×N，大于M，N的部分用零補(bǔ)充［8］。則觀察面的抽樣點(diǎn)數(shù)為MP×NQ，其中MP＝P×M，NQ＝Q×N。

頻移插值法通過將A（m，n；0）在頻率域平移P×Q次后，再疊加得到A～（p，q；0），相當(dāng)于對(duì)A（m，n；0）進(jìn)行了P×Q次插值。頻率間隔分別為Δfx＝，設(shè)x和y方向的平移量分別為αΔfx和βΔfx（其中α＝0，1，…，P－1；β＝0，1，…，Q－1），則：

式中，∑表示將M×N大小的矩陣A線性疊加得到MP×NQ大小矩陣。根據(jù)傅里葉變換的性質(zhì)，在頻域的平移相當(dāng)于在空域乘以一個(gè)線性相位因子［9］：

令：

經(jīng)傳播距離z0后，衍射光波的頻譜為：

觀察面上的衍射光場(chǎng)分布為：

式中，IFFT為快速傅里葉逆變換（inverse fast Fourier transform，IFFT），Hα，β和G（m，n；z0）的大小為M×N，u（p，q；z0）的大小為MP×NQ，因此通過對(duì)M×N大小的矩陣作P×Q次快速傅里葉變換，疊加后可以得到觀察面上的衍射光場(chǎng)［10-11］，衍射圖像的大小變?yōu)镸P×NQ。而補(bǔ)零法是直接對(duì)MP×NQ的矩陣進(jìn)行兩次傅里葉變換得到衍射圖像，因此（10）式較補(bǔ)零法計(jì)算量大大減小，能夠在有限離散的觀察面上得到更完整的衍射分布。

利用頻移插值法計(jì)算時(shí)，有限離散點(diǎn)數(shù)的觀察面上可以得到的衍射光波頻譜寬度為：

頻移插值法能夠增加確定尺寸的觀察面上衍射光波的頻譜寬度，因此，能夠更加精確地計(jì)算出衍射光波的光場(chǎng)分布。

3 數(shù)值模擬

為驗(yàn)證上述分析，編寫計(jì)算程序，根據(jù)（10）式進(jìn)行模擬計(jì)算。光波長(zhǎng)λ＝632.8nm，衍射物體為128pixel×128pixel的透光方孔，對(duì)應(yīng)衍射面寬度S＝2.5mm，抽樣點(diǎn)數(shù)M＝N＝256，平移次數(shù)P＝Q＝2，觀察面尺寸10mm×10mm。改變衍射距離z0的大小，得到不同距離上衍射光場(chǎng)的振幅及相位分布如圖3所示。

Fig.3 Amplitude and phase distribution in different propagating distance

圖3 a、圖3b和圖3c中分別給出衍射距離為0.3m，0.6m和1m時(shí)衍射光場(chǎng)的振幅及相位分布，用于計(jì)算觀察面上光場(chǎng)分布的MP×NQ總數(shù)達(dá)到512pixel×512pixel。通過頻移插值法對(duì)256pixel× 256pixel大小的矩陣進(jìn)行8次傅里葉變換得到衍射光場(chǎng)分布，而補(bǔ)零算法需要進(jìn)行2次512pixel× 512pixel大小的傅里葉變換，由（12）式可知，頻移插值算法能夠在有限離散點(diǎn)數(shù)的觀察面上得到更多頻率成分的衍射光波。因此，頻移插值法計(jì)算量較小，計(jì)算更為精確（根據(jù)數(shù)值模擬的條件編寫計(jì)算程序，利用內(nèi)存為4G的計(jì)算機(jī)進(jìn)行實(shí)驗(yàn)計(jì)算，補(bǔ)零算法的計(jì)算時(shí)間約為0.058453s，頻移插值法約為0.023571s，計(jì)算時(shí)間縮短近1.5倍）。

4 結(jié) 論

基于光傳播的角譜理論，推導(dǎo)出遠(yuǎn)距離光波衍射的數(shù)值計(jì)算方法，該方法能夠較好地滿足耐奎斯特抽樣定理。通過頻移插值，能夠在有限抽樣的觀察面上獲得更多頻率成分的衍射光波，可以更加準(zhǔn)確地計(jì)算出衍射光場(chǎng)的空間分布情況。利用頻移插值算法能夠?qū)⒋缶仃嚨母道锶~變換計(jì)算轉(zhuǎn)化成多個(gè)小矩陣的傅里葉變換計(jì)算，大大減小了計(jì)算量、縮短了計(jì)算時(shí)間。

［1］ WOLF B.Principles of optics［M］.7th ed.Beijing：Publishing House of Electronics Industry，2005：342（in Chinese）.

［2］ GOODMAN JW.Introduction to fourier optics［M］.3rd ed.Beijing：Publishing House of Electronics Industry，2011：46-61（in Chinese）.

［3］ LI JCh，XIONG B H.Information optics theory and calculation［M］.Beijing：Science Press，2009：93-137（in Chinese）.

［4］ ZHANG BW，MA Y，ZHANG PP，et al.Simulation of Gaussian laser standing wave based on diffraction by straight edge［J］.Laser Technology，2012，36（6）：810-813（in Chinese）.

［5］ LIJCh.FFT computation of angular spectrum diffraction formula and its application in wavefront reconstruction of digital holography［J］.Acta Optica Sinica，2009，29（5）：1163-1167（in Chinese）.

［6］ MUFFOLETTO R P，TYLER JM，TOHLINE JE.Shifted Fresnel diffraction for computational holography［J］.Optics Express，2007，15（9）：5631-5640.

［7］ KOZACKI T，F(xiàn)ALAGGIS K，KUJAWINSKA M.Computation of diffracted fields for the case of high numerical aperture using the angular spectrum method［J］.Applied Optics，2012，51（29）：7080-7088.

［8］ SU X Y，LIJT.Information optics［M］.Beijing：Science Press，1999：34-47（in Chinese）.

［9］ LOBAZ P.Reference calculation of light propagation between parallel planes of different sizes and sampling rates［J］.Optics Express，2010，19（1）：32-39.

［10］ LIU Y J.Diffraction field of orthogonal cosine-gratings［J］.Laser Technology，2013，37（1）：114-117（in Chinese）.

［11］ CHANG Sh，SANG Zh W，ZHUANG L.Calculation and analysis of spherical aberration on Gaussian image plane based on diffraction theory［J］.Laser Technology，2012，36（3）：405-407（in Chinese）.

Diffraction computation with frequency shift interpolation based on theory of angular spectrum

MIAO Zhengxiang，LI Chongguang，ZHANG Zhongheng
（Faculty of Science，Kunming University of Science and Technology，Kunming 650500，China）

In the discrete numerical calculation of diffraction distribution，the total number of sampling points following sampling theorem increases with the augment of propagating distance.Although the traditional zero-padding method resolves this problem，the calculation load increases inevitably and PC’s memory can not afford.A novel frequency shift interpolation，i.e.，utilizing frequency spectrum calculation after shift in Fourier domain，was proposed.In the method，more sufficient spectrum components can be got without the increase of sampling number.A complete diffraction distribution is accomplished by splicing the light field obtained by spectrum of each shifted interpolation.The results show that the proposed method successfully evades the great sampling number caused by larger diffracting distance.Comparing with traditional zero-padding method，calculation load is decreased notably and the requirement to memory in numerical calculation is depressed.

diffraction；gratings；angular spectrum；shift in frequency domain；interpolation

O436.1

A

10.7510／jgjs.issn.1001-3806.2014.01.019

1001-3806（2014）01-0087-04

云南省應(yīng)用基礎(chǔ)研究基金資助項(xiàng)目（2008F04OM）

繆正祥（1986-），男，碩士研究生，主要從事信息處理方面的研究。

*通訊聯(lián)系人。E-mail：2544248546＠qq.com

2013-04-23；

2013-05-28