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

プログラムの動作

deconvolution.py は、畳み込みによって鈍化（ぼやけ）した測定データから、元の鋭い信号を復元する逆畳み込み処理を行うPythonプログラムです。主に、実験で得られたスペクトルデータや波形データが、測定系の応答関数（フィルター関数）によって広げられた場合に、その応答関数の影響を取り除き、本来の物理的な形状を明らかにする目的で使用されます。

主な機能

逆畳み込み手法の提供:
- Jacobi法: 反復法による逆畳み込み。
- Gauss-Seidel法: Jacobi法を改良した反復法。
- Smoothing Penalty法: 正則化を伴う回帰による逆畳み込み。
- Ridge回帰: 正則化項を追加した線形回帰による逆畳み込み。
- FFT (高速フーリエ変換) による逆畳み込み。
- scipy.signal.deconvolve 関数による逆畳み込み。
畳み込み処理の確認: numpy.convolve を使用して、フィルター関数による畳み込み処理を確認する機能。
フィルター関数の選択: ガウス関数またはローレンツ関数をフィルターカーネルとして選択可能。半値幅 (Wa) と範囲 (Grange) を指定して調整できます。
データ前処理:
- 入力データ範囲の指定。
- ゼロパディングや線形拡張によるデータ範囲の調整（特にFFTや畳み込み時に重要）。
- 単純移動平均または多項式近似による平滑化。
リアルタイム処理状況の可視化: Jacobi法やGauss-Seidel法などの反復計算中に、現在の結果と入力データをグラフでリアルタイムに表示し、収束状況を確認できます。
結果の出力: 逆畳み込みされたデータ、入力データ、および畳み込みフィルター関数をグラフとして表示し、処理結果をExcelファイル（.xlsx）に保存します。

解決する課題

このプログラムは、以下のような課題を解決することを目的としています。

測定機器の分解能限界: 物理的な測定器には常に分解能の限界があり、真の信号が機器の応答関数によって広げられてしまいます。このプログラムは、その広がりを取り除き、より高分解能な情報を引き出すことを試みます。
ノイズへの頑健性: 逆畳み込みはノイズに非常に敏感な処理ですが、Jacobi法やGauss-Seidel法での平滑化処理、Ridge回帰やSmoothing Penalty法での正則化を通じて、ノイズの影響を抑制しつつ安定した解を得ることを目指します。
信号の歪み補正: 複雑な信号が特定のフィルターによって歪められた場合でも、そのフィルター特性が既知であれば、逆畳み込みによって元の信号特性を推定できます。

原理

畳み込みと逆畳み込み

畳み込み (convolution) は、一つの関数がもう一つの関数によって「ぼやける」様子を記述する数学的操作です。測定では、本来の信号 $I(t)$ が測定装置の応答関数 $H(t)$ によって畳み込まれ、測定信号 $Y_{raw}(t)$ が得られます。

\[ Y_{raw}(t) = (I * H)(t) = \int_{-\infty}^{\infty} I(\tau) H(t - \tau) d\tau \]

離散データの場合、これは以下の和で表されます。

\[ Y_{raw}[i] = \sum_{j} I[j] H[i-j] \]

逆畳み込み (deconvolution) は、測定信号 $Y_{raw}$ と応答関数 $H$ が与えられたときに、元の信号 $I$ を推定するプロセスです。

フィルター関数 (`convolve_func`)

プログラムでは、応答関数 $H$ としてガウス関数またはローレンツ関数が使用されます。これらの関数は、Wa (半値幅に相当) と Grange (フィルターの範囲) パラメータで定義されます。

ガウス関数

\[ G(x, x_0, whalf) = A \cdot \exp\left( -\left(\frac{x - x_0}{a}\right)^2 \right) \]

ここで、$a = \frac{whalf}{\sqrt{\ln 2}}$ であり、$A$ は規格化定数です。

ローレンツ関数

\[ L(x, x_0, whalf) = A \cdot \frac{1}{1 + \left(\frac{x - x_0}{whalf}\right)^2} \]

