python - 这些频谱波段以前是靠肉眼判断的，怎么用程序来判断呢？

Question

运算符(operator)过去常常检查光谱，知道每个峰的位置和宽度，并判断光谱属于哪个部分。在新的方式中，图像由相机捕获到屏幕上。并且必须以编程方式计算每个波段的宽度。

旧系统:分光镜->人眼新系统:分光镜->相机->程序

什么是计算每个带的宽度的好方法，给定它们的近似 X 轴位置；考虑到这个任务过去可以用眼睛完美地完成，现在必须由程序来完成？

对不起，如果我缺少细节，但它们很少。

---

生成上一个图的程序列表；我希望它是相关的:

import Image
from scipy import *
from scipy.optimize import leastsq

# Load the picture with PIL, process if needed
pic         = asarray(Image.open("spectrum.jpg"))

# Average the pixel values along vertical axis
pic_avg     = pic.mean(axis=2)
projection  = pic_avg.sum(axis=0)

# Set the min value to zero for a nice fit
projection /= projection.mean()
projection -= projection.min()

#print projection

# Fit function, two gaussians, adjust as needed
def fitfunc(p,x):
    return p[0]*exp(-(x-p[1])**2/(2.0*p[2]**2)) + \
        p[3]*exp(-(x-p[4])**2/(2.0*p[5]**2))
errfunc = lambda p, x, y: fitfunc(p,x)-y

# Use scipy to fit, p0 is inital guess
p0 = array([0,20,1,0,75,10])
X  = xrange(len(projection))
p1, success = leastsq(errfunc, p0, args=(X,projection))
Y = fitfunc(p1,X)

# Output the result
print "Mean values at: ", p1[1], p1[4]

# Plot the result
from pylab import *
#subplot(211)
#imshow(pic)
#subplot(223)
#plot(projection)
#subplot(224)
#plot(X,Y,'r',lw=5)
#show()

subplot(311)
imshow(pic)
subplot(312)
plot(projection)
subplot(313)
plot(X,Y,'r',lw=5)
show()

Answer 1

给定一个近似起点，您可以使用一个简单的算法来找到最接近该点的局部最大值。您的拟合代码可能已经这样做了(我不确定您是否成功使用了它)。

下面是一些代码，演示了从用户给定的起点进行简单的峰值查找:

#!/usr/bin/env python
from __future__ import division
import numpy as np
from matplotlib import pyplot as plt

# Sample data with two peaks: small one at t=0.4, large one at t=0.8
ts = np.arange(0, 1, 0.01)
xs = np.exp(-((ts-0.4)/0.1)**2) + 2*np.exp(-((ts-0.8)/0.1)**2)

# Say we have an approximate starting point of 0.35
start_point = 0.35

# Nearest index in "ts" to this starting point is...
start_index = np.argmin(np.abs(ts - start_point))

# Find the local maxima in our data by looking for a sign change in
# the first difference
# From http://stackoverflow.com/a/9667121/188535
maxes = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1

# Find which of these peaks is closest to our starting point
index_of_peak = maxes[np.argmin(np.abs(maxes - start_index))]

print "Peak centre at: %.3f" % ts[index_of_peak]

# Quick plot showing the results: blue line is data, green dot is
# starting point, red dot is peak location
plt.plot(ts, xs, '-b')
plt.plot(ts[start_index], xs[start_index], 'og')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.show()

此方法仅在从您的起点完全平稳攀登山峰时才有效。如果这需要对噪音更有弹性，我没有使用它，但是 PyDSTool似乎它可能有帮助。这SciPy post详细介绍了如何使用它来检测嘈杂数据集中的一维峰。

假设此时您已经找到了峰的中心。现在是宽度:您可以使用多种方法，但最简单的方法可能是“半峰全宽”(FWHM)。同样，这很简单，因此很脆弱。它会因接近的双峰或嘈杂的数据而中断。

FWHM 顾名思义:您可以找到峰值宽度一半的最大值。这是执行此操作的一些代码(它只是从上面继续):

# FWHM...
half_max = xs[index_of_peak]/2

# This finds where in the data we cross over the halfway point to our peak. Note
# that this is global, so we need an extra step to refine these results to find
# the closest crossovers to our peak.

# Same sign-change-in-first-diff technique as above
hm_left_indices = (np.diff(np.sign(np.diff(np.abs(xs[:index_of_peak] - half_max)))) > 0).nonzero()[0] + 1
# Add "index_of_peak" to result because we cut off the left side of the data!
hm_right_indices = (np.diff(np.sign(np.diff(np.abs(xs[index_of_peak:] - half_max)))) > 0).nonzero()[0] + 1 + index_of_peak

# Find closest half-max index to peak
hm_left_index = hm_left_indices[np.argmin(np.abs(hm_left_indices - index_of_peak))]
hm_right_index = hm_right_indices[np.argmin(np.abs(hm_right_indices - index_of_peak))]

# And the width is...    
fwhm = ts[hm_right_index] - ts[hm_left_index]

print "Width: %.3f" % fwhm

# Plot to illustrate FWHM: blue line is data, red circle is peak, red line
# shows FWHM
plt.plot(ts, xs, '-b')
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or')
plt.plot(
    [ts[hm_left_index], ts[hm_right_index]],
    [xs[hm_left_index], xs[hm_right_index]], '-r')
plt.show()

它不一定是最大一半的全宽——正如一位评论者指出的那样，您可以尝试找出您的运算符(operator)的正常峰值检测阈值在哪里，然后将其转向转化为该过程的这一步的算法。

一种更稳健的方法可能是将高斯曲线(或您自己的模型)拟合到以峰值为中心的数据子集——例如，从一侧的局部最小值到另一侧的局部最小值——并使用该曲线的参数之一(例如西格玛)来计算宽度。

我意识到这是很多代码，但我故意避免分解出索引查找函数以更多地“展示我的工作”，当然绘图函数只是为了演示。

希望这至少能给你一个好的起点，让你想出更适合你的特定场景的东西。