python - 从两个未归一化的光谱创建复合光谱

Question

在 Python 中，我试图规范化两个数组，然后取它们重叠区域的平均值来创建一个新的复合数组。

为此，我认为我必须:

1. 找到重叠区域，
2. 插入重叠的 y 值，
3. 遍历找到最佳拟合的归一化常数，然后
4. 将各个部分粘贴在一起形成我的新曲线

有了一些半随机值，这就是它的样子:

此代码适用于 y 值相距不太远的小型数据集，但当 Y1 和 Y2 之间存在数量级时(显然是由于迭代)，Python 会崩溃。这是代码:

X1o = [x for x in X1 if x > X2[0]]
X2o = [x for x in X2 if x < X1[-1]]
Y1o = [y for y in Y1[(len(Y1)-len(X1o)):]]
Y2o = [y for y in Y2[:len(X2o)]]
Y2o = list(interp(X1o,X2o,Y2o))

c = abs(min(Y1o)-max(Y2o))
Y2test = [y2+c for y2 in Y2o]
Y2s = []
d = 0.01*min(Y2test)
while min(Y2test) < max(Y1o):
  Y2test = [y+d for y in Y2test]
  Y2s.append(Y2test)
  plot(X1o,Y2test,c='k',alpha=0.5)

idx = min(map(lambda i: (u.squaredError(Y1o, i), i, Y2s.index(i)), Y2s))[-1]         
Yavg = [(y1+y2)/2 for y1,y2 in zip(Y1o,Y2s[idx])]
diff = Y2s[idx][0]-Y2o[0]

X = [x for x in X1 if x < X2[0]] + X1o + [x for x in X2 if x > X1[-1]]
Y = [y for x,y in zip(X1,Y1) if x < X2[0]] + Yavg + [y+diff for x,y in zip(X2,Y2) if x > X1[-1]]

我确实需要使用具有数千个数据点和 y 值之间最多 20 个数量级分布的恒星光谱来执行此操作。

如有任何建议，我们将不胜感激!

Answer 1

您的代码将从 numpy 和使用较少的 python 列表中受益匪浅，这些列表效率低下，尤其是您的代码行 Y2s.append(Y2test)。当您的 while 周期太长时，您只会追加到一个很长的列表，这既慢又低效。

话虽如此，代码的瓶颈在于最小化。您目前正在使用 python 列表进行暴力破解。您将从使用 scipy.optimize 中的一个中获益匪浅功能。

以下是我将要做什么的一些广泛建议:

1. 找到两个光谱的 x 坐标极值，将两者插值到公共(public) x 值的网格。
2. 使用 scipy.optimize.fmin 的风格为您进行最小化并计算最佳归一化。
3. 将部分归一化光谱插值到公共(public)网格

下面是一些使用 fmin 的示例代码(未经测试):

import numpy as np
import scipy.optimize as opt

# y1 = interpolated values for one of the spectra
# y2 = interpolated values for the other spectra, normalise this one 

def errfunc(p, a1, a2):
    return np.sum(a1 - a2 * p)

p0 = 1.  # initial guess
norm_factor = opt.fmin(errfunc, p0, args=(y1, y2))

这应该为您提供最合适的 norm_factor。