Estas bandas de espectro solían juzgarse a ojo, ¿cómo hacerlo programáticamente?

Question

operadores utilizados para examinar el espectro, conociendo la ubicación y anchura de cada pico y juzgar la pieza del espectro pertenece. De la nueva manera, la cámara captura una imagen en una pantalla. Y el ancho de cada banda debe calcularse programáticamente.

sistema antiguo: espectroscopio -> ojo humano Nuevo sistema: espectroscopio -> cámara -> programa de

¿Qué es un buen método para calcular el ancho de cada banda, dadas sus posiciones aproximadas del eje X; dado que esta tarea solía realizarse perfectamente a simple vista, y ahora debe realizarse por programa?

Lo siento si me faltan detalles, pero son escasos.

---

Lista de programas que generó el gráfico anterior; Espero que sea relevante:

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

Dado un punto de partida aproximada, se puede usar un algoritmo simple que encuentra un máximo local más cercano a este punto. Su código de ajuste ya puede estar haciendo eso (no estaba seguro de si lo estaba usando con éxito o no).

Aquí hay algo de código que demuestra hallazgo pico simple a partir de un determinado punto de partida por el usuario:

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

Este método sólo funciona si el ascenso hasta el pico es perfectamente lisa desde su punto de partida. Si esto necesita ser más resistente al ruido, no lo he usado, pero PyDSTool parece que podría ayudar. Este SciPy post detalla cómo usarlo para detectar picos 1D en un conjunto de datos ruidoso.

Asuma que en este punto ha encontrado el centro del pico. Ahora para el ancho: hay varios métodos que puede usar, pero el más fácil es, probablemente, el "ancho completo a la mitad de máximo" (FWHM). De nuevo, esto es simple y, por lo tanto, frágil. Se romperá por dobles picos cercanos, o por datos ruidosos.

El FWHM es exactamente lo que su nombre sugiere: usted encuentra el ancho del pico si está a medio camino del máximo. Aquí hay algo de código que hace que (sólo continúa en la de arriba):

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

No tiene por qué ser la anchura total a medio máximo - como un comentarista señala, se puede tratar de averiguar donde es el umbral normal de los operadores para la detección de picos, y conviértalo en un algoritmo para este paso del proceso.

Una forma más robusta podría ser ajustar una curva gaussiana (o su propio modelo) a un subconjunto de los datos centrados en el pico, por ejemplo, desde un mínimo local en un lado hasta un mínimo local en el otro y utilice uno de los parámetros de esa curva (por ejemplo, sigma) para calcular el ancho.

Me doy cuenta de que esto es una gran cantidad de código, pero he evitado deliberadamente las funciones de búsqueda de índices para "mostrar mi trabajo" un poco más, y por supuesto las funciones de trazado están ahí para demostrar.

Afortunadamente, esto le brinda al menos un buen punto de partida para encontrar algo más adecuado para su conjunto en particular.

Answer 2

El mejor método podría ser comparar estadísticamente un grupo de métodos con resultados humanos.

Tomaría una gran variedad de datos y una gran variedad de estimaciones de medición (anchos en varios umbrales, área por encima de varios umbrales, diferentes métodos de selección de umbrales, segundos momentos, ajustes de curvas polinomiales de varios grados, coincidencia de patrones, etc.) y comparar estas estimaciones con mediciones humanas del mismo conjunto de datos. Elija el método de estimación que se correlacione mejor con los resultados humanos expertos. O tal vez recoger varios métodos, el mejor para cada una de las varias alturas, por diversas separaciones de otros picos, y etc

Answer 3

tarde a la fiesta, pero para cualquiera que venga a través de esta cuestión en el futuro ...

Los datos del movimiento de los ojos son muy similares a esto; Basaría un enfoque en el utilizado por Nystrom + Holmqvist, 2010. Alise los datos usando un filtro Savitsky-Golay (scipy.signal.savgol_filter en scipy v0.14 +) para eliminar parte del ruido de bajo nivel mientras se mantienen intactos los picos grandes. Los autores recomiendan usar un orden de 2 y un tamaño de ventana de aproximadamente el doble del ancho del pico más pequeño que desea ser capaz de detectar. Puede encontrar dónde están las bandas eliminando arbitrariamente todos los valores por encima de cierto valor y (configúrelos en numpy.nan). Luego tome la (nan) media y (nan) desviación estándar del resto, y elimine todos los valores mayores que la media + [parámetro] * std (creo que usan 6 en el papel). Itere hasta que no elimine ningún punto de datos, pero dependiendo de sus datos, es posible que ciertos valores de [parámetro] no se estabilicen. Luego use numpy.isnan() para encontrar eventos vs no eventos, y numpy.diff() para encontrar el inicio y el final de cada evento (valores de -1 y 1 respectivamente). Para obtener puntos de inicio y finalización aún más precisos, puede escanear los datos hacia atrás desde cada inicio y reenvío desde cada extremo para encontrar el mínimo local más cercano que tenga un valor menor que la media + [otro parámetro] * estándar (creo que utilizan 3 en el papel). Entonces solo necesita contar los puntos de datos entre cada inicio y final.

Esto no funcionará para ese doble pico; tendrías que hacer una extrapolación para eso.