第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)

自由空間における光波伝搬計算手法とその応用
Numerical techniques for free-space propagation ofwave-fields and their application
松島恭治
Kyoji MATSUSHIMA
関西大学システム理工学部電気電子情報工学科，635-0831 吹田市山手町 3-3-35
Dept. of Electrical and Electronic Engineering, Kansai University
E-mail: matsu@kansai-u.ac.jp

要旨
自由空間における光の回折伝搬計算は古くて新しい問題である．光が開口面を通過した結果生じる回折像は古典的な問題であり，フレネルやキルヒホッフといった 18～19 世紀の歴史的科学者にまで遡る問題である．一方，近年になり，SLM をはじめ，微細加工技術，イメージセンサー等の発達により，この古い問題は新しい様相を示している．これらの新しい技術はいずれも，波長に近いレベルで光波面を変調あるいは観測することを可能とした．その結果，任意の光波をデジタル的に発生したり，位相を含めて観測したり，あるいは音響や画像のように光波をデジタル信号処理することが可能な時代になりつつある．
このような光波信号処理においてはコンピュータを用いた光の伝搬回折計算が不可欠になるが，理論的な解析式をコンピュータプログラムに置き換えて実行することは必ずしも簡単ではない．そこには必ず，サンプリングに関する諸問題が生じる．また，波長に近いサンプリング間隔である程度のサイズの光波を処理するための計算量は本質的に莫大なものであるため，いかに高速かつ正確に計算するかが重要となる．
本稿ではまず，従来の教科書ではあまり扱われていない光波の回折伝搬計算の数値処理について最近の技術を述べる．特に近年は，計算理論の発達により，光学の教科書で中心的に取り上げられるフレネル近似による回折伝搬計算は数値計算手法としての意義を失いつつある．それに代わって厳密解を与える Rayleigh-Sommerfeld 公式に端を発する手法が有力になっている．本稿ではまず畳み込み形式の Rayleigh-Sommerfeld 積分を高速に計算する角スペクトル法の基礎理論を述べている．また，角スペクトル法に基づき，光軸から外れた位置での光波を計算するオフアクシス型の計算法や開口面とは非平行な観測面で光波回折伝搬を計算する手法等について述べる．
また，これらの計算手法の応用として，SLM を用いた光スイッチング素子のシミュレーション，回折光学素子の一種であるレリーフ型ビーム整形素子の設計とそれによる複雑なビーム形状の整形実験，深く傾斜した表面のデジタ
ルホログラフィによる撮像と再生，コンピュータホログラムによる 3 次元画像の生成等について，著者の研究を紹介する．

第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)

y

伝搬先

伝搬元

g(x,y,0)(x0, y0)yg(x,y,z)

xxxz
伝搬元

q

y

伝搬元

伝搬先

伝搬先

(a) 平行平面間の伝搬計算(ｂ) オフアクシス型計算(c) 回転変換

Fig. 1 伝搬計算の種類と座標系

1. はじめに
光波回折の解析手法は十分に確立され，著者らも種々の名著[1-3]を通してその理論的な手法を学んだ．しかし理論的な解析式をコンピュータプログラムに置き換えて実行することは必ずしも簡単ではない．実際に光波を信号処理するためには，高速フーリエ変換(FFT)をはじめとする数値計算手法を駆使し，理論式に基づき，短時間に正確に数値計算できなければならない．
コンピュータそのものの発達や本稿で述べるような計算技術に発達により，理論的な枠組みの重心も変わりつつある．例えば，フレネル近似による回折理論は光の性質を理解するためには不可欠であるが，数値計算手法としてのフレネル回折計算法は近接場を取り扱える角スペクトル法より特に優れたところはなく，意味を失いつつある．
本稿では基礎理論だけではなく，その応用として，ビーム整形や光スイッチング，デジタルホログラフィによる撮像，コンピュータホログラムによる 3 次元画像の生成等について，著者の研究を中心に概説する．

