[Python] 関数フィッティングなしで半値幅を求める

2024年2月24日

この記事で分かること

pythonで「フィッティングではなく数値関係のみで半値幅を出したい人」に向けた記事です。
不定形のピークを評価するのに利用できます

前書き

ピーク形状を評価したいものの、フィットすべき関数形がわからないとき、とりあえず半値幅を評価したいときがあります。

最大と最小強度を用いて機械的に半値幅を計算する方法を紹介します。

コード

# %%
import numpy as np
import scipy
import pandas as pd
import matplotlib.pyplot as plt

def gaussian(x, amp, cen, wid):
    return (amp / (np.sqrt(2*np.pi) * wid)) * np.exp(-(x-cen)**2 / (2*wid**2))

x = np.linspace(0, 10.0, 201)
y = gaussian(x, 20, 5.0, 1.0) +np.random.normal(scale=0.1, size=x.size)
sr = pd.Series(y, index = x)

#%%
x0_max = sr.idxmax()
y0_max = sr.max()
y0_min = sr.min()
#%%
x1_min = sr.loc[:sr.idxmax()].idxmin()
y1_min = sr.loc[:sr.idxmax()].min()

x2_min = sr.loc[sr.idxmax():].idxmin()
y2_min = sr.loc[sr.idxmax():].min()

# %%
y1_half = (y0_max - y1_min)/2
y2_half = (y0_max - y2_min)/2

# %%
f1 = scipy.interpolate.interp1d(x,y-y1_half) 
f2 = scipy.interpolate.interp1d(x,y-y2_half) 

# %%
x1 =scipy.optimize.brentq(f1, x.min(),x0_max)
x2 =scipy.optimize.brentq(f2, x0_max,x.max())
fwhm = x2-x1
print(fwhm)

# %%
fig,ax = plt.subplots()
ax.plot(x, y)
ax.scatter(x0_max, y0_max, color ="r")
ax.scatter(x1_min, y1_min, color ="r")
ax.scatter(x2_min, y2_min, color ="r")
ax.scatter(x1, y1_half, color ="b")
ax.scatter(x2, y2_half, color ="b")
ax.plot([x1,x1], [y0_min, y0_max], linestyle="--",color ="b")
ax.plot([x2,x2], [y0_min, y0_max], linestyle="--", color ="b")
ax.plot([x1,x2], [f1(x1),f2(x2)], "-x", color ="r")
plt.show()

プロファイルを最大値で左側と右側に分け、それぞれで最小値を求め、半値となる座標x1, x2を求めます。

半値幅はx2-x1として機械的に求まります。