第 7 回関西大学先端科学技術シンポジウム講演集, 47-51(2003.1)

波面展開法を用いた光波伝搬計算アルゴリズム—自由空間における光波伝搬のシミュレーションとその応用—
松島恭治
関西大学工学部電気工学科，関西大学 HRC

A novel method for calculating complex-amplitude distribution of the light wave diffracted by a tilted
aperture as well as the light wave diffracting onto tilted screen is presented. Since this method is based on
rotational transformation of angular spectrum of plane waves and is analytically rigorous, the method
does not involve any restriction on the propagation distance. Furthermore, the computation time of thismethod is comparable with conventional algorithms of diffraction between parallel planes, because only
twice FFT and an interpolation in Fourier spectrum are necessary for its numerical calculation. The basic
theory of the method and its application to digitally synthetic holograms is reported.

1. はじめに
近年の微細加工技術の発展に伴い，光の回折効果を利用したビーム成形素子や分岐素子が作製可能となってきた．これらの回折光学素子の設計・最適化では回折光波を正確に数値計算することが必要である．
また位相シフト干渉露光法を用いたデジタルホログ
ラフィの技術により，単に結像した被写体の光強度分布のみならず，その被写体からの位相を含んだ光波複素振幅そのものを電子的に記録することが可能となって来た．この場合，記録された光波複素振幅を数値的に逆回折伝播することにより，被写体の像が得られる．
また逆に，架空の物体からの光波を数値合成し，ホログラム干渉縞をデジタル的に描画することによって作 製 さ れ る デ ジ タ ル 合 成 ホ ロ グ ラ ム(Computer-Generated Hologram, CGH とも呼ばれる)では，点光源を仮想物体の表面に配置する従来の手法から脱却し，小さなパッチで構成された表面モデル 3

は，フレネル・キルヒホッフ積分やフレネル近似，フランフォーファ近似，フレネル・ゾンマーフェルド法，波面展開法などで数値計算することが可能である．特に波面展開法は，フレネル近似やフランフォーファ近似における伝播距離の制約が無く，2 回のフーリエ変換と若干の付加演算で実現できるため，数値計算向きの方法と言える．しかしながら波面展開法も含めて，従来から知られた回折光波計算法に共通の制約事項は，開口面と回折光観測面が平行でなければならないことである．傾いた開口面からの回折計算は，フレネル領域[1]やフランフォーファ領域[2]ではすでに報告されている．しかしながら，通常のフレネル回折計算・フラウンフォーファー回折計算と同様に，これらは伝搬距離に関する制約を抱えている．
そこで，本件研究では伝播距離の制約が無い波面展開法において，開口面あるいは観測面を任意に傾けて回折光波を高速に計算する方法を考案した[3-6]．

次元物体の表面そのものからの光波を数値合成する方法が，主に計算時間の問題から求められている．自2. 配置と座標系

由空間の光波伝搬・回折の理論は古くから多数知ら

Fig. 1 に示したように，平面 S 上に複素開口関数れているが，この様な応用に対しては，これらの従来g ( x, y ) を定義し，S と平行でない観測平面 R で回折の計算方法・理論は無力である．

光波の複素振幅 f ( xˆ , yˆ ) を求める．図中(a)は観測面 R

等方性で均一な媒質中における光の回折波の伝播が傾いた場合，また(b)は開口面 S が傾いた場合に相第 7 回関西大学先端科学技術シンポジウム講演集, 47-51(2003.1)

(a)

(b)
yy SŷSxŷxx̂

ϕzx̂

z=dẑ

g ( x, y )g ( x, y )RRf ( xˆ, yˆ )f ( xˆ, yˆ )Fig. 1 開口面と観測面の座標系の定義
当するが，両者の違いは平面の回転変換と並進伝播

従って， ( xˆ , yˆ , zˆ ) 座標系から ( x, y, z ) 座標系への回転計算の順序にすぎない．従ってこの様な回折シミュレ変換行列をーションでは光波複素振幅の回転変換が重要となる．

 a1
