python - 光谱数据的基线校正

Question

我正在研究拉曼光谱，它通常有一个基线与我感兴趣的实际信息叠加。因此，我想估计基线贡献。为此，我实现了 this question 中的解决方案。

我确实喜欢那里描述的解决方案，并且给出的代码在我的数据上运行良好。计算数据的典型结果如下所示，其中红色和橙色线是基线估计:Typical result of baseline estimation with calculated data

问题是:我经常在 pandas DataFrame 中收集数千个光谱，每一行代表一个光谱。我当前的解决方案是使用 for 循环一次迭代一个频谱的数据。然而，这使得该过程相当缓慢。由于我对 python 相当陌生，而且由于 numpy/pandas/scipy 的帮助，我已经习惯了几乎不必使用 for 循环，因此我正在寻找一种解决方案，使得也可以省略这个 for 循环。然而，使用的稀疏矩阵函数似乎仅限于二维，但我可能需要三个，而且我还无法想到另一个解决方案。有人有想法吗？

当前代码如下所示:

import numpy as np
import pandas as pd
from scipy.signal import gaussian
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse.linalg import spsolve

def baseline_correction(raman_spectra,lam,p,niter=10):
    #according to "Asymmetric Least Squares Smoothing" by P. Eilers and H. Boelens
    number_of_spectra = raman_spectra.index.size
    baseline_data = pd.DataFrame(np.zeros((len(raman_spectra.index),len(raman_spectra.columns))),columns=raman_spectra.columns)

    for ii in np.arange(number_of_spectra):
        curr_dataset = raman_spectra.iloc[ii,:]

        #this is the code for the fitting procedure        
        L = len(curr_dataset)
        w = np.ones(L)
        D = sparse.diags([1,-2,1],[0,-1,-2], shape=(L,L-2))

        for jj in range(int(niter)):
            W = sparse.spdiags(w,0,L,L)
            Z = W + lam * D.dot(D.transpose())
            z = spsolve(Z,w*curr_dataset.astype(np.float64))
            w = p * (curr_dataset > z) + (1-p) * (curr_dataset < z)
        #end of fitting procedure

        baseline_data.iloc[ii,:] = z
    return baseline_data

#the following four lines calculate two sample spectra
wavenumbers = np.linspace(500,2000,100)
intensities1 = 500*gaussian(100,2) + 0.0002*wavenumbers**2
intensities2 = 100*gaussian(100,5) + 0.0001*wavenumbers**2
raman_spectra = pd.DataFrame((intensities1,intensities2),columns=wavenumbers)
#end of smaple spectra calculataion

baseline_data = baseline_correction(raman_spectra,200,0.01)

#the rest is just for plotting the data
plt.figure(1)
plt.plot(wavenumbers,raman_spectra.iloc[0])
plt.plot(wavenumbers,baseline_data.iloc[0])
plt.plot(wavenumbers,raman_spectra.iloc[1])
plt.plot(wavenumbers,baseline_data.iloc[1])

Answer 1

根据Christian K.的建议，我研究了背景估计的SNIP算法，详细信息可以参见例如here 。这是我的 python 代码:

import numpy as np
import pandas as pd
from scipy.signal import gaussian
import matplotlib.pyplot as plt

def baseline_correction(raman_spectra,niter):

    assert(isinstance(raman_spectra, pd.DataFrame)), 'Input must be pandas DataFrame'

    spectrum_points = len(raman_spectra.columns)
    raman_spectra_transformed = np.log(np.log(np.sqrt(raman_spectra +1)+1)+1)

    working_spectra = np.zeros(raman_spectra.shape)

    for pp in np.arange(1,niter+1):
        r1 = raman_spectra_transformed.iloc[:,pp:spectrum_points-pp]
        r2 = (np.roll(raman_spectra_transformed,-pp,axis=1)[:,pp:spectrum_points-pp] + np.roll(raman_spectra_transformed,pp,axis=1)[:,pp:spectrum_points-pp])/2
        working_spectra = np.minimum(r1,r2)
        raman_spectra_transformed.iloc[:,pp:spectrum_points-pp] = working_spectra

    baseline = (np.exp(np.exp(raman_spectra_transformed)-1)-1)**2 -1
    return baseline

wavenumbers = np.linspace(500,2000,1000)
intensities1 = gaussian(1000,20) + 0.000002*wavenumbers**2
intensities2 = gaussian(1000,50) + 0.000001*wavenumbers**2
raman_spectra = pd.DataFrame((intensities1,intensities2),columns=np.around(wavenumbers,decimals=1))

iterations = 100
baseline_data = baseline_correction(raman_spectra,iterations)


#the rest is just for plotting the data
plt.figure(1)
plt.plot(wavenumbers,raman_spectra.iloc[0])
plt.plot(wavenumbers,baseline_data.iloc[0])
plt.plot(wavenumbers,raman_spectra.iloc[1])
plt.plot(wavenumbers,baseline_data.iloc[1])

它确实有效，并且看起来与基于非对称最小二乘平滑的算法同样可靠。它也更快。通过 100 次迭代，拟合 73 个真实的测量光谱大约需要 1.5 秒，结果总体良好，而相比之下，大约需要 1.5 秒。 2.2 对于非对称最小二乘平滑，所以是一种改进。

更好的是:所需的计算时间对于 3267 个真实光谱，使用 SNIP 算法只需 11.7 秒，而对于非对称最小二乘平滑则为 1 分 28 秒。这可能是由于 SNIP 算法没有使用任何 for 循环一次遍历每个频谱的结果。

A typical result of the SNIP algorithm with calculated examples is shown here 。

我对这个结果非常满意，感谢所有贡献者的支持!

更新:感谢 this question 中的 sascha ，我找到了一种使用不对称最小二乘平滑而不使用 for 循环迭代每个频谱的方法，基线校正函数如下所示:

def baseline_correction4(raman_spectra,lam,p,niter=10):
    #according to "Asymmetric Least Squares Smoothing" by P. Eilers and H. Boelens
    number_of_spectra = raman_spectra.index.size

    #this is the code for the fitting procedure        
    L = len(raman_spectra.columns)
    w = np.ones(raman_spectra.shape[0]*raman_spectra.shape[1])

    D = sparse.block_diag(np.tile(sparse.diags([1,-2,1],[0,-1,-2],shape=(L,L-2)),number_of_spectra),format='csr')

    raman_spectra_flattened = raman_spectra.values.ravel()

    for jj in range(int(niter)):
        W = sparse.diags(w,format='csr')
        Z = W + lam * D.dot(D.transpose())
        z = spsolve(Z,w*raman_spectra_flattened,permc_spec='NATURAL')
        w = p * (raman_spectra_flattened > z) + (1-p) * (raman_spectra_flattened < z)
    #end of fitting procedure

    baseline_data = pd.DataFrame(z.reshape(number_of_spectra,-1),index=raman_spectra.index,columns=raman_spectra.columns)
    return baseline_data

这种方法基于将所有稀疏矩阵组合成一个 block 对角稀疏矩阵。这样，无论您有多少个光谱，您都只需调用 spsolve 一次。这导致在 593 毫秒内对 73 个真实光谱进行基线校正(比 SNIP 快)，在 32.8 秒内对 3267 个真实光谱进行基线校正(比 SNIP 慢)。我希望这对将来的人有用。