如何从巨蟒的上升尾谱中减去基线？

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)

Original_spectrum

​

​

在高x值附近有一个清晰的尾。假设尾巴是一个工件，需要移除。我尝试过使用P. Eilers、rubberband方法和桃子包使用的解决方案，但是这些解决方案最终会减去尾部，在低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方法：

# 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)

Rubber_plot

​

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算法，这里。其中一种算法是非对称重加权惩罚最小二乘平滑算法(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算法是我能得到的最接近的算法，至少对于像这样的光谱是这样的。我们将看到该算法的鲁棒性如何适合不同形状的光谱。当然，我总是乐于接受建议或改进！

Answer 2

看看python中的RamPy库，它提出了各种基线减法算法。这包括样条，ARPLS，ALS，多项式函数，等等。它还提供了各种其他功能，如重采样、规范化和峰值拟合的例子。

在你的情况下，一个简单的样条函数拟合前后的峰值应该很容易做的工作。看看这个例子木星笔记本。