基于Dini展開的高階Hankel變換及其在光束傳輸中的應(yīng)用*

游開明 林燕玲 王友文 陳列尊 戴志平 陸世專

(湖南衡陽師范學(xué)院物理與電子信息科學(xué)系,衡陽 421002)

(2013年2月2日收到;2013年3月30日收到修改稿)

1 引言

由于Hankel變換(Hankel transform,HT)在物理學(xué)的很多分支(如光學(xué)、聲學(xué)、電磁學(xué)、分子生物學(xué)等)有著廣泛的應(yīng)用[1-4],引起了很多科學(xué)工作者的關(guān)注,并對它的數(shù)值計算算法進行了研究.1977年,Siegman[5]首先對徑向坐標(biāo)抽樣點用指數(shù)形式量化,將HT化為離散形式,用快速傅里葉變換(FFT)進行計算,得出了準(zhǔn)快速Hankel變換算法.此后,出現(xiàn)了很多用于物理學(xué)各個方面的HT算法,如文獻[6—18].文獻[9—11]把HT看作是一個二變量函數(shù)的簡單積分的傅氏變換,先用一般方法計算簡單積分,再利用FFT得到HT.2003年,Markham和Conchello[12]總結(jié)了計算HT的6種方法,并把它們用于計算振蕩函數(shù).文獻[13]和[14]分別用正弦和余弦變換進行快速HT及非均勻快速HT等.2007年,Cerjan[15]把變換函數(shù)在單位區(qū)間上用Zernike多項式表示,再根據(jù)Zernike多項式基礎(chǔ)集與Bessel函數(shù)展開的對偶關(guān)系進行HT,其展開的系數(shù)采用Simpson規(guī)則或Gauss-Kronrod求積法計算,使HT變換在精度上有所提高.1998年Yu等[16]對徑向坐標(biāo)按零階Bessel函數(shù)零點分布的比例進行離散,將變換函數(shù)做Fourier-Bessel級數(shù)展開,實現(xiàn)了準(zhǔn)離散HT(quasi-discrete HT,QDHT),2004年Guizar-Sicairos和Gutie′rrez-Vega[17]又把它推廣到p階HT(the p th-order quasi-discrete HT,p QDHT),這類方法利用變換矩陣使計算變得簡單,因而速度也較快,由于能量守恒,逆變換精度很高,不依賴輸入函數(shù),能適應(yīng)反復(fù)變換,具有很高的實用價值.遺憾的是基于Fourier-Bessel級數(shù)展開HT方法處于條件收斂,不能計算空間域和空間頻率域端點處的變換值,同時其0階HT也不能計算對稱軸上點的變換值,不能直接用于研究傳輸問題中對稱軸上的特性.文獻[18]給出了一種基于Dini級數(shù)展開的0階準(zhǔn)離散HT算法,由于Dini級數(shù)展開比Fourier-Bessel級數(shù)展開收斂更快,精度更高,對于0階HT還能計算對稱軸上的值,所以作者把它應(yīng)用到了解決光束非線性傳輸問題.然而物理學(xué)中更一般的傳輸問題是高階的,如高階Bessel光束傳輸問題等[19-21],又基于Dini級數(shù)展開0階HT與高階HT畢竟還有一定區(qū)別,因此本文開展基于Dini級數(shù)展開的高階p(p＞0)準(zhǔn)離散HT(the p th-order quasi-discrete Hankel transform based on Dini series expansion,p DQDHT)算法的研究,并給出它在光束傳輸中的應(yīng)用.

2 算法的推導(dǎo)

很多有限區(qū)域上圓對稱物理量分布問題,數(shù)值方法上使用HT更有效.歸一化坐標(biāo)的整數(shù)p(p＞0)階Hankel變換和逆變換分別表示為[7]

對于光束傳輸,式中r表示歸一化空間徑向坐標(biāo),ρ表示歸一化空間頻率徑向坐標(biāo),Jp表示第一類p階Bessel函數(shù),b和β分別表示空間變量和空間頻率變量的截斷半徑,S表示帶寬乘積bβ的2π倍,即

