python - 如何从 python 中尾部上升的频谱中减去基线？

Question

我有一个要从中减去基线的光谱。光谱数据为:

1.484043000000000001e+00    1.121043091000000004e+03
1.472555999999999976e+00    1.140899658000000045e+03
1.461239999999999872e+00    1.135047851999999921e+03
1.450093000000000076e+00    1.153286499000000049e+03
1.439112000000000169e+00    1.158624877999999853e+03
1.428292000000000117e+00    1.249718872000000147e+03
1.417629999999999946e+00    1.491854857999999922e+03
1.407121999999999984e+00    2.524922362999999677e+03
1.396767000000000092e+00    4.102439940999999635e+03
1.386559000000000097e+00    4.013319579999999860e+03
1.376497999999999999e+00    3.128252441000000090e+03
1.366578000000000070e+00    2.633181152000000111e+03
1.356797999999999949e+00    2.340077147999999852e+03
1.347154999999999880e+00    2.099404540999999881e+03
1.337645999999999891e+00    2.012083983999999873e+03
1.328268000000000004e+00    2.052154540999999881e+03
1.319018999999999942e+00    2.061067871000000196e+03
1.309895999999999949e+00    2.205770507999999609e+03
1.300896999999999970e+00    2.199266602000000148e+03
1.292019000000000029e+00    2.317792235999999775e+03
1.283260000000000067e+00    2.357031494000000293e+03
1.274618000000000029e+00    2.434981689000000188e+03
1.266089999999999938e+00    2.540746337999999923e+03
1.257675000000000098e+00    2.605709472999999889e+03
1.249370000000000092e+00    2.667244141000000127e+03
1.241172999999999860e+00    2.800522704999999860e+03

我只从实际数据文件中提取了每 20 个数据点，但保留了一般形状。

import matplotlib.pyplot as plt
share = the_above_array
plt.plot(share)

原始光谱

在高 x 值周围有明显的尾部。假设尾部是人工制品，需要去除。我尝试过使用 ALS algorithm 的解决方案作者:P. Eilers，一个 rubberband approach , 和 peakutils包，但这些最终会减去尾部并在低 x 值周围产生上升或不产生合适的基线。

ALS 算法，在这个例子中我使用了 lam=1E6 和 p=0.001；这些是我能够手动找到的最佳参数:

# ALS approach
from scipy import sparse
from scipy.sparse.linalg import spsolve
def baseline_als(y, lam, p, niter=10):
    L = len(y)
    D = sparse.csc_matrix(np.diff(np.eye(L), 2))
    w = np.ones(L)
    for i in range(niter):
      W = sparse.spdiags(w, 0, L, L)
      Z = W + lam * D.dot(D.transpose())
      z = spsolve(Z, w*y)
      w = p * (y > z) + (1-p) * (y < z)
    return z

baseline = baseline_als(share[:,1], 1E6, 0.001)
baseline_subtracted = share[:,1] - baseline
plt.plot(baseline_subtracted)

ALS_plot

橡皮筋方法:

# rubberband approach
from scipy.spatial import ConvexHull
def rubberband(x, y):
    # Find the convex hull
    v = ConvexHull(share).vertices
    # Rotate convex hull vertices until they start from the lowest one
    v = np.roll(v, v.argmax())
    # Leave only the ascending part
    v = v[:v.argmax()]
    # Create baseline using linear interpolation between vertices
    return np.interp(x, x[v], y[v])

baseline_rubber = rubberband(share[:,0], share[:,1])
intensity_rubber = share[:,1] - baseline_rubber
plt.plot(intensity_rubber)

橡胶图

peakutils 包:

# peakutils approach
import peakutils
baseline_peakutils = peakutils.baseline(share[:,1])
intensity_peakutils = share[:,1] - baseline_peakutils
plt.plot(intensity_peakutils)

Peakutils_plot

除了屏蔽低 x 值数据之外，是否有任何建议可以构建基线并减去尾部而不增加低 x 值？

Answer 1

我找到了一组类似的ALS算法here .其中一种算法，非对称重新加权惩罚最小二乘平滑 (arpls)，比 als 的拟合效果稍好。

# arpls approach
from scipy.linalg import cholesky
def arpls(y, lam=1e4, ratio=0.05, itermax=100):
    r"""
    Baseline correction using asymmetrically
    reweighted penalized least squares smoothing
    Sung-June Baek, Aaron Park, Young-Jin Ahna and Jaebum Choo,
    Analyst, 2015, 140, 250 (2015)
    """
    N = len(y)
    D = sparse.eye(N, format='csc')
    D = D[1:] - D[:-1]  # numpy.diff( ,2) does not work with sparse matrix. This is a workaround.
    D = D[1:] - D[:-1]
    H = lam * D.T * D
    w = np.ones(N)
    for i in range(itermax):
        W = sparse.diags(w, 0, shape=(N, N))
        WH = sparse.csc_matrix(W + H)
        C = sparse.csc_matrix(cholesky(WH.todense()))
        z = spsolve(C, spsolve(C.T, w * y))
        d = y - z
        dn = d[d < 0]
        m = np.mean(dn)
        s = np.std(dn)
        wt = 1. / (1 + np.exp(2 * (d - (2 * s - m)) / s))
        if np.linalg.norm(w - wt) / np.linalg.norm(w) < ratio:
            break
        w = wt
    return z

baseline = baseline_als(share[:,1], 1E6, 0.001)
baseline_subtracted = share[:,1] - baseline
plt.plot(baseline_subtracted, 'r', label='als')

baseline_arpls = arpls(share[:,1], 1e5, 0.1)
intensity_arpls = share[:,1] - baseline_arpls
plt.plot(intensity_arpls, label='arpls')

plt.legend()

ARPLS 图

幸运的是，当使用来自整个频谱的数据时，这种改进变得更好:

请注意，两种算法的参数都不同。目前，我认为 arpls 算法是我所能得到的最接近的算法，至少对于看起来像这样的光谱来说是这样。我们将看到该算法如何稳健地拟合具有不同形状的光谱。当然，我总是乐于接受建议或改进!