「EXAFS振動をフーリエ変換する」というフレーズはX線吸収分光をやっていると腐る程見かけるが、実際にはArtemisやらLarchやらのソフトで流し込んで、実際に気にするところはパラメータのフィッティングだけということがほとんどな気がする。
ここでは「そもそもフーリエ変換するとはどういうことか?」について触れて行きたいと思う。
EXAFSについては、田渕先生の以下のスライドが非常に分かりやすい。
http://titan.nusr.nagoya-u.ac.jp/Tabuchi/BL5S1/lib/exe/fetch.php?media=nssr-long-201406.pdf
EXAFS振動を超簡略化して、一つのsin関数(かつ余計な項(位相シフト)を含まない)で与えられるとする。これは最も理想的なEXAFS振動である。
「理想的な」という意味は、フーリエ変換(波数 距離)すると、原子間距離(求めたい量)のところで(原理的に)厳密に鋭いピークを持つ。(位相シフトなどゴチャゴチャ入ると、からピーク位置が少しズレる)
("2"が付いているのは距離を電子が「行って返って来る」ことに依る)
これをいきなりフーリエ変換する前に、通常のやり方(信号解析(時間 周波数))を復習して、比較出来るようにする。
周波数の信号がsin関数で次のように与えられているとする。
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()
- に対応したピークがキチンと得られている。
- sin関数のフーリエ変換なので、負の振動数成分と対で現れる( cos関数も同様)。
- 規格化は。はメッシュ数(今の場合、)。
- そのため、正の振動成分のみをプロットする場合、二倍しないと一見規格化が保たれていないように見える。
- フーリエ変換後のの範囲は](今の場合、、[ ])。
- これに対し、(今の場合、)。
- したがって、幾ら細かく()データを取っても、元データのデータ範囲が変わらないとフーリエ変換後のメッシュは細かくならない。
単純にフーリエ変換すると、固有振動数の位置にピークが得られたわけだが、これはつまり横軸が振動数()に変換されており、角振動数ではない。
に変換するには、以下の2パターンが考えられる。
1. の位置にピークが出て来ると解釈する。
よりと思える。
2. 横軸にを掛けて、に変えてしまう。()
さて、ここでEXAFSの話に戻る。だったので、この形はに近い。
やりたいことは、「をフーリエ変換して、にピークが出る動径分布関数を得る」ことである。なので、
- ではなく、でフーリエ変換する。
- はではなくに対応しているので、横軸にを掛ける。
これを踏まえて、から動径分布関数を得る。
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()
の付近にピークが得られているのがわかる。しかし、綺麗なピークとは言い難い。
これはの範囲を広げると改善される。
]にすれば、
しかし、現実的にはの実験データを得ることは難しい。
そこで、「窓関数が掛かっている」と見做して、ゼロ値のデータを足したものをフーリエ変換する。
ピーク強度が小さくなり、周りに小さなピーク構造が現れたが、その代わりにデータ点が増えて滑らかな関数を表すようになった。
これが、Artemis等のソフトで行われている有限フーリエ変換である。窓関数を掛けるのはデータの一部「のみ」でフーリエ変換するのではなく、追加でゼロデータを足しまくって無理矢理滑らかなピークにすることの宣言でもある。
これによって、一見、十分なデータ数があるように見えてしまい、「fittingパラメータを増やすためには広いの範囲と十分大きなが必要」と言われても全然しっくり来なくなってしまっているのである。
実態は、汚いピーク(少ないデータ点)を与えたフーリエ変換であり、それに対してfittingをすると思うと、「もっとたくさん(フーリエ変換後のピーク近傍の)データ点が欲しい!」と思えるのではないだろうか?