はじめに
レーザーを用いた測定ではスペックルと呼ばれる干渉縞が観測されることあります。ボツボツとした形状が特徴で測定ではノイズとなることが多いです。
この記事ではスペックルがどのように現れるかをPythonを使って、可視化してみます。フーリエ変換を用いるのですが、曖昧になりがちな縦軸の大きさについても解説します。
スペックル(speckle)とは？
レーザーのような可干渉性の高い光を粗面に照射した際、反射・散乱光が干渉し合って生じるランダムな明暗の斑点模様のことを言います。光学測定のノイズになったり画像の劣化につながったりします。どんな風に見えるの？って人はWikipediaに画像が載っているので参照ください。
https://en.wikipedia.org/wiki/Speckle_(interference)#Speckle_pattern
スペックルコントラスト
スペックルのボツボツの強さを表す量としてスペックルコントラストCがあります。
定義は強度の標準偏差と平均の比率で表されます。
後述のように位相が完全にランダムな時は、「完全発達スペックル」と呼ばれC=1になります。
強度平均と標準偏差が同じということですから、影響の大きさが分かります。
シミュレーション方法
計算条件
・波長:532nm
・計算領域:0.1 mm x 0.1 mmの矩形領域に照明光を照射
・照明光強度均一
計算コード
計算領域の表面には微細な凹凸があるため、各点における光の電場の位相はランダムと仮定します。
試料面の光電場をフーリエ変換すれば、遠方での光の電場を求めることができます。フラウンフォーファー回折ですね。
import numpy as np
import matplotlib.pyplot as plt
# 計算条件
N = 1024 # 分割数
L = 100e-6 # 計算領域の大きさ (0.1 mm)
wavelength = 532e-9 # 波長 (532 nm)
# ランダム位相（粗面の凹凸を模擬）
phase = np.random.rand(N, N) * 2 * np.pi
# 試料面の複素電場（強度は均一と仮定）
field = np.exp(1j * phase)
# フーリエ変換（遠方回折 / フラウンフォーファー回折）
# FFTのシフトを行い、強度（絶対値の2乗）を計算
fft_field = np.fft.fftshift(np.fft.fft2(field))
speckle = np.abs(fft_field)**2
# スペックルコントラストの計算 (標準偏差 / 平均値)
contrast = np.std(speckle) / np.mean(speckle)
# プロット
plt.figure(figsize=(6, 5))
im = plt.imshow(speckle, cmap='gray')
plt.title(f"Speckle Pattern\nContrast: {contrast:.2f}", fontsize=14)
plt.colorbar(im, label='Intensity')
plt.tight_layout()
plt.show()
スペックルらしいボツボツした干渉縞が得られました。スペックルコントラストが確かに1になっています。黒の単色に見える人はズームしてみてください。
また、光の強度自体は等方的なようです。試料の物性や照明条件の詳細を取り込んでいないため、純粋なスペックルの影響のみが計算されています。このシリーズのどこかでその辺りも扱えたらと思います。
この画像の縦軸の単位は何でしょうか？今はFFTしただけなのでただのインデックスになっています。
実空間上の粗面をFFTしたため、波数に変換されます。波数は光の伝播方向に対応します。
絶対値をそのまま表示すると分かりにくいので、通常は波長532nmの波数で規格化して再度グラフ化します。
k = 1 / wavelength#532nmの波数[cycle/m]
kx = np.fft.fftshift(np.fft.fftfreq(N, d=L/N))
kx /= k
ky = np.copy(kx)
KX,KY = np.meshgrid(kx,ky)
fig, ax = plt.subplots(figsize=(6,6))
c = ax.contourf(KX, KY, speckle, levels=50, cmap="gray")
ax.set_aspect("equal")
ax.set_xlabel("KX")
ax.set_ylabel("KY")
fig.colorbar(c, ax=ax, label="Intensity",shrink=0.7)
plt.tight_layout()
plt.savefig("k-value-speckle.png",dpi=300)
plt.show()
縦軸と横軸を付けることが出来ました。ただし、波数の絶対値が1を超えています。これはどういう意味でしょうか？
波数制限
FFTのようなデジタル信号処理においては波長 $\lambda $ の光が持つ最大の波数は $k_{max} = \frac{1}{\lambda}$ です。これを超える領域は「物体波が伝播しない領域（非伝播波・エバネッセント波）」、つまりレンズの開口数でいう $\text{NA} > 1$ の領域に対応するため、実際の光学系では観測されません。
そのため、実際の計算ではこの最大波数を超える外側の領域をマスクします。
最終的に下図のようになります。以前散乱光計算した際の天球分布と同じ形で求められました。この形で計算できるのであればスペックルを結像したらどう見えるか？も計算できます。像の劣化具合も議論できそうです。
波数に関する注意事項
波数には定数の$2\pi$がつく、つかないの2種類あります。
・物理の教科書（sin/cosの引数など）： $k = \frac{2\pi}{\lambda}$ のように、定数 $2\pi$ がつく（単位: [rad/m]）
・デジタル信号処理（NumPyのFFTなど）： デフォルトでは $2\pi$ がつかないサンプリング周波数ベース（単位: [cycle/m]）
感覚的な説明ですが、sinやcosの引数の場合は $2\pi$がつく、FFTで計算される波数は $2\pi$がつきません。
この点を意識しておかないとスケールが狂ってしまいます。
筆者は物理専攻なので、$2\pi$がつくのが普通だと思っていたのでとても混乱しました。
まとめ
粗面でランダムに位相が付加されるとして、FFTを使ってスペックルを計算しました。
・FFT後の縦軸横軸の単位は波数であり、$k_{max} = \frac{1}{\lambda}$で規格化するとよい。
・波数には定数$2\pi$付くつかないの2種類があり、FFT計算では付かない。
・波数が$k_{max} = \frac{1}{\lambda}$を超えたところは実際には伝播しない領域であり、マスクするとよいことが多い。
疑問
位相は0 – $2\pi$の範囲でランダムと仮定しましたが、これは現実の状況を表わしているのでしょうか？ ランダム位相は試料表面の高さの不均一さに由来するものです。例えば高さのばらつきが波長よりもずっと小さいと位相のばらつきは $2\pi$よりもずっと小さいはずです。また、そのような時スペックルコントラストはどうなるだろう？と色々疑問が湧いてきます。
そこで、次回は実際の物質の高さを使って計算してみます。
次回記事
参考文献
論文の引用元や海外コードの解説で必ず登場するグッドマンの不朽の名著。英語の専門用語や理論的なニュアンスをそのまま正確に理解したい、本格志向の技術者におすすめの一冊です。
回折計算（FFT）をコードに落とし込む際の物理的・数学的な根拠を最も丁寧に解説してくれる不朽の名著です。光学シミュレーションで迷ったときの道標として、手元に置きたい一冊。
コメント