ここで、$A$ は規格化定数です。

Jacobi法とGauss-Seidel法

これらは、線形方程式系 $Y_{raw} = H \cdot I$ を反復的に解くための手法です。プログラムでは、以下の更新式に基づいて $I$ の近似値 $y$ を逐次的に改善します。

\[ y_{new}[i] = y_{old}[i] + \frac{1}{H_{ii}} \left( Y_{raw}[i] - \frac{1}{S_g[i]}\sum_{j} H_{ij} y_{old}[j] \right) \cdot dump \]

ここで、$H_{ij}$ はフィルター関数 Hij(xstep, Wa, Grange, i, j) であり、$S_g[i]$ はフィルター関数 $H_{ij}$ の行方向の総和です。dump は緩和係数であり、収束を調整します。

Jacobi法: $y_{new}[i]$ を計算する際に、右辺の $\sum_{j} H_{ij} y_{old}[j]$ にはすべて前回の反復で得られた $y_{old}$ の値が使用されます。
Gauss-Seidel法: $y_{new}[i]$ を計算する際に、$\sum_{j} H_{ij} y[j]$ の部分で、$j < i$ の場合は既に更新された $y[j]$ の値が、$j \geq i$ の場合は前回の反復で得られた $y_{old}[j]$ の値が使用されます。これにより、Jacobi法よりも速い収束が期待されます。

両手法とも、各反復後に SmoothingByPolynomialFit 関数による平滑化が適用され、ノイズによる発散を防ぎます。また、zero_correction が有効な場合、負の値は0にクリップされます。

Ridge回帰とSmoothing Penalty法

これらは、線形最小二乗回帰に正則化項を導入することで、過学習を防ぎ、逆畳み込みの安定性を高める手法です。

基本的な考え方

信号 $Y_{raw}$ とカーネル $H$ から元の信号 $I$ を推定する際、以下の目的関数を最小化します。

\[ \text{min}_I \left\| Y_{raw} - H \cdot I \right\|^2 \]

ここで、$I$ は元の信号のベクトル、$H$ は畳み込み操作を表す行列です。

Ridge回帰

Ridge回帰では、過学習を防ぐためにL2正則化項が追加されます。目的関数は以下のようになります。

\[ \text{min}_I \left( \left\| Y_{raw} - H \cdot I \right\|^2 + \alpha \left\| I \right\|^2 \right) \]

ここで $\alpha \ge 0$ は正則化パラメータです。プログラムでは、この目的関数の最小化に対応する正規方程式を解きます。具体的には、Sij 行列の対角成分に alpha を加算することで、正則化が行われます。

Smoothing Penalty法

Smoothing Penalty法はRidge回帰と類似していますが、信号の「滑らかさ」を促進する正則化項を追加します。これにより、ノイズによる急激な変動を抑える効果があります。プログラムでは、Sij 行列の非対角成分に alpha を減算する形で平滑化ペナルティを導入しています。

FFTによる逆畳み込み

周波数領域では、畳み込みは単純な乗算になります。この特性を利用して逆畳み込みを行います。

測定信号 $Y_{raw}$ と応答関数 $H$ を高速フーリエ変換 (FFT) します: $$ \mathcal{F}(Y_{raw}) \quad \text{および} \quad \mathcal{F}(H) $$
周波数領域で除算を行います: $$ \mathcal{F}(I) = \frac{\mathcal{F}(Y_{raw})}{\mathcal{F}(H)} $$
結果を逆高速フーリエ変換 (IFFT) して、元の信号 $I$ を得ます: $$ I = \mathcal{F}^{-1}\left( \frac{\mathcal{F}(Y_{raw})}{\mathcal{F}(H)} \right) $$

この手法は計算効率が高いですが、ノイズに対して非常に敏感であり、周波数領域で $\mathcal{F}(H)$ がゼロに近い成分を持つ場合に不安定になることがあります。プログラムでは、FFTのためにデータ長を2のべき乗に拡張しています。

