はじめに
前回までに3回にわたってスペックル(speckle)で生じる光の分布の計算方法を詳細に解説しました。
実際にスペックルの影響を考える際にはここまで詳細に扱う必要性に迫られることは少なく、スペックルの統計的な性質を押さえておけばスケーリング則で見積もりを立てられることが多いと思います。
スペックルの統計的性質に絞り、スペックルコントラスト低減させる場合の限界についても議論します。
前回記事
ある一点における強度分布のヒストグラム
以下では完全発達スペックルに絞って解説します。
完全発達スペックルでは強度分布が指数分布になることが知られています。
指数分布は下式で表せます。
$$f(I) = \frac{1}{\mu} \exp\left(-\frac{I}{\mu}\right)$$
計算内容は単純で各点における位相がランダムに付加されるとして、フーリエ変換することで遠方での強度分布(フラウンフォーファー回折)を計算出来ます。
import numpy as np
import matplotlib.pyplot as plt
N = 1024
# ランダム位相
phase = np.random.rand(N, N) * 2 * np.pi
# 複素振幅
field = np.exp(1j * phase)
# フーリエ変換（遠方回折）
speckle=np.abs(np.fft.fftshift(np.fft.fft2(field)))**2
I = speckle.flatten()
I = I / np.mean(I)
plt.hist(I, bins=100, density=True, alpha=0.6, label="simulation")
x = np.linspace(0,7,200)
plt.plot(x, np.exp(-x), label="exponential distribution", linewidth=2)
plt.xlim(-0.01,7.1)
plt.title("Single speckle intensity distribution",fontsize=16)
plt.xlabel("Normalized intensity",fontsize=16)
plt.ylabel("Probability density",fontsize=16)
plt.legend()
plt.show()
オレンジ色の指数分布フィットがよく合っていることが分かります。
指数分布は、”0″が一番確率密度が高くなる分布です。これは、観測時には何も検出できない可能性が一番高いという意味で、少し変な感じがします。
上の状況設定は、ある1点に飛んでくる光を議論しているからです。実際の光学系ではある程度の広さをレンズなどを使って集光することが行われます。この場合は強度分布は指数分布から外れていきます。次節で計算してみましょう。
開口上で積分した時のスペックル統計
スペックルを1枚生成した後、空間的に近接した領域（3×3画素）を切り出して平均することで、開口による空間積算を模擬しました。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
N = 2048
def generate_speckle(N):
phase = np.random.rand(N, N) * 2*np.pi
field = np.exp(1j * phase)
speckle = np.abs(np.fft.fftshift(np.fft.fft2(field)))**2
return speckle
speckle = generate_speckle(N)
M = 3**2
block_size = int(np.sqrt(M))
samples = []
for _ in range(5000):
x = np.random.randint(0, N - block_size)
y = np.random.randint(0, N - block_size)
block = speckle[x:x+block_size, y:y+block_size]
samples.append(np.mean(block))
I_avg = np.array(samples) /np.mean(samples)
plt.hist(I_avg, bins=100, density=True, alpha=0.6, label="simulation")
x = np.linspace(0,3,200)
plt.plot(x, gamma.pdf(x, M, scale=1/M), linewidth=2, label="Gamma")
plt.title(f"Averaged speckle (M={M})")
plt.xlabel("Normalized intensity")
plt.ylabel("Probability density")
plt.legend()
plt.show()
ヒストグラムは右側に尾を引いていることが分かります。オレンジ色で示したガンマ分布でよくフィットできています。
指数分布はガンマ分布（shape=1）の特別な場合であり、独立な指数分布の和はガンマ分布（shape=M）になります。
積算する開口を広げていったらどうなるでしょうか？M=25の場合を計算してみます。
左右対称な分布となっています。
またばらつきもM=9に比べて小さくなっています。なぜでしょうか？
積算領域を広げると平均化されるサンプル数が増えるため、分布は対称に近づきます。
各サンプルが独立であると仮定すると中心極限定理により正規分布に収束するのです。
ここまでの結果をまとめると、積算する開口を広げていくと強度分布は指数分布→ガンマ分布→正規分布とクロスオーバーがおこるのです。
実際にスペックルの強度分布が本当に独立かという話は残りますが、空間相関があるときは完全ランダムよりもゆっくりと正規分布に近づいていくことになります。
ガンマ関数について
ガンマ分布の確率密度関数は、一般に以下の式で表されます。
$$f(x; k, \theta) = \frac{x^{k-1} e^{-\frac{x}{\theta}}}{\theta^k \Gamma(k)}$$
ここで、パラメータには以下のような意味と流儀（表現方法）があります。
- $k$（形状母数）： 分布の形を決めるパラメータ。本記事のコードにおけるスペックルモード数 $M$ に対応します。
- $\theta$（尺度母数）： 分布全体のスケールを調整するパラメータ。※分野や教科書によっては、逆数である $\beta = \frac{1}{\theta}$ （逆尺度母数や速度母数）が使われることもあります。
本記事のシミュレーションでは、平均値が $1$ になるよう $\theta = \frac{1}{M}$ （つまり $\theta = \frac{1}{k}$ ）として規格化しています。
したがって、「開口を広げて積算サンプル数 $M$ を増やすこと」は、ガンマ分布の形状母数 $k$ を大きくしていく操作そのものに対応しており、これが指数分布から正規分布へのクロスオーバーを引き起こします。
ガンマ分布は、特殊な場合として指数分布を包含し、形状母数を大きくしていくことで正規分布へと漸近する確率分布です。 個人的には、非対称な分布をフィットしたいときにガンマ分布を使いますが、こうして指数分布から正規分布への「橋渡し」をする分布と捉えると、自然なことだと感じました。
平均と分散の制約
もう少しガンマ関数と正規分布の関係を考察してみましょう。
ガンマ関数の平均と分散は下式で表されます。
平均：$\mu = k\theta$
分散：$\sigma^2 = k\theta^2$
分散を平均の2乗で割ってみると下式が得られます。
$$\frac{\sigma^2}{\mu^2} = \frac{1}{k} \quad (\text{または } \frac{1}{M})$$
つまり、「分散の2乗を平均の2乗で割った値は、常に形状母数の逆数、つまりスペックルモード数Mの逆数に等しい」という強い制約があることが分かります。
そして、この制約はkを無限に大きくした正規分布でも成立します。
スペックルの強さを表す指標として、スペックルコントラストがあるのでした。
$C = \sigma / \mu$
先ほど確認したように、綺麗に $\frac{1}{\sqrt{M}}$ になります。
つまり、「開口を広げて積算サンプル数を $M$ にすると、コントラストは元の $\frac{1}{\sqrt{M}}$ に低下（改善）するが、そこが限界である」といえるのです。
まとめ
本記事では、完全発達スペックルを対象として、その強度の統計的性質について解説してきました。ここまでのポイントは以下の通りです。
- 指数分布への従属：ある1点における完全発達スペックルの強度は指数分布に従い、最も観測されやすい強度が「0」になるという揺らぎの激しい状態を示します。
- ガンマ分布とモード数 $M$：レンズ等の開口を広げて空間的に積算するサンプル数（モード数）を $M$ とすると、強度分布はガンマ分布へと変化します。コード上の $M$ は、そのまま統計学におけるガンマ分布の形状母数 $k$ に対応しています。
- 正規分布へのクロスオーバー：積算サンプル数 $M$ を大きくしていくと、中心極限定理により分布は指数分布からガンマ分布を経て、最終的に左右対称な正規分布へとクロスオーバーします。
最初に「なんとなくの点数」として置いた $M$ というパラメータが、実はスペックルの揺らぎの統計的性質や、光学系におけるコントラストの限界を支配する非常に重要なファクターであることが分かりました。
光学系の設計やノイズの見積もりを行う上でも、この「 $M$ を変えることで分布とコントラストがどう変化するか」というスケーリング則は強力な指針になります。
参考文献
論文の引用元や海外コードの解説で必ず登場するグッドマンの不朽の名著。今回省略した理論の詳細についても解説されています。また、英語の専門用語や理論的なニュアンスをそのまま正確に理解したい、本格志向の技術者におすすめの一冊です。
本記事のように、Pythonを用いて物理現象や統計モデルを実際にコードで検証・シミュレーションしたいエンジニアや研究者に最適な実用書です。
コメント