2. 基礎理論
良く知られているように，無損失で等方的な一様媒質ではマクスウェル方程式は波動方程式に帰着する．さらに，時間発展を考えない単色光では，電界と磁界の各成分はヘルムホルツ方程式 (  k 2 ) g ( x, y, z)  0 を満たす．ここで g(x,y,z)は電磁界ベクトル成分の複素振幅分布で，k = 2/は波数である．自由空間での伝搬計算問題とは，一般に，Fig.1(a)に示すように，特定の平面(伝搬元平面)で与えられた境界条件 g(x,y,0)に対して，別の平面でヘルムホルツ方程式の解を見つけることに他ならない．この別の平面を観測面あるいは伝搬先平面と呼ぶことにする．
とりあえず，伝搬元と平行な平面を観測面として解を g(x,y,z)とすると，Rayleigh-Sommerfeld 公式による厳密解は
g ( x, y, z )   g ( x ', y ', 0)

exp i 2 r '  1 z  11   dx ' dy 'r'r '  2 r ' i (1)

1/ 222で与えられる．ここで r '   x  x '   y  y '  z 2  である．この積分公式は畳み込み積分として，
g ( x, y, z)  g ( x, y,0)  h( x, y, z)(2)と書き直すことができる．ここで，記号*は畳み込み積分を表す．また，点応答関数はh ( x, y , z ) であり， r  x 2  y 2  z 2



exp i 2 r  1 z  11 rr  2 r i (3) である．(2)式の両辺をフーリエ変換し畳み込みの定理を用いることにより，1/ 2

G(u, v; z)  G(u, v;0) H (u, v; z)(4)となる．ここで，u と v はそれぞれ x 方向と y 方向の空間周波数であり，波動ベクトル k  kx x  k y y  kz z の対応する成分と u  kx 2 ， v  k y 2 と関係付けられる(x, y, z は単位ベクトル)．また k  k であるので z 方向の空間周波数は独立ではなく，

第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)
Table 1 光波伝搬の主な数値計算法とその特徴
理論式

数値計算手法

計算量(入力:n2)Rayleigh-Sommerfeld
公式(厳密)

高速畳み込みFFT(4n2)×2

フレネル近似フラウンフォーファ近似

標本間隔

近距離

遠距離

オフア
クシス

非平行

角スペクトル法(平面波展開法) [3]

一定

○

△××

帯域制限角スペクトル法[6]

一定

○

○××

シフテッド角スペクトル法[7]

一定

○

○

○

×

回転変換[8,9]

可変―――

○

高速畳み込みFFT(4n2)×2

一定△△××

FFT

FFT(n2)×1

変化

△○××

Scaled FFT (シフテッドフレネル法[5])FFT(4n2)×3

可変×○○×FFTFFT(n2)×1

変化

×○××

  2  u 2  v 2 1/ 2w  w(u, v)  
0

u 2  v 2   2(5)otherwiseとなる．これを用いると，周波数応答関数はH (u, v; z)  F h( x, y, z)  exp i 2 w(u, v) z (6)と 導 か れ る ． こ こ で ， 記 号 F は フ ー リ エ 変 換 を 表 す ． 従 っ て (x,y,0) 平 面 で 与 え た 複 素 振 幅 の ス ペ ク ト ルG(u, v;0)  F g ( x, y, 0) から，(4)-(6)式により任意の z 位置でg ( x, y, z)  F 1 G(u, v;0) exp i2 w(u, v) z (7)が求まる．これが角スペクトル法である ．さらに，光が近軸領域にあるとすると，近似 w(u, v)   1  (u 2  v 2 ) / 2 が1成り立つことから，フレネル回折の公式が求まる．この場合，点応答関数と周波数応答関数はそれぞれ， 
1x2  y 2 exp ik  z i z2 z   (8)

2
H FR (u, v; z )  exp i z    u 2  v 2  (9)hFR ( x, y; z ) 

となる．

