DPS アルゴリズムを理解するための習作 -フーリエ変換による周期の検出-

DPS アルゴリズムでは、逆空間での散乱ベクトルを一次元に射影し、その長さの周期性をフーリエ変換によって検出する。

例えば a* 軸長が L の結晶で散乱ベクトルを a* 方向に射影したとすると、(1, k, l) という反射は全て長さ L に射影されるし、(2, k, l) という反射は 2Lに……という具合になる。したがって射影ベクトルの長さのヒストグラムを作ると、ちょうど L, 2L, 3L, ... の部分にのみピークが表れる。このような理想的な場合には、フーリエ変換をするまでもなく、最初のピークの位置から逆空間での基底の長さ L をすぐに求められる。

現実には、

  • 射影方向は離散的にサンプリングするので厳密ではない
  • 結晶のモザイク性・ビームのスペクトル幅・ビームの発散などの影響により逆格子点が有限の大きさを持つ
  • 振動写真なので逆格子点の中心がゴニオの回転角方向について厳密には求められない

といった問題があり、射影は L, 2L, 3L, ... 付近に幅を持って分布する。また、指数付けには強い反射のみが利用されるため、ヒストグラム上での各ピークの高さは一定ではない。このような状況では、ヒストグラムの最初のピークを選ぶとか、ヒストグラムの最大のピークを選ぶよりも、フーリエ変換を使ったほうが周期性が明確に分かる。

それを確かめるのが以下のコード。

L <- 1 / 40 # recoprocal vector length
N <- 1024 # number of grid points
MAX_IND <- 32 # maximal index

ls <- c()
for (i in 1:MAX_IND) {
#  ls <- c(ls, rep(i * L, 20)) # ideal
#  ls <- c(ls, rnorm(20, i * L, L / 10)) # noise in length
#  ls <- c(ls, rnorm(rpois(1, 20), i * L, L / 10)) # noise in length and number of vectors
  ls <- c(ls, rnorm(rpois(1, floor(40 * (1 - i / MAX_IND))), i * L, L / 10)) # fewer points in higher resolution
}

fft_in <- rep(0, N)
for (x in ls) {
  ind <- ceiling(x * N) # R vector is 1-indexed
  if (ind < 1) next
  fft_in[ind] <- fft_in[ind] + 1
}
par(mfrow=c(2, 1))
plot(seq(0, by=1 / N, length.out=N), fft_in, type='l', main='FFT in', xlab='reciprocal vector length (1 / angstrom)', ylab='count')
fft_out <- fft(fft_in)
plot(Mod(fft_out), type='l', main='FFT out', xlim=c(0, N / 2), xlab='real basis length (angstrom)', ylab='FFT out')
abline(v=1 + (1:10) / L, col=2) # fft_out[1] is DC component

最初に、ノイズがない場合を計算すると以下のとおり。

次に、射影した長さにノイズがあり、高分解能になるほど指数付けに利用できる点の数も減る場合。

ヒストグラムは汚くなっているが、フーリエ変換後にはきれいに格子長のピークが出ているのが分かる。