oversampling / density modification による位相回復のデモ

R で oversampling / density modification による位相回復のデモを作った。

実験的に得られた F とランダムな位相から始める。これを実空間に戻して、マスクをかけてからフーリエ変換する。得られた位相(hopefully improved) と実験的 F を組み合わせて、反復する。

(クリックして最初から再生)

正解の画像は R でフーリエ変換を行い、位相の重要性を確認する - biochem_fan's note 参照。ソースは以下の通り。

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)
scale <- N / 8 # [1, N] を [-4, 4] にマップする
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 - 4
      y <- k / scale - 4
      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)

mask <- matrix(0, N, N)
for (j in 1:N) {
  for (k in 1:N) {
    x <- j / scale - 4
    y <- k / scale - 4
    dist2 <- x ** 2 + y ** 2
    if (dist2 < 6.5)
      mask[j, k] <- 1
  }
}
image(mask, asp=1)
image(mask * real_space, asp=1)


Friedel_mate <- function(x, N) {(N - (x - 1)) %% N + 1}

symmetrize_phase <- function(mat) {
  N <- dim(mat)[1]
  
  phases <- mat
  for (h in 1:floor(N / 2)) {
    for (k in 1:floor(N / 2)) {
      phases[h, k] <- -phases[Friedel_mate(h, N), Friedel_mate(k, N)]
    }
  }
}

# ffmpeg -r 30 -i oversampling-%04d.png -b 1000k -s 200x200  out.avi
#  or
# convert -delay 3 -loop 1 oversampling-0[0-5]*.png oversampling.gif

setwd("/tmp")
png("oversampling-%04d.png", width=200, height=200)
par(mai=c(0, 0, 0, 0))

set.seed(2)
reci_space <- fft(real_space) / N / N
correct_mod <- Mod(reci_space)
cur_phase <- symmetrize_phase(matrix(runif(N * N, max=2*pi), N))

for (i in 0:3000) {
  cur_reci <- matrix(complex(modulus=correct_mod, argument=cur_phase), N)
  back_transform <- fft(cur_reci, inverse=T)
  
  if (i %% 10 == 0) {
    image(Re(back_transform), asp=1)
    print(i)
  }
  cur_real <- back_transform * mask
  cur_phase <- Arg(fft(cur_real) / N / N)
#  cur_phase <- symmetrize_phase(cur_phase) # 実数の拘束はなくても動くようだ
}
dev.off()

実空間のCCを計算して経過をプロットしてみるのも面白いだろう。