Skip to main content
 首页 » 编程设计

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

2025年05月04日7落叶无声

在 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 个数量级分布的恒星光谱来执行此操作。

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

请您参考如下方法:

您的代码将从 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。