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

introduzca la descripción de la imagen aquí

Los operadores solían examinar el espectro, conociendo la ubicación y el ancho de cada pico y juzgar la pieza a la que pertenece el espectro. De la nueva forma, la imagen es capturada por una cámara a una pantalla. Y el ancho de cada banda debe ser computado programáticamente.

Sistema antiguo: espectroscopio -> ojo humano Nuevo sistema: espectroscopio -> cámara -> progtwig

¿Cuál 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 ser realizada a simple vista, y ahora debe ser realizada por progtwig?

Lo siento si me faltan detalles, pero son escasos.


    Listado de progtwigs 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(-(xp[1])**2/(2.0*p[2]**2)) + \ p[3]*exp(-(xp[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() 

    Dado un punto de partida aproximado, podría usar un algoritmo simple que encuentre un máximo local más cercano a este punto. Es posible que su código de ajuste ya esté haciendo eso (no estaba seguro de si lo estaba usando con éxito o no).

    Aquí hay algo de código que demuestra la búsqueda de picos simples desde un punto de inicio dado 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 solo funcionará si el ascenso hasta el pico es perfectamente suave desde su punto de partida. Si esto necesita ser más resistente al ruido, no lo he usado, pero parece que PyDSTool podría ayudar. Esta publicación de SciPy detalla cómo usarla para detectar picos 1D en un conjunto de datos ruidosos.

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

    El FWHM es exactamente lo que su nombre sugiere: encuentra el ancho del pico donde está a mitad de camino al máximo. Aquí hay un código que hace eso (simplemente continúa desde 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 que ser el ancho total a la mitad del máximo - como lo señala un comentarista, puede intentar averiguar cuál es el umbral normal de los operadores para la detección de picos, y convertirlo 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 alrededor del pico (por ejemplo, desde un mínimo local en un lado a un mínimo local en el otro) y usar uno de los parámetros de esa curva (ej. sigma) para calcular el ancho.

    Me doy cuenta de que esto es un montón de código, pero he evitado deliberadamente factorizar 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í solo para demostrar.

    Esperemos que esto le brinde al menos un buen punto de partida para encontrar algo más adecuado para su conjunto en particular.

    Tarde a la fiesta, pero para cualquiera que se encuentre con esta pregunta en el futuro …

    Los datos del movimiento ocular se parecen mucho a esto; Me basaría un enfoque en el utilizado por Nystrom + Holmqvist, 2010 . Alise los datos con un filtro Savitsky-Golay ( scipy.signal.savgol_filter en scipy v0.14 +) para eliminar parte del ruido de bajo nivel mientras mantiene intactos los picos grandes: los autores recomiendan usar un orden de 2 y un Tamaño de la ventana de aproximadamente el doble del ancho del pico más pequeño que desea detectar. Puede encontrar dónde están las bandas eliminando arbitrariamente todos los valores por encima de un determinado valor y ( numpy.nan en numpy.nan ). Luego tome la media (nan) y la desviación estándar (nan) del rest, y elimine todos los valores mayores que la media + [parámetro] * std (creo que usan 6 en el documento). Iterar hasta que no esté eliminando ningún punto de datos, pero dependiendo de sus datos, ciertos valores de [parámetro] pueden no estabilizarse. Luego use numpy.isnan() para encontrar eventos frente a 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 hacia adelante desde cada extremo para encontrar el mínimo local más cercano que tenga un valor más pequeño que la media + [otro parámetro] * std (creo que usan 3 en el papel). Entonces solo necesitas contar los puntos de datos entre cada inicio y final.

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

    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 diversos grados, coincidencia de patrones, etc.) y los compararía Estimaciones a medidas humanas del mismo conjunto de datos. Elija el método de estimación que se correlaciona mejor con los resultados humanos expertos. O tal vez elija varios métodos, el mejor para cada una de varias alturas, para varias separaciones de otros picos, etc.