Optics Photonics Japan 2009 講演予稿集, pp.284–285 (25aF9), 新潟 (2009.11.25)

自由空間における光波伝搬シミュレーションのためのシフト角スペクトル法
The shifted angular spectrum method for simulation of freespace propagation of wave fields○松島恭治○Kyoji Matsushima

関西大学 電気電子情報工学科
Department of Electrical and Electronic Engineering, Kansai University
E-mail: matsu@kansai-u.ac.jp
A novel method for off-axis numerical propagation of optical wave fields is proposed. The shiftedFresnel method has been presented for off-axis propagation, but the method suffers from a strong
aliasing error in cases of short distance propagation. The proposed new method is based on the angular
spectrum method that features a rigorous calculation of Rayleigh-Sommerfeld equation without anyapproximation, and has the advantage of producing no aliasing error unlike the shifted-Fresnel method.
１．はじめに近年のディスプレイ用計算機合成ホログラムやディジタルホログラフィの進展により，自由空間における伝搬計算を正確かつ高い自由度で行う数値計算技術が重要となりつつある．イメージセンサの大型化・高精細化やホログラムの大型化により，特に伝搬元の光波分布のサンプリング窓からずれた位置で計算結果を高速高精度で得る技術が重要になりつつある．この様な軸外れの計算を行うためには，従来，伝搬元のサンプリング窓をゼロ拡張する手法が用いられていたが，伝搬先のサンプリング窓のシフトが大きい場合，サンプリング数が巨大になる問題があった
1)．一方近年，一度のフーリエ変換で行うフレネル回折計算を Scaled FFT と呼ばれる手法で可変ピッチ化したシフトフレネル法 2)が提案されており，軸外れ位置での回折像を計算できることから計算機のメモリに入りきらないような大規模な光波複素振幅分布の伝搬計算にも用いられている 3-5)．しかし，シフトフレネル法には伝搬距離が短い場合に強度のエイリアシング誤差が生じ，最悪の場合には全く解が得られないという問題点がある．
数値的な回折伝搬計算には種々の方式が知られているが， 2 回の FFT による高速畳み込みを用いて計算を行う角スペクトル法 6)はフレネル条件に拘束されず，原理的に伝搬距離の制約がないため数値計算用としてはもっとも有利である．著者は昨年，角スペクトル法において伝搬距離が長い場合に計算結果に生じる高周波ノイズ様の誤差を，帯域を制限することにより除去する手法を提案した 7)．本研究ではこの手法を発展させ，軸外れ位置での回折像が計算できるシフト角スペクトル法を提案する．この手法には伝搬距離の制約がほとんどなく，シフトフレネ
ル法で生じるようなエイリアシング誤差も起きないため，軸外れの回折計算が必要な問題に幅広く応用できると思われる．

G (u , v;0)  F  g ( x, y; 0)   g ( x, y;0) exp[ i 2 (ux  vy )]dxdy,H (u , v; z0 )  exp[i 2 wz0 ],w  ( 2  u 2  v 2 )1/ 2である．F{}はフーリエ変換を表し，求められる回折波 g (x, y; z0)=F1{G(u, v; z0)}は波動方程式の解となっている．この手法は角スペクトル法と呼ばれ，フレネル回折式は上式の放物面近似式になっている．
この角スペクトル法の式に Fig.1 に示すように座標変換を行う．すなわち， x̂  x  x0 ， ŷ  y  y0 とすることにより，gˆ ( xˆ , yˆ ; z0 )  g ( xˆ  x0 , yˆ  y0 , z0 )  g ( x, y;0)h(xˆ  x0  x, yˆ  y0  y; z0 ) dxdy

(x, y)( xˆ , yˆ )ẑz0z

Source sampling
window
Fig.1 Theoretical model for off-axis
numerical propagation and definition ofcoordinates systems.となるため，畳み込みの定理を用いて，gˆ ( xˆ , yˆ ; z0 )  F 1{Gˆ (u , v; z0 )},Destination
sampling
window

x0, y0Sx Sy２．シフト角スペクトル法
伝搬元である(x, y, 0)平面上の光波複素振幅分布 g(x, y; 0)を距離 z
伝搬したときの回折像は，Rayleigh-Sommerfeld の回折式の畳み込み形式で g(x, y; z) = g(x, y; 0)*h(x, y; z0)として与えられ，畳み込みの定理を用いて G(u, v; z0) = G(u, v; 0)H(u, v; z0)と表される．ここで，Gˆ (u , v; z0 )  G (u , v,0) H (u , v; z0 ) exp[i 2 ( x0u  y0 v)]

