色々と科学データを扱うとき、バックグラウンド強度の見積もりをすることは多い。たいていの場合は付属のソフトがついていて、バックグラウンドの計算をして減算してくれる機能がついている。
しかしそんな場合でも勝手に自動でやってくれるならいいのだけど、意外と１個１個マニュアルでやらなければならないこともある。１０個くらいのデータならマニュアルでやっても良いのだけど、１００個から１０００個とかになってくるとプログラムを走らせてまとめて見積もりをさせたくなる。
まとめてファイルを読み込む方法はいろいろあるけれど、下記の投稿などを参照していただきたい。
内部リンク：pythonのtkinterを使ってユーザー入力を取得する方法
ということでこの投稿ではまとめて読み込んだデータを、まとめて処理するのに便利なバックグラウンド評価や減算方法についていくつか書いて行きたい。私は主にX線散乱に使っているけれど、ラマンなど他のスペクトルデータでもなんでも良い。
数式を使ったバックグラウンド計算
以下データはnumpyのndarrayに入った散布データということで。０列にX軸の値、そして１列に対応する強度データなどが入っているとする。下の例ではpwdという名前になっているが、なんでこんな名前をつけたかは忘れた。ndarrayへのデータの収納などに興味があれば下記投稿などで。
内部リンク：pythonでのテキストデータの入出力方法いろいろ
さてどうバックグラウンドを処理するかということだけど、一番都合が良いのはバックグラウンドが理由があるタイプのもので数式に従うはずのケース。理論式をそのまま使っても良いし、フィッティングをしても良い。
さてそれをどうpythonで引くかってことだけど、普通に式を作ってデータポイントごとに引くだけで良い。
まずは意外とよく使うコンスタントな最小値をバックグラウンドとして引く場合。これは最小値をnumpyのnanminで見つけて引くだけで良い。範囲指定（この例の場合XからY）の最小値や平均値を拾ってもいいし、どこかの定点の値を恣意的に引っ張ってしまっても大した問題はないだろう。
#XからYの最小値を探して各データポイントから引く const = np.nanmin(pwd[X:Y,1]) pwd[:,1] = pwd[:,1]-const
続いて以下はデータをPower Lawでフィッティングをしてから引く場合の例。この例ではpower lawだけど、gaussianでもpolynomialでも好きな関数を使ってフィットすれば良い。フィッティングにはlmfitを使っている[1]。
内部リンク：pythonのlmfitを使ったカーブフィッティング
from lmfit.models import PowerLawModel
amp, exp = fit_PL(pwd)
def fit_PL(pwd):
#データファイルの４４から７０と５７４から５８４のエリアを使ってフィットしたかったのでnumpyのconcatenateで結合。astypeで変換してたのはデータファイル読み込みの都合だったと思う。たいていの場合なくて良いはず。
y = np.concatenate((pwd[44:70,1],pwd[574:584,1]),axis=0)
y = y.astype(float)
x = np.concatenate((pwd[44:70,0],pwd[574:584,0]),axis=0)
x = x.astype(float)
#lmfitのデフォルトのPowerLawModelを使ってフィット
PL_mod = PowerLawModel(prefix ='PL_')
pars = PL_mod.make_params()
pars['PL_amplitude'].set(10000)
pars['PL_exponent'].set(-1.2, min=-4, max = -1)
res = PL_mod.fit(y, pars, x=x)
print (res.fit_report())
#Power Lawのパラメータを取得
amp = res.best_values.get('PL_amplitude')
exp = res.best_values.get('PL_exponent')
#プロットするならmatplotlibなどで
plt.plot(pwd[:,0],pwd[:,1])
plt.plot(pwd[:,0],amp*(pwd[:,0]**exp))
#Power Lawパラメータをリターン
return amp, exp
この例の場合だと、データの最初と最後を使ってpower lawにフィットしたかったので、numpyのconcatenateを使ってそのようなデータセットを作ってからフィットした。
そしてamplitudeとexponentの二つのパラメータを関数から返す。
pwd[:,1] = pwd[:,1] - amp*(pwd[:,0]**ex)
この二つのパラメータを使って各データポイントの値を計算し、引き算をすればPower Lawの引き算がおしまい。こんな感じで式を使う場合は、簡単に引くことができる。
バックグラウンドがシンプルな場合のベースライン補正的なバックグラウンド演算
上記のような計算できるバックグラウンドならば良いのだけど、計算式によらないいろいろ混ざったバックグラウンドが引きたいってことも結構多い。以下はそんな時に使えるバックグラウンドの見積もり方について。こういった手法って結構いろいろ転がっていて、pythonだと割と簡単に導入できることが多い。
まずはバックグラウンドがどんどん下がって行くや、どんどんと上がっていくのようなシンプルなケースに使えるバックグラウンドの引き方。アイディアは少し古い１９９７年のクレイのX線研究の論文から[2]。たぶんオープンアクセス。
私がやったのは論文の手法と全く同じではないのかな。少しわかりにくかった場所があったのでそこらは適当に書いたので。結果はほぼ一緒になっていると思うけど。
大雑把にいうと、ある地点から直線を引く。直線の傾きを少しずつ変えていってデータのカーブと接したところでバックグラウンドラインとする。でその接点からまた直線でのバックグラウンド評価をやり直すということを繰り返す。基本的にはピーク以外のところがフラットになる操作。
装置的要因や、アモルファス的なバックグラウンドなどを考えて引くものではない。どんな実験したんだっていうような安定感のないデータ群を、とりあえずフラットなベースラインに揃えたい場合に使う感じ。
アモルファスなどの非ピーク情報を得たい場合にはあまりよろしくない手法。根元をがっつり削っちゃうので。
あとはスペクトルデータのベースライン補正などにも使えるかな。
#まずはスタート地点（init_ang）を指定する。conv_lineは値を行番号に変換する定義関数。意外と便利というかよく使う。
init_ang = 9.5
init_ang = conv_line(pwd,init_ang)
lipbkg = fit_lip(pwd,init_ang)
#バックグラウンド演算の関数fit_lip()。lipbkgにバックグラウンドのデータを詰めていく。
def fit_lip(pwd,init_ang):
pwd = pwd.astype(float)
plt.plot(pwd[:,0],pwd[:,1])
step = 1
init_loop = init_ang-1
lipbkg = pwd[:,1]
#バックグラウンドを決める直線を引き直すためのループ。回数はデータ次第だけどある程度大きければ良いはず、適当。
for j in range(100):
init_loop = init_loop+step+1
#直線の開始地点がデータポイントの行数を超えたらループ終了。この場合５８２だったけどもちろんデータ次第。
if init_loop > 582:
print ('break')
break
#バックグラウンドの直線とデータのカーブが接しているかチェックするループ。直線は単純に一次関数で定義して、傾きを１０ずつ変えている。データ次第で刻みは適宜調整する必要がある。
for i in range(40):
slope = -(i*10)
sec = pwd[init_loop,1] - slope*pwd[init_loop,0]
linear = slope*pwd[init_loop:,0]+sec
#plt.plot(pwd[init_loop:,0],linear,alpha=0.2)
diff = pwd[init_loop:,1]-linear
step = 1
#生データと直線の差分の最小値がプラスになったところで直線とカーブが接していると判断。argpartitionで何個目のデータポイントで接しているのかIndexを抜き出して、バックグラウンドの強度を得るのと次のループでのスタート地点を決めるのに利用。
if np.nanmin(diff[1:]) > 0:
step = np.argpartition(diff[1:],0)[0]
print('i,j',i,j)
lipbkg[init_loop:init_loop+step+2] = linear[0:step+2]
break
plt.plot(pwd[:,0],lipbkg[:])
plt.show()
return lipbkg
#データの値と行番号を変換する関数。差分の一番小さくなるIndexをargpartitionで抜き出すだけ。よく使うので比較的速度の速いargpartitionで。
def conv_line(plot,num):
num = float(num)
dif = np.absolute(plot-num)
mini = np.argpartition(dif,0,axis=0)
return mini[0][0]
というわけで、ある点から引いた直線がカーブと接するたびにその地点からバックグラウンド評価をやり直すことで、直線的なベースラインを得る手法でした。
評価範囲を限定しながら、直線の角度をもっとフレクシブルにすればデコボコしたバックグラウンドにも使えるだろう。しかしこのプログラムを書いた時は一様に減少するバックグラウンドだったのでこれだけで大丈夫だった。forループたくさんだし早くもないけれど、困るほど遅くもなかったと思う。
ピークvs滑らかなアモルファスカーブのスムージングを利用した分離法
こちらはシャープなピークとそうじゃないものをざっくばらんに分けたいという場合に有用。元々はX線のバックグラウンドで開発された手法だけれども[3]、スペクトロスコピーのベースラインの平滑化とかにも使える。
必ずしも正確にバックグラウンドを表現できる手法とはいえないかもしないが、いくつか利点がある。
比較的早い。プログラムの途中に組み込みやすい。評価手法次第だがある程度の定量性がある。演算の回数やスムージングの条件次第でブロードピークの除去から、スペクトルのベースラインを平滑化したい時まで利用できる。比較的複雑なバックグラウンドにも対応できる。などなど。
アイディアとしてはピークを含んだデータに対してスムージングを行うと、ピークは潰れるけどピークじゃないところは大体そのままになる。そしてスムージング後のデータとオリジナルデータを各データポイントごとに比べて値の小さな方を残す。これをどんどんと繰り返すことで、ピークが速やかに潰れていき最終的にバックグラウンドが残ることになる。弱点は大体どんな評価法でもそうだけどピークのオーバーラップ。
というわけで以下はスクリプト。オリジナルの論文とは違う論文のものを参考に、savgol_filterを利用したもの[4]。
#下に書いたbkg_estimate関数を使って最初にバックグラウンドを潰す。
bkg_init = bkg_estimate(Imin, Iave, pwd)
#bkg_iteration関数を一度回して１回目のバックグラウンドを得る。
BackG = bkg_iteration(bkg_init,pwd,lin1,lin2,win)
#後は好きなだけスムージングサイクルを回せば良い。
last = 599
for j in range(1,last+1):
BackG = bkg_iteration(BackG,pwd,lin1,lin2,win)
#バックグラウンドの初期値を決める関数。のちのバックグラウンド減算の回数を減らせる効果があるのだけど、今のコンピュータだと誤差範囲。オリジナルの論文が書かれた２０００年には効果があったのだろう。
def bkg_estimate(arr1, arr2, arr3):
Imin = arr1
Iave = arr2
Int = copy.deepcopy(arr3[:,1])
plt.plot(arr3[:,0],Int)
for i in range(len(arr3[:,1])):
if arr3[i,1] < (Iave +2*(Iave-Imin)):
Int[i] = arr3[i,1]
else:
Int[i] = Iave +2*(Iave-Imin)
plt.plot(arr3[:,0],Int)
return Int
#ピークをスムージングでつぶしながら、バックグラウンドを計算する関数。
#引数はarr2は諸事情でプロット用に入れていたので、普通はなくても良い。それからlin1,lin2もデータ範囲の指定に使ってるだけなのでなくても良い。winは結構この引き方に致命的に効くパラメータなので、メインの関数から書き直せるようにしているので引数になっている。もちろん直接関数の中に書けば省略できる。
def bkg_iteration(arr1,arr2,lin1,lin2,win):
bkg = arr1
area = 0
#scipyのsavgol_filterを使ってピークをスムージングで潰す。
Nbkg = savgol_filter(bkg,win+1,1,mode='nearest')
#あるデータ範囲でスムージング前後の小さい方の値をバックグランドとして残すことでピークがどんどん潰れて行く。
for i in range(lin1,lin2):
if (bkg[i]-Nbkg[i])<0:
bkg[i] = bkg[i]
else:
bkg[i] = Nbkg[i]
plt.plot(arr2[:,0], bkg)
return bkg
状況によっては千回とか回すので、forループの中でいちいち保存するよりも自分自身（BackG）を書き換えた方が良いのかなとプログラムを書いたのだった。
しかしこういう書き方は絵面があまりよろしくない。どういう風に書くと美しいのだろうか？なんかPython的な書き方もありそうなものだけど。
いずれにしろこの手法の一番の問題点はいつスムージングを終わらせるかというところ。
特にバックグラウンドを使うわけでなければ目で見て好きなところで止めれば良いけれど、バックグラウンドを使いたい場合はある程度定量性のある止め方を考えないといけない。しかしこの投稿の本題ではないので深くは突っ込まない。
参照文献
[１] lmfit, https://dx.doi.org/10.5281/zenodo.11813
[２] Clays and Clay Minerals, 45, 132-146, 1997.
[３] Journal of applied chemistry. 33, 977-979, 2000.
[４] Carbohydrate polymers. 78, 543-548, 2009
関連記事
D