T = a4
a7

3. フーリエスペクトル領域における座標系の回転
今，Fig.1 において d = 0 とし，回転変換のみを考える． g ( x, y ) のフーリエスペクトルを
G (u , v) = F {g ( x, y )}= ∫∫ g ( x, y ) exp[ −i 2π (ux + vy )]dxdy(1)g ( x, y ) = F

a5a8a3 a6 a9 (4)で与えるとき，平面 S における空間周波数 (u, v) は，平面 R での周波数 (uˆ , vˆ) と，
u = a1uˆ + a2 vˆ + a3 wˆ (uˆ , vˆ)v = a4uˆ + a5 vˆ + a6 wˆ (uˆ , vˆ)(5)の関係にある．そのため，平面 R でのスペクトルはとし，逆変換を−1a2

F (uˆ , vˆ) = G (a1uˆ + a2 vˆ + a3 wˆ (uˆ , vˆ), a4uˆ + a5 vˆ + a6 wˆ (uˆ , vˆ)) (6){G (u, v)}= ∫∫ G (u , v) exp[ −i 2π (ux + vy )]dudv(2)とする時，(2)式は平面波 G (u, v) exp[−i 2π (ux + vy )] を平面 S 上で積分していると解釈できるため，平面波の波動ベクトル k =2π[k]で与えられる．このスペクトルを再積分することにより，
平面 R 上における複素振幅が次のように求まる．
f ( xˆ , yˆ ) = ∫∫ F (uˆ , vˆ) J (uˆ , vˆ) exp[i 2π (uˆxˆ + vˆyˆ )]duˆdvˆ
∞{= F −1 F (uˆ, vˆ) J (uˆ , vˆ)}(7)

k z と空間周波数の間

ここで，再積分時の積分要素の変換が非線形であるには， k = 2π [u v w(u , v)] の関係があることがわかる．

ため，回転変換の前後で光波のエネルギーを保存すλxky

こ こ で ， k = 2π / λ ， k x2 + k y2 + k z2 = 1 で あ り ， ま たw(u , v) ≡ λ−2 − u 2 − v 2 である．同様に f ( xˆ , yˆ ) のスペるためにヤコビアン
J (uˆ , vˆ) = a1a5 − a2 a4 + (a2 a6 − a3 a5 ) uˆ wˆ (uˆ , vˆ)+ (a3 a4 − a1a6 ) vˆ wˆ (uˆ , vˆ)(8)が必要である．

クトルは，F (uˆ , vˆ) = F { f ( xˆ , yˆ )}= ∫∫ g ( xˆ , yˆ ) exp[ −i 2π (uˆxˆ + vˆyˆ )]dxˆdyˆ以上の表式(6)〜(8)の導出では，全く近似を用いて(3)で あ り ， ( xˆ , yˆ , zˆ ) 座 標 系 で の 波 動 ベ ク ト ル はkˆ = 2π [uˆ vˆ wˆ (uˆ , vˆ)] ， wˆ (uˆ , vˆ) ≡ λ−2 − uˆ 2 − vˆ 2 である．

いない．しかしながら，数値計算で離散フーリエ変換を用いる際には，(6)式のスペクトルの回転変換において，サンプリング点の分布が不均一になるため，スペクトルの補間が必要となり，このスペクトル補間の方法が計第 7 回関西大学先端科学技術シンポジウム講演集, 47-51(2003.1)
ング間隔は A と同じである．Sample C は A と同様のガ算精度決定することになる．

ウスビームであるが，焦点距離 15[mm]のレンズを通過した場合と同様の位相を与えており，レンズから距離4. 一軸回転の精度と計算時間
傾いた観測面での回折光波の一例として，Fig.1(a)の配置で， y 軸を中心に回転した観測面での回折光波を求めた．この場合，回転変換行列は，
cos ϕT =  0 sin ϕ0 − sin ϕ 10 0 cosϕ 

