nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

EXAFS振動のフーリエ変換について

「EXAFS振動をフーリエ変換する」というフレーズはX線吸収分光をやっていると腐る程見かけるが、実際にはArtemisやらLarchやらのソフトで流し込んで、実際に気にするところはパラメータのフィッティングだけということがほとんどな気がする。
ここでは「そもそもフーリエ変換するとはどういうことか?」について触れて行きたいと思う。

EXAFSについては、田渕先生の以下のスライドが非常に分かりやすい。
http://titan.nusr.nagoya-u.ac.jp/Tabuchi/BL5S1/lib/exe/fetch.php?media=nssr-long-201406.pdf

EXAFS振動 \chi( k )を超簡略化して、一つのsin関数(かつ余計な項(位相シフト)を含まない)で与えられるとする。これは最も理想的なEXAFS振動である。
「理想的な」という意味は、フーリエ変換(波数 k \rightarrow 距離 R)すると、原子間距離(求めたい量) R_0のところで(原理的に)厳密に鋭いピークを持つ。(位相シフトなどゴチャゴチャ入ると、 R_0からピーク位置が少しズレる)

 \displaystyle
\chi( k ) = \sin{( 2 k R_0 )}

("2"が付いているのは距離 R_0 を電子が「行って返って来る」ことに依る)

これをいきなりフーリエ変換する前に、通常のやり方(信号解析(時間 t \rightarrow 周波数 f))を復習して、比較出来るようにする。
周波数 f_0の信号 s(t)がsin関数で次のように与えられているとする。

 \displaystyle
s( t ) = \sin{( 2 \pi f_0 t ) } = \sin{( \omega_0 t ) }

これをpythonを使ってフーリエ変換してみる。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

N = 600
delta_t = 1.0 / 600.0

min_t = 0.0
max_t = N * delta_t
t = np.linspace( min_t, max_t, N )

f0 = 10
s = np.sin( 2. * np.pi * f0 * t )

plt.plot( t, s, "-o" )
plt.savefig( "s.png" )
plt.close()

tilde_s = fft( s )
tilde_s = fftshift( tilde_s )
tilde_t = fftfreq( N, delta_t )
tilde_t = fftshift( tilde_t )

plt.plot( tilde_t, 1.0/N * np.abs( tilde_s ), "-o" )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

plt.plot( tilde_t, 2.0 * 1.0/N * np.abs( tilde_s ), "-o" )
plt.xlim( 0, 50 )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

f:id:koideforest:20180716015720p:plain
f:id:koideforest:20180716020903p:plain
f:id:koideforest:20180716020921p:plain

  •  f_0 = 10に対応したピークがキチンと得られている。
  • sin関数のフーリエ変換なので、負の振動数成分と対で現れる( cos関数も同様)。
  • 規格化は 1/N Nはメッシュ数(今の場合、 N = 600)。
  • そのため、正の振動成分のみをプロットする場合、二倍しないと一見規格化が保たれていないように見える。
  • フーリエ変換後の \tilde{t}の範囲は [ - 1 / 2 \Delta_t, \,  1 / 2 \Delta_t ](今の場合、 \Delta_t = 1/600、[ -300, 300 ])。
  • これに対し、 \Delta_{ \tilde{t} } = 1 / \Delta_t N (今の場合、 \Delta_{ \tilde{t} } = 1 )。
  • したがって、幾ら細かく( \Delta_t \ll 1 )データを取っても、元データのデータ範囲  N \Delta_tが変わらないとフーリエ変換後のメッシュは細かくならない。

単純にフーリエ変換すると、固有振動数 f_0の位置にピークが得られたわけだが、これはつまり横軸が振動数( \tilde{t} = f)に変換されており、角振動数 \omegaではない。
 \omegaに変換するには、以下の2パターンが考えられる。

1.  \omega_0 / 2 \pi の位置にピークが出て来ると解釈する。

 \displaystyle
s( t ) = \sin( \omega_0 t ) = \sin \left( 2 \pi \frac{\omega_0}{2\pi} t \right )
より f_0 = \omega_0 / 2 \pi と思える。

2. 横軸 \tilde{t} 2\piを掛けて、 \omegaに変えてしまう。( \omega = 2 \pi f


さて、ここでEXAFSの話に戻る。 \chi( k ) = \sin( 2 k R_0 ) = \sin( R_0 ( 2 k ) ) だったので、この形は \sin( \omega_0 t )に近い。
やりたいことは、「 \chi( k )フーリエ変換して、 R_0にピークが出る動径分布関数を得る」ことである。なので、

  1.  kではなく、 2kフーリエ変換する。
  2.  R_0 f_0ではなく \omega_0に対応しているので、横軸 \tilde{2k} 2\piを掛ける。

これを踏まえて、 \chi(k) から動径分布関数を得る。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

min_k = 2
max_k = 12
delta_k = 0.1
k = np.arange( min_k, max_k, delta_k )
N = len( k )

R0 = 3
k2 = 2. * k
delta_k2 = 2. * delta_k
chi = np.sin( R0 * k2 )
plt.plot( k, chi, "-o")
plt.savefig( "chi.png" )
plt.close()

tilde_chi = fft( chi )
tilde_chi = fftshift( tilde_chi )
tilde_k2 = fftfreq( N, delta_k2 )
tilde_k2 = fftshift( tilde_k2 )
R = 2. * np.pi * tilde_k2
plt.plot( R, 2.0 / N * np.abs( tilde_chi ), "-o" )
plt.xlim( 0, 10 )
plt.ylim( 0, 1 )
plt.grid()
plt.savefig( "tilde_chi.png" )
plt.close()

f:id:koideforest:20180716031948p:plain
f:id:koideforest:20180716032005p:plain

 R_0=3の付近にピークが得られているのがわかる。しかし、綺麗なピークとは言い難い。
これは kの範囲を広げると改善される。
 k = [ 2, 24 ]にすれば、
f:id:koideforest:20180716033511p:plain
f:id:koideforest:20180716033527p:plain

しかし、現実的には k>12の実験データを得ることは難しい。
そこで、「窓関数が掛かっている」と見做して、ゼロ値のデータを足したものをフーリエ変換する。
f:id:koideforest:20180716034843p:plain
f:id:koideforest:20180716034922p:plain

ピーク強度が小さくなり、周りに小さなピーク構造が現れたが、その代わりにデータ点が増えて滑らかな関数を表すようになった。
これが、Artemis等のソフトで行われている有限フーリエ変換である。窓関数を掛けるのはデータの一部「のみ」でフーリエ変換するのではなく、追加でゼロデータを足しまくって無理矢理滑らかなピークにすることの宣言でもある。
これによって、一見、十分なデータ数があるように見えてしまい、「fittingパラメータを増やすためには広い kの範囲と十分大きな Nが必要」と言われても全然しっくり来なくなってしまっているのである。

実態は、汚いピーク(少ないデータ点)を与えたフーリエ変換であり、それに対してfittingをすると思うと、「もっとたくさん(フーリエ変換後のピーク近傍の)データ点が欲しい!」と思えるのではないだろうか?