我目前正在处理一些拉曼光谱数据，并且正在尝试纠正由荧光倾斜引起的数据。请看下图：
我已经非常接近实现我想要的了。如您所见，我试图在我的所有数据中拟合多项式，而我实际上应该只在局部最小值处拟合多项式。
理想情况下，我希望有一个多项式拟合，从我的原始数据中减去它会导致如下结果：
是否有任何内置库已经这样做了？
如果没有，有什么简单的算法可以推荐给我吗？
【问题讨论】：
rfft() 转换您的信号并将低频部分设置为零来设计一个高路径滤波器。
标签： python numpy scipy signal-processing
我找到了我的问题的答案，只是分享给所有偶然发现这个问题的人。
P. Eilers 和 H. Boelens 在 2005 年提出了一种名为“Asymmetric Least Squares Smoothing”的算法。该论文是免费的，您可以在 google 上找到它。
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 xrange(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
【讨论】：
z 是基线？所以最后需要用那个数组减去光谱来校正基线？
以下代码适用于 Python 3.6。
这是改编自公认的正确答案，以避免密集矩阵diff 计算（这很容易导致内存问题）并使用range（不是xrange）
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
def baseline_als(y, lam, p, niter=10):
L = len(y)
D = sparse.diags([1,-2,1],[0,-1,-2], shape=(L,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
【讨论】：
最近，我需要使用这种方法。答案中的代码运行良好，但显然过度使用了内存。所以，这是我优化内存使用的版本。
def baseline_als_optimized(y, lam, p, niter=10):
L = len(y)
D = sparse.diags([1,-2,1],[0,-1,-2], shape=(L,L-2))
D = lam * D.dot(D.transpose()) # Precompute this term since it does not depend on `w`
w = np.ones(L)
W = sparse.spdiags(w, 0, L, L)
for i in range(niter):
W.setdiag(w) # Do not create a new matrix, just update diagonal values
Z = W + D
z = spsolve(Z, w*y)
w = p * (y > z) + (1-p) * (y < z)
return z
根据我下面的基准，它也快了大约 1.5 倍。
%%timeit -n 1000 -r 10 y = randn(1000)
baseline_als(y, 10000, 0.05) # function from @jpantina's answer
# 20.5 ms ± 382 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)
%%timeit -n 1000 -r 10 y = randn(1000)
baseline_als_optimized(y, 10000, 0.05)
# 13.3 ms ± 874 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)
注意 1： 原文说：
为了强调算法的基本简单性，迭代次数已固定为 10。在实际应用中，应检查权重是否有任何变化；如果没有，则已经达到收敛。
所以，这意味着停止迭代更正确的方法是检查||w_new - w|| < tolerance
注意 2： 另一个有用的引用（来自@glycoaddict 的评论）给出了如何选择参数值的想法。
有两个参数：p 表示不对称性，λ 表示平滑度。两者都必须是 调整到手头的数据。我们发现通常 0.001 ≤ p ≤ 0.1 是一个不错的选择（对于具有正峰值的信号）和 102 ≤ λ ≤ 109，但可能会出现例外情况。在任何情况下，都应该在对 log λ 近似线性的网格上改变 λ。目测通常足以获得良好的参数值。
【讨论】：
有一个 python 库可用于基线校正/删除。它具有 Modpoly、IModploy 和 Zhang 拟合算法，当您将原始值输入为 python 列表或 pandas 系列并指定多项式次数时，可以返回基线校正结果。
将库安装为pip install BaselineRemoval。下面是一个例子
from BaselineRemoval import BaselineRemoval
input_array=[10,20,1.5,5,2,9,99,25,47]
polynomial_degree=2 #only needed for Modpoly and IModPoly algorithm
baseObj=BaselineRemoval(input_array)
Modpoly_output=baseObj.ModPoly(polynomial_degree)
Imodpoly_output=baseObj.IModPoly(polynomial_degree)
Zhangfit_output=baseObj.ZhangFit()
print('Original input:',input_array)
print('Modpoly base corrected values:',Modpoly_output)
print('IModPoly base corrected values:',Imodpoly_output)
print('ZhangFit base corrected values:',Zhangfit_output)
Original input: [10, 20, 1.5, 5, 2, 9, 99, 25, 47]
Modpoly base corrected values: [-1.98455800e-04 1.61793368e+01 1.08455179e+00 5.21544654e+00
7.20210508e-02 2.15427531e+00 8.44622093e+01 -4.17691125e-03
8.75511661e+00]
IModPoly base corrected values: [-0.84912125 15.13786196 -0.11351367 3.89675187 -1.33134142 0.70220645
82.99739548 -1.44577432 7.37269705]
ZhangFit base corrected values: [ 8.49924691e+00 1.84994576e+01 -3.31739230e-04 3.49854060e+00
4.97412948e-01 7.49628529e+00 9.74951576e+01 2.34940300e+01
4.54929023e+01
【讨论】：
我在之前的评论中使用了glinka 引用的算法版本，这是对相对recent paper 中发布的惩罚加权线性平方方法的改进。我用Rustam Guliev's 代码构建了这个：
from scipy import sparse
from scipy.sparse import linalg
import numpy as np
from numpy.linalg import norm
def baseline_arPLS(y, ratio=1e-6, lam=100, niter=10, full_output=False):
L = len(y)
diag = np.ones(L - 2)
D = sparse.spdiags([diag, -2*diag, diag], [0, -1, -2], L, L - 2)
H = lam * D.dot(D.T) # The transposes are flipped w.r.t the Algorithm on pg. 252
w = np.ones(L)
W = sparse.spdiags(w, 0, L, L)
crit = 1
count = 0
while crit > ratio:
z = linalg.spsolve(W + H, W * y)
d = y - z
dn = d[d < 0]
m = np.mean(dn)
s = np.std(dn)
w_new = 1 / (1 + np.exp(2 * (d - (2*s - m))/s))
crit = norm(w_new - w) / norm(w)
w = w_new
W.setdiag(w) # Do not create a new matrix, just update diagonal values
count += 1
if count > niter:
print('Maximum number of iterations exceeded')
break
if full_output:
info = {'num_iter': count, 'stop_criterion': crit}
return z, d, info
else:
return z
为了测试算法，我创建了一个类似于论文图 3 所示的光谱，首先生成一个由多个高斯峰组成的模拟光谱：
def spectra_model(x):
coeff = np.array([100, 200, 100])
mean = np.array([300, 750, 800])
stdv = np.array([15, 30, 15])
terms = []
for ind in range(len(coeff)):
term = coeff[ind] * np.exp(-((x - mean[ind]) / stdv[ind])**2)
terms.append(term)
spectra = sum(terms)
return spectra
x_vals = np.arange(1, 1001)
spectra_sim = spectra_model(x_vals)
然后，我使用直接取自论文的 4 个点创建了一个三阶插值多项式：
from scipy.interpolate import CubicSpline
x_poly = np.array([0, 250, 700, 1000])
y_poly = np.array([200, 180, 230, 200])
poly = CubicSpline(x_poly, y_poly)
baseline = poly(x_vals)
noise = np.random.randn(len(x_vals)) * 0.1
spectra_base = spectra_sim + baseline + noise
最后，我使用基线校正算法从改变的光谱中减去基线 (spectra_base)：
_, spectra_arPLS, info = baseline_arPLS(spectra_base, lam=1e4, niter=10,
full_output=True)
结果是（作为参考，我与Rustam Guliev's的纯ALS实现进行了比较，使用lam = 1e4和p = 0.001）：
【讨论】：
我知道这是一个老问题，但几个月前我偶然发现了它，并使用辣味.sparse 例程实现了等效的答案。
# Baseline removal
def baseline_als(y, lam, p, niter=10):
s = len(y)
# assemble difference matrix
D0 = sparse.eye( s )
d1 = [numpy.ones( s-1 ) * -2]
D1 = sparse.diags( d1, [-1] )
d2 = [ numpy.ones( s-2 ) * 1]
D2 = sparse.diags( d2, [-2] )
D = D0 + D2 + D1
w = np.ones( s )
for i in range( niter ):
W = sparse.diags( [w], [0] )
Z = W + lam*D.dot( D.transpose() )
z = spsolve( Z, w*y )
w = p * (y > z) + (1-p) * (y < z)
return z
干杯，
佩德罗。
【讨论】：