Эти полосы спектра раньше судили по глазу, как это сделать программно?

enter image description here

Операторы использовали для изучения спектра, зная местоположение иwidth каждой вершины и оцените часть, к которой принадлежит спектр. По-новому, изображение захватывается камерой на экран. И ширина каждой полосы должна быть рассчитана программно.

Старая система: спектроскоп - & gt; человеческий глаз Новая система: спектроскоп - & gt; камера - & gt; программа

Что является хорошим методом дляcompute the width of each band, учитывая их приблизительное положение по оси 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()
 Jesvin Jose26 мая 2012 г., 11:22
@ детально, я не совсем информирован, но положение вершин было бы известно заранее. И да, меня беспокоит, находятся ли два известных пика слишком близко друг к другу.
 Blender26 мая 2012 г., 09:56
Найти порог, при котором человеческий глаз не может отличить цвет от фона. Это должно быть довольно постоянным. Затем просто определите пороговые значения, чтобы те, которые превышают это пороговое значение, были "видны". человеческим глазом и просто скопировать эти точки данных и найти ширину каждого кластера.
 Jesvin Jose26 мая 2012 г., 11:01
@ Блендер, не могли бы вы уточнить вашу точку зрения?
 detly26 мая 2012 г., 11:05
Является ли уровень шума достаточно плоским? Как вы определяете (на глаз), что является пиком, а что слишком низким? Например, является ли этот всплеск около 660 × 670 дублетом или не будет засчитан? Как насчет 740?
 Lennart Regebro26 мая 2012 г., 10:17
Может быть, взять производную кривой, чтобы легко найти пики и склоны? Разве ширина полос довольно постоянна? Это амплитуды и положения, которые меняются, насколько я понимаю.

Ответы на вопрос(3)

результатами, полученными человеком.

Вы должны взять большое разнообразие данных и большое разнообразие оценок измерений (ширины при различных порогах, площадь выше различных порогов, разные методы выбора порогов, 2-е моменты, подборы полиномиальных кривых различной степени, сопоставление с образцом и т. Д.) И сравнить их оценки человеческого измерения одного и того же набора данных. Выберите метод оценки, который лучше всего соотносится с результатами эксперта. Или, может быть, выбрать несколько методов, наилучший для каждой из различных высот, для различных отделений от других пиков и т. Д.

но для любого, кто столкнется с этим вопросом в будущем ...

Данные о движении глаз выглядят очень похоже на это; Я основываю подход, который использовалНистром + Холмквист, 2010, Сгладить данные с помощью фильтра Савицкого-Голея (scipy.signal.savgol_filter в scipy v0.14 +), чтобы избавиться от некоторых низкоуровневых шумов при сохранении больших пиков без изменений - авторы рекомендуют использовать порядок 2 и размер окна примерно вдвое больше ширины наименьшего пика, которым вы хотите быть в состоянии обнаружить. Вы можете найти, где находятся полосы, путем произвольного удаления всех значений выше определенного значения y (установите их наnumpy.nan). Затем возьмите (nan) среднее и (nan) стандартное отклонение от остатка и удалите все значения, превышающие среднее значение + [параметр] * std (я думаю, что они используют 6 в статье). Повторяйте, пока вы не удаляете какие-либо точки данных, но в зависимости от ваших данных некоторые значения [параметра] могут не стабилизироваться. Тогда используйтеnumpy.isnan() чтобы найти события против не-событий, иnumpy.diff() найти начало и конец каждого события (значения -1 и 1 соответственно). Чтобы получить еще более точные начальную и конечную точки, вы можете сканировать данные в обратном направлении от каждого начала и вперед от каждого конца, чтобы найти ближайший локальный минимум, значение которого меньше среднего + [другой параметр] * std (я думаю, что они используют 3 в газете). Тогда вам просто нужно посчитать точки данных между каждым началом и концом.

Это не сработает для этого двойного пика; Вы должны сделать некоторую экстраполяцию для этого.

Решение Вопроса

вы можете использовать простой алгоритм, который находит локальные максимумы, наиболее близкие к этой точке. Ваш подходящий код может уже делать это (я не был уверен, успешно ли вы его использовали или нет).

Вот некоторый код, который демонстрирует простое нахождение пика из заданной пользователем начальной точки:

#!/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()

Этот метод будет работать, только если восхождение на вершинуperfectly smooth с вашей отправной точки. Если это должно быть более устойчивым к шуму, я не использовал его, ноPyDSTool Кажется, это может помочь. этоSciPy сообщение подробно, как использовать его для обнаружения пиков 1D в наборе данных с шумом.

Поэтому предположим, что в этот момент вы нашли центр пика. Теперь для ширины: есть несколько методов, которые вы могли бы использовать, но самый простой - это, вероятно, «полная ширина при половине максимума»; (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()

Это не обязательно должна быть полная ширина вhalf максимум & # x2014; как указывает один из комментаторов, вы можете попытаться выяснить, где находятся ваши операторы & apos; нормальный порог для обнаружения пика есть, и превратить его в алгоритм для этого шага процесса.

Более надежным способом может быть подгонка кривой Гаусса (или вашей собственной модели) к подмножеству данных, центрированных вокруг пика & # x2014; скажем, от локальных минимумов с одной стороны до локальных минимумов с другой & # x2014; и использовать один из параметров этой кривой (например, сигма) для расчета ширины.

Я понимаю, что это много кода, но я сознательно избегал выделять функции поиска индексов, чтобы "показать мою работу" немного больше, и, конечно, функции построения графиков предназначены только для демонстрации.

Надеюсь, это даст вам хотя бы хорошую отправную точку для того, чтобы придумать что-то более подходящее для вашего конкретного набора.

 08 июн. 2012 г., 07:20
@aitchnyu Ах, я понимаю тебя. У меня нет времени, чтобы сделать это самому, но вы могли бы реализовать то, что смотрит на первое различие в диапазоне пика. Оно должно быть строго положительным слева от пика и строго отрицательным справа. По определению, он равен нулю прямо на пике, хотя, возможно, вам следует также учитывать нулевые значения с обеих сторон.
 Jesvin Jose07 июн. 2012 г., 08:08
Не могли бы вы включить модификацию, которая выдает исключение, если & quot; base & quot; ниже порога (0,5 пика в этом случае)?
 Jesvin Jose08 июн. 2012 г., 07:18
Да, и я хочу, чтобы алгоритм генерировал исключение, если он встречает плохо сформированный пик: пик, который не падает монотонно при пиковом значении до 0,5 (см. Дублет на моем графике). А мне нужноonly ширина, амплитуда может измен тьс из-за свойств прибора. «Уровень шума» должен быть идеально тихим, соответствующим интенсивности черного.
 Jesvin Jose08 июн. 2012 г., 07:30
Такhm_left_index:index_of_peak первого различия должно быть все положительным иindex_of_peak:hm_right_index должно быть все отрицательно, верно?
 08 июн. 2012 г., 04:36
@aitchnyu Хорошо, здесь я сделал несколько упрощений: я только что сказал 0,5 от общей высоты, но FWHM лучше измерять как половину высоты пика.above the base (или уровень шума). Таким образом, не имеет смысла, чтобы основание было ниже 0,5 от пика.

Ваш ответ на вопрос