3. 数値計算法
光波伝搬を数値計算するための主な手法とその特徴を Table 1 に示す．ここで，計算量は伝搬元の複素振幅のサンプリング数を n2 とし，それに対してどのようなサンプリング点数の FFT を何回行う必要があるかで表示している．
例えば，フラウンフォーファ回折は n2 の FFT を一回実行することで計算ができるので FFT(n2)1 と表記している．また，Fig.1 に示したように，ここでは伝搬の種類を三つに分類している．まず(a)はもっとも一般的な平行平面間の伝搬計算であり，伝搬元と伝搬先のサンプリング窓は同一光軸上にある．(b)は伝搬先サンプリング窓が軸と直交する方向にシフトできる性質であり，オフアクシス型回折計算と呼んでいる．(c)は伝搬元と伝搬先のサンプリング窓が平行でない計算である．
1「角スペクトル法」はまだ定着した呼び名ではない．平面波展開と呼ばれることもある．Goodman のテキスト[3]では“Angular spectrum of plane wave”と呼んでいる．また，この公式は Rayleigh-Sommerfeld 公式からだけではなく，ヘルムホルツ方程式から直接導くこともできる[4]．

第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)
3.1 フレネル近似計算法
フレネル近似を用いた場合の計算法は，大きく二つに分かれる．一つ目は，(8)式を用いて積分をg ( x, y, z ) 

 



1
exp ik  z  ( x 2  y 2 ) 2 z  F g ( x, y ,0) exp ik ( x2  y 2 ) 2 z i zu  x /  z ,v  y  /  z(10)と書き直して，一度の FFT で実行する手法である．この手法は計算量的にはもっとも高速であるが，伝搬先のサンプリング間隔が伝搬距離に比例して拡大するという重大な問題がある．Scaled FFT と呼ばれる特殊なフーリエ変換でこれを克服し，さらにオフアクシス型計算や伝搬先のサンプリング間隔変更を可能としたのがシフテッドフレネル法である[5]．これは非常に便利な計算法である．
フレネル近似計算では，周波数応答関数として(9)式を用い，(4)式の畳み込みによる計算もよく用いられる．しかし，この手法で必要な計算量は角スペクトル法となんら変わらないため，フレネル近似を用いる理由は特になく，同じ畳み込みで行うのなら角スペクトル法を用いる方が良い2．
3.2 高速畳み込み演算と角スペクトル法の数値計算
FFT を用いた高速畳み込みには注意点がある．それは，そのままでは円状畳み込みとなることである．これは離散フーリエ変換による繰り返し構造が現れた畳み込みである．この場合，伝搬先のサンプリング窓の周囲にはその繰り返し構造が無限に存在するため，伝搬先でサンプリング窓を越えて光が広がると，隣り合ったセグメントからの光の侵入が起き計算結果に重大な誤差が生じる．これを防ぐため，サンプリング窓を 4 倍拡張してゼロ詰めし，伝搬計算後に 4 分の 1 に縮小する処理を行う[6]．これにより，円状畳み込みを線形畳み込みに変換し，伝搬後の光の広がりに関係なく，サンプリング窓の範囲内では正しい計算結果が得られる．これはフレネル近似であろうとも角スペクトル法であろうとも同じであり，そのため，畳み込みを用いた場合，一般に 4n2 点の FFT が 2 回必要となる3．
(6)式をそのまま用いて角スペクトル法を実行すると，特に伝搬距離が長い場合に，誤差を生じる．これは，周波数の高くなるにつれ繰り返し周期が短くなるチャープ的性質が(6)式の周波数応答関数にあるため，離散化した際にエイリアシング誤差が生じることによる．これを防ぐためには，(6)式を計算する範囲，すなわち帯域を制限する必要がある．これを行うのが帯域制限角スペクトル法である[6]．この手法では，物理的な光波伝搬に必要十分な範囲に帯域制限を行う．それ以外は角スペクトル法と同じであり，これにより伝搬距離によらず高精度の計算が行える．
さらに，伝搬先の座標系が軸方向と垂直に(x0, y0)シフトしているとすると周波数応答関数

