2014-07-09 14 views
6

Ho uno spettro TOF e vorrei implementare un algoritmo utilizzando python (numpy) che trova tutti i massimi dello spettro e restituisce i corrispondenti valori x.
Ho cercato online e ho trovato l'algoritmo riportato di seguito.
trova la posizione dei picchi in uno spettro numpy

L'assunto è che vicino al massimo la differenza tra il valore precedente e il valore al massimo è maggiore di un numero DELTA. Il problema è che il mio spettro è composto da punti equamente distribuiti, anche vicino al massimo, così che DELTA non viene mai superato e la funzione peakdet restituisce una matrice vuota.

Avete qualche idea su come superare questo problema? Gradirei davvero i commenti per capire meglio il codice dato che sono abbastanza nuovo in Python.

Grazie!

import sys 
from numpy import NaN, Inf, arange, isscalar, asarray, array 

def peakdet(v, delta, x = None): 
    maxtab = [] 
    mintab = [] 

    if x is None: 
     x = arange(len(v)) 
     v = asarray(v) 

    if len(v) != len(x): 
     sys.exit('Input vectors v and x must have same length') 
    if not isscalar(delta): 
     sys.exit('Input argument delta must be a scalar') 
    if delta <= 0: 
     sys.exit('Input argument delta must be positive') 

    mn, mx = Inf, -Inf 
    mnpos, mxpos = NaN, NaN 

    lookformax = True 

    for i in arange(len(v)): 
     this = v[i] 
     if this > mx: 
     mx = this 
     mxpos = x[i] 
     if this < mn: 
     mn = this 
     mnpos = x[i] 

     if lookformax: 
     if this < mx-delta: 
      maxtab.append((mxpos, mx)) 
      mn = this 
      mnpos = x[i] 
      lookformax = False 
     else: 
     if this > mn+delta: 
      mintab.append((mnpos, mn)) 
      mx = this 
      mxpos = x[i] 
      lookformax = True 

return array(maxtab), array(mintab) 

Di seguito viene mostrata parte dello spettro. In realtà ho più picchi di quelli mostrati qui.

enter image description here

+0

Correggere il seguente: questa> mn + Delta e questo (mn + delta) e questo <(mx -delta) – TommasoF

+0

Il codice era senza parentesi. Tuttavia, anche con loro non cambia molto. Ho ancora un array vuoto. – diegus

+0

Non puoi semplicemente usare la convolvola invece e cercare tutti gli zero-crossing con un primo kernel derivativo adatto? – deinonychusaur

risposta

10

Questo, penso, potrebbe funzionare come punto di partenza. Io non sono un esperto di elaborazione del segnale, ma ho provato questo su un segnale generato Y che sembra piuttosto come la tua e uno con molto più rumore:

from scipy.signal import convolve 
import numpy as np 
from matplotlib import pyplot as plt 
#Obtaining derivative 
kernel = [1, 0, -1] 
dY = convolve(Y, kernel, 'valid') 

#Checking for sign-flipping 
S = np.sign(dY) 
ddS = convolve(S, kernel, 'valid') 

#These candidates are basically all negative slope positions 
#Add one since using 'valid' shrinks the arrays 
candidates = np.where(dY < 0)[0] + (len(kernel) - 1) 

#Here they are filtered on actually being the final such position in a run of 
#negative slopes 
peaks = sorted(set(candidates).intersection(np.where(ddS == 2)[0] + 1)) 

plt.plot(Y) 

#If you need a simple filter on peak size you could use: 
alpha = -0.0025 
peaks = np.array(peaks)[Y[peaks] < alpha] 

plt.scatter(peaks, Y[peaks], marker='x', color='g', s=40) 

I risultati dei campioni: Smooth curve Per quello rumoroso, ho filtrato picchi con alpha: Rough curve

Se il alpha ha bisogno di più sofisticato si potrebbe provare a impostare dinamicamente alpha dalle vette scoperti utilizzando ad esempio le ipotesi su di loro sono un misto gaussiano (il mio preferito è la soglia Otsu, esiste in cv e skimage) o una sorta di clustering (k-means potrebbe funzionare).

