ft_exafs_scan.py - EXAFSデータ変換ツール
プログラムの動作
ft_exafs_scan.py は、EXAFS (Extended X-ray Absorption Fine Structure) 測定で得られたk空間データに対して、フーリエ変換 (FFT: Fast Fourier Transform) または最大エントロピ法 (MEM: Maximum Entropy Method) を適用し、動径分布関数 (R空間) を計算・解析するためのPythonスクリプトです。
プログラムの目的: EXAFSデータをk空間からR空間へ変換し、原子間距離情報（動径分布関数）を抽出すること。特に、EXAFS解析特有のR軸補正、窓関数処理、k重み付けをサポートし、複数のパラメータ設定でのスキャン解析機能を提供することで、解析の効率化と最適化を支援します。
主な機能:
データ入力:
.xlsx,.csv,.txt,.dat形式のファイルからk軸と信号データを読み込みます。EXAFS R軸補正: 通常のフーリエ変換の周波数軸をEXAFS解析に特化したR軸（実空間）へ変換します。
FFT計算: 高速フーリエ変換を実行し、R空間における実部、虚部、振幅、パワースペクトルを出力します。
MEM計算: 最大エントロピ法を用いてR空間スペクトル（パワースペクトル密度: PSD）を計算します。MEMの自己回帰 (AR) モデル次数は、AIC (Akaike Information Criterion) / BIC (Bayesian Information Criterion) 基準または固定値で選択可能です。
データ前処理:
k重み付け: 入力信号に \(k^{\text{korder}}\) を乗算してk重み付けを適用します。
Hanning窓関数: 指定されたk範囲とテーパー幅でEXAFSスタイルのHanning窓関数を適用し、変換範囲を限定します。
平均値除去: k重み付け後に信号の平均値を差し引くオプションを提供します。
ゼロパディング: FFT/MEM計算のNFFT長を制御し、R空間の分解能を調整します。
パラメータスキャン: 特定のパラメータ（MEM次数、Hanning窓範囲/幅、k重み付け、NFFT長）を段階的に変化させ、複数の変換結果を重ね合わせて比較する機能を提供します。
正規化: 変換結果の曲線を最大値で正規化し、形状の比較を容易にします。
結果の出力:
計算結果の詳細をExcelファイル (
.xlsx) に保存します。入力データと変換結果のグラフをMatplotlibで表示し、画像ファイルとして保存する機能も持ちます。
解決する課題:
EXAFS解析において、R空間への変換は窓関数の設定、k重み付け、変換手法（FFT/MEM）の選択、MEMにおけるARモデル次数の決定など、多くのパラメータに依存します。ft_exafs_scan.py はこれらの複雑な設定を一元的に管理し、特にパラメータスキャン機能を通じて、最適な解析条件を見つけるプロセスを効率化し、結果の解釈を支援します。
原理
EXAFS R軸変換の原理
EXAFS振動関数 \(\chi(k)\) は、光電子の散乱経路に依存し、\(k\) 空間で振動する関数として表されます。一般に \(\chi(k)\) は、シングル散乱経路に対して \(A(k) \sin(2kR + \phi(k))\) の形でモデル化されます。ここで \(R\) は原子間距離、\(\phi(k)\) は位相シフトです。
通常のフーリエ変換 (FT) では、時間軸データ \(f(t)\) を角周波数軸 \(\omega\) へ変換する際に、\(\omega = 2\pi \cdot \text{frequency}\) の関係があります。しかし、EXAFS解析においては、フーリエ変換後の軸を物理的な原子間距離 \(R\) と対応させる慣習があります。
ft_exafs_scan.py では、numpy.fft.rfftfreq(nfft, d=dk) によって得られる周波数 freq (単位は \(A\)) に対して、以下の式を用いてR軸を定義しています。
これは、\(\chi(k)\) のフーリエ変換が、通常のFTで得られる周波数成分の2倍の物理的な距離に対応するというEXAFSの慣例に基づいています。具体的には、\(\sin(2kR + \phi(k))\) のような項のフーリエ変換のピークは、通常の周波数軸で \(2R\) に対応する位置に現れます。しかし、EXAFSのR空間プロットでは直接原子間距離 \(R\) を表示するため、周波数軸を \(\pi\) 倍することで、この変換を実現しています。
Hanning窓関数の原理
フーリエ変換やMEMスペクトル計算を行う際、データが無限に連続していると仮定されます。しかし実際の測定データは有限な範囲 \([k_{min}, k_{max}]\) でしか得られません。この有限なデータに変換を適用すると、データの開始点と終了点での不連続性が人工的なサイドローブ（ギブス現象）を引き起こし、R空間スペクトルの品質を低下させる可能性があります。
Hanning窓関数は、この不連続性の影響を軽減するために使用される重み付け関数です。データの両端で値を0に近づけ、中央で1となる滑らかな形状をしています。ft_exafs_scan.py では、EXAFS解析で一般的に用いられる以下のスタイルのHanning窓関数を実装しています。
与えられたk範囲 \([k_1, k_2]\) とテーパー幅 \(w\) (hwidth) に対して、窓関数 \(W(k)\) は以下のように定義されます。
\(k < k_1\) または \(k > k_2\) の場合: \(W(k) = 0\)
\(k_1 \le k < k_1+w\) の場合（左側のテーパー領域）: $\(W(k) = 0.5 \left( 1 - \cos \left( \pi \frac{k - k_1}{w} \right) \right)\)$
\(k_1+w \le k \le k_2-w\) の場合（平坦領域）: \(W(k) = 1.0\)
\(k_2-w < k \le k_2\) の場合（右側のテーパー領域）: $\(W(k) = 0.5 \left( 1 + \cos \left( \pi \frac{k - (k_2 - w)}{w} \right) \right)\)$ この窓関数を信号に乗じることで、有限範囲のデータでも変換によるアーティファクトを抑制し、より信頼性の高いR空間スペクトルを得ることができます。
FFT (高速フーリエ変換) の原理
FFTは、離散フーリエ変換 (DFT) を高速に計算するアルゴリズムです。ft_exafs_scan.py では numpy.fft.rfft を使用しています。これは実数入力信号に対するDFTを効率的に計算し、その結果は複素数で得られます。
FFTの結果は、積分スケールを模倣するためにサンプリング間隔 dk で乗算されます。これにより、FTの振幅がデータ中の信号強度に比例するようになります。
ここで、signal_used は窓関数が適用された信号、n はゼロパディング後のデータ点数 (nfft)、dk はk軸のサンプリング間隔です。
FFT結果からは、実部、虚部、振幅 (np.abs(FT))、パワースペクトル (np.abs(FT)**2) が導出されます。
MEM (最大エントロピ法) の原理
MEMは、データからパワースペクトル密度 (PSD) を推定するための強力な手法の一つです。特にデータ点数が少ない場合や、FFTでは検出が難しい鋭いピークを持つスペクトル成分を高い分解能で検出するのに適しています。MEMは自己回帰 (AR) モデルに基づいており、データが過去の値の線形結合とノイズで表されると仮定します。
ここで \(X_t\) は時刻 \(t\) での信号、\(\alpha_i\) はAR係数、\(p\) はARモデルの次数、\(\epsilon_t\) はホワイトノイズです。
ft_exafs_scan.py では以下のライブラリとアルゴリズムを使用しています。
spectrum.pburg: Burg法を用いてARモデルの係数を推定し、そのモデルからパワースペクトル密度を計算します。Burg法は自己相関関数を直接計算せずにAR係数を推定するため、短いデータセットでも比較的安定した結果が得られます。ARモデル次数選択: ARモデルの次数 \(p\) はスペクトル推定の品質に大きく影響します。適切な次数を選択するために、以下の情報量基準が用いられます。
AIC (Akaike Information Criterion): モデルの適合度と複雑さのバランスを取る基準で、AIC値が最小となる次数が最適とされます。
BIC (Bayesian Information Criterion): AICと同様にモデル選択基準ですが、より複雑なモデルにペナルティを課すため、一般にAICよりも低い次数を選択する傾向があります。
PACF (Partial Autocorrelation Function): 偏自己相関関数は、特定のラグにおける自己相関から、間のラグの影響を取り除いたもので、ARモデルの次数を視覚的に判断する手がかりとなります。PACFが急激にゼロに落ちるラグがARモデルの次数として適切であると考えられます。
スクリプトは、デフォルトで preset_order に基づくか、--order-method で指定されたAICまたはBIC基準によってAR次数を自動的に選択します。
必要な非標準ライブラリとインストール方法
ft_exafs_scan.py を実行するには、以下の非標準ライブラリが必要です。
numpy: 数値計算pandas: データフレーム操作、Excelファイル入出力matplotlib: グラフ描画spectrum: 最大エントロピ法 (MEM) スペクトル計算statsmodels: 時系列分析、ARモデル次数選択 (AIC/BIC、PACF)
これらのライブラリは pip コマンドでインストールできます。PandasがExcelファイルを扱うために openpyxl も必要になる場合があります。
pip install numpy pandas matplotlib spectrum statsmodels openpyxl
必要な入力ファイル
プログラムは、k軸データと対応する信号データを含むファイルを入力として受け取ります。
サポートされるファイル形式:
Excelファイル (
.xlsx)CSVファイル (
.csv)テキストファイル (
.txt,.dat)
データ構造:
データは2次元のテーブル形式である必要があります。
デフォルトでは、0列目 がk軸の値 (
k [A^-1])、1列目 が信号の値 (signal, 例: \(\chi(k)\), \(k^2\chi(k)\) など) であると期待されます。--xcolおよび--ycolオプションを使用して、異なる列インデックス（0-based）を指定できます。Excelファイルの場合、
--sheetオプションでシート名またはインデックスを指定できます。
データの要件:
k軸は単調増加である必要があります（単調でない場合は警告が表示され、自動的にソートされます）。
k軸のサンプリング間隔
dkはほぼ均一であると仮定されます（不均一な場合は警告が表示されます）。変換を行うのに十分なデータ点数が必要です（最低4点）。
生成される出力ファイル
プログラムは、解析結果をExcelファイルとオプションでグラフ画像ファイルとして出力します。
Excel出力ファイル
出力ファイル名は、入力ファイル名、変換モード、k重み付け次数、およびスキャンパラメータに基づいて自動生成されます（例: data_exafs_fft_korder=0_scan_hwidth.xlsx）。--outfile オプションで任意のファイル名を指定することも可能です。
生成されるExcelファイルには、以下のシートが含まれます。
Inputシート:k_A^-1: 入力k軸データ。signal_raw: 読み込まれた生の信号データ。signal_k^{korder:g}_weighted: k重み付け (korder) が適用された信号。signal_centered: k重み付け後、オプションで平均値が除去された信号。hanning_window: 適用されたHanning窓関数の値。signal_used: 最終的にFFT/MEMに使用された信号（重み付け、中心化、窓関数適用後）。
FFTシート (モードが "fft" または "both" の場合):R_A: EXAFS R軸 [Å]。FT_real: FFTの実部。FT_imag: FFTの虚部。FT_amplitude: FFTの振幅 (|FT|)。FT_power: FFTのパワースペクトル (|FT|^2)。FT_real_norm,FT_imag_norm,FT_amplitude_norm,FT_power_norm:--normalize maxオプションが指定された場合の正規化された値。
MEMシート (モードが "mem" または "both" の場合):Frequency_MEM_cycles_per_sample: MEMで計算された周波数 (cycles/sample)。R_A: EXAFS R軸 [Å]。MEM_PSD: MEMパワースペクトル密度。MEM_PSD_norm:--normalize maxオプションが指定された場合の正規化されたPSD。
ScanSummaryシート (パラメータスキャンが実行された場合):scan_parameter: スキャンされたパラメータの名前。scan_value: 各スキャンステップでのパラメータ値。label: プロットに使用されるラベル。mode: 変換モード。y_name: プロットされたy軸のデータ種別。order: MEMのAR次数 (MEMモードの場合)。npadding: NFFT長。korder: k重み付け次数。hrange: Hanning窓のk範囲。hwidth: Hanning窓のテーパー幅。peak_R_A: ピークR値 [Å]。peak_y: ピークのy値。normalize: 正規化方法。peak_R_A_norm: 正規化されたピークR値 [Å]。peak_y_norm: 正規化されたピークのy値。
ScanCurvesシート (パラメータスキャンが実行された場合):R_A: 各スキャン曲線に対応するR軸データ。[label]: 各スキャンステップで計算された正規化済みYデータ (例:order=10,hwidth=0.5など)。raw_[label]:--normalizeがnoneでない場合、正規化前の生Yデータ。異なる
npaddingの値が使用されR軸が異なる場合、R_A_[label]とY_[label]のペアが列として追加されます。
Summaryシート:item: 設定項目や計算結果の要約項目。value: 各項目の値（入力ファイル名、モード、データ点数、krange、hwidth、AR次数、ピーク情報など）。
グラフ画像ファイル
--savefig [ファイル名] オプションが指定された場合、Matplotlibによって生成されたグラフは指定されたファイル名で画像ファイル（例: .png, .jpg）として保存されます。このグラフには、入力信号、窓関数、およびR空間変換結果（FFT振幅/パワー/実虚部、またはMEM PSD）が表示されます。パラメータスキャンが実行された場合は、すべてのスキャン結果が同じグラフ上に重ねてプロットされます。
コマンドラインでの使用例 (Usage)
usage: ft_exafs_scan.py [-h] [--mode {fft,mem,both}] [--sheet SHEET] [--xcol XCOL] [--ycol YCOL] [--npadding NPADDING]
[--korder KORDER] [--hrange KSTART KEND] [--hwidth HWIDTH] [--remove-mean] [--max-order MAX_ORDER]
[--order-method {bic,aic,preset}] [--order ORDER]
[--scan_order START STEP N] [--scan_hmin START STEP N] [--scan_hmax START STEP N]
[--scan_hwidth START STEP N] [--scan_korder START STEP N] [--scan_npadding START STEP N]
[--yscale {linear,log}] [--fft-plot {amplitude,power,realimag}] [--normalize {none,max}]
[--xlim XMIN XMAX] [--outfile OUTFILE] [--savefig SAVEFIG] [--no-show]
infile [preset_order]
EXAFS-oriented FFT/MEM transform. Input column 0 is k [A^-1], column 1 is chi(k)-like signal.
positional arguments:
infile Input file: .xlsx, .csv, .txt, or .dat
preset_order Fallback/fixed AR model order for MEM. Default: 8.
options:
-h, --help show this help message and exit
--mode {fft,mem,both}
Transform mode: fft, mem, or both. Default: fft.
--sheet SHEET Excel sheet name or index. Default: 0.
--xcol XCOL k-axis column index, 0-based. Default: 0.
--ycol YCOL signal column index, 0-based. Default: 1.
--npadding NPADDING FFT/MEM NFFT length. If smaller than data length, data length is used. Default: 1024.
--korder KORDER Apply k^korder to input signal before transform. If input is already k^2 chi(k), use 0. Default: 0.0.
--hrange KSTART KEND Hanning window k range [A^-1], e.g. --hrange 2 12. Default: None.
--hwidth HWIDTH Hanning taper width [A^-1]. Default: 0.5.
--remove-mean Subtract mean after k-weighting and before windowing.
--max-order MAX_ORDER
Maximum AR order tested by AIC/BIC. Default: 40.
--order-method {bic,aic,preset}
AR order selection method for MEM. Default: preset.
--order ORDER Fixed AR order for MEM. Overrides positional preset_order and --order-method.
--scan_order START STEP N
Scan MEM order. Example: --scan_order 10 10 5 gives 10,20,30,40,50.
--scan_hmin START STEP N
Scan lower bound of --hrange.
--scan_hmax START STEP N
Scan upper bound of --hrange.
--scan_hwidth START STEP N
Scan Hanning taper width.
--scan_korder START STEP N
Scan k-weight order.
--scan_npadding START STEP N
Scan zero-padding NFFT length.
--yscale {linear,log}
Y-axis scale for plot. Default: linear.
--fft-plot {amplitude,power,realimag}
FFT plot type. Default: amplitude.
--normalize {none,max}
Normalize each plotted/saved comparison curve. 'max' divides each curve by max(abs(y)). Default: none.
--xlim XMIN XMAX R-axis plot range, e.g. --xlim 0 6. Default: None.
--outfile OUTFILE Output Excel file. Default: auto-generated.
--savefig SAVEFIG Save figure filename, e.g. result.png. Default: "".
--no-show Do not show matplotlib window.
コマンドラインでの具体的な使用例
以下の例では、data.xlsx という架空のEXAFSデータファイル（0列目にk、1列目に\(\chi(k)\)）を使用することを想定しています。
通常のFFT計算 入力がすでに \(k^2\chi(k)\) の場合、
--korder 0を指定します。Hanning窓はk=2から12の範囲に、テーパー幅0.5 \(A^{-1}\) で適用されます。python ft_exafs_scan.py data.xlsx --mode fft --korder 0 --hrange 2 12 --hwidth 0.5 --savefig fft_result.png
実行結果:
data_exafs_fft_korder=0.xlsxというExcelファイルが生成され、FFTの振幅、実部、虚部、パワーがシートに保存されます。また、fft_result.pngに結果のグラフが保存され、ウィンドウが表示されます。MEM計算 (固定次数) MEMモードでARモデル次数を20に固定して計算します。
python ft_exafs_scan.py data.xlsx 20 --mode mem --korder 0 --hrange 2 12 --hwidth 0.5 --outfile mem_order20.xlsx
実行結果:
mem_order20.xlsxというExcelファイルが生成され、MEMのPSDがシートに保存されます。グラフにはMEM PSDが表示されます。MEM次数スキャン MEMのARモデル次数を10から始まりステップ10で5回（10, 20, 30, 40, 50）スキャンします。結果は最大値で正規化され、Excelファイルとグラフに重ねて表示されます。
python ft_exafs_scan.py data.xlsx --mode mem --scan_order 10 10 5 --korder 0 --hrange 2 12 --hwidth 0.5 --normalize max --savefig mem_order_scan.png
実行結果:
data_exafs_mem_order=XX_korder=0_scan_order.xlsx(XXは選択されたMEM次数) というExcelファイルが生成され、ScanSummaryとScanCurvesシートにスキャン結果がまとめられます。mem_order_scan.pngには正規化されたMEM PSDが次数ごとに重ねてプロットされます。Hanning窓の幅スキャン FFTモードでHanning窓のテーパー幅を0.2から始まりステップ0.2で5回（0.2, 0.4, 0.6, 0.8, 1.0）スキャンします。
python ft_exafs_scan.py data.xlsx --mode fft --scan_hwidth 0.2 0.2 5 --korder 0 --hrange 2 12 --normalize max --savefig hwidth_scan.png --no-show
実行結果: ExcelファイルにはHanning窓幅スキャン結果が保存されます。
hwidth_scan.pngには正規化されたFFT振幅が各テーパー幅ごとに重ねてプロットされます。--no-showオプションにより、グラフウィンドウは表示されずにファイルに保存されるのみとなります。k重み付けとR軸範囲指定 \(\chi(k)\) データに対して \(k^2\) 重み付けを行い、R軸のプロット範囲を0から6 \(A\) に設定します。
python ft_exafs_scan.py data.xlsx --mode fft --korder 2 --hrange 3 10 --hwidth 0.8 --xlim 0 6
実行結果:
data_exafs_fft_korder=2.xlsxが生成され、プロットのR軸範囲が0-6 \(A\) に限定されます。両モードでの比較プロットとAIC次数選択 FFTとMEMの両方を計算し、MEMの次数はAIC基準で自動選択します。
python ft_exafs_scan.py data.xlsx --mode both --order-method aic --korder 1 --hrange 2.5 13 --hwidth 0.7 --savefig both_aic.png
実行結果:
data_exafs_fft_mem_order=XX_korder=1.xlsx(XXはAICで選択された次数) が生成されます。グラフにはFFT振幅とMEM PSDの両方が表示されます。