對于第二類齊次邊界條件,第一類p階Bessel函數(shù)滿足如下關(guān)系[18]

式中 αn(n=1,2,···)是 p 階 Bessel函數(shù)導(dǎo)數(shù) J′p(x)的正實根,δnm為克羅內(nèi)克符號,m=n時為1,否則為0.

這里考慮p＞0的情況,將方程(1),(2)中的變換函數(shù) f(br)和g(βρ)展開為Jp的級數(shù),分別有

式中的系數(shù) fn和gn由方程(4)確定,再結(jié)合方程(1)和(2),分別有

上述級數(shù)展開稱為Dini級數(shù)展開[22].

將歸一化徑向坐標(biāo)變量r和ρ分別量化為

綜合方程(5)到(8),得到p階離散Hankel變換對為

這里N為數(shù)值計算點數(shù).為了減少方程(10)和(11)右邊求和符中數(shù)乘次數(shù),分別定義向量G(m)和F(n)為

于是方程(10)和(11)簡化為如下形式

其中Cnm為N×N變換矩陣C的矩陣元,表示為

根據(jù)方程(14)—(16),容易看出矩陣C的一些特性：由方程(16)知C是實對稱矩陣,滿足C=CT;將方程(14)代入(15)或反過來代換,發(fā)現(xiàn)F(n)和G(m)均為自身變換,說明正、逆變換毫無誤差時還需要C是正交矩陣,滿足CCT=I(I表示單位矩陣),即C的行列式det[C]=±1;由方程(16)還知C是單參數(shù)S的函數(shù),C=C(S),S不同,矩陣C中元素的值不同,所以如何選擇S以使C成為正交矩陣對提高該數(shù)值算法的精度非常重要.

根據(jù)方程(1),(2)和(9)應(yīng)取S=αN,然而我們通過數(shù)值驗證,得知矩陣C并不是很好的正交矩陣.若取S=jN+1(其中 jN+1為Jp(x)的第N+1個零點,滿足αN＜ jN+1＜αN+1),矩陣C非常接近正交矩陣,這是因為當(dāng)方程(1)和(2)中的r=1和ρ=1時,其中的變換核函數(shù)為Jp(S)=0,可以使變換后的g(ρ)和 f(b)的值為0,所以本算法一般的數(shù)值計算中取S=jN+1.若要進一步提高精度,可利用矩陣C的性質(zhì)對 jN+1進行修正來得到S,其修正表達式為

利用矩陣C不僅能使在一次HT中減少N(3N-4)次復(fù)數(shù)或?qū)崝?shù)乘法運算,若將C存儲起來,還可避免下一次變換時對數(shù)據(jù)的重復(fù)計算,這對分步傳輸問題非常重要,可以極大地提高程序的執(zhí)行速度.

這樣p DQDHT算法為：先由方程(13)將 f(br)轉(zhuǎn)換為F(n),乘以矩陣C后再按方程(12)將G(m)轉(zhuǎn)換為g(βρ).逆變換時只要將方程(12)和(13)右邊的ρ和b交換即可.

3 算法的測試與驗證

我們用C++語言編程對p DQDHT算法的正確性與精度進行測試,根據(jù)C為實對稱矩陣,對數(shù)據(jù)的存儲采用三角矩陣形式,可以分別節(jié)省N(N-1)/2個雙精度型存儲單元的存儲空間,這用C++語言編程是非常容易實現(xiàn)的.對于Bessel函數(shù)及其導(dǎo)數(shù)的零點很容易先用公式估算出來[23],然后用牛頓迭代法分別解方程Jp(x)=0和 J′p(x)=0.5(Jp-1(x)-Jp+1(x))=0 很快得到精確值.注意,在本文后面的敘述中r和ρ表示實際坐標(biāo),前后的關(guān)系分別是r=br和ρ=βρ,同時取S=jN+1.