E per riferimento, questo ho usato per generare il segnale:

Y = np.zeros(1000) 

def peaker(Y, alpha=0.01, df=2, loc=-0.005, size=-.0015, threshold=0.001, decay=0.5): 
    peaking = False 
    for i, v in enumerate(Y): 
     if not peaking: 
      peaking = np.random.random() < alpha 
      if peaking: 
       Y[i] = loc + size * np.random.chisquare(df=2) 
       continue 
     elif Y[i - 1] < threshold: 
      peaking = False 

     if i > 0: 
      Y[i] = Y[i - 1] * decay 

peaker(Y) 

EDIT: Supporto per degradanti base-line

ho simulato una base-line oblique in questo modo:

Z = np.log2(np.arange(Y.size) + 100) * 0.001 
Y = Y + Z[::-1] - Z[-1] 

Poi da rilevare con un alfa fissa (nota che ho cambiato firmare alpha):

from scipy.signal import medfilt 

alpha = 0.0025 
Ybase = medfilt(Y, 51) # 51 should be large in comparison to your peak X-axis lengths and an odd number. 
peaks = np.array(peaks)[Ybase[peaks] - Y[peaks] > alpha] 

risultato i seguenti risultati (linea di base viene tracciata come linea nera tratteggiata): Degrading base-line

EDIT 2: semplificazione e commenti

I semplificare la codice per utilizzare un kernel per entrambi gli convolve come commentato @skymandr. Questo ha anche rimosso il numero magico nel regolare il restringimento in modo che qualsiasi dimensione del kernel dovesse fare.

Per la scelta del "valid" come opzione per convolve. Si avrebbe probabilmente funzionato altrettanto bene con "same", ma ho scelto "valid" quindi non ho dovuto pensare alle edge-condizioni e se l'algoritmo in grado di rilevare i picchi spurios lì.

+0

Grazie mille! Lo userò come punto di partenza. Il problema con i dati reali è che la linea di base non è costante, cioè non è sempre intorno allo zero, ma scende al tempo maggiore rendendo difficile l'uso di alfa (insieme con il rumore anche picchi di bassa intensità sono tagliati). Cercherò di impostare dinamicamente alfa .. – diegus

+0

Dal momento che la maggior parte del vostro segnale è la vostra linea di base (anche se cadere fuori), si potrebbe usare un 'signal.medfilt' con la grande dimensione del kernel e impostare alfa per quanto ci si aspetterebbe per il' x 'di ogni singolo picco. – deinonychusaur

+1

Dal kernel = [1, -1] trova tecnicamente il (denormalizzati) derivato in posizioni tra i datapoints, piuttosto che ai datapoints, vi suggerirei di usare kernel2 sia per trovare il derivato e per il segno flipping. Ciò eviterà un pregiudizio che, anche se piccolo, potrebbe essere difficile eseguire il debug. Per i dati altamente risolti e ben educati come questo, il risultato dovrebbe essere altrettanto buono. – skymandr

0

Trovare un minimo o un massimo non è così semplice, perché non esiste una definizione universale per "massimo locale".

Il tuo codice sembra cercare un miximum e poi accettarla come massimo se il segnale cade dopo il massimo al di sotto del valore di un certo delta massimo meno. Dopo di ciò inizia a cercare un minimo con criteri simili. Non importa se i tuoi dati diminuiscono o aumentano lentamente, poiché il massimo viene registrato quando viene raggiunto e aggiunto alla lista dei massimi una volta che il livello si abbassa al di sotto della soglia di isteresi.

Questo è un possibile modo per trovare minimi locali e massimi, ma ha diversi difetti. Uno di questi è che il metodo non è simmetrico, cioè se gli stessi dati vengono eseguiti all'indietro, i risultati non sono necessariamente gli stessi.

Purtroppo, non posso aiutare molto di più, perché il metodo corretto dipende molto dai dati che si stanno guardando, dalla sua forma e dalla sua rumorosità. Se si dispone di alcuni campioni, potremmo essere in grado di fornire alcuni suggerimenti.

Problemi correlati