10[mm]伝搬した位置での複素振幅である．これらの計算では，スペクトルの補間には Cubic 補間法を用いた．
この結果の精度を評価するため，Fig.3 に示した様(9)であるので，平面 R でのスペクトルとヤコビアンは，に，通常の波面展開法で多数の平面上で複素振幅を計算することにより，厳密な 1 次元分布 f ( xˆ ,0) を計算し，本手法による計算結果との SN 比，z = d でのスペクトル Gd (u, v) に対して

(10)

SourcePlaneとなる．これを(7)式に代入して逆フーリエ変換すること

g ( x, y )F (uˆ , vˆ) = Gd (uˆ cos ϕ − wˆ (uˆ, vˆ) sin ϕ , vˆ)J (uˆ , vˆ) = cos ϕ + uˆ sin ϕ wˆ (uˆ , vˆ)xx̂

ϕにより，傾いた観測面上での複素振幅を求めることがz

できる．
z=dFig.2 は，この方法で計算した，傾いた観測面上のf rig ( xˆ, yˆ )光波複素振幅である．Sample A は波長 633[nm]のガウスビームの断面を傾いた観測面上で見たものであり，Reference
Plane

左列が振幅，右列がグレイスケールで表わした位相である．サンプリング数 256×256，サンプリング間隔 20µm

Fig. 3 精度評価のための厳密解の計算である．一方 Sample B は矩形開口から 30[mm]の位置での回折光波であり，波長，サンプリング数，サンプリ

Sample A

Sample B

Sample C

ϕ = 0°

ϕ = 30°ϕ = 60°ŷx̂Fig. 2傾いた観測面上で求められた複素振幅．各例について左列は振幅，右列は位相をグレイスケールで描いている

第 7 回関西大学先端科学技術シンポジウム講演集, 47-51(2003.1)

36.9

Inverse FFT

40

Sample A

Propagation
Rotation(Interpolation)

Sample C

Time [sec]Sample B

3.1
202040
SNR [dB]

60

9.6

8.5
10012.1

FFT30

1st2nd3rd
Cubic 8

0.1

0.5

1282

2562

80

2.0

12.1

0Fig. 4 S/N 比で表わした計算精度

5122

10242

20482

Fig. 6 Pentium III 1GHz による計算時間．

∫ f ( xˆ,0) dxˆ
2

SNR =(11)∫ f ( xˆ,0) − f rig ( xˆ,0) dxˆ
2分の 2 を FFT および逆 FFT が占めており，フーリエ領域における回転変換(スペクトル補間)が占める割合はを様々な補間法で求めた結果を Fig. 4 に示す．この図26%に過ぎない．回転変換以外の時間は，通常の並で“n-th”は n 次の反復線形補間法，“Cubic8”は 8 点の進的な光波伝搬計算でも必要な時間であり，回転変Cubic 補間法を表わす．この結果からわかるとおり，補換にともなう時間の増加はそれほど大きなものでない

間法によって多少の差異はあるが，概ね 50dB 以上のことがわかる．

SN 比が得られている．
Cubic 補間法を用い，回転角を変えて測定した SN

5. デジタル合成ホログラムの物体光波合成への応用

比を Fig. 5 に示す．本手法で計算した複素振幅分布

前節では主に，観測面の回転を扱ったが，複素開の SN 比はいずれの角度でも 50dB 以上であることがわ

口関数を定義している平面 S の回転は，デジタル合成かった．

ホログラムにおける表面モデル 3 次元物体からの物体2 次元の場合について，これらの一軸回転した観測光波の合成に有用である．Fig. 7 に物体光波合成へ面上への光波計算の計算速度を測定した結果を Fig.の応用する場合の配置と座標系を示す．ホログラムは6 に示す．これは Pentium III 1GHz による計算時間で( X , Y ) 平面上に位置するとし，ホログラムの前方向にある．この結果からわかるとおり，計算時間のちょうど 3Z 軸を取る．一方，小平面により構成されている物体の80

