今回は弊社DSPやMCAに搭載されているスペクトル解析の機能の一つであるpeak search analysisのガウスフィッテングについて解説します。

　テクノエーピーのpeak search analysisに搭載されているガウスフィッテングは非線形最小二乗法を用いて適合関数をデータへの当てはめを行う手法ですが、アルゴリズムはLevenberg-Marquardt 法を利用しています。高速な収束速度と安定性を兼ね備えており、十分な精度で最適値を求めることができます。パラメータの初期値が必要ですが、ピークサーチから得られた情報より自動で算出するようになっており、ユーザーが入力する必要はありません。簡単ではありますがアルゴリズムを解説します。

　まず、スペクトルデータy_iとモデル関数f(x_i;a)ガウス関数＋1次式があります。

　Eq.(1)の係数（パラメータ）をaと統一します。

　１．モデル関数f(x_i;a)の独立変数x_iを与え、パラメータの初期値a⁽⁰⁾をセットします。
　２．データy_iとモデルf(x_i;a)の適合具合を評価するχ²（重み付き残差二乗の和）
　　　Eq.(3)を定義して、χ²(a⁽⁰⁾)を計算します。この値をχ₀²とします。

　３．f(x_i;a)を各パラメータで偏微分したm×mの対称行列α_jkを作ります。

　４．残差に偏微分をかけたm 行 1 列の行列β_jを作ります。

　５．対角成分α_jjに強調させる係数(1+λ)を乗じたのがα'です。

　６．連立方程式Eq.(7)を解いてδaを求め、χ2(a⁽⁰⁾+δa)を計算します。

　７．χ²>χ₀²ならば、λを10倍にして、５．に戻ってやり直します。
　８．χ²-χ₀²<εχ²ならば計算終了です。εは収束条件です。
　９．χ²<χ₀²ならば、a_j=a_j⁰,χ₀²=χ²とし、λを1/10にして、３．に戻ります。

　Levenberg-Marquardt 法はフィッテングの最初の段階（a_j⁽⁰⁾が真の階に近くない時）にはλを大きく採り、対角成分を強調してχ²の勾配を急降下させます。この場合、α_jjが大きいので変化量δαは少ないです。χ²が順調に真の解に近づくにつれλ→0と減っている場合は、Newton法に近づけて大胆に変化量δαを動かして一気に収束させます。

　(Fig.1 ガウスフィッテングコード一部)

　ガウスフィッテングなどの非線形最小二乗法の演算ライブラリはpython、rootなど数値計算が得意なインタプリタ言語には無償で準備されているようなのでそれを使えばよいと思ったのですが、弊社のようにアプリケーションを開発している場合、DLLとして組み込む必要があります。peakfit、MATLABで有償で対応しているようです。改造や調整もできないことを考えると自作した方がよいということで、開発することにしました。ピークサーチと同様にC言語で記述してDLL化（Fig.1）することで高速化を実現、リアルタイム演算を可能にしております。

(Fig.2 統計が少ないスペクトル)

　先ずガウスフィッテングの良いところは、統計が足りてなくて手計算だと全然まともな値にならない場合（Fig.3）でも、それなりの数値（Fig.2）を出してくれる所です。再現性が高いので長時間測定しなくても安定した数値を得ることができます。

(Fig.3 main画面の手計算)

　次にフルスケールを10MeVとか非常に広く取ったとき(Fig.4)で、低エネルギー（この場合は121.8keV）を着目した場合、ヒストグラムのチャネルを16kchと細かく取ったとしてもチャネルの分解能が足りなくなる場合があります。(Fig5) net(cps)rawは黒点からネット計数率を算出した値であり、ガウスフィッテングで算出した値net(cps)fitより小さい数字になってしまっています。理由はチャネルの分解能が足りないからでしょう。このような場合でもガウスフィッテングは有効です。

(Fig.4 10MeVフルスケールγ線スペクトル)


(Fig.5 チャネルの分解能が足りないスペクトル)

　アプリケーション上ではピークサーチから得られた結果より自動的に初期値を求めガウスフィッテングを行いますので、これと言った使い方はないのですが、リアルタイムに情報が更新しますので便利かと思います。

(Fig.6 2.6MeVのピーク)

　Fig.6は流通しているTIG溶接用のトリウム入りタングステン電極棒から発生している自然放射線トリウム系列からのTl-208の2.6MeVγ線のピークです。高エネルギーガンマ線の試験をするときに便利です。

　Fig.7はフルスケール3MeVで35分測定したスペクトルです。青枠に注目して下さい。

（Fig.7 3MeVフルスケールγ線スペクトル）

　Fig.8は青枠を拡大しているスペクトルです。1分間計測しました。443keVと463keVはまだ埋もれている感じですが、

（Fig.8 埋もれ気味のピーク）

　きちんとガウスフィッテングがかかってくれます。（Fig.9）

　
(Fig.9 フィットがかかったスペクトル1min計測)

　1minスペクトル(Fig.9)と10minスペクトル(Fig.10)で3つのピークについて比較します。net(cps)fit、FWHM(keV)fitは誤差の範囲内だと思います。

　（Fig.10 フィットがかかったスペクトル10min計測)）

　それでは、より良い製品が作れるように社員一同全力で頑張ります。応援の程どうぞよろしくお願いします。

参考文献
[1] Willaim H. Press, NUMERICAL RECIPES in C, 技術評論社, 1993年