Optics Photonics Japan 2008 講演予稿集, pp.324–325, つくば(2008.11.05)

高速畳み込みを用いた伝搬回折計算における帯域制限と可変ピッチ角スペクトル法の検討
Discussion on band-limited numerical diffraction by using
fast convolution and variable pitch angular spectrum methodin free-space propagation○松島恭治○Kyoji Matsushima

関西大学 電気電子情報工学科
Department of Electrical and Electronic Engineering, Kansai University
E-mail: matsu@kansai-u.ac.jp
Angular spectrum method, which gives wave fields satisfying Helmholtz equation, is the best method
in calculation of diffraction. However, numerical errors sometimes occur in the method, especially forlong distance propagation. The errors are attributed in aliasing error of the frequency response function.
A technique of band limit for avoiding the aliasing error and method for variable pitch in calculation
are proposed and discussed.
１．はじめに近年のディジタルホログラフィと計算機ハードウェアの進展により，位相を含む光波分布を大規模に取得し，回転変換 1,2)や視点移動 3)等のディジタル信号処理を行うことが可能になりつつある．この様なディジタル信号処理においては，自由空間における伝搬計算を正確かつ高い自由度で行う数値計算技術が重要となる．光の回折伝搬公式としてはフレネルの回折式がよく知られている．しかし，2 回の FFT による高速畳み込みを用いて回折計算をする場合には，計算時間にほとんど差がないため，回折積分式を直接計算する角スペクトル法 4)が有利であり，しばしばこちらが用いられる．一方近年，1 回の FFT で行うフレネ
ル回折計算を Scaled FFT と呼ばれる手法 5)で可変ピッチ化したシフトフレネル法 6)が提案されている．
角スペクトル法はフレネル近似条件に拘束されないため，本来フレネル回折計算より精度が高い．しかしながら，伝搬距離が長い場合には計算結果に高周波ノイズ様の誤差が生じ，しばしばシフトフレネル法よりも計算精度が落ちる．そこで本報告では，帯域制限を行うことにより，角スペクトル法において長い伝搬距離でも高精度に計算する手法，および Scaled FFT を用いて可変ピッチ化する手法について検討する．
２．角スペクトル法と周波数応答関数の帯域制限
いま(x, y, 0)平面上で光波の複素振幅分布 g(x, y; 0)が与えられる時，伝搬距離 z の回折積分は g(x, y; z) =g∗h(x, y; z)で与えられ，畳み込み定理を用いて G(u, v; z) = G(u, v; 0)H(u, v; z)と表される．ここで，G (u , v; 0) = F { g ( x, y; 0)} = ∫∫ g ( x, y;0) exp[−i 2π (ux + vy )]dxdy,

H (u , v; z ) = exp[i 2π wz ],w = (λ −2 − u 2 − v 2 )1/ 2

である．F{}はフーリエ変換を表し，回折波 g (x, y; z)= F−1{G(u, v; z)}は完全な波動方程式の解となっている．
この手法は角スペクトル法と呼ばれ，フレネル回折式は上式の放物面近似式になっている．
数値計算時には離散化して FFT と高速畳み込みを用いるが，この時問題となるのは周波数応答 H(u, v; z)がチャープ関数となっていることである．そのため，特に z が大きい場合には離散化した H(u, v; z)に大きなエイリアシング誤差が発生する．H(u, v; z)の局所的な信号周波数は， ∂ ( wz ) / ∂u と ∂ ( wz ) / ∂v で表される．
そこで，簡単化のため 1 次元で考え，周波数 u のサンプリング間隔を Δu とすると，エイリアシング誤差を避けるためには， 2 ∂ ( wz ) / ∂u ≤ Δu −1 が必要である．このことから，周波数応答関数はu ≤ ulimit = ⎡⎣ (2Δuz ) 2 + 1⎤⎦−1/ 2

