繆正祥，李重光，張中恒

（昆明理工大學理學院，昆明650500）

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

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

（昆明理工大學理學院，昆明650500）

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

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

引 言

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

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

1 空間頻率分析

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

Fig.1 Schematic of wave propagation

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

Fig.2 Spatial frequency changing with spatial position

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

2 頻移插值算法

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

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

衍射光波的頻譜為：

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

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

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

令：

經傳播距離z0后，衍射光波的頻譜為：

觀察面上的衍射光場分布為：

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

利用頻移插值法計算時，有限離散點數的觀察面上可以得到的衍射光波頻譜寬度為：

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

3 數值模擬

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

Fig.3 Amplitude and phase distribution in different propagating distance

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

4 結 論

基于光傳播的角譜理論，推導出遠距離光波衍射的數值計算方法，該方法能夠較好地滿足耐奎斯特抽樣定理。通過頻移插值，能夠在有限抽樣的觀察面上獲得更多頻率成分的衍射光波，可以更加準確地計算出衍射光場的空間分布情況。利用頻移插值算法能夠將大矩陣的傅里葉變換計算轉化成多個小矩陣的傅里葉變換計算，大大減小了計算量、縮短了計算時間。

［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，FALAGGIS 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

云南省應用基礎研究基金資助項目（2008F04OM）

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

*通訊聯系人。E-mail：2544248546＠qq.com

2013-04-23；

2013-05-28