R でフーリエ変換を行い、位相の重要性を確認する

R でフーリエ変換を行うには、stat パッケージ(つまり、デフォルトで読み込まれる環境)に含まれる fft を使う。ついでに、いわゆる Fourier duck の実験をしてみた。

N <- 128 # 実空間の点の数

# トルエンっぽい平面分子。座標は手計算。原子距離は適当。
atoms <- matrix(c(0, 1, -1.73 / 2, 0.5, -1.73 / 2, -0.5, 0, -1, 1.73 / 2, -0.5, 1.73 / 2, 0.5, 1.73, 1), nrow=2)
plot(atoms[1,], atoms[2, ], asp=1) # 確認

real_space <- matrix(0, N, N) # 電子密度格納用 N x N
scale <- N / 5 # [1, N] を [-2, 2] にマップする
for (i in 1:ncol(atoms)) {
  ax <- atoms[1, i]
  ay <- atoms[2, i]
  for (j in 1:N) {
    for (k in 1:N) {
      x <- j / scale - 2
      y <- k / scale - 2
      dist2 <- (x - ax) ** 2 + (y - ay) ** 2
      real_space[j, k] <- real_space[j, k] + 10 * exp(-dist2 / 0.2)
      # Gaussian 型の電子密度。係数は適当
    }
  }
}
image(real_space, asp=1)
reci_space <- fft(real_space) / N / N # この割り算が必要
image(Mod(reci_space), col=Arg(reci_space) %% pi, asp=1)

back_transform <- fft(reci_space, inverse=T) # 元に戻す
all.equal(real_space, Re(back_transform)) # OK
sum(abs(Im(back_transform))) # ほぼ 0 → 実数

atoms2 <- atoms[, 1:6] # メチルを外して、ベンゼンにする
real_space2 <- matrix(0, N, N)
scale <- N / 5 # [1, N] を [-2, 2] にマップする
for (i in 1:ncol(atoms2)) {
  ax <- atoms2[1, i]
  ay <- atoms2[2, i]
  print(ax**2 + ay**2)
  for (j in 1:N) {
    for (k in 1:N) {
      x <- j / scale - 2
      y <- k / scale - 2
      dist2 <- (x - ax) ** 2 + (y - ay) ** 2
      real_space2[j, k] <- real_space2[j, k] + 10 * exp(-dist2 / 0.2)
    }
  }
}
image(real_space2, asp=1)
reci_space2 <- fft(real_space2) / N / N

# トルエン由来の絶対値と、ベンゼン由来の位相を組み合わせて逆変換
combined <- matrix(complex(modulus=Mod(reci_space), argument=Arg(reci_space2)), N, N)
# その逆
combined2 <- matrix(complex(modulus=Mod(reci_space2), argument=Arg(reci_space)), N, N)

par(mfrow=c(2, 2), mar=c(3,3,3,3))
image(Re(fft(reci_space, inverse=T)), asp=1, main='Toluene') 
image(Re(fft(reci_space2, inverse=T)), asp=1, main='Benzene') 
image(Re(fft(combined, inverse=T)), asp=1, main='Abs=Tol, Arg=Bz')
image(Re(fft(combined2, inverse=T)), asp=1, main='Abs=Bz, Arg=Tol')
par(mfrow=c(1, 1))

実行結果はこちら。

# 2014/3/2 以前のコードはバグっていたので修正。さらに結果画像を掲載。