S q ( xq , y q )Sample A

Y60

SNR [dB]

Sample C

Sample B

yq

4020zqxq
XZ00204060
Rotation Angle [deg]Fig. 5 計算精度の回転角度依存性

80Hologram
Fig. 7 複素開口関数の回転により物体光波を数値合成するための原理と座標系

第 7 回関西大学先端科学技術シンポジウム講演集, 47-51(2003.1)(-1,-1,-1)(a)(+1,-1,-1)(b)(+2,-4,-1)(c)

Fig 8 表面特性関数(複素開口関数)を回転する方法により数値合成したホログラムの光学再生像．
物体はフラットシェーディングされており，図中の矢印と数字は照明光ベクトルを表わしている．
q 番目の平面を ( x q , y q ) 平面内に含む傾いた座標系

それほど大きなものではない．

( x q , y q , z q ) を定義する．さらに，この ( x q , y q ) 平面内

観測面の回転についてはシミュレーションによりそで表面特性関数 s q ( x q , y q ) を定義する．この関数は，の精度を，また開口面の回転については実際にホログ物体を構成する小平面の形状，散乱性，質感，テクスラムを作製し，その効果を確認した．

チャなどの情報を含む複素関数であり，複素開口関数
g ( x, y ) に相当している．

この表面特性関数からの回折光波 hq ( X , Y ) をホログラム上で求め，物体を構成するすべての平面について総和，H ( X , Y ) = ∑ hq ( X , Y )qを求めることにより，物体光波を合成する．
このような考え方に基づいて数値合成した物体光波に対して作製したホログラムの光学再生像の例を Fig.
8 に示す．ホログラムはいずれも 8192×4096 ピクセルで，ピクセル間隔 2µm×4µm で作製したものである．従来から用いられている点光源モデルで同様の物体光波を合成した場合数 10 時間を要するのに対して，本手法を用いるといずれのホログラムも 30 分以下で数値合成が可能であった[7]．

6. まとめ従来困難であった傾いた平面への光波回折や，傾いた開口からの光波回折を計算する手法を提案した．この手法は波面展開法に基づいてフーリエ空間で座標回転する方法であり，伝搬距離に関する制約が全くない．また，回転変換に伴う計算時間の増加は，通常の平行な平面間の回折・伝搬計算に対して約 26%であり，参考文献[1] D. Leseberg and C. Frére: Appl. Opt. 27,(1988)3020.
[2] H. J. Rabel, N. Bolognini, and E. E. Sicre: Opt.
Acta 32, (1985)1309.
[3] 松島恭治，H. Schimmel，F. Wyrowski: 波面展開法を用いた表面モデル CGH の高速計算法，3 次元画像カンファレンス 2001，13-16(2001).
[4] 松島恭治, H. Schimmel, F. Wyrowski : 波面展開法による傾いた平面への光波回折シミュレーシ
ョン, Optics Japan 2001, pp.229-230 (2001).
[5] K. Matsushima, H. Schimmel, F. Wyrowski: New
Creation Algorithm for Digitally Synthesized
Holograms in Surface Model by Diffraction fromTilted Planes, Practical Holography XVI andHolographic Materials VIII, SPIE Proc. #4659, SanJose, 53-60(2002).
[6] K. Matsushima and H. Schmmel and F. Wyrowski:
Fast Calculation Method for Optical Diffraction on
Tilted Planes Using the Angular Spectrum of Plane
Waves, (to be published).
[7] K. Matsushima, A. Kondoh: Wave Optical
Algorithm for Creating Digitally Synthetic
Holograms of Three-Dimensional Surface Objects,
Practical Holography XVII and Holographic
Materials IX, SPIE Proc., Santa Clara, (2003) (inpress).