deconvolution.py 技術ドキュメント
プログラムの動作
deconvolution.py は、測定データに畳み込まれた応答関数（主にガウス関数）を逆畳み込みすることによって、元のスペクトル形状を復元するためのPythonプログラムです。主に光電子分光（PES）などのスペクトルデータ解析を想定しています。
主な機能:
多様な逆畳み込みアルゴリズムのサポート:
Jacobi法: 繰り返し計算による逆畳み込み。収束性や安定性を制御するためのダンピング係数や平滑化機能を持つ。
Gauss-Seidel法: Jacobi法と同様の繰り返し計算だが、更新された値を即座に利用することで収束が速い可能性がある。
FFT (高速フーリエ変換) 法: 周波数領域での除算による効率的な逆畳み込み。
SciPyの
deconvolve関数: Pythonの科学計算ライブラリSciPyに実装されている逆畳み込み機能を利用。
畳み込みシミュレーション:
numpy.convolveを使用して、元のデータにフィルター関数を畳み込むことで、逆畳み込みの効果を確認できる。ガウス関数フィルターの生成: 半値幅（Wa）と有効範囲（Grange）を指定してガウス応答関数を生成する。
データの前処理:
平滑化: 移動平均または多項式フィット（Savitzky-Golayフィルターに類似）によるノイズ低減。
データ拡張: FFTやSciPyの
deconvolveで発生しがちな端点効果を抑制するため、データの両端をゼロや線形補間で拡張する。
結果の可視化:
matplotlibを用いて、元のデータ、畳み込み後のデータ、逆畳み込み後のデータをリアルタイムまたは最終的にグラフ表示する。CSV形式での入出力: 測定データの読み込みと結果の保存にCSV形式を使用する。
解決する課題:
測定データ（特に分光データ）は、測定装置の分解能や自然幅によって、本来のシャープなピークがブロードニング（広がり）を起こして観測されます。このプログラムは、ブロードニングを引き起こす応答関数（フィルター関数）が既知である場合に、観測データからその応答関数の影響を取り除き、より高分解能な元のスペクトル形状を推定することを目的としています。これにより、ピーク位置の精密な決定や、隠れた構造の検出が可能になります。
原理
このプログラムでは、いくつかの異なる逆畳み込みアルゴリズムが実装されています。基本的な畳み込みの概念から始め、各アルゴリズムの原理を説明します。
1. 畳み込みの基本
観測されるデータ \(I_{\text{conv}}(x)\) は、本来のデータ \(I_{\text{raw}}(x)\) と装置の応答関数（またはフィルター関数）\(H(x)\) の畳み込みとして表されます。
離散データの場合、これは以下のように近似されます。
逆畳み込みは、この \(I_{\text{conv},i}\) と \(H_{i-j}\) から \(I_{\text{raw},j}\) を推定する問題です。
2. ガウス関数フィルター
応答関数 \(H(x)\) として、このプログラムではガウス関数が使用されます。半値幅 Wa のガウス関数 \(G(x, x_0, W_a)\) は以下のように定義されます。
まず、プログラムでは半値幅 \(W_a\) を用いて正規化定数 \(A\) とスケール因子 \(a\) を計算します。 \(A = \frac{1}{W_a} \sqrt{\frac{\ln 2}{\pi}}\) \(a = \frac{W_a}{\sqrt{\ln 2}}\)
これらを用いて、ガウス関数は次の形式で計算されます。 $\(G(x, x_0, W_a) = A \cdot \exp\left(-\left(\frac{x - x_0}{a}\right)^2\right)\)$
ここで、\(x_0\) はピーク中心、\(W_a\) は半値幅です。プログラム内では、フィルター関数の中心 \(x_0\) は0として扱われ、\(H_{ij}\) が \((j-i) \cdot \text{xstep}\) におけるガウス関数として定義されます。
3. Jacobi法とGauss-Seidel法による繰り返し逆畳み込み
これらの方法は、上記の離散畳み込み方程式を \(I_{\text{raw}}\) について解くための繰り返しアルゴリズムです。プログラムでは、南茂夫 編著「科学計測のための波形データ処理」の繰り返し法をベースに、以下の更新式を用いています。
観測データ \(y_{\text{raw},i}\) と、現在の推定値 \(y_{\text{old},j}\) が与えられたとき、次の推定値 \(y_{\text{new},i}\) は以下のように計算されます。
まず、現在の推定値 \(y_{\text{old}}\) を用いて、畳み込み後の期待値 \(Hx_i\) を計算します。プログラムの Hij 関数は \(H_{j-i}\) に対応します。
$\(Hx_i = \sum_j H_{j-i} y_{\text{old},j}\)$
また、フィルター関数の総和 \(Sg_i\) を計算します。 $\(Sg_i = \sum_j H_{j-i}\)$
そして、中心要素 \(h = H_{0}\) を用いて、新しい推定値 \(y_{\text{new},i}\) を導出します。dump はダンピング係数です。
Jacobi法 (
deconvolute_jacobi): 新しい推定値 \(y_{\text{new},i}\) の計算に、常に一つ前の反復の推定値 \(y_{\text{old},j}\) を用います。つまり、ある反復内の全ての \(i\) について \(y_{\text{new},i}\) を計算し終えてから、その結果を次の反復の \(y_{\text{old},j}\) に更新します。Gauss-Seidel法 (
deconvolute_gauss_seidel): Jacobi法と異なり、新しい推定値 \(y_{\text{new},i}\) の計算において、既に同じ反復内で更新された \(y_j\) (\(j < i\)) を直ちに利用します。これにより、一般にJacobi法よりも収束が速いとされています。Gauss-Seidel法では、\(Hx_i\) の計算が以下のように変わります。 $\(Hx_i = \sum_{j < i} H_{j-i} y_{j} + \sum_{j \ge i} H_{j-i} y_{\text{old},j}\)\( ここで、\)y_{j}\( は現在の反復で既に計算された新しい値、\)y_{\text{old},j}$ は前の反復の値です。
これらの繰り返し法では、各反復後に平滑化処理（ノイズ抑制）やゼロ補正（負の値の除去）が適用されることがあります。
4. FFTによる逆畳み込み (deconvolute_fft)
畳み込み定理によれば、空間領域での畳み込みは、フーリエ変換された周波数領域での乗算に相当します。
したがって、逆畳み込みは周波数領域での除算によって行えます。
元のデータをFFTし、応答関数をFFTし、周波数領域で除算を行った後、逆フーリエ変換（IFFT）することで \(I_{\text{raw}}(x)\) を得ます。この方法では、データの長さが2のべき乗である必要があり、ゼロパディングによるデータ拡張が行われます。また、\(\hat{H}(\nu)\) が0に近い値をとる場合、除算が不安定になるため、ノイズの影響を受けやすいという欠点があります。
5. SciPyの scipy.signal.deconvolve (deconvolute_deconvolve)
SciPyライブラリの deconvolve 関数は、観測データと応答関数が与えられたときに、応答関数がどのようなデータに畳み込まれたかを推定します。これは最小二乗法などの数値的手法に基づいて実装されています。
6. データの前処理
平滑化 (Smoothing):
SmoothingBySimpleAverageは単純移動平均を、SmoothingByPolynomialFitはSavitzky-Golayフィルターに類似した多項式フィットに基づく重み付き移動平均を適用してノイズを低減します。データ拡張 (Extend Smooth): FFTや
scipy.signal.deconvolveを使用する際に、データの両端で発生する不連続性や端点効果を抑えるために、データの先頭にゼロや線形補間されたデータを追加します。
必要な非標準ライブラリとインストール方法
このプログラムは、以下の非標準ライブラリに依存しています。
NumPy: 数値計算を効率的に行うためのライブラリ。配列操作、フーリエ変換などに使用されます。
SciPy: 科学技術計算のためのライブラリ。信号処理モジュール（
scipy.signal.deconvolve）を使用します。Matplotlib: グラフ描画のためのライブラリ。結果の可視化に使用されます。
これらのライブラリは、Pythonのパッケージインストーラー pip を使用してインストールできます。
pip install numpy scipy matplotlib
必要な入力ファイル
プログラムは、CSV形式またはタブ区切り（.TXT 拡張子）のテキストファイルを読み込みます。
ファイル名: コマンドライン引数で指定します。デフォルトは
pes.csvです。ファイル形式:
CSV (.csv): カンマ (
,) 区切り。TXT (.txt): タブ (
\t) 区切り。
データ構造:
1行目はヘッダー行として扱われます（例:
Energy,Intensity）。2行目以降はデータ行であり、各行は少なくとも2つの数値カラムを持つ必要があります。
最初のカラムがX軸データ（例: エネルギー）、2番目のカラムがY軸データ（例: 強度）として扱われます。
.TXTファイルの場合、X軸データは読み込み時に符号が反転されます（PESデータにおける結合エネルギーの慣例に対応）。
データ範囲: コマンドライン引数
xminとxmaxで指定された範囲内のデータのみが読み込まれ、処理されます。
例 (pes.csv):
Energy,Intensity
-5.0,10.2
-4.9,11.5
-4.8,13.1
...
1.9,8.7
2.0,9.0
生成される出力ファイル
プログラムは、処理された逆畳み込み結果を新しいCSVファイルとして保存します。
ファイル名: 入力ファイル名に基づいて自動的に生成されます。
入力ファイルが
input.csvの場合、出力ファイルはinput-deconvoluted.csvとなります。入力ファイルが
data.txtの場合、出力ファイルはdata-deconvoluted.csvとなります。
ファイル内容: 出力CSVファイルには以下の3つのカラムが含まれます。
x: 処理後のX軸データ（エネルギーなど）。
y(input): 逆畳み込み処理の入力として使用されたY軸データ。前処理（平滑化、データ拡張、畳み込みなど）が適用された後のデータです。
y(deconvoluted): 逆畳み込み処理によって得られたY軸データ。
例 (input-deconvoluted.csv):
x,y(input),y(deconvoluted)
-4.5,10.0,9.5
-4.4,10.5,11.0
-4.3,11.2,13.2
...
コマンドラインでの使用例 (Usage)
deconvolution.py は、実行モードによって必要な引数が異なります。
1. 畳み込み (convolve) または FFT (fft) モードの場合:
python deconvolution.py file mode convmode smoothmode xmin xmax Wa Grange kzero klin
file: 入力データファイル名 (例:pes.csv)mode: 実行モード (convolveまたはfft)convmode:convolveまたはfftモードでの畳み込みの挙動 ('',full,same)''(空): デフォルトのfullと同等（NumPyのデフォルト）full: 出力は \((N+M-1)\) の長さになり、全ての畳み込み結果を含む。same: 出力は入力と同じ長さ \(N\) になるようにクリッピングされる。
smoothmode: 前処理の平滑化およびデータ拡張の方法。none,convolve,extend,averageを+で結合した文字列。none: 何も適用しない。convolve: 畳み込みによる平滑化（別途フィルターが適用されるわけではない。データ拡張の前にconvolve関数を通すことを意味する）extend: データ端をゼロまたは線形補間で拡張。average: 移動平均による平滑化。例:
convolve+extend,average+extend
xmin: 解析範囲のX軸最小値 (float)xmax: 解析範囲のX軸最大値 (float)Wa: ガウスフィルターの半値幅 (float, eV)Grange: ガウスフィルターを考慮するWaの倍率 (float)kzero: データ拡張時に先頭に追加するゼロデータの幅をnwの倍数で指定 (nwはガウスフィルター半分の幅) (float)klin: データ拡張時に線形補間する幅をnwの倍数で指定 (float)
2. Jacobi法 (jacobi) または Gauss-Seidel法 (gs) モードの場合:
python deconvolution.py file mode xmin xmax Wa dump nmaxiter eps nsmooth zeroc
file: 入力データファイル名 (例:pes.csv)mode: 実行モード (gsまたはjacobi)gs: Gauss-Seidel法jacobi: Jacobi法
xmin: 解析範囲のX軸最小値 (float)xmax: 解析範囲のX軸最大値 (float)Wa: ガウスフィルターの半値幅 (float, eV)dump: ダンピング係数 (float)。収束を制御する。nmaxiter: 最大反復回数 (int)。eps: 収束判定に使用する相対誤差の閾値 (float)。nsmooth: 各反復後に適用する平滑化の窓幅 (int)。zeroc: ゼロ補正 (0: 無効,1: 有効)。負の値を0にクリップする。
コマンドラインでの具体的な使用例
例1: Gauss-Seidel法による逆畳み込み
pes.csv という入力ファイルに対して、Gauss-Seidel法を用いて逆畳み込みを実行します。
python deconvolution.py pes.csv gs -6.0 2.0 0.12 1.0 300 1.0e-4 5 0
pes.csv: 入力データファイル。gs: Gauss-Seidel法で実行。-6.0: 解析範囲のX軸最小値。2.0: 解析範囲のX軸最大値。0.12: ガウスフィルターの半値幅を \(0.12\) eV に設定。1.0: ダンピング係数を \(1.0\) に設定。300: 最大反復回数を \(300\) 回に設定。1.0e-4: 収束判定の相対誤差を \(1.0 \times 10^{-4}\) に設定。5: 各反復後に窓幅 \(5\) のSavitzky-Golayフィルターで平滑化。0: 負の値のゼロ補正を無効化。
実行結果の説明:
プログラムはまず、入力データの概要、フィルターパラメータ、Gauss-Seidel法の反復条件を表示します。
その後、matplotlib を用いたグラフウィンドウが開き、リアルタイムで逆畳み込みの反復計算の進捗が表示されます。上段のグラフには、元の生データと現在の逆畳み込み結果が重ねて表示され、下段にはガウスフィルターの形状が表示されます。
各反復ステップごとに、誤差と収束状況がコンソールに表示されます。
収束条件を満たすか、最大反復回数に達すると計算が終了し、最終的なグラフが表示された後、結果が pes-deconvoluted.csv というファイルに保存されます。このCSVファイルには、X軸データ、逆畳み込みへの入力データ、および最終的な逆畳み込み結果が含まれます。
例2: FFT法による逆畳み込みとデータ拡張・平滑化
pes.csv という入力ファイルに対して、FFT法を用いて逆畳み込みを実行し、前処理としてデータ拡張と移動平均による平滑化を適用します。
python deconvolution.py pes.csv fft full average+extend -4.5 2.0 0.12 2.0 5 5
pes.csv: 入力データファイル。fft: FFT法で実行。full:numpy.convolveのモードをfullに設定（FFTでは内部的にデータ長が調整されるため、ここでのfullは主に前処理のconvolve関数に影響）。average+extend: 前処理として、移動平均 (average) とデータ拡張 (extend) を適用。-4.5: 解析範囲のX軸最小値。2.0: 解析範囲のX軸最大値。0.12: ガウスフィルターの半値幅を \(0.12\) eV に設定。2.0: ガウスフィルターを半値幅の \(2.0\) 倍の範囲で考慮。5: データ拡張時に先頭にnw*5個のゼロを追加。5: データ拡張時にnw*5個の線形補間を適用。
実行結果の説明:
プログラムは、入力データの概要、フィルターパラメータ、FFT法固有のデータ長変更（2のべき乗への調整）、およびデータ拡張の詳細を表示します。
前処理として、入力データに対して移動平均による平滑化と、FFT処理での端点効果を軽減するためのデータ拡張が行われます。
その後、拡張・平滑化されたデータとガウスフィルター波形に対してFFTベースの逆畳み込みが実行されます。
計算終了後、matplotlib のグラフウィンドウに3つのプロットが表示されます。上段には入力として使用された（畳み込まれた）データと逆畳み込み結果、中段には元の生データと畳み込み入力データ、下段にはFFT処理に用いられたフィルター波形が表示されます。
最後に、結果は pes-deconvoluted.csv というファイルに保存されます。