peaksearch.py 技術ドキュメント

プログラムの動作

peaksearch.py は、X線回折（XRD）パターンなどのスペクトルデータから、自動的にピークを検出し、その結果を視覚化することを目的としたPythonプログラムです。データファイルから指定されたX軸とY軸のデータ範囲を読み込み、Savitzky-Golayフィルターを用いてデータを平滑化し、その導関数を計算します。これらの導関数と、強度や形状に関する閾値に基づいて、スペクトル中のピークを識別します。検出されたピークは、元のデータ、平滑化されたデータ、そして推定されたガウス関数フィッティングとともにグラフ上に表示されます。

主な機能は以下の通りです。

データ読み込み: Excel (.xlsx) や CSV などの表形式データファイルから、X軸とY軸に対応するデータを読み込みます。
データ範囲指定: X軸の最小値と最大値を指定して、データの範囲をフィルタリングできます。
平滑化と導関数計算: Savitzky-Golayフィルターを用いて、データセットを平滑化し、さらに1次、2次、3次導関数を計算します。
ピーク検出: 3次導関数のゼロ交差をピーク候補として識別し、平滑化されたデータの強度、1次導関数の絶対値と平滑化されたY値の比、2次導関数の符号を用いて、候補をフィルタリングします。
ピーク幅の推定: 検出されたピークに対し、3次導関数のゼロ交差から半値半幅 (HWHM) を推定します。
結果の可視化: Matplotlib を使用して、元のデータ、平滑化されたデータ、検出されたピークとその推定ガウス関数フィッティングをプロットします。デバッグモードでは、各導関数も同時にプロットされます。
ログ出力: プログラムの実行パラメータ、ピーク検出プロセスに関する情報が標準出力およびログファイルに出力されます。

原理

peaksearch.py におけるピーク検出は、主にSavitzky-Golayフィルターによるデータの平滑化と導関数計算、およびそれらの導関数に基づいたゼロ交差検出、そして閾値によるフィルタリングを組み合わせて行われます。

Savitzky-Golay (S-G) フィルター

S-Gフィルターは、データ平滑化や導関数計算に用いられる手法です。移動平均フィルターと異なり、ノイズ除去と同時にデータのピーク形状や幅を維持する特性があります。このフィルターは、指定されたウィンドウサイズ内のデータ点に対し、最小二乗法を用いて指定された次数の多項式をフィッティングします。そして、ウィンドウ中央のデータ点を、フィッティングされた多項式の値（またはその導関数の値）で置き換えます。このプロセスをデータ全体にわたって繰り返します。プログラムでは、このフィルターを用いて入力データ $y(x)$ を平滑化し、$y_{smooth}(x)$ を得ます。さらに、$y_{smooth}(x)$ の1次導関数 $y'(x)$、2次導関数 $y''(x)$、3次導関数 $y'''(x)$ もこのフィルターで計算されます。

ピーク検出アルゴリズム