λ −1で打ち切られなければならないことがわかる．
この様な帯域制限の効果を調べるため，簡単な 1 次元分布 g(x; 0)を用いて数値実験を行った．この分布は間隔 Δx = 2λ でサンプリングされている．ここでλ は波長である．またそのサンプリング範囲は S x = N x Δxであり，ここでは N x = 1024 とした．なお，FFT を用いてこの様な高速畳み込み演算を行う際には，円状畳み込みを直線状畳み込みに変換するため，g(x; 0)のサンプリング範囲を 2 倍(2 次元では縦横 4 倍)に拡張し，g(x; z)の計算後にそれを 2 分の 1 に切り捨てる処理をしなければならない．そのため， Δu = (2S x ) −1 となる．
伝搬距離 z を S x /10 ≤ z ≤ 100 S x で変化した際のいた．また厳密数値解は回折積分を数値積分することにより得た．この結果から，帯域制限をしない通常の角スペクトル法では z ≥ 7 S x では S/N 比40

SNR [dB]

計算精度を Fig.1 に示す．ここで，計算精度の評
価は厳密数値解に対する複素振幅の S/N 比 1)で行い，g(x; 0)としては幅 Wx = S x / 2 の単純な開口を用

Shifted Fresnel
Ordinary angular spectrum
Band-limited angular spectrum

50302010が下がっているのに対して，帯域制限をした場合には伝搬距離全体に高い S/N 比を保持していることがわかる．

00.1

1

100z [Sx]

３．回折像計算に必要な帯域と可変ピッチの検討
H (u; z ) の 帯 域 を 制 限 す る こ と は す な わ ち

Fig.1 Precision in calculation of diffraction.

G (u;0) の帯域を制限することになるため，それが別の形で悪影響を及ぼさないかが疑問である．この点を調べるため次のようなモデルで必要な帯域を見積もった．すなわち，幅 Wx の開口の下端に点光源があるとして，その点光源から出た光のうち，スクリーン上の計算領域 Sx の上端を掠める平面波を考え，その空間周波数がスクリーン上で回折波を正しく求めるために必要な最大周波数の指標を与えると考えた．これによると必要な帯域 uneed は，−1/ 2

uneed = ⎡{2 z /(Wx + S x )} + 1⎤ λ −1⎣⎦となる．z = 50[Sx]の場合について帯域制限する周波数を変えて S/N 比を測定した結果を Fig.3 に示す．やはり，uneed 以下まで帯域制限すると回折像は劣化するが，それ以上であれば帯域制限して G (u;0) を打ち
210Aperture

ScreenSxWxzFig.2 Model for required
frequency band.
切っても問題ないことがわかる．
Fig.2 と上式からもわかるとおり，z が増加すると uneed の値は減少し，回折像計算に必要なスペクトルの範囲は減少する．Scaled FFT を用いると同じサンプリング数の g ( x;0) からより狭い範囲のスペクトルを計算することができる．すなわち Δu を可変にすることができるため，uneedを超える範囲であればスペクトルの計算範囲を変えることができる．サンプリング数は変わらないため，このような狭小化したスペクトルを逆 FFT することにより求まる回折像のサンプリング間隔は開口側のサンプリング間隔とは異なった値となる．すなわち，Scaled FFT を用いてスペクトルの計算範囲を変えることにより回折像のサンプリング間隔をある範囲内で可変にすることができる．詳細は発表時に述べる．

50u needu limit

40

SNR [dB]

４．まとめ
解析的には波動方程式の完全な解を与える
角スペクトル法が，周波数応答関数のエイリアシングにより数値的には計算誤差を発生することを示し，その対策として帯域制限する方法を示した．また，それを用いて可変ピッチ計算する方法を検討した．

3020
文献1) K. Matsushima et al.: J. Opt. Soc. Am A20,
1755 (2003).
102) K. Matsushima: Appl. Opt. 47, D110 (2008).
00.01
0.02
0.03
0.04
0.05
0.06
3) T. Nakatsuji, K. Matsushima: Appl. Opt. 47,-1limiting frequency, u [λ ]D136 (2008).
4) J. W. Goodman: Introduction to Fourier
Optics, 2nd ed. (McGraw-Hill, 1996),
Fig.3 Dependence of precision on limiting frequency.
chap. 3.10.
5) D. H. Vailey, et al.: SIAM Rev. 33, 389(1991).
6) R. P. Muffoletto, et al.: Opt. Express 15, 5631(2007).