首先我們用連續(xù)函數(shù) f(r)=r2exp(-πr2)作為輸入函數(shù)對p DQDHT算法進行測試,它在[0,∞)上的2階HT的精確表達式為

以 g(ρ)作為基準(zhǔn),用 g?(ρ)表示 p DQDHT 算法的變換值,在空間頻率域定義絕對誤差為|g-g?|,并分別用Max E表示最大絕對誤差,Min E表示最小絕對誤差,Mean E表示平均誤差|g-g?|/N,取b=β=(S/2π)1/2,2階HT的計算結(jié)果如表1所示.

表1 用p DQDHT計算f(r)的誤差隨抽樣點數(shù)N的變化

可見,對于連續(xù)函數(shù)當(dāng)抽樣點數(shù)N很小時p DQDHT算法就有很高的精度,進一步說明Dini展開收斂很快.

圖1表示函數(shù) f(r)的精確值與經(jīng)過40次2階p DQDHT后的結(jié)果比較,計算中空間域截斷半徑b=4,抽樣點數(shù)N=100.圖中的實線表示精確值,圓點線表示2階HT的數(shù)值結(jié)果,它們幾乎完全重合.在PIV 1.5 GHz的CPU,512 MB內(nèi)存的PC機上執(zhí)行程序,p DQDHT算法經(jīng)過40次變換共花費的時間為0.03 s.可見,p DQDHT算法經(jīng)過多次變換仍保持極高的精度,而且程序執(zhí)行速度快.

圖1 函數(shù) f(r)的精確值與經(jīng)過40次p DQDHT后的結(jié)果比較,其中p=2,b=4,N=100

其次,我們用top-hat函數(shù)作為輸入函數(shù),它的定義是 f(r)=rp[H(r)-H(r-b)],這里H(r)為Heaviside階躍函數(shù),其p階HT的精確表達式為

分別計算經(jīng)過一次p DQDHT和p QDHT與(19)式隨空間頻率ρ變化絕對誤差的對數(shù),再進行誤差比較,它們的動態(tài)比較結(jié)果如圖2所示.計算中空間域截斷半徑b=5,抽樣點數(shù)N=200,空間頻率域截斷半徑β在20.025附近.圖2是p=3時p DQDHT和p QDHT兩種算法的結(jié)果比較,圖中整個空間頻率范圍內(nèi) p DQDHT算法得到的曲線(實線)均位于p QDHT算法得到的曲線(圓點線)的下方,表示p DQDHT算法比p QDHT算法具有更高的精度.圖2還表現(xiàn)出基于Dini級數(shù)展開的p DQDHT算法與基于Fourier-Bessel級數(shù)展開的p QDHT算法的不同特性,前者誤差較大的地方主要出現(xiàn)在低頻端較小的區(qū)域內(nèi),在其他的地方保持較高的精度,這對再進行一次HT變換到空間域時減少Gibbs現(xiàn)象是非常有益的;后者在低頻端精度較高,隨著頻率的增大,誤差越來越大,特別是在最高頻率處誤差最大,因此p DQDHT算法更具有優(yōu)勢.

圖2 “top-hat”函數(shù)經(jīng)過一次 p DQDHT和 p QDHT后與精確值絕對誤差的對數(shù)隨空間頻率半徑變化的比較,其中p=3,b=5,N=200

4 p DQDHT在光束傳輸中的應(yīng)用

為了研究p階準(zhǔn)離散Hankel變換在光束傳輸中的應(yīng)用,我們考慮圓對稱4階Bessel光束f(r)=J4(ktr)exp(i4φ+i kzz)通過凸透鏡后的聚焦演變.讓光束入射到半徑R=4 mm,焦距 f=0.5 m,位于平面z=0的凸透鏡上,并使光束在透鏡光軸上對稱分布,波數(shù)k與橫向波數(shù)kt和縱向波數(shù)kz滿足關(guān)系 k2=(2π/λ)2=kt2+kz2.

光束通過透鏡后為非傍軸傳輸,數(shù)值方法上采用分步傳輸方法,所以運用p DQDHT解決類似于上述光束傳輸問題的計算步驟可為：

