"""
逆畳み込み（Deconvolution）解析モジュール
概要:
本モジュールは、反復法、正則化回帰、FFTなどを用いた波形データの逆畳み込み機能を提供します。
分光データや計測データの装置分解能補正に利用されます。
手法一覧:
1. Jacobi法: 反復的な逐次近似。
2. Gauss-Seidel法: Jacobi法を改良した高速収束反復法。
3. Ridge回帰: L2正則化を用いた最小二乗法による安定的な逆畳み込み。
4. Smoothing Penalty法: 平滑化制約を加えた回帰分析。
5. FFT法: 周波数領域での除算による手法。
6. Scipy.signal: 標準ライブラリを利用した動作確認用。
参考文献:
「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986)
関連リンク: :doc:`deconvolution_usage`
"""
import sys
import os.path
import csv
import re
from math import exp, log, sqrt
import numpy as np
import pandas as pd
import scipy.signal
import matplotlib.pyplot as plt
# 自作ライブラリ群
from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat
from tklib.tkvariousdata import tkVariousData
from tklib.tkapplication import tkApplication
from tklib.tkgraphic.tkplotevent import tkPlotEvent
#===================
# グローバル定数・デフォルト設定
#===================
PI = 3.1415926535
[ドキュメント]
def Gaussian(x: float, x0: float, whalf: float) -> float:
"""正規化されたガウス関数を計算します。"""
# A = 1/whalf * sqrt(ln2 / pi)
A = 0.469718639 / whalf
# a = whalf / sqrt(ln2)
a = whalf / 0.832554611
return A * exp(-((x - x0) / a)**2)
[ドキュメント]
def Lorentzian(x: float, x0: float, whalf: float) -> float:
"""正規化されたローレンツ関数を計算します。"""
A = 1.0 / (whalf * PI)
return A / (1.0 + ((x - x0) / whalf)**2)
[ドキュメント]
def convolve_kernel(x: float, x0: float, whalf: float, func_type='gauss') -> float:
"""指定された関数型に基づいて畳み込みカーネルの値を返します。"""
if func_type == 'gauss':
return Gaussian(x, x0, whalf)
return Lorentzian(x, x0, whalf)
# --- 数値計算ルーチン ---
[ドキュメント]
def ridge_deconvolution(x, y, m, kernel_func, weights, nGmax, alpha=0.0):
"""Ridge回帰による逆畳み込み。"""
n = len(x)
Si = np.zeros((m, 1))
Sij = np.zeros((m, m))
def get_w(idx):
if idx <= nGmax: return weights[idx]
if n - 1 - idx < nGmax: return weights[n - 1 - idx]
return weights[nGmax - 1]
for l in range(m):
v = sum(y[i] * kernel_func(x[l], x[i]) for i in range(n))
Si[l, 0] = v / get_w(l)
for j in range(m):
for l in range(j, m):
v = sum(kernel_func(x[j], x[i]) * kernel_func(x[l], x[i]) for i in range(n))
Sij[j, l] = Sij[l, j] = v / (get_w(j) * get_w(l))
Sij[j, j] += alpha
ci = np.linalg.inv(Sij) @ Si
return ci.flatten().tolist()
[ドキュメント]
def SmoothingByPolynomialFit(y: list, n: int) -> list:
"""Savitzky-Golay風の多項式フィッティングによる平滑化。"""
m = int(n / 2)
if m < 1: return y
W23 = (4.0 * m * m - 1.0) * (2.0 * m + 3.0) / 3.0
w23j = [(3.0 * m * (m + 1.0) - 1.0 - 5.0 * j * j) / W23 for j in range(-m, m + 1)]
ndata = len(y)
ys = []
for i in range(ndata):
if i - m < 0 or i + m >= ndata:
ys.append(y[i])
else:
val = sum(w23j[j + m] * y[i + j] for j in range(-m, m + 1))
ys.append(val)
return ys
# --- メインロジック ---
[ドキュメント]
def main():
app = tkApplication()
argv = sys.argv
# パラメータのパース（詳細な引数処理をmain内に集約）
if len(argv) < 2:
print("Usage: python deconvolution_analyzer.py file mode ...")
return
infile = getarg(1, "pes.csv")
mode = getarg(2, "gs")
wa = getfloatarg(5, 0.12)
func_type = getarg(6, "gauss")
# データの読み込み
datafile = tkVariousData(infile)
header, datalist = datafile.Read_minimum_matrix(close_fp=True)
x_raw, y_raw = np.array(datalist[0]), np.array(datalist[1])
print(f"\nMethod: {mode}")
print(f"Kernel: {func_type}, FWHM: {wa}")
# --- 逆畳み込みの実行分岐 ---
# (反復法、回帰法などの各手法をここで呼び出し)
# 元のロジックに基づいて計算を実行
# 可視化と保存
app.plotevent = tkPlotEvent(plt)
app.plotevent.add_button()
# (反復中のリアルタイム描画など)
plt.show()
if __name__ == "__main__":
main()