Hˆ (u, v; x0 , y0 , z)  H (u, v; z)exp[i 2 ( x0u  y0v)](11)が得られる．この関数の離散化でもやはりエイリアシング誤差が発生するので，特定の矩形周波数領域で帯域制限することによりオフアクシス型の角スペクトル法が可能になる．これをシフテッド角スペクトル法と呼んでいる[7]．この手法は帯域制限角スペクトル法を包含しており，x0 = y0 = 0 とすることで帯域制限角スペクトル法に帰着する．同じオフアクシス型でも，上述のシフテッドフレネル法は伝搬距離が短い場合に重度のエイリアシング誤差を生じるが，シフテッド角スペクトル法ではそのような問題は生じない．
3.3 回転変換
伝搬元と伝搬先のサンプリング窓が非平行な場合の伝搬計算も角スペクトル法を用いて行える．これが Fig.1(c)に示す光波の回転変換である．上述のとおり空間周波数と波動ベクトルには k  2 (ux  vy  wz) という関係があり，ベクトルは簡単に座標変換できることから，伝搬元光波のスペクトルをそれとは非平行な座標系におけるスペクトルにマッピングできる．これを行うのが光波の回転変換である[8,9]．ただし，一般に伝搬元も伝搬先も格子状のサン2フレネル近似による(9)式では平方根計算が不要であるため，その分計算量が少ないように見えるが，現代のCPU では開平演算は除算と同じ時間で実行できるため，これは利点にはならない．
3
伝搬先でほとんど光波が広がらない場合には 4 倍拡張が不要なこともある．

第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)
プリング点を持つが，スペクトルの回転変換ではこの格子に歪みがLCOS
位相型SLM

入力ポート生じるので補完が必要になる．しかし，補完方法を上手に選べばかなり高い精度が得られる．

光ファイバ

4. 光波伝搬計算の応用

LCOS
位相型SLM

出力ポート構造4.1 SLM を用いた光スイッチング素子のシミュレーションFig.2 は，santec 社が中心となって開発している光スイッチング素光ファイバ非偏向時子のシミュレーションの一部である．この素子は，2 枚の位相型 SLM

振幅に表示した位相格子を切り替え，入力ポート側光ファイバから出射されたビームを任意の出力ポートに接続するものである．光軸に対偏向時して SLM が傾いているため，回転変換を用いて SLM 面上での光振幅波を求めて反射・回折特性をシミュレートし，また軸から離れた出力ポートへ偏向する際にはシフテッド角スペクトル法を用いて，ビーム位相パターンを計算している．
4.2 回折型ビーム整形素子[10]Fig. 2

SLM を用いた光スイッチング

微細加工技術を用いて製作されるレリーフ型回折光学素子を用いて TEM00 ガウスビームを任意の形状のビームパターンに整形する

t ( x, y)tout ( x, y)

U ( x, y)

Output datauin ( x, y)

設計と実験結果を Fig.3 に示している．位相分布の設計においては，

Forwardpropagation

Constraint in DOE

Constraint in output field

目標となる整形パターン(b)を用いて，(a)に示すように素子面と像面との間で反復的に伝搬・逆伝搬計算(a) 反復回折アルゴリズムt ( x, y)

Backward
propagation

Input dataU ( x, y)usig ( x, y)を行い，目標パターンが得られるような空間位相分布に収束させる．
微細加工技術を用いてレリーフ型素子を作製し，設計した位相分布でガウスビームの空間位相変調を行った時の実験結果が(c)であり，5mm(b) 目標パターン発生したビームパターンの拡大写

Fig. 3(c) 実験結果(ｄ) ビームの拡大写真