と変形することができる．従って，周波数応答関数 Hˆ (u , v; z0 )  H (u , v; z0 ) exp[i 2 ( x0 u  y0 v)] を用いることにより，伝搬元の座標系から，軸外れの伝搬先の座標系(シフトした座標系)での回折波を求めることができる．
しかし，この新しい周波数応答関数をサンプリング間隔 u ， v で離散化して数値計算に用いた場合，特定の領域以外では周波数応答関数自体のエイリアシングが生じて意味のある結果が得られない．そのため，エイリアシングが生じない範囲で入力光波の空間周波数帯域を制限しなければならない．この制限帯域は，1 次元では()()ulimit u  ulimitS x  x0()()ulimit  u  ulimit   S x  x0  S x ，()()ulimit u  ulimit
x0   S x()ただし ulimit




2
 x 1z02  10

2u


1/ 2 1となる．ここで Sx = (2u)は円状離散畳み込みを直線畳み込みに変換する拡張処理を行った場合のサンプリング領域の大きさである．2 次元での帯域制限領域は，厳密には複雑な楕円領域の組み合わせとなるが，たいていの場合 1 次元の結果を水平垂直にそのまま組み合わせても良い近似となっている．
Shifted Fresnel

Shifted angular spectrum

SNR [dB]

３．1 次元での計算例と計算誤差(this work)
提案する伝搬計算法を検証するため，
z0 = 50Sx
サンプリング点数 Nx = 1024 の 1 次元で数値計算を行った．この計算ではサンプ
Aliasing error
リング間隔x = 2，  = 532[nm]とした．
z0 = 5Sx
従ってサンプリング窓の大きさは Sx =Nxx である．
Sx (Source sampling size)
Fig.2 は伝搬元サンプリング窓の中央
Fig.
2One-dimensional
amplitude
distributionof the wave field diffractedに幅 Sx /10 の矩形開口がある場合の回
by a rectangular aperture with a width of Sx/10.
折波の計算結果である．回折波は伝搬元サンプリング窓からシフト量を変えて 4 回計算し，それを40貼りあわせて 4Sx の領域の計算結果を得ている．また左右対称な結果となるため，ここでは x 軸の正の側には提案法
30
Shifted angular spectrumで，また負の側には既存のシフトフレネル法で計算した結(this work)果を示している．伝搬距離が z0 = 50Sx と長い場合はどちら20もほぼ同じ結果となっているが，z0 = 5Sx ではエイリアシング
誤差と思われる誤ったピークがシフトフレネル法で現れて10いる．それに対して，提案法では一切そのよう問題は生じShifted Fresnelない．
伝搬距離を変え，回折積分によって得た厳密解との00.1
110100S/N 比として両手法の誤差を比較した結果を Fig.3 に示す．
Propagationdistance,z[S]0x
シフトフレネル法では伝搬距離 7Sx 以下程度で急激に S/N
Fig.3 Precision in calculation of diffraction.
比が悪化するのに対して，提案手法では全領域で良好な結果を保っていることがわかる．
４．おわりにサンプリング窓の拡張を行わずに軸から外れた位置での回折波を求める新しい伝搬計算法として，シフト角スペクトル法を提案した．同様の手法としてシフトフレネル法が提案されているが，シフトフレネル法には伝搬距離が短い場合にエイリアシングによる多重像が生じるという大きな問題がある．提案手法ではそのような問題は生じず，また計算時間も従来の角スペクトル法とほぼ同じである．伝搬距離による制約はほぼ全くないため，提案手法はCGH やデジタルホログラフィに有用であり，実際にこの手法を用いて 40 億サンプリング規模の大規模な光波伝搬計算が行われている 8)．
文献1) T. Nakatsuji and K. Matsushima: Appl. Opt. 47, D136 (2008).
2) R. P. Muffoletto, J. M. Tyler and J. E. Tohline: Opt. Express 15, 5631(2007).
3) K. Matsushima and S. Nakahara: SPIE Proc. #7233, 72330E(2009).
4) 松島, 中原: 3 次元画像コンファレンス 2009 講演論文集, 61(2009).
5) K. Matsushima and S. Nakahara: Appl. Opt. 48, No. 34, (2009) in press.
6) J. W. Goodman: Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996), chap. 3.10.
7) 松島: Optics Photonics Japan 2008 講演予稿集, 324-325(2008).
8) 松島, 中原: Optics Photonics Japan 2009 講演予稿集, (2009).