1)選取b,p和N,計算程序中所需的Bessel函數(shù)及其導(dǎo)數(shù)的零點,將Bessel函數(shù)導(dǎo)數(shù)的零點存儲于一維數(shù)組中;

2)計算矩陣C,并使用指針數(shù)組按三角形式存儲起來;

4)計算z=0處初始光束 f(r)=J4(ktr)的值,乘以光束從透鏡前表面傍軸傳輸?shù)胶蟊砻娴淖饔靡蜃觘xp[i kr2/(2f)]后按復(fù)數(shù)形式存儲于一維數(shù)組中(在C和C++中偶地址存實部,奇地址存虛部);

5)選取傳輸步數(shù)M,計算傳輸步長Δz=z/M,再進行分步傳輸計算：

for j=0,···,M

對 f(r)進行一次p DQDHT得空間頻率域上的g(ρ);

再進行一次p DQDHT得空間域上的 f(r),存儲光強徑向分布數(shù)據(jù),end.

計算中取光的波長λ=632.8 nm,kt=19858.32 m-1,徑向截斷半徑b=4 mm,徑向取樣點數(shù)N=256,縱向傳輸距離z=0.75 m,傳輸步數(shù)M=300.按照這些參數(shù)利用p DQDHT算法編程計算,得到會聚Bessel光束光強I(r,z)的傳播演變?nèi)鐖D3所示.

圖3 4階會聚Bessel光束通過透鏡后光強在焦平面附近的傳播演變

從圖3看出,會聚Bessel光束通過凸透鏡后在焦平面附近三處會聚,會聚的峰值并不出現(xiàn)在對稱軸上,而是在以對稱軸為中心半徑不同的三個圓環(huán)上,這是高階貝塞爾光束為空心光束的原因.進一步用數(shù)值方法求三個會聚峰值的半徑和位置,得到透鏡焦平面前、焦平面上和焦平面后三個峰值光環(huán)的半徑分別為r=0.062647 mm,0.996897 mm和0.110658 mm,會聚平面的位置分別為z=0.38 m,0.5 m和0.72 m,它們的截面光強分布如圖4所示.這些結(jié)果與FFT算法得到的結(jié)果是一致的,特別是在焦平面上的結(jié)果還與文獻[20]用幾何方法得到的會聚圓環(huán)的半徑r=fkt/kz很好地相符.

圖4 4階Bessel光束通過透鏡后不同會聚處的截面光強分布 (a)z=0.38 m;(b)z=0.5 m;(c)z=0.72 m

5 結(jié)論

基于Dini級數(shù)展開,導(dǎo)出了高階Hankel變換的離散表達式,稱為p DQDHT算法,能實現(xiàn)復(fù)數(shù)或?qū)崝?shù)形式的HT.用向量表示,可把它們分別寫成為輸出向量是變換矩陣與輸入向量的乘積,計算變得更簡單,計算速度更快.討論了變換矩陣C的性質(zhì),以及使它成為正交矩陣與兩個截斷半徑乘積S的關(guān)系,定性分析了S的取值,得到S=jN+1時C接近于正交矩陣.為了進一步提高該算法的精度,定量給出了對S進行修正的表達式,同時用數(shù)值方法給出了其中的可選參數(shù)k與抽樣點數(shù)N的關(guān)系.

在算法的實施上,事先計算出Bessel函數(shù)導(dǎo)數(shù)的零點進行靜態(tài)存儲,計算變換矩陣C采用三角形式存儲,以免下一次變換對數(shù)據(jù)重復(fù)計算,進一步提高了計算速度.算法中的正變換和逆變換雖然用了兩種表達式,實際上它們的區(qū)別只是截斷半徑不同,所以正變換和逆變換可以采用相同的程序段,只是在每次變換后將兩個截斷半徑的值進行交換即可,使程序變得簡單.