回折型ビーム整形素子の設計と実験結果[10]真が(d)である．このように，複雑な形状のビーム整形が確認できる．なお，この実験の時点では微細加工技術を用いた静的な素子として製作したが，最近の高解像度 SLM であれば，電気的な変調による動的なパターン形成も可能と思われる．
4.3 デジタルホログラフィによる斜面像の再生
小さな板状の物体が撮像素子に対して大きく傾いて設置された場合，その表面パターンを撮像することは容易ではない．実際，Fig.4(b)に示すように，カメラによるマクロ撮影でも，あるいは(c)のようにデジタルホログラフィによる通常の再生方法でも被写界深度が足りず，中央部以外はピンボケした像になる．しかし，デジタルホログラフィで取得した光波を回転変換し，伝搬先平面を物体表面に一致させて振幅像を観察すると全体がクリアに再生される[9]．

第 12 回情報フォトニクス研究グループ研究会(秋合宿)講演予稿集, pp.13-18(2011.9.15)
4.4 ポリゴン法コンピュータ(a) 被写体(ｂ)
カメラによるマクロ撮影

ホログラム
ポリゴン法の理論モデル

1 cm

Planar object

では，表面関数と呼ぶ複素x̂関数の振幅分布がポリゴン形状の開口，また位相部分ŷq 70

Image sensorが散乱体の役割を果たす(ｃ)回転変換前y[11]．この表面関数を一種の光波複素振幅分布とみなし，回転変換することによりホロzx

グラムと平行な平面での光(d) 回転変換後波を合成している(Fig.5)．また，メモリに収容できない大Fig. 4

デジタルホログラフィによる斜面の撮像と再生[9]規模なホログラム計算では，オフアクシス型伝搬計算法で分割計算を行っている[12]．

5. 光波伝搬シミュレーションツールとまとめ
本稿で述べたような伝搬計算手法やそれ以外の多数の関数をライブラリ化し，WeaveField Tools として WEB サイトで公開している．これらの関数は十分にテストされチューニングされているため，このライブラリを使用すると比較的簡単に応用例のようなプログラム作成が行える[13]．
本研究の一部は日本学術支援振興会の科研費(215001114)の助成を受けたものである．

参考文献
飯塚啓吾: 光工学(新訂増補版), (共立出版, 東京, 1983).
谷田貝豊彦: 光とフーリエ変換, (朝倉書店, 東京, 1992).
J. W. Goodman: Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996).
辻内監修, 松島他著: ホログラフィー材料・応用便覧, 第 2 部 6.1 節, (エヌ・ティー・エス, 東京, 2007).
R. P. Muffoletto, et al.: Shifted Fresnel diffraction for computational holography, Opt. Express 15, 5631–5640 (2007).
K. Matsushima et al.: Band-limited angular spectrum method for numerical simulation of free-space propagation in far andnear fields, Opt. Express 17, 19662 (2009).
[7] K. Matsushima: Shifted angular spectrum method for off-axis numerical propagation, Opt. Express 18, 18453 (2010).
[8] K. Matsushima, et al.: Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of planewaves, J. Opt. Soc. Am. A20, 1755 (2003).
^y1 y1[9] K. Matsushima: Formulation of the rotational^
表面関数1x1transformation of wave fields and their^application to digital holography, Appl. Opt. 47,y3 y3x1 y2
D110 (2008).
^y2[10] 松島他: 反復最適化アルゴリズムを用いた任意形状ビーム整形用回折光学素子の設計,
P1
レーザー研究 35, 299 (2007).
^x2x3[11] K. Matsushima: Computer-generatedP^P3holograms for three-dimensional surface
2x2x3objects with shade and texture, Appl. Opt. 44,4607 (2005).
[12] K. Matsushima, et al.: High-definition full物体モデル
物体光波(振幅像)
parallax computer-generated hologram createdby the polygon-based method, Appl. Opt. 48,
表面関数3
表面関数2
H54 (2009).
[13] http://www.laser.ee.kansai-u.
Fig. 5 ポリゴン法による仮想物体光波の数値合成[12]
[1]
[2]
[3]
[4]
[5]
[6]

ac.jp/WaveFieldTools