データの前処理:
- 入力データ $y(x)$ をSavitzky-Golayフィルターで平滑化し、$y_{smooth}(x)$ を計算します。
- 同様に、1次導関数 $y'(x)$、2次導関数 $y''(x)$、3次導関数 $y'''(x)$ を計算します。
ピーク候補の特定:
- ピークの候補は、3次導関数 $y'''(x)$ のゼロ交差点として識別されます。これは、ピークの変曲点近傍で3次導関数が符号を変える特性を利用しています。
ピークのフィルタリング: 検出されたピーク候補に対し、以下の基準でフィルタリングを行います。
- 強度閾値: ピーク中心での平滑化された強度 $y_{smooth}(x_{peak})$ が、設定された閾値 $threshold$ ($cparams.threshold$) より小さい場合、その候補は無視されます。これは、弱すぎるピークやノイズを除外するためです。 $$y_{smooth}(x_{peak}) < threshold$$
- 傾き閾値: ピーク中心での1次導関数の絶対値と平滑化された強度の比 $|y'(x_{peak}) / y_{smooth}(x_{peak})|$ が、設定された閾値 $diff1\_ratio\_th$ ($cparams.ydiff1\_threshold$ に基づいて計算) より大きい場合、その候補は無視されます。これは、非常に急峻な変化やノイズによる偽陽性を排除するためです。 $$\frac{|y'(x_{peak})|}{y_{smooth}(x_{peak})} > diff1\_ratio\_th$$
- 凹凸判定: ピーク中心での2次導関数 $y''(x_{peak})$ が正の場合 ($y''(x_{peak}) > 0$)、それはピークではなく谷（最小値）であるため、無視されます。ピークでは2次導関数は負またはゼロに近いはずです。
ピーク幅の推定: フィルタリングを通過したピークに対して、半値半幅 (HWHM) に相当する幅 _w を推定します。これは、ピーク中心から左右の最も近い3次導関数のゼロ交差点までの距離として計算されます。

ガウス関数による表示

検出されたピークは、グラフ上でガウス関数としてプロットされます。ガウス関数の一般形は以下の通りです。

\[f(x) = I_0 \exp\left(-\frac{(x - x_0)^2}{a^2}\right)\]

ここで、$x_0$ はピーク中心の位置、$I_0$ はピークの最大強度、$a$ はピークの幅を決定するパラメータです。半値半幅 (HWHM) を $\Delta x_{1/2}$ とすると、$a$ との関係は次のようになります。

\[\begin{split} \frac{I_0}{2} = I_0 \exp\left(-\frac{(\Delta x_{1/2})^2}{a^2}\right) \\ \frac{1}{2} = \exp\left(-\frac{(\Delta x_{1/2})^2}{a^2}\right) \\ -\ln 2 = -\frac{(\Delta x_{1/2})^2}{a^2} \\ (\Delta x_{1/2})^2 = a^2 \ln 2 \\ \Delta x_{1/2} = a \sqrt{\ln 2} \\ a = \frac{\Delta x_{1/2}}{\sqrt{\ln 2}} \end{split}\]

プログラムでは、推定されたHWHMが _w として渡され、$\sqrt{\ln 2} \approx 0.832554611$ を用いて $a = \_w / 0.832554611$ として計算されます。これにより、検出されたピークの位置と幅に対応するガウス関数が生成され、可視化されます。

必要な非標準ライブラリとインストール方法

peaksearch.py の実行には、以下の非標準ライブラリが必要です。

numpy: 数値計算の基盤ライブラリ。
pandas: データ構造とデータ解析ツールを提供するライブラリ。特にExcelやCSVファイルの読み込みに使用されます。
scipy: 科学技術計算ライブラリ。特にSavitzky-Golayフィルター (scipy.signal.savgol_filter) の実装に利用されます。
matplotlib: グラフ描画ライブラリ。検出されたピークの可視化に使用されます。

これらのライブラリは、Pythonのパッケージマネージャー pip を使用してインストールできます。

pip install numpy pandas scipy matplotlib

また、peaksearch.py は、tklib というカスタムライブラリに依存しています。これは一般的な pip でインストールできるライブラリではないため、peaksearch.py と同じディレクトリ階層に tklib ディレクトリを配置するか、tklib ディレクトリがPythonのパス (PYTHONPATH) に含まれるように設定する必要があります。tklib ディレクトリには、tkutils.py, tkparams.py, tkvariousdata.py, tkapplication.py などのモジュールが含まれている必要があります。

必要な入力ファイル

peaksearch.py は、表形式のデータを含むファイルを入力として期待します。

ファイル形式: Excelファイル (.xlsx) または CSVファイル (.csv) など、pandas ライブラリがサポートする形式が利用可能です。
データ構造: ファイルはヘッダー行を持つ複数の列で構成されている必要があります。X軸データとY軸データに対応する列名をコマンドライン引数で指定します。

例 (ExcelファイルまたはCSVファイルの内容):

Angle,Intensity,Sample_ID
0,120.5,A
1,135.2,A
2,150.8,A
3,142.1,A
4,110.0,A
...
0,250.3,B
1,300.7,B
2,280.1,B
...

上記の例では、Angle がX軸、Intensity がY軸として使用できます。

生成される出力ファイル

peaksearch.py は、プログラムの実行ログをファイルとして出力します。

ファイル名: 入力ファイル名に基づいた名前が自動的に生成されます。例えば、入力ファイルが xrd-narrow.xlsx の場合、出力ログファイルは xrd-narrow-searched.xlsx となります。
内容: このファイルには、標準出力に表示されるすべての情報がリダイレクトされて保存されます。具体的には、以下の情報が含まれます。
- プログラムの実行パラメータ（モード、入力ファイル、X/Y軸名、X範囲、フィルターパラメータ、閾値など）。
- データの読み込み状況（データ点数など）。
- 各ピーク候補に対するフィルタリングの詳細（強度閾値、傾き閾値、凹凸判定の結果など）。
- 最終的に検出されたピークのインデックス、X値、推定される幅。

注意点として、この出力ファイルには検出されたピークのデータ自体（座標や強度、幅などをまとめた表形式のデータ）は直接含まれず、あくまで実行ログのテキストデータが保存されます。

コマンドラインでの使用例 (Usage)

peaksearch.py は、コマンドライン引数によって動作を制御します。基本的な実行コマンドと引数の説明は以下の通りです。

python peaksearch.py [mode] [infile] [xmin] [xmax] [xlabel] [ylabel] [threshold] [ydiff1_threshold] [norder] [nsmooth]

[mode] (文字列, デフォルト: 'peak search'):
- 'peak search' の場合、検出されたピークと平滑化された入力データのみを含むグラフが表示されます。
- 他の文字列（例: 'debug'）を指定した場合、平滑化されたデータと、その1次、2次、3次導関数がそれぞれ個別のサブプロットとして表示され、ピーク検出プロセスの詳細を視覚的に確認できます。
[infile] (文字列, デフォルト: 'xrd-narrow.xlsx'):
- ピーク検索を行う入力データファイルへのパス。Excel (.xlsx) または CSV (.csv) 形式をサポートします。
[xmin], [xmax] (浮動小数点数, デフォルト: -1e100, 1e100):
- X軸データの範囲を制限するための最小値と最大値。指定しない場合、または '*' や空文字列を指定した場合は、データ全体が使用されます。
[xlabel], [ylabel] (文字列, デフォルト: '0', '1'):
- 入力ファイル内でX軸およびY軸データに対応する列名。
[threshold] (浮動小数点数, デフォルト: 100.0):
- ピークと認識される平滑化されたY値の最小強度閾値。この値より小さいピーク候補は無視されます。
[ydiff1_threshold] (浮動小数点数, デフォルト: 1.0e-2):
- 1次導関数の絶対値と平滑化されたY値の比に関する閾値。この比がこの値より大きいピーク候補は、急峻すぎる変化やノイズとして無視される傾向があります。
[norder] (整数, デフォルト: 4):
- Savitzky-Golayフィルターでデータにフィッティングする多項式の次数。
[nsmooth] (整数, デフォルト: 11):
- Savitzky-Golayフィルターのウィンドウサイズ（フィルターを適用するデータ点の数）。この値は奇数である必要があります。偶数を指定した場合、プログラム内で自動的に1加算され、奇数に調整されます。norder より大きく設定する必要があります。

コマンドラインでの具体的な使用例

ここでは、peaksearch.py の具体的な使用例を2つ示します。入力ファイルとして xrd-narrow.xlsx が存在し、その中に 2Theta と Intensity という列があることを前提とします。

例1: 基本的なピーク検索

python peaksearch.py "peak search" "xrd-narrow.xlsx" "*" "*" "2Theta" "Intensity" 100.0 0.01 4 11

引数の説明:

"peak search": 通常のピーク検索モードで実行し、メインのプロットのみを表示します。
"xrd-narrow.xlsx": 入力データファイルとして xrd-narrow.xlsx を使用します。
"*": X軸の最小値および最大値を指定せず、全データ範囲を使用します。
"2Theta": X軸のデータとして、入力ファイル内の 2Theta 列を使用します。
"Intensity": Y軸のデータとして、入力ファイル内の Intensity 列を使用します。
100.0: 平滑化後のピーク強度が 100.0 未満のものはピークとして検出しません。
0.01: 1次導関数の絶対値とY値の比が 0.01 より大きいピークは、急峻すぎるとして除外されやすくなります。
4: Savitzky-Golayフィルターに4次の多項式を適用します。
11: Savitzky-Golayフィルターのウィンドウサイズを11点とします。

実行結果の説明:

このコマンドを実行すると、xrd-narrow.xlsx から 2Theta と Intensity のデータを読み込み、指定されたSavitzky-Golayフィルターのパラメータ（次数4、ウィンドウサイズ11）で平滑化処理を行います。その後、平滑化されたY値が100.0以上で、かつ1次導関数とY値の比が0.01以下の、3次導関数のゼロ交差をピークとして識別します。検出されたピークは、元のデータと平滑化されたデータ、さらに検出されたピークに対応するガウス関数フィッティングとともに、一つのMatplotlibグラフウィンドウに表示されます。プログラムの実行パラメータとピーク検出に関する詳細は、標準出力と xrd-narrow-searched.xlsx という名前で生成されるログファイルに記録されます。

例2: デバッグモードでの実行 (詳細な導関数プロット)

python peaksearch.py "debug_mode" "xrd-narrow.xlsx" "20.0" "30.0" "2Theta" "Intensity" 50.0 0.05 4 15

引数の説明:

"debug_mode": デバッグモードで実行します。これにより、平滑化されたデータに加えて、その1次、2次、3次導関数も個別のサブプロットとして表示されます。
"xrd-narrow.xlsx": 例1と同じ入力ファイルを使用します。
"20.0", "30.0": X軸のデータを 20.0 から 30.0 の範囲に限定します。
"2Theta", "Intensity": 例1と同じX, Y軸を使用します。
50.0: ピーク強度閾値を 50.0 に設定します。例1よりも低い強度のピークも検出対象となります。
0.05: 傾き閾値を 0.05 に設定します。例1よりも高い比率のピークも受け入れるため、より急峻なピークも検出されやすくなります。
4: Savitzky-Golayフィルターの多項式次数は例1と同じ 4 です。
15: Savitzky-Golayフィルターのウィンドウサイズを 15 点とします。例1より平滑化が強めにかかります。

実行結果の説明:

このコマンドは、xrd-narrow.xlsx からX軸範囲を 20.0 から 30.0 に絞ってデータを読み込みます。Savitzky-Golayフィルターは次数4、ウィンドウサイズ15で適用され、ピーク強度閾値50.0、傾き閾値0.05でピークを検出します。"debug_mode" が指定されているため、Matplotlibのグラフウィンドウには、元のデータと平滑化されたデータ（検出されたピークとガウス関数フィッティングを含む）に加えて、1次導関数、2次導関数、3次導関数がそれぞれ独立したサブプロットとして表示されます。これにより、各導関数の変化がピーク検出にどのように寄与しているかを詳細に分析できます。プログラムのログも同様に標準出力と xrd-narrow-searched.xlsx に出力されます。