`scipy.signal.deconvolve`

scipy.signal ライブラリの deconvolve 関数は、フーリエ変換を使用しない逆畳み込みを実行します。これは、より一般的な信号処理の状況に対応できる可能性があります。

平滑化関数

プログラムには、以下の平滑化関数が実装されています。

SmoothingBySimpleAverage(y, n): $n$ 点移動平均による平滑化。
SmoothingByPolynomialFit(y, n): $n$ 点の範囲で多項式近似を行い、中心点の値をその近似値で置き換える平滑化。Savitzky-Golayフィルターに類似した挙動を示します。

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

このプログラムは、以下のPythonライブラリに依存しています。

numpy: 数値計算を効率的に行うためのライブラリ。
pandas: データ解析と操作のためのライブラリ。
scipy: 科学技術計算のためのライブラリ（信号処理機能を使用）。
matplotlib: グラフ描画のためのライブラリ。

これらのライブラリは、通常 pip コマンドでインストールできます。

pip install numpy pandas scipy matplotlib

また、プログラムは tklib というカスタムライブラリも使用しています。

tklib: tklib.tkutils, tklib.tksci.tksci, tklib.tkvariousdata, tklib.tkapplication, tklib.tkgraphic.tkplotevent モジュールが含まれています。このライブラリは標準のPythonパッケージインデックスには存在しないため、別途入手してPythonのパスが通っている場所に配置するか、deconvolution.py と同じディレクトリに配置する必要があります。本ドキュメントでは、tklib のインストール方法は提供できません。

必要な入力ファイル

プログラムは、逆畳み込みを行うためのデータを含むCSVファイルを必要とします。

ファイル形式: CSV (Comma Separated Values) 形式。
データ構造: ヘッダー行を含み、X軸データとY軸データがそれぞれ独立した列に記述されている必要があります。
ファイル名: コマンドライン引数で指定します（デフォルトは pes.csv）。
X軸/Y軸の指定: コマンドライン引数で、データファイル内のどの列をX軸、Y軸として使用するかを、ヘッダー名または列インデックス（0から始まる整数）で指定します。

例: pes.csv

Energy(eV),Intensity(a.u.),Background
-10.0,10,0
-9.9,15,0
-9.8,20,0
...
0.0,100,0
...
9.9,15,0
10.0,10,0

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

プログラムの実行後、以下のファイルが生成されます。

結果のExcelファイル:
- ファイル名: 入力ファイル名に -deconvoluted.xlsx が追加された形式で保存されます。例: pes.csv → pes-deconvoluted.xlsx
- 内容: 以下の3つの列を含むデータシートが保存されます。
  - x: 処理されたX軸データ。
  - y(input): 逆畳み込み処理に入力されたY軸データ（前処理後のデータ）。
  - y(deconvoluted): 逆畳み込みによって得られたY軸データ。
- テンプレート: tkVariousData().to_excel メソッドが呼び出されており、内部的に StandardGraph.xlsm というExcelテンプレートファイルを使用する可能性があります。ただし、このテンプレートファイルの具体的な内容はコードからは不明です。
ログファイル:
- ファイル名: 入力ファイル名と同じ名前で、.log 拡張子が追加されたファイル（例: pes.log）。
- 内容: プログラムの標準出力（実行時のパラメータ、途中経過のメッセージ、エラーなど）が記録されます。

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

deconvolution.py は、実行モードに応じて異なる引数を取ります。

Jacobi法、Gauss-Seidel法、Smoothing Penalty法、Ridge回帰モード

usage: python deconvolution.py file mode xmin xmax Wa func_type dump nmaxiter eps nsmooth zeroc xlable ylabel
       mode = [gs|jacobi|sp|ridge]   gs=Gauss-Seidel method   jacobi=Jacobi method   sp=Smoothing penalty method
       zeroc = [0|1]  zero correction after a Jacobi/Gauss-Seidel cycle]