通過兩個不同的輸入函數(shù)對p DQDHT算法的測試和光束分步傳輸?shù)膽?yīng)用實例,結(jié)果表明：p DQDHT算法經(jīng)過多次變換而不丟失精度,特別適合分步傳輸?shù)膱龊?這是很多Hankel變換算法做不到的;在整個變換區(qū)域內(nèi),p DQDHT算法都具有很高的精度,特別是在空間頻率域的高頻端精度更高,在空間域近軸端精度更高,這正符合一般實際應(yīng)用的需要;程序執(zhí)行的速度與一般的快速算法相當(dāng).

我們相信,借助本文p DQDHT算法在光束分步傳輸中的應(yīng)用方法,p DQDHT算法在解決物理學(xué)各個分支類似的問題中必將發(fā)揮重要作用.

[1]Du G H 1988 Acta Phys.Sin.37 769(in Chinese)[杜功煥1988物理學(xué)報37 769]

[2]Yan CC,Xue G G,Liu C,Gao SM 2005 Acta Phys.Sin.54 3058(in Chinese)[閏長春,薛國剛,劉誠,高淑梅2005物理學(xué)報54 3058]

[3]Yang X J,Zang WP,Tian JG,Liu ZB,Zhou WY,Zhang CP,Zhang G Y 2005 Acta Phys.Sin.54 2735(in Chinese)[楊新江,臧維平,田建國,劉智波,周文遠,張春平,張光寅2005物理學(xué)報54 2735]

[4]Chen G B,Wang H N,Yao JJ,Han Z Y,Yang SW 2009 Acta Phys.Sin.58 1608(in Chinese)[陳桂波,汪宏年,姚敬金,韓子夜,楊守文2009物理學(xué)報58 1608]

[5]Siegman E 1977 Opt.Lett.1 13

[6]Agrawal G P 1981 Opt.Lett.6 171

[7]Magni V,Cerulle G,De Silvestri S 1992 J.Opt.Soc.Am.A 9 2031

[8]Agnesi A,Reali GC,Patrini G,Tomaselli A 1993 J.Opt.Soc.Am.A 10 1872

[9]Ferrari JA 1995 J.Opt.Soc.Am.A 12 1812

[10]Ferrari JA,Perciante D,Dubra A 1999 J.Opt.Soc.Am.A 16 2581

[11]Perciante CD 2004 J.Opt.Soc.Am.A 2 1911

[12]Markham J,Conchello A J2003 J.Opt.Soc.Am.A 20 621

[13]Knockaert L 2000 IEEETrans.Signal Process.48 1695

[14]Liu Q H,Zhang ZQ 1999 Appl.Opt.38 6705

[15]Cerjan C 2007 J.Opt.Soc.Am.A 24 1609

[16]Yu L,Huang M,Chen M,Chen W,Huang W,Zhu Z 1998 Opt.Lett.23 409

[17]Guizar-Sicairos M,Guti′errez-Vega JC 2004 J.Opt.Soc.Am.A 21 53

[18]You K M,Wen SC,Chen L Z,Wang Y W,Hu Y H 2009 Chin.Phys.B 18 3893

[19]Wang Z,Gao C Q,Xin JT 2012 Acta Phys.Sin.61 124209(in Chi-nese)[王錚,高春清,辛璟燾2012物理學(xué)報61 124209]

[20]Gutie’rrez-Vega JC,Rodr′?guez-Masegosa R,Cha’vez-Cerda S 2003 Pure Appl.Opt.5 276

[21]Zhang Q A,Wu FT,Zheng W T,Pu JX 2011 Sci.Sin.Phys.Mech.Astron.41 1131(in Chinese)[張前安,吳逢鐵,鄭維濤,蒲繼雄2011中國科學(xué)：物理學(xué)·力學(xué)·天文學(xué)41 1131]

[22]Arfken G B,Weber H J 2001 Mathematical Methods for Physicists(California：Harcourt-Academic)

[23]Wolfram S 1991 Mathematica：A System for Doing Mathematics by Computer(Massachusetts：Addison-Wesley)