file: 入力CSVファイルのパス。
mode: 実行モード (gs, jacobi, sp, ridge のいずれか)。
xmin: 解析範囲のX軸の最小値。
xmax: 解析範囲のX軸の最大値。
Wa: フィルター関数の半値幅に相当するパラメータ。
func_type: フィルター関数の種類 (gauss または lorentz)。
dump: 緩和係数（Jacobi/Gauss-Seidel法）または正則化パラメータ alpha (Ridge/Smoothing Penalty法)。
nmaxiter: 最大反復回数（Jacobi/Gauss-Seidel法）。
eps: 収束判定の許容誤差（Jacobi/Gauss-Seidel法）。
nsmooth: 各反復後の平滑化回数（Jacobi/Gauss-Seidel法およびSmoothing Penalty/Ridge回帰の平滑化）。
zeroc: ゼロ補正を行うかどうか（0: 行わない, 1: 行う）。
xlable: X軸のデータラベルまたはインデックス。
ylabel: Y軸のデータラベルまたはインデックス。

Convolve、FFTモード

usage: python deconvolution.py file mode convmode smoothmode xmin xmax Wa func_type Grange kzero klin
       mode = [convolve|fft]
       convmode = [''|full|same]
       smoothmode = combination of [none|convolve|extend|average]

file: 入力CSVファイルのパス。
mode: 実行モード (convolve または fft のいずれか)。
convmode: 畳み込みの出力モード (full, same, または空文字列)。
smoothmode: 前処理としての平滑化モード（none, convolve, extend, average の組み合わせ。例: convolve+extend）。
xmin: 解析範囲のX軸の最小値。
xmax: 解析範囲のX軸の最大値。
Wa: フィルター関数の半値幅に相当するパラメータ。
func_type: フィルター関数の種類 (gauss または lorentz)。
Grange: フィルター関数の範囲を Wa の倍数で指定。
kzero: ゼロパディングの幅をフィルター幅の倍数で指定。
klin: 線形拡張の幅をフィルター幅の倍数で指定。

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

例1: Gauss-Seidel法による逆畳み込み

この例では、pes.csv ファイルに対してGauss-Seidel法を用いて逆畳み込みを行います。

python deconvolution.py pes.csv gs -6.0 2.0 0.12 gauss 1.0 300 1.0e-4 5 0 0 1

引数の説明:

pes.csv: 入力ファイル。
gs: Gauss-Seidel法で実行。
-6.0: X軸の最小値。
2.0: X軸の最大値。
0.12: ガウスフィルターの半値幅 $W_a$。
gauss: ガウス関数をフィルターとして使用。
1.0: 緩和係数 dump。
300: 最大反復回数 nmaxiter。
1.0e-4: 収束判定の許容誤差 eps。
5: 各反復後の平滑化点数 nsmooth。
0: ゼロ補正は行わない。
0: X軸は0列目のデータを使用。
1: Y軸は1列目のデータを使用。

実行結果: プログラムが起動し、リアルタイムで3つのグラフが表示されます。上段のグラフでは、元の入力データ（raw/initial）と、Gauss-Seidel法によって繰り返し更新されていく逆畳み込み結果（updated）がプロットされます。反復が進行するにつれて、updated プロットは徐々に収束していきます。中段には元の入力データ、下段には使用されたフィルター関数が表示されます。指定された最大反復回数 (nmaxiter) に達するか、または相対誤差が eps より小さくなった場合に計算が停止します。最終的な逆畳み込み結果と元のデータ、および処理に使用されたフィルター関数がグラフに表示され、pes-deconvoluted.xlsx というファイルに x, y(input), y(deconvoluted) のデータが保存されます。

例2: FFTによる逆畳み込みとデータ前処理

この例では、SnSe-DOS.csv ファイルに対してFFT逆畳み込みを行います。畳み込みモードは full、平滑化モードは convolve+extend を指定しています。

python deconvolution.py SnSe-DOS.csv fft full convolve+extend -4.5 2.0 0.12 gauss 2.0 5 5